Mercurial > octave-nkf

changeset 7918:78eef61f75d5

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add matlab V4 sparse matrixload/save
author David Bateman <dbateman@free.fr>
date Wed, 09 Jul 2008 22:35:06 +0200
parents fa8f13a05687
children 9d080df0c843
files src/ChangeLog src/ls-mat4.cc
diffstat 2 files changed, 173 insertions(+), 40 deletions(-) [+]
line wrap: on
line diff
--- a/src/ChangeLog	Thu Jul 10 10:57:55 2008 -0400
+++ b/src/ChangeLog	Wed Jul 09 22:35:06 2008 +0200
@@ -1,3 +1,8 @@
+2008-07-10  David Bateman  <dbateman@free.fr>
+
+	* ls-mat4.cc (read_mat_binary_data, save_mat_binary_data): Add
+	loading and saving of sparse matrices.
+	
 2008-07-10  Michael Goffioul  <michael.goffioul@gmail.com>
 
 	* Makefile.in: Add OPENGL_LIBS to liboctinterp link command. Add
--- a/src/ls-mat4.cc	Thu Jul 10 10:57:55 2008 -0400
+++ b/src/ls-mat4.cc	Wed Jul 09 22:35:06 2008 +0200
@@ -61,6 +61,7 @@
 #include "variables.h"
 #include "version.h"
 #include "dMatrix.h"
+#include "dSparse.h"
 
 #include "ls-mat4.h"
 
@@ -301,12 +302,6 @@
       return retval;
     }
 
-  if (type != 0 && type != 1)
-    {
-      error ("load: can't read sparse matrices");
-      return retval;
-    }
-
   if (imag && type == 1)
     {
       error ("load: encountered complex matrix with string flag set!");
@@ -335,42 +330,101 @@
 	nc = tmp;
       }
 
-      re.resize (nr, nc);
-
-      read_mat_binary_data (is, re.fortran_vec (), prec, dlen, swap, flt_fmt);
+    if (type == 2)
+      {
+	if (nc == 4)
+	  {
+	    octave_idx_type nr_new, nc_new;
+	    Array<Complex> data (dim_vector (1, nr - 1));
+	    Array<octave_idx_type> c (dim_vector (1, nr - 1));
+	    Array<octave_idx_type> r (dim_vector (1, nr - 1));
+	    OCTAVE_LOCAL_BUFFER (double, dtmp, nr);
+	    OCTAVE_LOCAL_BUFFER (double, ctmp, nr);
 
-      if (! is || error_state)
-	{
-	  error ("load: reading matrix data for `%s'", name);
-	  goto data_read_error;
-	}
+	    read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
+	    for (octave_idx_type i = 0; i < nr - 1; i++)
+	      r.xelem(i) = dtmp[i] - 1;
+	    nr_new = dtmp[nr - 1];
+	    read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
+	    for (octave_idx_type i = 0; i < nr - 1; i++)
+	      c.xelem(i) = dtmp[i] - 1;
+	    nc_new = dtmp[nr - 1];
+	    read_mat_binary_data (is, dtmp, prec, nr - 1, swap, flt_fmt);
+	    read_mat_binary_data (is, ctmp, prec, 1, swap, flt_fmt);
+	    read_mat_binary_data (is, ctmp, prec, nr - 1, swap, flt_fmt);
 
-      if (imag)
-	{
-	  Matrix im (nr, nc);
+	    for (octave_idx_type i = 0; i < nr - 1; i++)
+	      data.xelem(i) = Complex (dtmp[i], ctmp[i]);
+	    read_mat_binary_data (is, ctmp, prec, 1, swap, flt_fmt);
+
+	    SparseComplexMatrix smc = SparseComplexMatrix (data, r, c, 
+							   nr_new, nc_new);
 
-	  read_mat_binary_data (is, im.fortran_vec (), prec, dlen, swap,
-				flt_fmt);
+	    tc = order ? smc.transpose () : smc;
+	  }
+	else
+	  {
+	    octave_idx_type nr_new, nc_new;
+	    Array<double> data (dim_vector (1, nr - 1));
+	    Array<octave_idx_type> c (dim_vector (1, nr - 1));
+	    Array<octave_idx_type> r (dim_vector (1, nr - 1));
+	    OCTAVE_LOCAL_BUFFER (double, dtmp, nr);
 
-	  if (! is || error_state)
-	    {
-	      error ("load: reading imaginary matrix data for `%s'", name);
-	      goto data_read_error;
-	    }
+	    read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
+	    for (octave_idx_type i = 0; i < nr - 1; i++)
+	      r.xelem(i) = dtmp[i] - 1;
+	    nr_new = dtmp[nr - 1];
+	    read_mat_binary_data (is, dtmp, prec, nr, swap, flt_fmt);
+	    for (octave_idx_type i = 0; i < nr - 1; i++)
+	      c.xelem(i) = dtmp[i] - 1;
+	    nc_new = dtmp[nr - 1];
+	    read_mat_binary_data (is, data.fortran_vec (), prec, nr - 1, swap, flt_fmt);
+	    read_mat_binary_data (is, dtmp, prec, 1, swap, flt_fmt);
+
+	    SparseMatrix sm = SparseMatrix (data, r, c, nr_new, nc_new);
 
-	  ComplexMatrix ctmp (nr, nc);
+	    tc = order ? sm.transpose () : sm;
+	  }
+      }
+    else
+      {
+	re.resize (nr, nc);
+
+	read_mat_binary_data (is, re.fortran_vec (), prec, dlen, swap, flt_fmt);
+
+	if (! is || error_state)
+	  {
+	    error ("load: reading matrix data for `%s'", name);
+	    goto data_read_error;
+	  }
 
-	  for (octave_idx_type j = 0; j < nc; j++)
-	    for (octave_idx_type i = 0; i < nr; i++)
-	      ctmp (i, j) = Complex (re (i, j), im (i, j));
+	if (imag)
+	  {
+	    Matrix im (nr, nc);
+
+	    read_mat_binary_data (is, im.fortran_vec (), prec, dlen, swap,
+				  flt_fmt);
+
+	    if (! is || error_state)
+	      {
+		error ("load: reading imaginary matrix data for `%s'", name);
+		goto data_read_error;
+	      }
 
-	  tc = order ? ctmp.transpose () : ctmp;
-	}
-      else
-	tc = order ? re.transpose () : re;
+	    ComplexMatrix ctmp (nr, nc);
+
+	    for (octave_idx_type j = 0; j < nc; j++)
+	      for (octave_idx_type i = 0; i < nr; i++)
+		ctmp (i, j) = Complex (re (i, j), im (i, j));
 
-      if (type == 1)
-	tc = tc.convert_to_str (false, true, '\'');
+	    tc = order ? ctmp.transpose () : ctmp;
+	  }
+	else
+	  tc = order ? re.transpose () : re;
+
+	if (type == 1)
+	  tc = tc.convert_to_str (false, true, '\'');
+      }
 
       return retval;
     }
