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

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

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.