[pure-lang-svn] SF.net SVN: pure-lang:[785] pure/trunk
Status: Beta
Brought to you by:
agraef
Menu
▾
▴
[pure-lang-svn] SF.net SVN: pure-lang:[785] pure/trunk
|
From: <ag...@us...> - 2008-09-18 13:19:21
|
Revision: 785
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=785&view=rev
Author: agraef
Date: 2008-09-18 20:19:31 +0000 (Thu, 18 Sep 2008)
Log Message:
-----------
Add some basic matrix operations (matrix size and dimensions, indexing, slicing).
Modified Paths:
--------------
pure/trunk/ChangeLog
pure/trunk/lib/primitives.pure
pure/trunk/runtime.cc
pure/trunk/runtime.h
Modified: pure/trunk/ChangeLog
===================================================================
--- pure/trunk/ChangeLog 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/ChangeLog 2008-09-18 20:19:31 UTC (rev 785)
@@ -1,17 +1,36 @@
2008-09-18 Albert Graef <Dr....@t-...>
+ * lib/primitives.pure: Add definitions of basic matrix operations.
+ Currently only size, dimensions and indexing are supported.
+
+ Synopsis:
+ - #x total number of elements
+ - dim x number of rows and columns (as a pair)
+ - x!i ith matrix element in row-major order
+ - x!(i,j) jth element in ith row of the matrix
+
* expr.cc, interpreter.cc, runtime.cc, printer.cc: Add basic GSL
matrix support. GSL double, complex and integer matrices can be
created with the new {x,y;u,v} syntax, which works more or less
like Octave/MATLAB matrices, but using curly braces instead of
- brackets.
+ brackets. Moreover, various basic operations to handle this kind
+ of objects (conversions, determining sizes, indexing, slicing)
+ have been added to the runtime, including some public API
+ operations to create and inspect matrix objects.
+ Note that the {...} syntax can be used only on the right-hand side
+ of a definition, matrix patterns are not supported right now. As a
+ remedy, there are three new type tags, matrix, cmatrix and
+ imatrix, which can be used in patterns to match double, complex
+ and integer matrices, respectively.
+
GSL matrices are always homogeneous, i.e., they only contain
values from one numeric type. Integer matrices can be created from
any combination of Pure machine int and bigint values (the latter
are converted to machine ints automatically). Matrices with at
least one double or complex element become double and complex
- matrices, respectively.
+ matrices, respectively, casting the other matrix elements as
+ needed.
Complex matrices can be created from either pairs of double or
integer values, such as {(1,2),(3,4)}, or from Pure complex
@@ -20,18 +39,10 @@
math.pure has been loaded in which case complex values are printed
in rectangular format x+:y. Also note that, for performance
reasons, the expression printer doesn't use the __show__ function
- to print matrix elements, but of course it is possible to override
- the default print representation of matrix values as a whole.
+ to print individual matrix elements, but of course it is possible
+ to override the default print representation of matrix values as a
+ whole.
- Note that the {...} syntax can be used only on the right-hand side
- of a definition, matrix patterns are not supported right now. As a
- remedy, there are three new type tags, matrix, cmatrix and
- imatrix, which can be used in patterns to match double, complex
- and integer matrices, respectively.
-
- Operations to handle Pure matrix expressions have been added to
- the public runtime API.
-
2008-09-15 Albert Graef <Dr....@t-...>
* configure.ac: Bump version number.
Modified: pure/trunk/lib/primitives.pure
===================================================================
--- pure/trunk/lib/primitives.pure 2008-09-18 17:35:01 UTC (rev 784)
+++ pure/trunk/lib/primitives.pure 2008-09-18 20:19:31 UTC (rev 785)
@@ -42,6 +42,10 @@
stringp x = case x of _::string = 1; _ = 0 end;
pointerp x = case x of _::pointer = 1; _ = 0 end;
+matrixp x = case x of _::matrix = 1; _ = 0 end;
+cmatrixp x = case x of _::cmatrix = 1; _ = 0 end;
+imatrixp x = case x of _::imatrix = 1; _ = 0 end;
+
/* Predicates to check for function objects, global (unbound) variables,
function applications, proper lists, list nodes and tuples. */
@@ -383,6 +387,31 @@
x::pointer==y::pointer = bigint x == bigint y;
x::pointer!=y::pointer = bigint x != bigint y;
+/* Basic matrix operations. */
+
+private matrix_size matrix_dim;
+extern int matrix_size(expr *x), expr* matrix_dim(expr *x);
+
+#x::matrix | #x::cmatrix | #x::imatrix = matrix_size x;
+dim x::matrix | dim x::cmatrix | dim x::imatrix = matrix_dim x;
+
+private matrix_elem_at;
+extern expr* matrix_elem_at(expr* x, int i);
+
+x::matrix!i::int | x::cmatrix!i::int | x::imatrix!i::int
+ = matrix_elem_at x i if i>=0 && i<#x;
+ = throw out_of_bounds otherwise;
+
+private matrix_elem;
+extern expr* matrix_elem(expr* x, int i, int j);
+
+x::matrix!(i::int,j::int) | x::cmatrix!(i::int,j::int) |
+x::imatrix!(i::int,j::int)
+ = matrix_elem x i j
+ if (i>=0 && i<n && j>=0 && j<m
+ when n,m = dim x end);
+ = throw out_of_bounds otherwise;
+
/* 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 17:35:01 UTC (rev 784)
+++ pure/trunk/runtime.cc 2008-09-18 20:19:31 UTC (rev 785)
@@ -3611,6 +3611,184 @@
return interp.errmsg.c_str();
}
+extern "C"
+uint32_t matrix_size(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return m->size1*m->size2;
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_dim(pure_expr *x)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return pure_tuplel(2, pure_int(m->size1), pure_int(m->size2));
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+static inline pure_expr *make_complex(double a, double b)
+{
+ interpreter& interp = *interpreter::g_interp;
+ symbol *rect = interp.symtab.complex_rect_sym();
+ if (rect)
+ return pure_appl(pure_symbol(rect->f), 2, pure_double(a), pure_double(b));
+ else
+ return pure_tuplel(2, pure_double(a), pure_double(b));
+}
+
+extern "C"
+pure_expr *matrix_elem_at(pure_expr *x, uint32_t i)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ return pure_double(m->data[i]);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ return make_complex(m->data[2*i], m->data[2*i+1]);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ return pure_int(m->data[i]);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_elem(pure_expr *x, uint32_t i, uint32_t j)
+{
+#ifdef HAVE_GSL
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ size_t k = i*m->tda+j;
+ return pure_double(m->data[k]);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ size_t k = 2*(i*m->tda+j);
+ return make_complex(m->data[k], m->data[k+1]);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ size_t k = i*m->tda+j;
+ return pure_int(m->data[k]);
+ }
+ default:
+ return 0;
+ }
+#else
+ return 0;
+#endif
+}
+
+extern "C"
+pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
+ uint32_t i2, uint32_t j2)
+{
+#ifdef HAVE_GSL
+ void *p = 0;
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_double_matrix(create_double_matrix(n1, n2));
+ gsl_matrix_view v = gsl_matrix_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix *m1 = (gsl_matrix*)malloc(sizeof(gsl_matrix));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_complex_matrix(create_complex_matrix(n1, n2));
+ gsl_matrix_complex_view v =
+ gsl_matrix_complex_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_complex *m1 =
+ (gsl_matrix_complex*)malloc(sizeof(gsl_matrix_complex));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ size_t n1 = (i2>=i1)?(i2+1-i1):0, n2 = (j2>=j1)?(j2+1-j1):0;
+ if (n1 == 0 || n2 == 0) // empty matrix
+ return pure_int_matrix(create_int_matrix(n1, n2));
+ gsl_matrix_int_view v = gsl_matrix_int_submatrix(m, i1, j1, n1, n2);
+ // take a copy of the view matrix
+ gsl_matrix_int *m1 = (gsl_matrix_int*)malloc(sizeof(gsl_matrix_int));
+ assert(m1 && v.matrix.data);
+ *m1 = v.matrix;
+ p = m1;
+ break;
+ }
+ default:
+ return 0;
+ }
+ // create a new expression for the slice, update the reference counter for
+ // the underlying GSL matrix
+ pure_expr *y = new_expr();
+ y->tag = EXPR::MATRIX;
+ y->data.mat.p = p;
+ y->data.mat.refc = x->data.mat.refc;
+ *y->data.mat.refc++;
+ MEMDEBUG_NEW(y)
+ return y;
+#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 17:35:01 UTC (rev 784)
+++ pure/trunk/runtime.h 2008-09-18 20:19:31 UTC (rev 785)
@@ -602,6 +602,30 @@
const char *lasterr();
+/* Basic matrix operations. These work with all supported GSL matrix types.
+ matrix_size determines the number of elements in a matrix, matrix_dim the
+ number of rows and columns, which are returned as a pair (n,m). */
+
+uint32_t matrix_size(pure_expr *x);
+pure_expr *matrix_dim(pure_expr *x);
+
+/* Matrix elements can be retrieved either by a single index (using row-major
+ order), or by row and column index. All indices are zero-based. Indices
+ aren't range-checked, if this is needed you have to do it beforehand using
+ matrix_size or matrix_dim above. */
+
+pure_expr *matrix_elem_at(pure_expr *x, uint32_t i);
+pure_expr *matrix_elem(pure_expr *x, uint32_t i, uint32_t j);
+
+/* The following operation retrieves a slice a.k.a. submatrix of a matrix and
+ returns it as a matrix object. The new matrix object shares the underlying
+ storage with the original matrix (i.e., matrix elements are *not* copied)
+ and so this is a comparatively cheap operation. Indices are zero-based and
+ not checked. */
+
+pure_expr *matrix_slice(pure_expr *x, uint32_t i1, uint32_t j1,
+ uint32_t i2, uint32_t j2);
+
/* 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.
|
