Mercurial > hg > octave-lyh
changeset 5282:5bdc3f24cd5f
[project @ 2005-04-14 22:17:26 by dbateman]
author | dbateman |
---|---|
date | Thu, 14 Apr 2005 22:17:27 +0000 |
parents | f3266e7dbb99 |
children | bb224e6c26a7 |
files | PROJECTS doc/ChangeLog doc/interpreter/sparse.txi liboctave/ChangeLog liboctave/SparseCmplxLU.cc liboctave/SparseCmplxLU.h liboctave/SparseType.cc liboctave/SparseType.h liboctave/SparsedbleLU.cc liboctave/SparsedbleLU.h src/ChangeLog src/DLD-FUNCTIONS/luinc.cc src/DLD-FUNCTIONS/sparse.cc src/DLD-FUNCTIONS/splu.cc src/Makefile.in src/ov-bool-sparse.h src/ov-mapper.cc src/ov-re-sparse.cc src/ov-re-sparse.h |
diffstat | 19 files changed, 1020 insertions(+), 306 deletions(-) [+] |
line wrap: on
line diff
--- a/PROJECTS +++ b/PROJECTS @@ -116,22 +116,13 @@ * Once dmperm is implemented, use the technique to detect permuted triangular matrices. Test the permuted triangular matrix solver code - * Accelerate the copying of the data from a sparse matrix to a banded matrix - in the solvers, that takes a significant portion of the computation time - for narrow matrices. This is not obvious, due to the treatment of zero - elements in the matrix. Maybe current solution is optimal. - - * Perhaps split the overly long ::solve functions up, either by the type - of solver, or seperate factorization functions, so that they can be - reused in each of 4 different ::solve functions. - * Sparse inverse function, based on Andy's code from octave-forge. * Implement fourth argument to the sprand and sprandn that the leading brand implements. - * Mapper functions such as real, imag, abs, etc need to be treated, either - with a dispatch or by changing the mapper function code. + * Sparse logical indexing in idx_vector class so that something like + "a=sprandn(1e6,1e6,1e-6); a(a<1) = 0" won't cause a memory overflow. * Write the rest of the sparse docs
--- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,3 +1,7 @@ +2005-03-14 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Add luinc function. + 2005-03-09 John W. Eaton <jwe@octave.org> * Makefile.in (bin-dist): Delete target.
--- a/doc/interpreter/sparse.txi +++ b/doc/interpreter/sparse.txi @@ -872,7 +872,7 @@ @table @asis @item License -UMFPACK Version 4.3 (Jan. 16, 2004), Copyright @copyright{} 2004 by +UMFPACK Version 4.4 (Jan. 28, 2005), Copyright @copyright{} 2005 by Timothy A. Davis. All Rights Reserved. Your use or distribution of UMFPACK or any modified version of @@ -1002,6 +1002,8 @@ @emph{Not implemented} @item gmres @emph{Not implemented} +@item luinc +Produce the incomplete LU factorization of the sparse matrix A. @item lsqr @emph{Not implemented} @item minres @@ -1062,6 +1064,8 @@ sparse * issparse:: Return 1 if the value of the expression EXPR is a sparse matrix. +* luinc:: Produce the incomplete LU factorization of the sparse + A. * nnz:: returns number of non zero elements in SM See also: sparse * nonzeros:: Returns a vector of the non-zero values of the sparse matrix S @@ -1138,12 +1142,17 @@ @DOCSTRING(full) -@node issparse, nnz, full, Function Reference +@node issparse, luinc, full, Function Reference @subsubsection issparse @DOCSTRING(issparse) -@node nnz, nonzeros, issparse, Function Reference +@node luinc, nnz, issparse, Function Reference +@subsubsection luinc + +@DOCSTRING(luinc) + +@node nnz, nonzeros, luinc, Function Reference @subsubsection nnz @DOCSTRING(nnz)
--- a/liboctave/ChangeLog +++ b/liboctave/ChangeLog @@ -1,3 +1,13 @@ +2005-04-14 David Bateman <dbateman@free.fr> + + * SparseCmplxLU.cc: Add flags for incomplete factorization. + * SparsedbleLU.cc: Ditto. + * SparseCmplxLU.h: Definition. + * SparsedbleLU.h: ditto. + + * SparseType.cc (transpose): New function. + * SparseType.h (transpose): Definition. + 2005-04-11 John W. Eaton <jwe@octave.org> * lo-specfun.cc: Use F77_XFCN instead of F77_FUNC for calls to
--- a/liboctave/SparseCmplxLU.cc +++ b/liboctave/SparseCmplxLU.cc @@ -223,117 +223,113 @@ SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit, - double piv_thres, bool FixedQ) + double piv_thres, bool FixedQ, + double droptol, bool milu, bool udiag) { #ifdef HAVE_UMFPACK - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - // Setup the control parameters - Matrix Control (UMFPACK_CONTROL, 1); - double *control = Control.fortran_vec (); - umfpack_zi_defaults (control); - - double tmp = Voctave_sparse_controls.get_key ("spumoni"); - if (!xisnan (tmp)) - Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) - { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; - } - else - { - tmp = Voctave_sparse_controls.get_key ("piv_tol"); - if (!xisnan (tmp)) - { - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; - Control (UMFPACK_PIVOT_TOLERANCE) = tmp; - } - } - - // Set whether we are allowed to modify Q or not - if (FixedQ) - Control (UMFPACK_FIXQ) = 1.0; + if (milu) + (*current_liboctave_error_handler) + ("Modified incomplete LU not implemented"); else { - tmp = Voctave_sparse_controls.get_key ("autoamd"); - if (!xisnan (tmp)) - Control (UMFPACK_FIXQ) = tmp; - } + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); - // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - - umfpack_zi_report_control (control); + // Setup the control parameters + Matrix Control (UMFPACK_CONTROL, 1); + double *control = Control.fortran_vec (); + umfpack_zi_defaults (control); - const octave_idx_type *Ap = a.cidx (); - const octave_idx_type *Ai = a.ridx (); - const Complex *Ax = a.data (); - - umfpack_zi_report_matrix (nr, nc, Ap, Ai, X_CAST (const double *, Ax), NULL, - 1, control); + double tmp = Voctave_sparse_controls.get_key ("spumoni"); + if (!xisnan (tmp)) + Control (UMFPACK_PRL) = tmp; + if (piv_thres >= 0.) + { + piv_thres = (piv_thres > 1. ? 1. : piv_thres); + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; + Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + } + else + { + tmp = Voctave_sparse_controls.get_key ("piv_tol"); + if (!xisnan (tmp)) + { + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + } + } - void *Symbolic; - Matrix Info (1, UMFPACK_INFO); - double *info = Info.fortran_vec (); - int status; - - // Null loop so that qinit is imediately deallocated when not needed - do { - OCTAVE_LOCAL_BUFFER (int, qinit, nc); + if (droptol >= 0.) + Control (UMFPACK_DROPTOL) = droptol; - for (int i = 0; i < nc; i++) - qinit [i] = static_cast<int> (Qinit (i)); + // Set whether we are allowed to modify Q or not + if (FixedQ) + Control (UMFPACK_FIXQ) = 1.0; + else + { + tmp = Voctave_sparse_controls.get_key ("autoamd"); + if (!xisnan (tmp)) + Control (UMFPACK_FIXQ) = tmp; + } - status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, X_CAST (const double *, Ax), - NULL, qinit, &Symbolic, control, info); - } while (0); + // Turn-off UMFPACK scaling for LU + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - if (status < 0) - { - (*current_liboctave_error_handler) - ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); + umfpack_zi_report_control (control); + + const octave_idx_type *Ap = a.cidx (); + const octave_idx_type *Ai = a.ridx (); + const Complex *Ax = a.data (); - umfpack_zi_report_status (control, status); - umfpack_zi_report_info (control, info); + umfpack_zi_report_matrix (nr, nc, Ap, Ai, + X_CAST (const double *, Ax), NULL, + 1, control); + + void *Symbolic; + Matrix Info (1, UMFPACK_INFO); + double *info = Info.fortran_vec (); + int status; - umfpack_zi_free_symbolic (&Symbolic) ; - } - else - { - umfpack_zi_report_symbolic (Symbolic, control); + // Null loop so that qinit is imediately deallocated when not + // needed + do { + OCTAVE_LOCAL_BUFFER (int, qinit, nc); - void *Numeric; - status = umfpack_zi_numeric (Ap, Ai, X_CAST (const double *, Ax), NULL, - Symbolic, &Numeric, control, info) ; - umfpack_zi_free_symbolic (&Symbolic) ; + for (int i = 0; i < nc; i++) + qinit [i] = static_cast<int> (Qinit (i)); - cond = Info (UMFPACK_RCOND); + status = umfpack_zi_qsymbolic (nr, nc, Ap, Ai, + X_CAST (const double *, Ax), + NULL, qinit, &Symbolic, control, + info); + } while (0); if (status < 0) { (*current_liboctave_error_handler) - ("SparseComplexLU::SparseComplexLU numeric factorization failed"); + ("SparseComplexLU::SparseComplexLU symbolic factorization failed"); umfpack_zi_report_status (control, status); umfpack_zi_report_info (control, info); - umfpack_zi_free_numeric (&Numeric); + umfpack_zi_free_symbolic (&Symbolic) ; } else { - umfpack_zi_report_numeric (Numeric, control); + umfpack_zi_report_symbolic (Symbolic, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2, - &ignore3, Numeric) ; - + void *Numeric; + status = umfpack_zi_numeric (Ap, Ai, + X_CAST (const double *, Ax), NULL, + Symbolic, &Numeric, control, info) ; + umfpack_zi_free_symbolic (&Symbolic) ; + + cond = Info (UMFPACK_RCOND); + if (status < 0) { (*current_liboctave_error_handler) - ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); + ("SparseComplexLU::SparseComplexLU numeric factorization failed"); umfpack_zi_report_status (control, status); umfpack_zi_report_info (control, info); @@ -342,71 +338,102 @@ } else { - int n_inner = (nr < nc ? nr : nc); - - if (lnz < 1) - Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (1)); - else - Lfact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); - - octave_idx_type *Ltp = Lfact.cidx (); - octave_idx_type *Ltj = Lfact.ridx (); - Complex *Ltx = Lfact.data (); - - if (unz < 1) - Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (1)); - else - Ufact = SparseComplexMatrix (static_cast<octave_idx_type> (n_inner), nc, unz); + umfpack_zi_report_numeric (Numeric, control); - octave_idx_type *Up = Ufact.cidx (); - octave_idx_type *Uj = Ufact.ridx (); - Complex *Ux = Ufact.data (); - - P.resize (nr); - int *p = P.fortran_vec (); - - Q.resize (nc); - int *q = Q.fortran_vec (); - - int do_recip; - status = umfpack_zi_get_numeric (Ltp, Ltj, X_CAST (double *, Ltx), - NULL, Up, Uj, - X_CAST (double *, Ux), NULL, p, - q, NULL, NULL, &do_recip, - NULL, Numeric) ; - - umfpack_zi_free_numeric (&Numeric) ; - - if (status < 0 || do_recip) + int lnz, unz, ignore1, ignore2, ignore3; + status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, + &ignore2, &ignore3, Numeric); + + if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); umfpack_zi_report_status (control, status); + umfpack_zi_report_info (control, info); + + umfpack_zi_free_numeric (&Numeric); } else { - Lfact = Lfact.transpose (); + int n_inner = (nr < nc ? nr : nc); + + if (lnz < 1) + Lfact = SparseComplexMatrix + (static_cast<octave_idx_type> (n_inner), nr, + static_cast<octave_idx_type> (1)); + else + Lfact = SparseComplexMatrix + (static_cast<octave_idx_type> (n_inner), nr, + static_cast<octave_idx_type> (lnz)); + + octave_idx_type *Ltp = Lfact.cidx (); + octave_idx_type *Ltj = Lfact.ridx (); + Complex *Ltx = Lfact.data (); - umfpack_zi_report_matrix (nr, n_inner, Lfact.cidx (), - Lfact.ridx (), - X_CAST (double *, Lfact.data()), - NULL, 1, control); + if (unz < 1) + Ufact = SparseComplexMatrix + (static_cast<octave_idx_type> (n_inner), nc, + static_cast<octave_idx_type> (1)); + else + Ufact = SparseComplexMatrix + (static_cast<octave_idx_type> (n_inner), nc, unz); + + octave_idx_type *Up = Ufact.cidx (); + octave_idx_type *Uj = Ufact.ridx (); + Complex *Ux = Ufact.data (); + + P.resize (nr); + int *p = P.fortran_vec (); + + Q.resize (nc); + int *q = Q.fortran_vec (); - umfpack_zi_report_matrix (n_inner, nc, Ufact.cidx (), - Ufact.ridx (), - X_CAST (double *, Ufact.data()), - NULL, 1, control); - umfpack_zi_report_perm (nr, p, control); - umfpack_zi_report_perm (nc, q, control); + int do_recip; + status = + umfpack_zi_get_numeric (Ltp, Ltj, + X_CAST (double *, Ltx), + NULL, Up, Uj, + X_CAST (double *, Ux), + NULL, p, q, NULL, NULL, + &do_recip, NULL, Numeric) ; + + umfpack_zi_free_numeric (&Numeric) ; + + if (status < 0 || do_recip) + { + (*current_liboctave_error_handler) + ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); + + umfpack_zi_report_status (control, status); + } + else + { + Lfact = Lfact.transpose (); + + umfpack_zi_report_matrix (nr, n_inner, + Lfact.cidx (), + Lfact.ridx (), + X_CAST (double *, Lfact.data()), + NULL, 1, control); + + umfpack_zi_report_matrix (n_inner, nc, + Ufact.cidx (), + Ufact.ridx (), + X_CAST (double *, Ufact.data()), + NULL, 1, control); + umfpack_zi_report_perm (nr, p, control); + umfpack_zi_report_perm (nc, q, control); + } + + umfpack_zi_report_info (control, info); } - - umfpack_zi_report_info (control, info); } } + + if (udiag) + (*current_liboctave_error_handler) + ("Option udiag of incomplete LU not implemented"); } #else (*current_liboctave_error_handler) ("UMFPACK not installed");
--- a/liboctave/SparseCmplxLU.h +++ b/liboctave/SparseCmplxLU.h @@ -38,7 +38,9 @@ SparseComplexLU (const SparseComplexMatrix& a, double piv_thres = -1); SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit, - double piv_thres = -1, bool FixedQ = false); + double piv_thres = -1, bool FixedQ = false, + double droptol = -1., bool milu = false, + bool udiag = false); SparseComplexLU (const SparseComplexLU& a) : sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double> (a) { }
--- a/liboctave/SparseType.cc +++ b/liboctave/SparseType.cc @@ -688,6 +688,37 @@ typ = SparseType::Lower; } +SparseType +SparseType::transpose (void) const +{ + SparseType retval (*this); + if (typ == SparseType::Upper) + retval.typ = Lower; + else if (typ == SparseType::Permuted_Upper) + { + int *tmp = retval.row_perm; + retval.row_perm = retval.col_perm; + retval.col_perm = tmp; + retval.typ = Lower; + } + else if (typ == SparseType::Lower) + retval.typ = Upper; + else if (typ == SparseType::Permuted_Upper) + { + int *tmp = retval.row_perm; + retval.row_perm = retval.col_perm; + retval.col_perm = tmp; + retval.typ = Upper; + } + else if (typ == SparseType::Banded) + { + retval.upper_band = lower_band; + retval.lower_band = upper_band; + } + + return retval; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/SparseType.h +++ b/liboctave/SparseType.h @@ -132,6 +132,8 @@ void mark_as_unpermuted (void); + SparseType transpose (void) const; + private: void type (int new_typ) { typ = static_cast<matrix_type>(new_typ); }
--- a/liboctave/SparsedbleLU.cc +++ b/liboctave/SparsedbleLU.cc @@ -217,116 +217,108 @@ } SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, - double piv_thres, bool FixedQ) + double piv_thres, bool FixedQ, double droptol, + bool milu, bool udiag) { #ifdef HAVE_UMFPACK - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.cols (); - - // Setup the control parameters - Matrix Control (UMFPACK_CONTROL, 1); - double *control = Control.fortran_vec (); - umfpack_di_defaults (control); - - double tmp = Voctave_sparse_controls.get_key ("spumoni"); - if (!xisnan (tmp)) - Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) - { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; - } - else - { - tmp = Voctave_sparse_controls.get_key ("piv_tol"); - if (!xisnan (tmp)) - { - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; - Control (UMFPACK_PIVOT_TOLERANCE) = tmp; - } - } - - // Set whether we are allowed to modify Q or not - if (FixedQ) - Control (UMFPACK_FIXQ) = 1.0; + if (milu) + (*current_liboctave_error_handler) + ("Modified incomplete LU not implemented"); else { - tmp = Voctave_sparse_controls.get_key ("autoamd"); - if (!xisnan (tmp)) - Control (UMFPACK_FIXQ) = tmp; - } + octave_idx_type nr = a.rows (); + octave_idx_type nc = a.cols (); - // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - - umfpack_di_report_control (control); + // Setup the control parameters + Matrix Control (UMFPACK_CONTROL, 1); + double *control = Control.fortran_vec (); + umfpack_di_defaults (control); - const octave_idx_type *Ap = a.cidx (); - const octave_idx_type *Ai = a.ridx (); - const double *Ax = a.data (); - - umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + double tmp = Voctave_sparse_controls.get_key ("spumoni"); + if (!xisnan (tmp)) + Control (UMFPACK_PRL) = tmp; + if (piv_thres >= 0.) + { + piv_thres = (piv_thres > 1. ? 1. : piv_thres); + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; + Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + } + else + { + tmp = Voctave_sparse_controls.get_key ("piv_tol"); + if (!xisnan (tmp)) + { + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + } + } - void *Symbolic; - Matrix Info (1, UMFPACK_INFO); - double *info = Info.fortran_vec (); - int status; + if (droptol >= 0.) + Control (UMFPACK_DROPTOL) = droptol; - // Null loop so that qinit is imediately deallocated when not needed - do { - OCTAVE_LOCAL_BUFFER (int, qinit, nc); - for (int i = 0; i < nc; i++) - qinit [i] = static_cast<int> (Qinit (i)); + // Set whether we are allowed to modify Q or not + if (FixedQ) + Control (UMFPACK_FIXQ) = 1.0; + else + { + tmp = Voctave_sparse_controls.get_key ("autoamd"); + if (!xisnan (tmp)) + Control (UMFPACK_FIXQ) = tmp; + } - status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit, - &Symbolic, control, info); - } while (0); + // Turn-off UMFPACK scaling for LU + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; - if (status < 0) - { - (*current_liboctave_error_handler) - ("SparseLU::SparseLU symbolic factorization failed"); + umfpack_di_report_control (control); - umfpack_di_report_status (control, status); - umfpack_di_report_info (control, info); + const octave_idx_type *Ap = a.cidx (); + const octave_idx_type *Ai = a.ridx (); + const double *Ax = a.data (); + + umfpack_di_report_matrix (nr, nc, Ap, Ai, Ax, 1, control); + + void *Symbolic; + Matrix Info (1, UMFPACK_INFO); + double *info = Info.fortran_vec (); + int status; - umfpack_di_free_symbolic (&Symbolic) ; - } - else - { - umfpack_di_report_symbolic (Symbolic, control); + // Null loop so that qinit is imediately deallocated when not needed + do { + OCTAVE_LOCAL_BUFFER (int, qinit, nc); - void *Numeric; - status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, - control, info) ; - umfpack_di_free_symbolic (&Symbolic) ; + for (int i = 0; i < nc; i++) + qinit [i] = static_cast<int> (Qinit (i)); - cond = Info (UMFPACK_RCOND); + status = umfpack_di_qsymbolic (nr, nc, Ap, Ai, Ax, qinit, + &Symbolic, control, info); + } while (0); if (status < 0) { (*current_liboctave_error_handler) - ("SparseLU::SparseLU numeric factorization failed"); + ("SparseLU::SparseLU symbolic factorization failed"); umfpack_di_report_status (control, status); umfpack_di_report_info (control, info); - umfpack_di_free_numeric (&Numeric); + umfpack_di_free_symbolic (&Symbolic) ; } else { - umfpack_di_report_numeric (Numeric, control); + umfpack_di_report_symbolic (Symbolic, control); - int lnz, unz, ignore1, ignore2, ignore3; - status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, - &ignore3, Numeric) ; - + void *Numeric; + status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric, + control, info) ; + umfpack_di_free_symbolic (&Symbolic) ; + + cond = Info (UMFPACK_RCOND); + if (status < 0) { (*current_liboctave_error_handler) - ("SparseLU::SparseLU extracting LU factors failed"); + ("SparseLU::SparseLU numeric factorization failed"); umfpack_di_report_status (control, status); umfpack_di_report_info (control, info); @@ -335,67 +327,100 @@ } else { - int n_inner = (nr < nc ? nr : nc); - - if (lnz < 1) - Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (1)); - else - Lfact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nr, - static_cast<octave_idx_type> (lnz)); - - octave_idx_type *Ltp = Lfact.cidx (); - octave_idx_type *Ltj = Lfact.ridx (); - double *Ltx = Lfact.data (); - - if (unz < 1) - Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (1)); - else - Ufact = SparseMatrix (static_cast<octave_idx_type> (n_inner), nc, - static_cast<octave_idx_type> (unz)); + umfpack_di_report_numeric (Numeric, control); - octave_idx_type *Up = Ufact.cidx (); - octave_idx_type *Uj = Ufact.ridx (); - double *Ux = Ufact.data (); - - P.resize (nr); - int *p = P.fortran_vec (); - - Q.resize (nc); - int *q = Q.fortran_vec (); - - int do_recip; - status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj, - Ux, p, q, (double *) NULL, - &do_recip, (double *) NULL, - Numeric) ; - - umfpack_di_free_numeric (&Numeric) ; - - if (status < 0 || do_recip) + int lnz, unz, ignore1, ignore2, ignore3; + status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2, + &ignore3, Numeric) ; + + if (status < 0) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); umfpack_di_report_status (control, status); + umfpack_di_report_info (control, info); + + umfpack_di_free_numeric (&Numeric); } else { - Lfact = Lfact.transpose (); - umfpack_di_report_matrix (nr, n_inner, Lfact.cidx (), - Lfact.ridx (), Lfact.data (), - 1, control); - umfpack_di_report_matrix (n_inner, nc, Ufact.cidx (), - Ufact.ridx (), Ufact.data (), - 1, control); - umfpack_di_report_perm (nr, p, control); - umfpack_di_report_perm (nc, q, control); + int n_inner = (nr < nc ? nr : nc); + + if (lnz < 1) + Lfact = SparseMatrix + (static_cast<octave_idx_type> (n_inner), nr, + static_cast<octave_idx_type> (1)); + else + Lfact = SparseMatrix + (static_cast<octave_idx_type> (n_inner), nr, + static_cast<octave_idx_type> (lnz)); + + octave_idx_type *Ltp = Lfact.cidx (); + octave_idx_type *Ltj = Lfact.ridx (); + double *Ltx = Lfact.data (); + + if (unz < 1) + Ufact = SparseMatrix + (static_cast<octave_idx_type> (n_inner), nc, + static_cast<octave_idx_type> (1)); + else + Ufact = SparseMatrix + (static_cast<octave_idx_type> (n_inner), nc, + static_cast<octave_idx_type> (unz)); + + octave_idx_type *Up = Ufact.cidx (); + octave_idx_type *Uj = Ufact.ridx (); + double *Ux = Ufact.data (); + + P.resize (nr); + int *p = P.fortran_vec (); + + Q.resize (nc); + int *q = Q.fortran_vec (); + + int do_recip; + status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Uj, + Ux, p, q, + (double *) NULL, + &do_recip, + (double *) NULL, + Numeric) ; + + umfpack_di_free_numeric (&Numeric) ; + + if (status < 0 || do_recip) + { + (*current_liboctave_error_handler) + ("SparseLU::SparseLU extracting LU factors failed"); + + umfpack_di_report_status (control, status); + } + else + { + Lfact = Lfact.transpose (); + umfpack_di_report_matrix (nr, n_inner, + Lfact.cidx (), + Lfact.ridx (), + Lfact.data (), + 1, control); + umfpack_di_report_matrix (n_inner, nc, + Ufact.cidx (), + Ufact.ridx (), + Ufact.data (), + 1, control); + umfpack_di_report_perm (nr, p, control); + umfpack_di_report_perm (nc, q, control); + } + + umfpack_di_report_info (control, info); } - - umfpack_di_report_info (control, info); } } + + if (udiag) + (*current_liboctave_error_handler) + ("Option udiag of incomplete LU not implemented"); } #else (*current_liboctave_error_handler) ("UMFPACK not installed");
--- a/liboctave/SparsedbleLU.h +++ b/liboctave/SparsedbleLU.h @@ -36,7 +36,8 @@ SparseLU (const SparseMatrix& a, double piv_thres = -1.0); SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, - double piv_thres = -1.0, bool FixedQ = false); + double piv_thres = -1.0, bool FixedQ = false, + double droptol = -1., bool milu = false, bool udiag = false); SparseLU (const SparseLU& a) : sparse_base_lu <SparseMatrix, double, SparseMatrix, double> (a) { }
--- a/src/ChangeLog +++ b/src/ChangeLog @@ -1,3 +1,27 @@ +2005-04-14 David Bateman <dbateman@free.fr> + + * Makefile.in: Add luinc.cc to DLD_XSRC. + * DLD-FUNCTIONS/luinc.cc: New file for incomplete LU factorization. + + * ov-bool-sparse.h (index_vector): New function. + * ov-re-sparse.cc (index_vector): Ditto. + * ov-re-sparse.h (index_vector): Definition. + + * ov.mapper (any_element_less_than, any_element_greater_than): + New versions for SparseMatrix + (SPARSE_MAPPER_LOOP_2, SPARSE_MAPPER_LOOP_1, SPARSE_MAPPER_LOOP): + New macros. + (apply): Add special cases for sparse arguments to the mapper + functions + + * ov-re-sparse.cc (streamoff_array_value): Use octave_idx_type. + (convert_to_str_internal): New function. + * ov-re-sparse.h (convert_to_str_internal): Definition. + + * DLD-FUNCTIONS/sparse.cc (Fsparse): More care for nargin=2 case. + + * DLD-FUNCTIONS/splu.cc (Fsplu): Use octave_idx_type. + 2005-04-14 John W. Eaton <jwe@octave.org> * strfns.cc (Fchar): If arg is a dq string, return a dq string.
new file mode 100644 --- /dev/null +++ b/src/DLD-FUNCTIONS/luinc.cc @@ -0,0 +1,364 @@ +/* + +Copyright (C) 2005 David Bateman + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the Free +Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" +#include "oct-map.h" + +#include "SparseCmplxLU.h" +#include "SparsedbleLU.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" + +DEFUN_DLD (luinc, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, '0')\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, @var{droptol})\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} luinc (@var{a}, @var{opts})\n\ +@cindex LU decomposition\n\ +Produce the incomplete LU factorization of the sparse matrix @var{a}.\n\ +Two types of incomplete factorization are possible, and the type\n\ +is determined by the second argument to @dfn{luinc}.\n\ +\n\ +Called with a second argument of '0', the zero-level incomplete\n\ +LU factorization is produced. This creates a factorization of @var{a}\n\ +where the position of the non-zero arguments correspond to the same\n\ +positions as in the matrix @var{a}.\n\ +\n\ +Alternatively, the fill-in of the incomplete LU factorization can\n\ +be controlled through the variable @var{droptol} or the structure\n\ +@var{opts}. The UMFPACK multifrontal factorization code by Tim A.\n\ +Davis is used for the incomplete LU factorication, (availability\n\ +@url{http://www.cise.ufl.edu/research/sparse/umfpack/})\n\ +\n\ +@var{droptol} determines the values below which the values in the LU\n\ +factorization are dropped and replaced by zero. It must be a positive\n\ +scalar, and any values in the factorization whose absolute value are\n\ +less than this value are dropped, expect if leaving them increase the\n\ +sparsity of the matrix. Setting @var{droptol} to zero results in a\n\ +complete LU factorization which is the default.\n\ +\n\ +@var{opts} is a structure containing one or more of the fields\n\ +\n\ +@table @code\n\ +@item droptol\n\ +The drop tolerance as above. If @var{opts} only contains @code{droptol}\n\ +then this is equivalent to using the variable @var{droptol}.\n\ +\n\ +@item milu\n\ +A logical variable flagging whether to use the modified incomplete LU\n\ +factorization. In the case that @code{milu} is true, the dropped values\n\ +are subtract from the diagonal of the matrix U of the factorization.\n\ +The default is @code{false}.\n\ +\n\ +@item udiag\n\ +A logical variable that flags whether zero elements on the diagonal of U\n\ +should be replaced with @var{droptol} to attempt to avoid singular\n\ +factors. The default is @code{false}.\n\ +\n\ +@item thresh\n\ +Defines the pivot threshold in the interval [0,1]. Values outside that\n\ +range are ignored.\n\ +@end table\n\ +\n\ +All other fields in @var{opts} are ignored. The outputs from @dfn{luinc}\n\ +are the same as for @dfn{lu}.\n\ +@end deftypefn\n\ +@seealso{sparse, lu, cholinc}") +{ + int nargin = args.length (); + octave_value_list retval; + + if (nargin == 0) + print_usage ("luinc"); + else if (nargin != 2) + error ("luinc: incorrect number of arguments"); + else + { + bool zero_level = false; + bool milu = false; + bool udiag = false; + bool thresh = -1; + double droptol = -1.; + + if (args(1).is_string ()) + { + if (args(1).string_value () == "0") + zero_level = true; + else + error ("luinc: unrecognized string argument"); + } + else if (args(1).is_map ()) + { + Octave_map map = args(1).map_value (); + + if (map.contains ("droptol")) + droptol = map.contents ("droptol")(0).double_value (); + + if (map.contains ("milu")) + { + double tmp = map.contents ("milu")(0).double_value (); + + milu = (tmp == 0. ? false : true); + } + + if (map.contains ("udiag")) + { + double tmp = map.contents ("udiag")(0).double_value (); + + udiag = (tmp == 0. ? false : true); + } + + if (map.contains ("thresh")) + thresh = map.contents ("thresh")(0).double_value (); + } + else + droptol = args(1).double_value (); + + // XXX FIXME XXX Add code for zero-level factorization + if (zero_level) + error ("luinc: zero-level factorization not implemented"); + + if (!error_state) + { + if (args(0).class_name () == "sparse") + { + SparseMatrix sm = args(0).sparse_matrix_value (); + octave_idx_type sm_nc = sm.cols (); + ColumnVector Qinit (sm_nc); + + for (octave_idx_type i = 0; i < sm_nc; i++) + Qinit (i) = i; + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + case 2: + { + SparseLU fact (sm, Qinit, thresh, true, droptol, + milu, udiag); + + SparseMatrix P = fact.Pr (); + SparseMatrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + { + SparseLU fact (sm, Qinit, thresh, true, droptol, + milu, udiag); + + retval(2) = fact.Pr (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + + case 4: + default: + { + SparseLU fact (sm, Qinit, thresh, false, droptol, + milu, udiag); + + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + } + } + } + else if (args(0).class_name () == "sparse complex") + { + SparseComplexMatrix sm = + args(0).sparse_complex_matrix_value (); + octave_idx_type sm_nc = sm.cols (); + ColumnVector Qinit (sm_nc); + + for (octave_idx_type i = 0; i < sm_nc; i++) + Qinit (i) = i; + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + case 2: + { + SparseComplexLU fact (sm, Qinit, thresh, true, + droptol, milu, udiag); + + SparseMatrix P = fact.Pr (); + SparseComplexMatrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + { + SparseComplexLU fact (sm, Qinit, thresh, true, + droptol, milu, udiag); + + retval(2) = fact.Pr (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + + case 4: + default: + { + SparseComplexLU fact (sm, Qinit, thresh, false, + droptol, milu, udiag); + + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + } + } + } + else + error ("luinc: first argument must be sparse"); + } + } + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ + +/* + LUINC Sparse Incomplete LU factorization. + LUINC produces two different kinds of incomplete LU factorizations -- the + drop tolerance and the 0 level of fill-in factorizations. These factors + may be useful as preconditioners for a system of linear equations being + solved by an iterative method such as BICG (BiConjugate Gradients). + + LUINC(X,DROPTOL) performs the incomplete LU factorization of + X with drop tolerance DROPTOL. + + LUINC(X,OPTS) allows additional options to the incomplete LU + factorization. OPTS is a structure with up to four fields: + droptol --- the drop tolerance of incomplete LU + milu --- modified incomplete LU + udiag --- replace zeros on the diagonal of U + thresh --- the pivot threshold (see also LU) + + Only the fields of interest need to be set. + + droptol is a non-negative scalar used as the drop + tolerance for the incomplete LU factorization. This factorization + is computed in the same (column-oriented) manner as the + LU factorization except after each column of L and U has + been calculated, all entries in that column which are smaller + in magnitude than the local drop tolerance, which is + droptol * NORM of the column of X, are "dropped" from L or U. + The only exception to this dropping rule is the diagonal of the + upper triangular factor U which remains even if it is too small. + Note that entries of the lower triangular factor L are tested + before being scaled by the pivot. Setting droptol = 0 + produces the complete LU factorization, which is the default. + + milu stands for modified incomplete LU factorization. Its + value is either 0 (unmodified, the default) or 1 (modified). + Instead of discarding those entries from the newly-formed + column of the factors, they are subtracted from the diagonal + of the upper triangular factor, U. + + udiag is either 0 or 1. If it is 1, any zero diagonal entries + of the upper triangular factor U are replaced by the local drop + tolerance in an attempt to avoid a singular factor. The default + is 0. + + thresh is a pivot threshold in [0,1]. Pivoting occurs + when the diagonal entry in a column has magnitude less + than thresh times the magnitude of any sub-diagonal entry in + that column. thresh = 0 forces diagonal pivoting. thresh = 1 is + the default. + + Example: + + load west0479 + A = west0479; + nnz(A) + nnz(lu(A)) + nnz(luinc(A,1e-6)) + + This shows that A has 1887 nonzeros, its complete LU factorization + has 16777 nonzeros, and its incomplete LU factorization with a + drop tolerance of 1e-6 has 10311 nonzeros. + + + [L,U,P] = LUINC(X,'0') produces the incomplete LU factors of a sparse + matrix with 0 level of fill-in (i.e. no fill-in). L is unit lower + trianglar, U is upper triangular and P is a permutation matrix. U has the + same sparsity pattern as triu(P*X). L has the same sparsity pattern as + tril(P*X), except for 1's on the diagonal of L where P*X may be zero. Both + L and U may have a zero because of cancellation where P*X is nonzero. L*U + differs from P*X only outside of the sparsity pattern of P*X. + + [L,U] = LUINC(X,'0') produces upper triangular U and L is a permutation of + unit lower triangular matrix. Thus no comparison can be made between the + sparsity patterns of L,U and X, although nnz(L) + nnz(U) = nnz(X) + n. L*U + differs from X only outside of its sparsity pattern. + + LU = LUINC(X,'0') returns "L+U-I", where L is unit lower triangular, U is + upper triangular and the permutation information is lost. + + Example: + + load west0479 + A = west0479; + [L,U,P] = luinc(A,'0'); + isequal(spones(U),spones(triu(P*A))) + spones(L) ~= spones(tril(P*A)) + D = (L*U) .* spones(P*A) - P*A + + spones(L) differs from spones(tril(P*A)) at some positions on the + diagonal and at one position in L where cancellation zeroed out a + nonzero element of P*A. The entries of D are of the order of eps. + + LUINC works only for sparse matrices. + + See also LU, CHOLINC, BICG. + +*/
--- a/src/DLD-FUNCTIONS/sparse.cc +++ b/src/DLD-FUNCTIONS/sparse.cc @@ -139,7 +139,8 @@ if (nargin == 2 && ! args(0).is_scalar_type() && args(1).is_scalar_type()) mutate = (args(1).double_value() != 0.); - if (nargin == 1 || (nargin == 2 && mutate)) + if (nargin == 1 || (nargin == 2 && ! args(0).is_scalar_type() && + args(1).is_scalar_type())) { octave_value arg = args (0);
--- a/src/DLD-FUNCTIONS/splu.cc +++ b/src/DLD-FUNCTIONS/splu.cc @@ -82,8 +82,8 @@ octave_value arg = args(0); - int nr = arg.rows (); - int nc = arg.columns (); + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); int arg_is_empty = empty_arg ("splu", nr, nc); @@ -114,21 +114,21 @@ thres = tmp (0); else if (dv(0) == 1 || dv(1) == 1) { - int nel = tmp.numel (); + octave_idx_type nel = tmp.numel (); Qinit.resize (nel); - for (int i = 0; i < nel; i++) + for (octave_idx_type i = 0; i < nel; i++) Qinit (i) = tmp (i) - 1; have_Qinit = true; } else { - int t_nc = tmp.cols (); + octave_idx_type t_nc = tmp.cols (); if (tmp.nnz () != t_nc) error ("splu: Not a valid permutation matrix"); else { - for (int i = 0; i < t_nc + 1; i++) + for (octave_idx_type i = 0; i < t_nc + 1; i++) if (tmp.cidx(i) != i) { error ("splu: Not a valid permutation matrix"); @@ -138,7 +138,7 @@ if (!error_state) { - for (int i = 0; i < t_nc; i++) + for (octave_idx_type i = 0; i < t_nc; i++) if (tmp.data (i) != 1.) { error ("splu: Not a valid permutation matrix"); @@ -168,9 +168,9 @@ thres = tmp (0); else if (dv(0) == 1 || dv(1) == 1) { - int nel = tmp.numel (); + octave_idx_type nel = tmp.numel (); Qinit.resize (nel); - for (int i = 0; i < nel; i++) + for (octave_idx_type i = 0; i < nel; i++) Qinit (i) = tmp (i) - 1; have_Qinit = true; } @@ -178,13 +178,13 @@ { SparseMatrix tmp2 (tmp); - int t_nc = tmp2.cols (); + octave_idx_type t_nc = tmp2.cols (); if (tmp2.nnz () != t_nc) error ("splu: Not a valid permutation matrix"); else { - for (int i = 0; i < t_nc + 1; i++) + for (octave_idx_type i = 0; i < t_nc + 1; i++) if (tmp2.cidx(i) != i) { error ("splu: Not a valid permutation matrix"); @@ -194,7 +194,7 @@ if (!error_state) { - for (int i = 0; i < t_nc; i++) + for (octave_idx_type i = 0; i < t_nc; i++) if (tmp2.data (i) != 1.) { error ("splu: Not a valid permutation matrix"); @@ -219,9 +219,9 @@ if (nargout < 4 && ! have_Qinit) { - int m_nc = m.cols (); + octave_idx_type m_nc = m.cols (); Qinit.resize (m_nc); - for (int i = 0; i < m_nc; i++) + for (octave_idx_type i = 0; i < m_nc; i++) Qinit (i) = i; } @@ -291,9 +291,9 @@ if (nargout < 4 && ! have_Qinit) { - int m_nc = m.cols (); + octave_idx_type m_nc = m.cols (); Qinit.resize (m_nc); - for (int i = 0; i < m_nc; i++) + for (octave_idx_type i = 0; i < m_nc; i++) Qinit (i) = i; }
--- a/src/Makefile.in +++ b/src/Makefile.in @@ -50,7 +50,7 @@ lpsolve.cc lsode.cc lu.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc schur.cc sort.cc sparse.cc spdet.cc \ splu.cc spparms.cc sqrtm.cc svd.cc syl.cc time.cc gplot.l \ - __glpk__.cc + __glpk__.cc luinc.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC))
--- a/src/ov-bool-sparse.h +++ b/src/ov-bool-sparse.h @@ -82,9 +82,9 @@ octave_value *try_narrowing_conversion (void); -#if 0 - idx_vector index_vector (void) const { return idx_vector (matrix); } -#endif + // XXX FIXME XXX Adapt idx_vector to allow sparse logical indexing!! + idx_vector index_vector (void) const + { return idx_vector (bool_array_value ()); } bool is_bool_matrix (void) const { return true; }
--- a/src/ov-mapper.cc +++ b/src/ov-mapper.cc @@ -57,6 +57,25 @@ } static bool +any_element_less_than (const SparseMatrix& a, double val) +{ + octave_idx_type len = a.nonzero (); + + if (val > 0. && len != a.numel ()) + return true; + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (a.data(i) < val) + return true; + } + + return false; +} + +static bool any_element_greater_than (const NDArray& a, double val) { octave_idx_type len = a.length (); @@ -72,6 +91,25 @@ return false; } +static bool +any_element_greater_than (const SparseMatrix& a, double val) +{ + octave_idx_type len = a.nonzero (); + + if (val < 0. && len != a.numel ()) + return true; + + for (octave_idx_type i = 0; i < len; i++) + { + OCTAVE_QUIT; + + if (a.data(i) > val) + return true; + } + + return false; +} + // In most cases, we could use the map member function from the NDArray // classes, but as currently implemented, they don't allow us to // detect errors and abort properly. So use these macros to do the @@ -104,6 +142,72 @@ #define MAPPER_LOOP(T, F, M) \ MAPPER_LOOP_1 (T, F, M, ) +#define SPARSE_MAPPER_LOOP_2(T, ET, F, M, CONV, R) \ + do \ + { \ + ET f_zero = CONV (F (0.)); \ + \ + if (f_zero != 0.) \ + { \ + octave_idx_type nr = M.rows (); \ + octave_idx_type nc = M.cols (); \ + \ + T result (nr, nc, f_zero); \ + \ + for (octave_idx_type j = 0; j < nc; j++) \ + for (octave_idx_type i = M.cidx(j); i < M.cidx (j+1); i++) \ + { \ + OCTAVE_QUIT; \ + result.elem (M.ridx (i), j) = CONV (F (M.data(i))); \ + \ + if (error_state) \ + return retval; \ + } \ + \ + result.maybe_compress (true); \ + retval = R; \ + } \ + else \ + { \ + octave_idx_type nnz = M.nonzero (); \ + octave_idx_type nr = M.rows (); \ + octave_idx_type nc = M.cols (); \ + \ + T result (nr, nc, nnz); \ + ET zero = ET (0.); \ + octave_idx_type ii = 0; \ + result.cidx (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++) \ + { \ + ET val = CONV (F (M.data (i))); \ + if (val != zero) \ + { \ + result.data (ii) = val; \ + result.ridx (ii++) = M.ridx (i); \ + } \ + OCTAVE_QUIT; \ + \ + if (error_state) \ + return retval; \ + } \ + result.cidx (j+1) = ii; \ + } \ + \ + result.maybe_compress (false); \ + retval = R; \ + } \ + } \ + while (0) + +#define SPARSE_MAPPER_LOOP_1(T, ET, F, M, CONV) \ + SPARSE_MAPPER_LOOP_2 (T, ET, F, M, CONV, result) + +#define SPARSE_MAPPER_LOOP(T, ET, F, M) \ + SPARSE_MAPPER_LOOP_1 (T, ET, F, M, ) + octave_value octave_mapper::apply (const octave_value& arg) const { @@ -136,6 +240,32 @@ error ("%s: unable to handle real arguments", name().c_str ()); } + else if (arg.class_name () == "sparse") + { + const SparseMatrix m = arg.sparse_matrix_value (); + + if (error_state) + return retval; + + if (can_ret_cmplx_for_real + && (any_element_less_than (m, lower_limit) + || any_element_greater_than (m, upper_limit))) + { + if (c_c_map_fcn) + SPARSE_MAPPER_LOOP (SparseComplexMatrix, Complex, + c_c_map_fcn, m); + else + error ("%s: unable to handle real arguments", + name().c_str ()); + } + else if (d_d_map_fcn) + SPARSE_MAPPER_LOOP (SparseMatrix, double, d_d_map_fcn, m); + else if (d_b_map_fcn) + SPARSE_MAPPER_LOOP (SparseBoolMatrix, bool, d_b_map_fcn, m); + else + error ("%s: unable to handle real arguments", + name().c_str ()); + } else { NDArray m = arg.array_value (); @@ -178,6 +308,25 @@ error ("%s: unable to handle complex arguments", name().c_str ()); } + else if (arg.class_name () == "sparse") + { + SparseComplexMatrix cm = arg.sparse_complex_matrix_value (); + + if (error_state) + return retval; + + if (d_c_map_fcn) + SPARSE_MAPPER_LOOP (SparseMatrix, double, d_c_map_fcn, cm); + else if (c_c_map_fcn) + SPARSE_MAPPER_LOOP (SparseComplexMatrix, Complex, + c_c_map_fcn, cm); + else if (c_b_map_fcn) + SPARSE_MAPPER_LOOP (SparseBoolMatrix, bool, + c_b_map_fcn, cm); + else + error ("%s: unable to handle complex arguments", + name().c_str ()); + } else { ComplexNDArray cm = arg.complex_array_value ();
--- a/src/ov-re-sparse.cc +++ b/src/ov-re-sparse.cc @@ -47,6 +47,19 @@ DEFINE_OV_TYPEID_FUNCTIONS_AND_DATA (octave_sparse_matrix, "sparse matrix", "sparse"); +idx_vector +octave_sparse_matrix::index_vector (void) const +{ + if (matrix.numel () == matrix.nonzero ()) + return idx_vector (array_value ()); + else + { + std::string nm = type_name (); + error ("%s type invalid as index value", nm.c_str ()); + return idx_vector (); + } +} + octave_value * octave_sparse_matrix::try_narrowing_conversion (void) { @@ -146,11 +159,11 @@ octave_sparse_matrix::streamoff_array_value (void) const { streamoff_array retval (dims ()); - int nc = matrix.cols (); - int nr = matrix.rows (); + octave_idx_type nc = matrix.cols (); + octave_idx_type nr = matrix.rows (); - for (int j = 0; j < nc; j++) - for (int i = matrix.cidx(i); i < matrix.cidx(i+1); i++) + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = matrix.cidx(j); i < matrix.cidx(j+1); i++) { double d = matrix.data(i); @@ -169,6 +182,67 @@ return retval; } +octave_value +octave_sparse_matrix::convert_to_str_internal (bool, bool) const +{ + octave_value retval; + dim_vector dv = dims (); + octave_idx_type nel = dv.numel (); + + if (nel == 0) + { + char s = '\0'; + retval = octave_value (&s); + } + else + { + octave_idx_type nr = matrix.rows (); + octave_idx_type nc = matrix.cols (); + charNDArray chm (dv, static_cast<char> (0)); + + bool warned = false; + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = matrix.cidx(j); + i < matrix.cidx(j+1); i++) + { + OCTAVE_QUIT; + + double d = matrix.data (i); + + if (xisnan (d)) + { + ::error ("invalid conversion from NaN to character"); + return retval; + } + else + { + int ival = NINT (d); + + if (ival < 0 || ival > UCHAR_MAX) + { + // XXX FIXME XXX -- is there something + // better we could do? + + ival = 0; + + if (! warned) + { + ::warning ("range error for conversion to character value"); + warned = true; + } + } + + chm (matrix.ridx(i) + j * nr) = + static_cast<char> (ival); + } + } + retval = octave_value (chm, 1); + } + + return retval; +} + bool octave_sparse_matrix::save_binary (std::ostream& os, bool&save_as_floats) {
--- a/src/ov-re-sparse.h +++ b/src/ov-re-sparse.h @@ -81,9 +81,7 @@ octave_value *try_narrowing_conversion (void); -#if 0 - idx_vector index_vector (void) const { return idx_vector (matrix); } -#endif + idx_vector index_vector (void) const; bool is_real_matrix (void) const { return true; } @@ -114,6 +112,8 @@ streamoff_array streamoff_array_value (void) const; + octave_value convert_to_str_internal (bool pad, bool force) const; + #if 0 int write (octave_stream& os, int block_size, oct_data_conv::data_type output_type, int skip,