[pure-lang-svn] SF.net SVN: pure-lang:[828] pure/trunk

[<ag...@us...>]

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.