[pure-lang-svn] SF.net SVN: pure-lang:[829] pure/trunk
Status: Beta
Brought to you by:
agraef
|
From: <ag...@us...> - 2008-09-22 20:53:45
|
Revision: 829
http://pure-lang.svn.sourceforge.net/pure-lang/?rev=829&view=rev
Author: agraef
Date: 2008-09-22 20:53:37 +0000 (Mon, 22 Sep 2008)
Log Message:
-----------
Add (sub-,super-) diagonal operations.
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 19:20:24 UTC (rev 828)
+++ pure/trunk/lib/primitives.pure 2008-09-22 20:53:37 UTC (rev 829)
@@ -469,6 +469,16 @@
cols x::matrix = map (col x) (0..m-1) when _,m::int = dim x end;
+/* Extract (sub-,super-) diagonals from a matrix. Sub- and super-diagonals for
+ k=0 return the main diagonal. Indices for sub- and super-diagonals can also
+ be negative, in which case the corresponding super- or sub-diagonal is
+ returned instead. */
+
+private matrix_diag matrix_subdiag matrix_supdiag;
+extern expr* matrix_diag(expr* x) = diag;
+extern expr* matrix_subdiag(expr* x, int k) = subdiag;
+extern expr* matrix_supdiag(expr* x, int k) = supdiag;
+
/* 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. */
Modified: pure/trunk/runtime.cc
===================================================================
--- pure/trunk/runtime.cc 2008-09-22 19:20:24 UTC (rev 828)
+++ pure/trunk/runtime.cc 2008-09-22 20:53:37 UTC (rev 829)
@@ -4313,6 +4313,154 @@
}
extern "C"
+pure_expr *matrix_diag(pure_expr *x)
+{
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n = (n1<n2)?n1:n2;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
+pure_expr *matrix_subdiag(pure_expr *x, int32_t k)
+{
+ if (k<0) return matrix_supdiag(x, -k);
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*k0+2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k*m->tda;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
+pure_expr *matrix_supdiag(pure_expr *x, int32_t k)
+{
+ if (k<0) return matrix_subdiag(x, -k);
+ switch (x->tag) {
+ case EXPR::MATRIX: {
+ gsl_matrix_symbolic *m = (gsl_matrix_symbolic*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_symbolic *m1 = create_symbolic_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_symbolic_matrix(m1);
+ }
+#ifdef HAVE_GSL
+ case EXPR::DMATRIX: {
+ gsl_matrix *m = (gsl_matrix*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix *m1 = create_double_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_double_matrix(m1);
+ }
+ case EXPR::CMATRIX: {
+ gsl_matrix_complex *m = (gsl_matrix_complex*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_complex *m1 = create_complex_matrix(1, n);
+ for (size_t i = 0; i < n; i++) {
+ const size_t k = 2*k0+2*i*(m->tda+1);
+ m1->data[2*i] = m->data[k];
+ m1->data[2*i+1] = m->data[k+1];
+ }
+ return pure_complex_matrix(m1);
+ }
+ case EXPR::IMATRIX: {
+ gsl_matrix_int *m = (gsl_matrix_int*)x->data.mat.p;
+ const size_t n1 = m->size1, n2 = m->size2, n0 = (n1<n2)?n1:n2;
+ const size_t n = (n0>(size_t)k)?n0-k:0, k0 = k;
+ gsl_matrix_int *m1 = create_int_matrix(1, n);
+ for (size_t i = 0; i < n; i++)
+ m1->data[i] = m->data[k0+i*(m->tda+1)];
+ return pure_int_matrix(m1);
+ }
+#endif
+ default:
+ return 0;
+ }
+}
+
+extern "C"
pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m)
{
void *p = 0;
Modified: pure/trunk/runtime.h
===================================================================
--- pure/trunk/runtime.h 2008-09-22 19:20:24 UTC (rev 828)
+++ pure/trunk/runtime.h 2008-09-22 20:53:37 UTC (rev 829)
@@ -697,6 +697,12 @@
pure_expr *matrix_redim(pure_expr *x, int32_t n, int32_t m);
+/* Retrieve (sub-,super-)diagonals of a matrix, as a 1xn matrix. */
+
+pure_expr *matrix_diag(pure_expr *x);
+pure_expr *matrix_subdiag(pure_expr *x, int32_t k);
+pure_expr *matrix_supdiag(pure_expr *x, int32_t k);
+
/* 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.
|