[pure-lang-svn] SF.net SVN: pure-lang:[828] pure/trunk
Status: Beta
Brought to you by:
agraef
|
From: <ag...@us...> - 2008-09-22 19:20:32
|
Revision: 828
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=828&view=rev
Author: agraef
Date: 2008-09-22 19:20:24 +0000 (Mon, 22 Sep 2008)
Log Message:
-----------
Add operation to change the dimensions of a matrix without changing its size.
Modified Paths:
--------------
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/lib/primitives.pure 2008-09-22 19:20:24 UTC (rev 828)
@@ -447,7 +447,8 @@
nth x n = catch (cst {}) (row x n);
mth x m = catch (cst {}) (col x m);
end;
-x::matrix!!ns = colcatmap (nth x) ns with
+x::matrix!!ns = if packed x then rowvector x!!(1,ns)
+ else colcatmap (nth x) ns with
nth x n = catch (cst {}) {x!n};
end;
@@ -468,11 +469,28 @@
cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
-/* Extract a submatrix of a given size at a given offset. */
+/* Extract a submatrix of a given size at a given offset. The result shares
+ the underlying storage with the input matrix (i.e., matrix elements are
+ *not* copied) and so this is a comparatively cheap operation. */
submat x::matrix (i::int,j::int) (n::int,m::int)
= matrix_slice x i j (i+n-1) (j+m-1);
+/* Change the dimensions of a matrix without changing its size. The total
+ number of elements must match that of the input matrix. Reuses the
+ underlying storage of the input matrix if possible. */
+
+private matrix_redim;
+extern expr* matrix_redim(expr* x, int n, int m);
+
+redim x::matrix (n::int,m::int)
+ = matrix_redim x n m if n*m==#x;
+
+/* Convenience functions for converting a matrix to a row or column vector. */
+
+rowvector x::matrix = redim x (1,#x);
+colvector x::matrix = redim x (#x,1);
+
/* Construct matrices from lists of rows and columns. These take either
scalars or submatrices as inputs; corresponding dimensions must match.
rowcat combines submatrices vertically, like {x;y}; colcat combines them
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/runtime.cc 2008-09-22 19:20:24 UTC (rev 828)
@@ -4313,6 +4313,114 @@
}
extern "C"
+pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m)
+{
+ void *p = 0;
+ if (n<0 || m<0) return 0;
+ const size_t n1 = (size_t)n, n2 = (size_t)m;
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_symbolic *m1 =
+ (gsl_matrix_symbolic*)malloc(sizeof(gsl_matrix_symbolic));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_symbolic_matrix_dup(m);
+ if (y) {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix *m1 = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_double_matrix_dup(m);
+ if (y) {
+ gsl_matrix *m = (gsl_matrix*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_complex *m1 =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_complex_matrix_dup(m);
+ if (y) {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ if (n1*n2!=m->size1*m->size2) return 0;
+ if (m->tda == m->size2) {
+ // No copying necessary, just create a new view of this matrix.
+ gsl_matrix_int *m1 =
+ (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m1);
+ *m1 = *m;
+ m1->size1 = n1; m1->tda = m1->size2 = n2;
+ p = m1;
+ } else {
+ // Create a new matrix containing the same elements.
+ pure_expr *y = pure_int_matrix_dup(m);
+ if (y) {
+ gsl_matrix_int *m = (gsl_matrix_int*)y->data.mat.p;
+ m->size1 = n1; m->tda = m->size2 = n2;
+ }
+ return y;
+ }
+ break;
+ }
+#endif
+ default:
+ return 0;
+ }
+ pure_expr *y = new_expr();
+ y->tag = x->tag;
+ y->data.mat.p = p;
+ y->data.mat.refc = x->data.mat.refc;
+ (*y->data.mat.refc)++;
+ MEMDEBUG_NEW(y)
+ return y;
+}
+
+extern "C"
pure_expr *matrix_rows(pure_expr *xs)
{
size_t n;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-22 17:46:47 UTC (rev 827)
+++ pure/trunk/runtime.h 2008-09-22 19:20:24 UTC (rev 828)
@@ -692,6 +692,11 @@
pure_expr *matrix_slice(pure_expr *x, int32_t i1, int32_t j1,
int32_t i2, int32_t j2);
+/* Convert a matrix to a new matrix with same size but different dimensions.
+ Reuse the storage of the original matrix if its data is contiguous. */
+
+pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m);
+
/* Matrix construction. These work like the pure_matrix_rows/
pure_matrix_columns functions in the public API, but take their input from
a Pure list. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|