Mercurial > hg > octave-lojdl
changeset 7918:78eef61f75d5
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 +++ b/src/ChangeLog @@ -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 +++ b/src/ls-mat4.cc @@ -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 ();