[pure-lang-svn] SF.net SVN: pure-lang:[787] pure/trunk
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
[pure-lang-svn] SF.net SVN: pure-lang:[787] pure/trunk
|
From: <ag...@us...> - 2008-09-18 22:26:27
|
Revision: 787
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=787&view=rev
Author: agraef
Date: 2008-09-19 05:26:38 +0000 (Fri, 19 Sep 2008)
Log Message:
-----------
Add matrix transposition and conversion operations.
Modified Paths:
--------------
pure/trunk/lib/prelude.pure
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/lib/prelude.pure
===================================================================
--- pure/trunk/lib/prelude.pure 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/lib/prelude.pure 2008-09-19 05:26:38 UTC (rev 787)
@@ -58,6 +58,7 @@
infixl 6 + - or ; // addition, bitwise or
infixl 7 * / div mod and ; // multiplication, bitwise and
prefix 7 ~ ; // bitwise not
+postfix 7 ' ; // matrix transposition
infixr 8 ^ ; // exponentiation
prefix 8 # ; // size operator
infixl 9 ! !! ; // indexing, slicing
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/lib/primitives.pure 2008-09-19 05:26:38 UTC (rev 787)
@@ -428,6 +428,25 @@
when n::int,m::int = dim x end);
= throw out_of_bounds otherwise;
+private matrix_transpose;
+extern expr* matrix_transpose(expr *x);
+
+x::matrix' | x::cmatrix' | x::imatrix' = matrix_transpose x;
+
+private matrix_double matrix_complex matrix_int;
+extern expr* matrix_double(expr *x), expr* matrix_complex(expr *x),
+ expr* matrix_int(expr *x);
+
+dmatrix x::matrix | dmatrix x::imatrix = matrix_double x;
+imatrix x::matrix | imatrix x::imatrix = matrix_int x;
+cmatrix x::matrix | cmatrix x::cmatrix | cmatrix x::imatrix = matrix_complex x;
+
+private matrix_re matrix_im;
+extern expr* matrix_re(expr *x), expr* matrix_im(expr *x);
+
+re x::matrix | re x::cmatrix | re x::imatrix = matrix_re x;
+im x::matrix | im x::cmatrix | im x::imatrix = matrix_im x;
+
/* IEEE floating point infinities and NaNs. Place these after the definitions
of the built-in operators so that the double arithmetic works. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/runtime.cc 2008-09-19 05:26:38 UTC (rev 787)
@@ -3789,6 +3789,201 @@
#endif
}
+extern "C"
+pure_expr *matrix_transpose(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[j*m2->tda+i] = m1->data[i*m1->tda+j];
+ return pure_double_matrix(m2);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = 2*(j*m2->tda+i);
+ m2->data[l] = m1->data[k];
+ m2->data[l+1] = m1->data[k+1];
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(m, n);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[j*m2->tda+i] = m1->data[i*m1->tda+j];
+ return pure_int_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_double(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX:
+ return x;
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[i*m2->tda+j] = (double)m1->data[i*m1->tda+j];
+ return pure_double_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_complex(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m2->tda+j);
+ m2->data[k] = m1->data[i*m1->tda+j];
+ m2->data[k+1] = 0.0;
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_complex *m2 = create_complex_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m2->tda+j);
+ m2->data[k] = (double)m1->data[i*m1->tda+j];
+ m2->data[k+1] = 0.0;
+ }
+ return pure_complex_matrix(m2);
+ }
+ case EXPR::CMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_int(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++)
+ m2->data[i*m2->tda+j] = (int)m1->data[i*m1->tda+j];
+ return pure_int_matrix(m2);
+ }
+ case EXPR::IMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_re(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = i*m2->tda+j;
+ m2->data[l] = m1->data[k];
+ }
+ return pure_double_matrix(m2);
+ }
+ case EXPR::MATRIX:
+ case EXPR::IMATRIX:
+ return x;
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_im(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m1 = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ for (size_t i = 0; i < n; i++)
+ for (size_t j = 0; j < m; j++) {
+ size_t k = 2*(i*m1->tda+j), l = i*m2->tda+j;
+ m2->data[l] = m1->data[k+1];
+ }
+ return pure_double_matrix(m2);
+ }
+ case EXPR::MATRIX: {
+ gsl_matrix *m1 = (gsl_matrix*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix *m2 = create_double_matrix(n, m);
+ memset(m2->data, 0, n*m*sizeof(double));
+ return pure_double_matrix(m2);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m1 = (gsl_matrix_int*)x->data.mat.p;
+ size_t n = m1->size1, m = m1->size2;
+ gsl_matrix_int *m2 = create_int_matrix(n, m);
+ memset(m2->data, 0, n*m*sizeof(int));
+ return pure_int_matrix(m2);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
static uint32_t mpz_hash(const mpz_t z)
{
uint32_t h = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-18 20:35:54 UTC (rev 786)
+++ pure/trunk/runtime.h 2008-09-19 05:26:38 UTC (rev 787)
@@ -626,6 +626,27 @@
pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
uint32_t i2, uint32_t j2);
+/* Transpose a matrix. The resulting matrix has the rows of the original
+ matrix as its columns, and vice versa. */
+
+pure_expr *matrix_transpose(pure_expr *x);
+
+/* Convert an existing matrix to a double, complex or int matrix,
+ respectively. Any kind of matrix can be converted to a complex matrix, but
+ the input must be a double or integer matrix for the other conversions (see
+ matrix_re and matrix_im below to handle the complex->double case). */
+
+pure_expr *matrix_double(pure_expr *x);
+pure_expr *matrix_complex(pure_expr *x);
+pure_expr *matrix_int(pure_expr *x);
+
+/* Extract the real and imaginary parts of a matrix. If the input is a complex
+ matrix, the result is a double matrix. Otherwise the type of the result is
+ the same as that of the input matrix. */
+
+pure_expr *matrix_re(pure_expr *x);
+pure_expr *matrix_im(pure_expr *x);
+
/* Compute a 32 bit hash code of a Pure expression. This makes it possible to
use arbitary Pure values as keys in a hash table. */
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