@@ -389,7 +443,7 @@
 {
   int32_t mopt = 0;
 
-  mopt += tc.is_string () ? 1 : 0;
+  mopt += tc.is_sparse_type () ? 2 : tc.is_string () ? 1 : 0;
 
   oct_mach_info::float_format flt_fmt =
     oct_mach_info::native_float_format ();;
@@ -398,16 +452,34 @@
 
   os.write (reinterpret_cast<char *> (&mopt), 4);
   
+  octave_idx_type len;
   int32_t nr = tc.rows ();
-  os.write (reinterpret_cast<char *> (&nr), 4);
 
   int32_t nc = tc.columns ();
-  os.write (reinterpret_cast<char *> (&nc), 4);
+
+  if (tc.is_sparse_type ())
+    {
+      len = tc.nnz ();
+      uint32_t nnz = len + 1;
+      os.write (reinterpret_cast<char *> (&nnz), 4);
+
+      uint32_t iscmplx = tc.is_complex_type () ? 4 : 3;
+      os.write (reinterpret_cast<char *> (&iscmplx), 4);
 
-  octave_idx_type len = nr * nc;
+      uint32_t tmp = 0;
+      os.write (reinterpret_cast<char *> (&tmp), 4);
+    }
+  else
+    {
+      os.write (reinterpret_cast<char *> (&nr), 4);
+      os.write (reinterpret_cast<char *> (&nc), 4);
 
-  int32_t imag = tc.is_complex_type () ? 1 : 0;
-  os.write (reinterpret_cast<char *> (&imag), 4);
+      int32_t imag = tc.is_complex_type () ? 1 : 0;
+      os.write (reinterpret_cast<char *> (&imag), 4);
+
+      len = nr * nc;
+    }
+
 
   // LEN includes the terminating character, and the file is also
   // supposed to include it.
@@ -457,6 +529,62 @@
       double tmp = tc.double_value ();
       os.write (reinterpret_cast<char *> (&tmp), 8);
     }
+  else if (tc.is_sparse_type ())
+    {
+      double ds;
+      OCTAVE_LOCAL_BUFFER (double, dtmp, len);
+      if (tc.is_complex_matrix ())
+	{
+	  SparseComplexMatrix m = tc.sparse_complex_matrix_value ();
+
+	  for (octave_idx_type i = 0; i < len; i++)
+	    dtmp [i] = m.ridx(i) + 1;
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  ds = nr;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+
+	  octave_idx_type ii = 0;
+	  for (octave_idx_type j = 0; j < nc; j++)
+	    for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+	      dtmp[ii++] = j + 1;
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  ds = nc;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+
+	  for (octave_idx_type i = 0; i < len; i++)
+	    dtmp [i] = std::real (m.data(i));
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  ds = 0.;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+
+	  for (octave_idx_type i = 0; i < len; i++)
+	    dtmp [i] = std::imag (m.data(i));
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+	}
+      else
+	{
+	  SparseMatrix m = tc.sparse_matrix_value ();
+
+	  for (octave_idx_type i = 0; i < len; i++)
+	    dtmp [i] = m.ridx(i) + 1;
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  ds = nr;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+
+	  octave_idx_type ii = 0;
+	  for (octave_idx_type j = 0; j < nc; j++)
+	    for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++)
+	      dtmp[ii++] = j + 1;
+	  os.write (reinterpret_cast<const char *> (dtmp), 8 * len);
+	  ds = nc;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+
+	  os.write (reinterpret_cast<const char *> (m.data ()), 8 * len);
+	  ds = 0.;
+	  os.write (reinterpret_cast<const char *> (&ds), 8);
+	}
+    }
   else if (tc.is_real_matrix ())
     {
       Matrix m = tc.matrix_value ();