Mercurial > octave-nkf
changeset 7515:f3c00dc0912b
Eliminate the rest of the dispatched sparse functions
line wrap: on
line diff
--- a/doc/ChangeLog Fri Feb 22 13:47:38 2008 -0500 +++ b/doc/ChangeLog Fri Feb 22 15:50:51 2008 +0100 @@ -1,3 +1,9 @@ +2008-02-22 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Remove refernces to spdiag, spcumprod, + spcumsum, spprod, spsum, spsumsq, spchol, spchol2inv, spcholinv, + spinv and splu. + 2008-02-20 David Bateman <dbateman@free.fr> * interpreter/sparse.txi: Remove references to spmin, spmax,
--- a/doc/interpreter/sparse.txi Fri Feb 22 13:47:38 2008 -0500 +++ b/doc/interpreter/sparse.txi Fri Feb 22 15:50:51 2008 +0100 @@ -174,7 +174,7 @@ @table @asis @item Returned from a function There are many functions that directly return sparse matrices. These include -@dfn{speye}, @dfn{sprand}, @dfn{spdiag}, etc. +@dfn{speye}, @dfn{sprand}, @dfn{diag}, etc. @item Constructed from matrices or vectors The function @dfn{sparse} allows a sparse matrix to be constructed from three vectors representing the row, column and data. Alternatively, the @@ -202,22 +202,17 @@ elements of @var{d}. Other functions of interest that directly create sparse matrices, are -@dfn{spdiag} or its generalization @dfn{spdiags}, that can take the +@dfn{diag} or its generalization @dfn{spdiags}, that can take the definition of the diagonals of the matrix and create the sparse matrix that corresponds to this. For example @example -s = spdiag (sparse(randn(1,n)), -1); +s = diag (sparse(randn(1,n)), -1); @end example creates a sparse (@var{n}+1)-by-(@var{n}+1) sparse matrix with a single diagonal defined. -@DOCSTRING(spcumprod) - -@DOCSTRING(spcumsum) - -@DOCSTRING(spdiag) @DOCSTRING(spdiags) @@ -231,18 +226,12 @@ @DOCSTRING(spones) -@DOCSTRING(spprod) - @DOCSTRING(sprand) @DOCSTRING(sprandn) @DOCSTRING(sprandsym) -@DOCSTRING(spsum) - -@DOCSTRING(spsumsq) - The recommended way for the user to create a sparse matrix, is to create two vectors containing the row and column index of the data and a third vector of the same size containing the data to be stored. For example @@ -484,10 +473,7 @@ @dfn{dmperm}, @dfn{symamd}, @dfn{randperm}, @dfn{symrcm} @item Linear algebra: - @dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv}, - @dfn{spchol2inv}, @dfn{spinv}, - @dfn{splchol}, @dfn{splu}, @dfn{normest}, @dfn{condest}, - @dfn{sprank} + @dfn{matrix_type}, @dfn{normest}, @dfn{condest}, @dfn{sprank} @c @dfn{spaugment} @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now @@ -497,9 +483,7 @@ @c @dfn{lsqr}, @dfn{minres}, @dfn{qmr}, @dfn{symmlq} @item Miscellaneous: - @dfn{spparms}, @dfn{symbfact}, @dfn{spstats}, - @dfn{spprod}, @dfn{spcumsum}, @dfn{spsum}, - @dfn{spsumsq}, @dfn{spdiag} + @dfn{spparms}, @dfn{symbfact}, @dfn{spstats} @end table In addition all of the standard Octave mapper functions (ie. basic @@ -833,18 +817,6 @@ @DOCSTRING(condest) -@DOCSTRING(spchol) - -@DOCSTRING(spcholinv) - -@DOCSTRING(spchol2inv) - -@DOCSTRING(spinv) - -@DOCSTRING(splchol) - -@DOCSTRING(splu) - @DOCSTRING(spparms) @DOCSTRING(sprank)
--- a/liboctave/CSparse.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/CSparse.cc Fri Feb 22 15:50:51 2008 +0100 @@ -1042,7 +1042,7 @@ Qinit(i) = i; MatrixType tmp_typ (MatrixType::Upper); - SparseComplexLU fact (*this, Qinit, -1.0, false); + SparseComplexLU fact (*this, Qinit, Matrix (), false, false); rcond = fact.rcond(); double rcond2; SparseComplexMatrix InvL = fact.L().transpose().
--- a/liboctave/ChangeLog Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/ChangeLog Fri Feb 22 15:50:51 2008 +0100 @@ -1,3 +1,39 @@ +2008-02-22 David Bateman <dbateman@free.fr> + + * boolSparse.cc (SparseBoolMatrix SparseBoolMatrix::diag + (octave_idx_type) const): New method. + * boolSparse.h (SparseBoolMatrix SparseBoolMatrix::diag + (octave_idx_type) const): Declare it. + + * base-lu.h (lu_type Y (void) const): New method to return + factorization of xGETRF directly. + * sparse-base-lu.cc (template <class lu_type, class lu_elt_type, + class p_type, class p_elt_type> lu_type sparse_base_lu <lu_type, + lu_elt_type, p_type, p_elt_type> :: Y (void) const): New method + to simulate the retirn of xGETRF. + * sparse-base-lu.h (template <class lu_type, class lu_elt_type, + class p_type, class p_elt_type> lu_type sparse_base_lu <lu_type, + lu_elt_type, p_type, p_elt_type> :: Y (void) const): Declare it + (SparseMatrix R (void) const): Method to return scaling factors. + * SparsedbleLU.cc: Allow two element pivot thresholding and + scaling. + * SparseCmplxLU.cc: ditto. + * SparsedbleLU.h: Modify constructors to allow passing of two + element pivoting thresholds and flag for scaling + * SparseCmplxLU.h: ditto. + + * base-lu.cc (ColumnVector P_vec (void) const): New method to + return permutations as a vector. + * base-lu.h (ColumnVector P_vec (void) const): Declare it. + * sparse-base-lu.cc (ColumnVector Pr_vec (void) const): New method + return row permutations as a vector. + (ColumnVector Pc_vec (void) const): New method return column + permutations as a vector. + * sparse-base-lu.h (ColumnVector Pr_vec (void) const): Declare it. + (ColumnVector Pc_vec (void) const): Declare it. + + * oct-spparms.cc: Add sym_tol field. + 2008-02-20 David Bateman <dbateman@free.fr> * SparseComplexQR.cc (ComplexMatrix
--- a/liboctave/SparseCmplxLU.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/SparseCmplxLU.cc Fri Feb 22 15:50:51 2008 +0100 @@ -42,7 +42,7 @@ #include "oct-sparse.h" SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, - double piv_thres) + const Matrix& piv_thres, bool scale) { #ifdef HAVE_UMFPACK octave_idx_type nr = a.rows (); @@ -56,20 +56,24 @@ double tmp = octave_sparse_params::get_key ("spumoni"); if (!xisnan (tmp)) Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) + if (piv_thres.nelem() != 2) { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); + if (!xisnan (tmp)) + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } else { tmp = octave_sparse_params::get_key ("piv_tol"); if (!xisnan (tmp)) - { + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + + tmp = octave_sparse_params::get_key ("sym_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 @@ -78,7 +82,10 @@ Control (UMFPACK_FIXQ) = tmp; // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; + if (scale) + Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; + else + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; UMFPACK_ZNAME (report_control) (control); @@ -173,6 +180,15 @@ octave_idx_type *Uj = Ufact.ridx (); Complex *Ux = Ufact.data (); + Rfact = SparseMatrix (nr, nr, nr); + for (octave_idx_type i = 0; i < nr; i++) + { + Rfact.xridx (i) = i; + Rfact.xcidx (i) = i; + } + Rfact.xcidx (nr) = nr; + double *Rx = Rfact.data (); + P.resize (nr); octave_idx_type *p = P.fortran_vec (); @@ -185,11 +201,11 @@ NULL, Up, Uj, reinterpret_cast <double *> (Ux), NULL, p, q, NULL, NULL, - &do_recip, NULL, Numeric); + &do_recip, Rx, Numeric); UMFPACK_ZNAME (free_numeric) (&Numeric) ; - if (status < 0 || do_recip) + if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); @@ -200,6 +216,10 @@ { Lfact = Lfact.transpose (); + if (do_recip) + for (octave_idx_type i = 0; i < nr; i++) + Rx[i] = 1.0 / Rx[i]; + UMFPACK_ZNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (), reinterpret_cast<double *> (Lfact.data()), @@ -224,8 +244,9 @@ SparseComplexLU::SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit, - double piv_thres, bool FixedQ, - double droptol, bool milu, bool udiag) + const Matrix& piv_thres, bool scale, + bool FixedQ, double droptol, + bool milu, bool udiag) { #ifdef HAVE_UMFPACK if (milu) @@ -244,20 +265,24 @@ double tmp = octave_sparse_params::get_key ("spumoni"); if (!xisnan (tmp)) Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) + if (piv_thres.nelem() != 2) { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); + if (!xisnan (tmp)) + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } else { tmp = octave_sparse_params::get_key ("piv_tol"); if (!xisnan (tmp)) - { - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; - Control (UMFPACK_PIVOT_TOLERANCE) = tmp; - } + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + + tmp = octave_sparse_params::get_key ("sym_tol"); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } if (droptol >= 0.) @@ -274,7 +299,10 @@ } // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; + if (scale) + Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; + else + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; UMFPACK_ZNAME (report_control) (control); @@ -379,6 +407,15 @@ octave_idx_type *Uj = Ufact.ridx (); Complex *Ux = Ufact.data (); + Rfact = SparseMatrix (nr, nr, nr); + for (octave_idx_type i = 0; i < nr; i++) + { + Rfact.xridx (i) = i; + Rfact.xcidx (i) = i; + } + Rfact.xcidx (nr) = nr; + double *Rx = Rfact.data (); + P.resize (nr); octave_idx_type *p = P.fortran_vec (); @@ -392,22 +429,25 @@ NULL, Up, Uj, reinterpret_cast<double *> (Ux), NULL, p, q, NULL, NULL, - &do_recip, NULL, Numeric) ; + &do_recip, Rx, Numeric) ; UMFPACK_ZNAME (free_numeric) (&Numeric) ; - if (status < 0 || do_recip) + if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexLU::SparseComplexLU extracting LU factors failed"); - UMFPACK_ZNAME (report_status) (control, - status); + UMFPACK_ZNAME (report_status) (control, status); } else { Lfact = Lfact.transpose (); + if (do_recip) + for (octave_idx_type i = 0; i < nr; i++) + Rx[i] = 1.0 / Rx[i]; + UMFPACK_ZNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (),
--- a/liboctave/SparseCmplxLU.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/SparseCmplxLU.h Fri Feb 22 15:50:51 2008 +0100 @@ -38,10 +38,13 @@ SparseComplexLU (void) : sparse_base_lu <SparseComplexMatrix, Complex, SparseMatrix, double> () { } - SparseComplexLU (const SparseComplexMatrix& a, double piv_thres = -1); + SparseComplexLU (const SparseComplexMatrix& a, + const Matrix& piv_thres = Matrix (), + bool scale = false); SparseComplexLU (const SparseComplexMatrix& a, const ColumnVector& Qinit, - double piv_thres = -1, bool FixedQ = false, + const Matrix& piv_thres = Matrix (), + bool scale = false, bool FixedQ = false, double droptol = -1., bool milu = false, bool udiag = false);
--- a/liboctave/SparsedbleLU.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/SparsedbleLU.cc Fri Feb 22 15:50:51 2008 +0100 @@ -41,7 +41,7 @@ #include "oct-sparse.h" -SparseLU::SparseLU (const SparseMatrix& a, double piv_thres) +SparseLU::SparseLU (const SparseMatrix& a, const Matrix& piv_thres, bool scale) { #ifdef HAVE_UMFPACK octave_idx_type nr = a.rows (); @@ -56,20 +56,24 @@ if (!xisnan (tmp)) Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) + if (piv_thres.nelem() != 2) { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); + if (!xisnan (tmp)) + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } else { tmp = octave_sparse_params::get_key ("piv_tol"); if (!xisnan (tmp)) - { + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + + tmp = octave_sparse_params::get_key ("sym_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 @@ -77,8 +81,10 @@ if (!xisnan (tmp)) Control (UMFPACK_FIXQ) = tmp; - // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; + if (scale) + Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; + else + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; UMFPACK_DNAME (report_control) (control); @@ -167,6 +173,15 @@ octave_idx_type *Uj = Ufact.ridx (); double *Ux = Ufact.data (); + Rfact = SparseMatrix (nr, nr, nr); + for (octave_idx_type i = 0; i < nr; i++) + { + Rfact.xridx (i) = i; + Rfact.xcidx (i) = i; + } + Rfact.xcidx (nr) = nr; + double *Rx = Rfact.data (); + P.resize (nr); octave_idx_type *p = P.fortran_vec (); @@ -176,12 +191,12 @@ octave_idx_type do_recip; status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, NULL, - &do_recip, NULL, + &do_recip, Rx, Numeric) ; UMFPACK_DNAME (free_numeric) (&Numeric) ; - if (status < 0 || do_recip) + if (status < 0) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); @@ -192,6 +207,10 @@ { Lfact = Lfact.transpose (); + if (do_recip) + for (octave_idx_type i = 0; i < nr; i++) + Rx[i] = 1.0 / Rx[i]; + UMFPACK_DNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (), Lfact.data (), 1, control); @@ -212,8 +231,8 @@ } SparseLU::SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, - double piv_thres, bool FixedQ, double droptol, - bool milu, bool udiag) + const Matrix& piv_thres, bool scale, bool FixedQ, + double droptol, bool milu, bool udiag) { #ifdef HAVE_UMFPACK if (milu) @@ -232,20 +251,25 @@ double tmp = octave_sparse_params::get_key ("spumoni"); if (!xisnan (tmp)) Control (UMFPACK_PRL) = tmp; - if (piv_thres >= 0.) + + if (piv_thres.nelem() != 2) { - piv_thres = (piv_thres > 1. ? 1. : piv_thres); - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = piv_thres; - Control (UMFPACK_PIVOT_TOLERANCE) = piv_thres; + tmp = (piv_thres (0) > 1. ? 1. : piv_thres (0)); + if (!xisnan (tmp)) + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + tmp = (piv_thres (1) > 1. ? 1. : piv_thres (1)); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } else { tmp = octave_sparse_params::get_key ("piv_tol"); if (!xisnan (tmp)) - { - Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; - Control (UMFPACK_PIVOT_TOLERANCE) = tmp; - } + Control (UMFPACK_PIVOT_TOLERANCE) = tmp; + + tmp = octave_sparse_params::get_key ("sym_tol"); + if (!xisnan (tmp)) + Control (UMFPACK_SYM_PIVOT_TOLERANCE) = tmp; } if (droptol >= 0.) @@ -262,8 +286,10 @@ Control (UMFPACK_FIXQ) = tmp; } - // Turn-off UMFPACK scaling for LU - Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; + if (scale) + Control (UMFPACK_SCALE) = UMFPACK_SCALE_SUM; + else + Control (UMFPACK_SCALE) = UMFPACK_SCALE_NONE; UMFPACK_DNAME (report_control) (control); @@ -363,6 +389,15 @@ octave_idx_type *Uj = Ufact.ridx (); double *Ux = Ufact.data (); + Rfact = SparseMatrix (nr, nr, nr); + for (octave_idx_type i = 0; i < nr; i++) + { + Rfact.xridx (i) = i; + Rfact.xcidx (i) = i; + } + Rfact.xcidx (nr) = nr; + double *Rx = Rfact.data (); + P.resize (nr); octave_idx_type *p = P.fortran_vec (); @@ -373,11 +408,11 @@ status = UMFPACK_DNAME (get_numeric) (Ltp, Ltj, Ltx, Up, Uj, Ux, p, q, NULL, &do_recip, - NULL, Numeric) ; + Rx, Numeric) ; UMFPACK_DNAME (free_numeric) (&Numeric) ; - if (status < 0 || do_recip) + if (status < 0) { (*current_liboctave_error_handler) ("SparseLU::SparseLU extracting LU factors failed"); @@ -387,6 +422,11 @@ else { Lfact = Lfact.transpose (); + + if (do_recip) + for (octave_idx_type i = 0; i < nr; i++) + Rx[i] = 1.0 / Rx[i]; + UMFPACK_DNAME (report_matrix) (nr, n_inner, Lfact.cidx (), Lfact.ridx (),
--- a/liboctave/SparsedbleLU.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/SparsedbleLU.h Fri Feb 22 15:50:51 2008 +0100 @@ -36,11 +36,13 @@ SparseLU (void) : sparse_base_lu <SparseMatrix, double, SparseMatrix, double> () { } - SparseLU (const SparseMatrix& a, double piv_thres = -1.0); + SparseLU (const SparseMatrix& a, const Matrix& piv_thres = Matrix(), + bool scale = false); SparseLU (const SparseMatrix& a, const ColumnVector& Qinit, - double piv_thres = -1.0, bool FixedQ = false, - double droptol = -1., bool milu = false, bool udiag = false); + const Matrix& piv_thres = Matrix(), bool scale = false, + 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/liboctave/base-lu.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/base-lu.cc Fri Feb 22 15:50:51 2008 +0100 @@ -98,6 +98,37 @@ return p; } +template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> +ColumnVector +base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: P_vec (void) const +{ + octave_idx_type a_nr = a_fact.rows (); + + Array<octave_idx_type> pvt (a_nr); + + for (octave_idx_type i = 0; i < a_nr; i++) + pvt.xelem (i) = i; + + for (octave_idx_type i = 0; i < ipvt.length(); i++) + { + octave_idx_type k = ipvt.xelem (i); + + if (k != i) + { + octave_idx_type tmp = pvt.xelem (k); + pvt.xelem (k) = pvt.xelem (i); + pvt.xelem (i) = tmp; + } + } + + ColumnVector p (a_nr); + + for (octave_idx_type i = 0; i < a_nr; i++) + p.xelem (i) = static_cast<double> (pvt.xelem (i) + 1); + + return p; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/base-lu.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/base-lu.h Fri Feb 22 15:50:51 2008 +0100 @@ -24,6 +24,7 @@ #define octave_base_lu_h 1 #include "MArray.h" +#include "dColVector.h" template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> class @@ -51,8 +52,12 @@ lu_type U (void) const; + lu_type Y (void) const { return a_fact; } + p_type P (void) const; + ColumnVector P_vec (void) const; + protected: lu_type a_fact;
--- a/liboctave/boolSparse.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/boolSparse.cc Fri Feb 22 15:50:51 2008 +0100 @@ -140,6 +140,93 @@ SPARSE_ANY_OP (dim); } +SparseBoolMatrix +SparseBoolMatrix::diag (octave_idx_type k) const +{ + octave_idx_type nnr = rows (); + octave_idx_type nnc = cols (); + + if (k > 0) + nnc -= k; + else if (k < 0) + nnr += k; + + SparseBoolMatrix d; + + if (nnr > 0 && nnc > 0) + { + octave_idx_type ndiag = (nnr < nnc) ? nnr : nnc; + + // Count the number of non-zero elements + octave_idx_type nel = 0; + if (k > 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i, i+k) != 0.) + nel++; + } + else if ( k < 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i-k, i) != 0.) + nel++; + } + else + { + for (octave_idx_type i = 0; i < ndiag; i++) + if (elem (i, i) != 0.) + nel++; + } + + d = SparseBoolMatrix (ndiag, 1, nel); + d.xcidx (0) = 0; + d.xcidx (1) = nel; + + octave_idx_type ii = 0; + if (k > 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + bool tmp = elem (i, i+k); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + else if ( k < 0) + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + bool tmp = elem (i-k, i); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + else + { + for (octave_idx_type i = 0; i < ndiag; i++) + { + bool tmp = elem (i, i); + if (tmp != 0.) + { + d.xdata (ii) = tmp; + d.xridx (ii++) = i; + } + } + } + } + else + (*current_liboctave_error_handler) + ("diag: requested diagonal out of range"); + + return d; +} + boolMatrix SparseBoolMatrix::matrix_value (void) const {
--- a/liboctave/boolSparse.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/boolSparse.h Fri Feb 22 15:50:51 2008 +0100 @@ -88,6 +88,8 @@ SparseBoolMatrix concat (const SparseBoolMatrix& rb, const Array<octave_idx_type>& ra_idx); + SparseBoolMatrix diag (octave_idx_type k = 0) const; + boolMatrix matrix_value (void) const; SparseBoolMatrix squeeze (void) const;
--- a/liboctave/dSparse.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/dSparse.cc Fri Feb 22 15:50:51 2008 +0100 @@ -1116,7 +1116,7 @@ Qinit(i) = i; MatrixType tmp_typ (MatrixType::Upper); - SparseLU fact (*this, Qinit, -1.0, false); + SparseLU fact (*this, Qinit, Matrix(), false, false); rcond = fact.rcond(); double rcond2; SparseMatrix InvL = fact.L().transpose().tinverse(tmp_typ,
--- a/liboctave/oct-spparms.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/oct-spparms.cc Fri Feb 22 15:50:51 2008 +0100 @@ -111,35 +111,37 @@ void octave_sparse_params::do_defaults (void) { - params(0) = 0; // spumoni - params(1) = 1; // ths_rel - params(2) = 1; // ths_abs - params(3) = 0; // exact_d - params(4) = 3; // supernd - params(5) = 3; // rreduce - params(6) = 0.5; // wh_frac - params(7) = 1; // autommd - params(8) = 1; // autoamd - params(9) = 0.1; // piv_tol - params(10) = 0.5; // bandden - params(11) = 1; // umfpack + params(0) = 0; // spumoni + params(1) = 1; // ths_rel + params(2) = 1; // ths_abs + params(3) = 0; // exact_d + params(4) = 3; // supernd + params(5) = 3; // rreduce + params(6) = 0.5; // wh_frac + params(7) = 1; // autommd + params(8) = 1; // autoamd + params(9) = 0.1; // piv_tol + params(10) = 0.5; // bandden + params(11) = 1; // umfpack + params(12) = 0.001; // sym_tol } void octave_sparse_params::do_tight (void) { - params(0) = 0; // spumoni - params(1) = 1; // ths_rel - params(2) = 0; // ths_abs - params(3) = 1; // exact_d - params(4) = 1; // supernd - params(5) = 1; // rreduce - params(6) = 0.5; // wh_frac - params(7) = 1; // autommd - params(8) = 1; // autoamd - params(9) = 0.1; // piv_tol - params(10) = 0.5; // bandden - params(11) = 1; // umfpack + params(0) = 0; // spumoni + params(1) = 1; // ths_rel + params(2) = 0; // ths_abs + params(3) = 1; // exact_d + params(4) = 1; // supernd + params(5) = 1; // rreduce + params(6) = 0.5; // wh_frac + params(7) = 1; // autommd + params(8) = 1; // autoamd + params(9) = 0.1; // piv_tol + params(10) = 0.5; // bandden + params(11) = 1; // umfpack + params(12) = 0.001; // sym_tol } void @@ -157,6 +159,7 @@ keys(9) = "piv_tol"; keys(10) = "bandden"; keys(11) = "umfpack"; + keys(12) = "sym_tol"; } double
--- a/liboctave/oct-spparms.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/oct-spparms.h Fri Feb 22 15:50:51 2008 +0100 @@ -33,7 +33,7 @@ #include "dColVector.h" #include "dNDArray.h" -#define OCTAVE_SPARSE_CONTROLS_SIZE 12 +#define OCTAVE_SPARSE_CONTROLS_SIZE 13 class OCTAVE_API
--- a/liboctave/sparse-base-lu.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/sparse-base-lu.cc Fri Feb 22 15:50:51 2008 +0100 @@ -28,9 +28,45 @@ #include "sparse-base-lu.h" template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> +lu_type +sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Y (void) const +{ + octave_idx_type nr = Lfact.rows (); + octave_idx_type nc = Ufact.rows (); + octave_idx_type rcmin = (nr > nc ? nr : nc); + + lu_type Yout (nr, nc, Lfact.nnz() + Ufact.nnz()); + octave_idx_type ii = 0; + Yout.xcidx(0) = 0; + + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = Ufact.cidx (j); i < Ufact.cidx(j + 1); i++) + { + Yout.xridx (ii) = Ufact.ridx(i); + Yout.xdata (ii++) = Ufact.data(i); + } + if (j < rcmin) + { + // Note the +1 skips the 1.0 on the diagonal + for (octave_idx_type i = Lfact.cidx (j) + 1; + i < Lfact.cidx(j +1); i++) + { + Yout.xridx (ii) = Lfact.ridx(i); + Yout.xdata (ii++) = Lfact.data(i); + } + } + Yout.xcidx(j + 1) = ii; + } + + return Yout; +} + +template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> p_type sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr (void) const { + octave_idx_type nr = Lfact.rows (); p_type Pout (nr, nr, nr); @@ -47,6 +83,21 @@ } template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> +ColumnVector +sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pr_vec (void) const +{ + + octave_idx_type nr = Lfact.rows (); + + ColumnVector Pout (nr); + + for (octave_idx_type i = 0; i < nr; i++) + Pout.xelem (i) = static_cast<double> (P(i) + 1); + + return Pout; +} + +template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> p_type sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc (void) const { @@ -65,6 +116,21 @@ return Pout; } +template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> +ColumnVector +sparse_base_lu <lu_type, lu_elt_type, p_type, p_elt_type> :: Pc_vec (void) const +{ + + octave_idx_type nr = Lfact.rows (); + + ColumnVector Pout (nr); + + for (octave_idx_type i = 0; i < nr; i++) + Pout.xelem (i) = static_cast<double> (Q(i) + 1); + + return Pout; +} + /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/liboctave/sparse-base-lu.h Fri Feb 22 13:47:38 2008 -0500 +++ b/liboctave/sparse-base-lu.h Fri Feb 22 15:50:51 2008 +0100 @@ -26,6 +26,7 @@ #define octave_sparse_base_lu_h 1 #include "MArray.h" +#include "dSparse.h" template <class lu_type, class lu_elt_type, class p_type, class p_elt_type> class @@ -57,10 +58,18 @@ lu_type U (void) const { return Ufact; } + SparseMatrix R (void) const { return Rfact; } + + lu_type Y (void) const; + p_type Pc (void) const; p_type Pr (void) const; + ColumnVector Pc_vec (void) const; + + ColumnVector Pr_vec (void) const; + const octave_idx_type * row_perm (void) const { return P.fortran_vec (); } const octave_idx_type * col_perm (void) const { return Q.fortran_vec (); } @@ -71,6 +80,7 @@ lu_type Lfact; lu_type Ufact; + SparseMatrix Rfact; double cond;
--- a/scripts/ChangeLog Fri Feb 22 13:47:38 2008 -0500 +++ b/scripts/ChangeLog Fri Feb 22 15:50:51 2008 +0100 @@ -1,3 +1,8 @@ +2008-02-22 David Bateman <dbateman@free.fr> + + * sparse/pcg.m, sparse/spdiags, spstats.m: Remove references to + spdiag. + 2008-02-22 John W. Eaton <jwe@octave.org> * miscellaneous/fullfile.m: Improve handling of empty args and
--- a/scripts/sparse/pcg.m Fri Feb 22 13:47:38 2008 -0500 +++ b/scripts/sparse/pcg.m Fri Feb 22 15:50:51 2008 +0100 @@ -123,7 +123,7 @@ ## @example ## @group ## N = 10; -## A = spdiag ([1:N]); +## A = diag (sparse([1:N])); ## b = rand (N, 1); ## [L, U, P, Q] = luinc (A,1.e-3); ## @end group
--- a/scripts/sparse/spdiags.m Fri Feb 22 13:47:38 2008 -0500 +++ b/scripts/sparse/spdiags.m Fri Feb 22 15:50:51 2008 +0100 @@ -21,7 +21,7 @@ ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{a}, @var{c}) ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{a}) ## @deftypefnx {Function File} {@var{b} =} spdiags (@var{v}, @var{c}, @var{m}, @var{n}) -## A generalization of the function @code{spdiag}. Called with a single +## A generalization of the function @code{diag}. Called with a single ## input argument, the non-zero diagonals @var{c} of @var{A} are extracted. ## With two arguments the diagonals to extract are given by the vector ## @var{c}. @@ -56,8 +56,8 @@ if (nargin == 1 || nargin == 2) ## extract nonzero diagonals of v into A,c + [nr, nc] = size (v); [i, j, v] = find (v); - [nr, nc] = size (v); if (nargin == 1) ## c contains the active diagonals
--- a/scripts/sparse/spstats.m Fri Feb 22 13:47:38 2008 -0500 +++ b/scripts/sparse/spstats.m Fri Feb 22 15:50:51 2008 +0100 @@ -45,14 +45,14 @@ endif [n, m] = size (S); - count = spsum (sparse (i, j, 1, n, m)); + count = sum (sparse (i, j, 1, n, m)); if (nargout > 1) - mean = spsum (S) ./ count; + mean = sum (S) ./ count; endif if (nargout > 2) ## FIXME Variance with count = 0 or 1? diff = S - sparse (i, j, mean(j), n, m); - var = spsum (diff .* diff) ./ (count - 1); + var = sum (diff .* diff) ./ (count - 1); endif endfunction
--- a/src/ChangeLog Fri Feb 22 13:47:38 2008 -0500 +++ b/src/ChangeLog Fri Feb 22 15:50:51 2008 +0100 @@ -1,3 +1,31 @@ +2008-02-21 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/chol.c (Fchol, Fcholinv, Fchol2iinv): Treat + sparse matrices. Add the "lower" and "vector" flags. + * DLD-FUNCTIONS/inv.c (Finv): Treat sparse matrices. + * DLD-FUNCTION/sparse.cc (static bool is_sparse (const + octave_value&)): Remove and use arg.is_sparse_type () instead. + (Fspcumprod, Fspcumsum, Fspprod, spsum, spsumsq, spdiag): Remove. + * DLD-FUNCTIONS/splu.cc: Remove. + * DLD-FUNCTIONS/lu.cc: Treat sparse matrices, Add "vector" flag. + * DLD-FUNCTIONS/spchol.cc: Move to symbfact.cc + * DLD-FUNCTIONS/symbfact.cc: Remove cholesky functions Fspchol, + Fspcholinv, Fspchol2inv. + * DLD-FUNCTIONS/luinc.cc: Modify for new sparse LU + constructors. Ass the 'vector' flag. + * DLD-FUNCTIONS/spparms.cc: Add the sum_tol flag. + + * Makifile.in (DLD_XSRC): Remove spchol.cc, splu.cc and add + symbfact.cc + * data.cc (NATIVE_REDUCTION): Treat sparse matrices + (DATA_REDUCTION): Ditto. + (template <class T> static octave_value make_spdiag (const T&, + octave_idx_type)): New template function for sparse diag + function. Instantiate it. + (static octave_value make_diag (const octave_value&, + octave_idx_type): Use make_spdiag for sparse matrice. + (Fatan2): Compatibility fixes for mixed full/sparse matrices. + 2008-02-21 John W. Eaton <jwe@octave.org> * DLD-FUNCTIONS/fsolve.cc (fsolve_user_jacobian):
--- a/src/DLD-FUNCTIONS/chol.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/chol.cc Fri Feb 22 15:50:51 2008 +0100 @@ -27,7 +27,13 @@ #include "CmplxCHOL.h" #include "dbleCHOL.h" +#include "SparseCmplxCHOL.h" +#include "SparsedbleCHOL.h" +#include "oct-spparms.h" +#include "sparse-util.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" #include "defun-dld.h" #include "error.h" #include "gripes.h" @@ -36,7 +42,11 @@ DEFUN_DLD (chol, args, nargout, "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} chol (@var{a})\n\ +@deftypefn {Loadable Function} {@var{r}} = chol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{r}, @var{p}]} = chol (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}]} = chol (@var{s})\n\ +@deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}]} = chol (@var{s}, 'vector')\n\ +@deftypefnx {Loadable Function} {[@var{l}, @dots{}]} = chol (@dots{}, 'lower')\n\ @cindex Cholesky factorization\n\ Compute the Cholesky factor, @var{r}, of the symmetric positive definite\n\ matrix @var{a}, where\n\ @@ -48,72 +58,216 @@ @ifinfo\n\ \n\ @example\n\ -r' * r = a.\n\ +@var{r}' * @var{r} = @var{a}.\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +Called with one output argument @code{chol} fails if @var{a} or @var{s} is\n\ +not positive definite. With two or more output arguments @var{p} flags\n\ +whether the matrix was positive definite and @code{chol} does not fail. A\n\ +zero value indicated that the matrix was positive definite and the @var{r}\n\ +gives the factorization, annd @var{p} will have a positive value otherwise.\n\ +\n\ +If called with 3 outputs then a sparsity preserving row/column permutation\n\ +is applied to @var{a} prior to the factorization. That is @var{r}\n\ +is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ +@iftex\n\ +@tex\n\ +$ R^T R = Q^T A Q$.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +@var{r}' * @var{r} = @var{q}' * @var{a} * @var{q}.\n\ @end example\n\ @end ifinfo\n\ +\n\ +The sparsity preserving permutation is generally returned as a matrix.\n\ +However, given the flag 'vector', @var{q} will be returned as a vector\n\ +such that\n\ +@iftex\n\ +@tex\n\ +$ R^T R = A (Q, Q)$.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +@var{r}' * @var{r} = a (@var{q}, @var{q}).\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +Called with either a sparse or full matrix and uing the 'lower' flag,\n\ +@code{chol} returns the lower triangular factorization such that\n\ +@iftex\n\ +@tex\n\ +$ L L^T = A $.\n\ +@end tex\n\ +@end iftex\n\ +@ifinfo\n\ +\n\ +@example\n\ +@var{l} * @var{l}' = @var{a}.\n\ +@end example\n\ +@end ifinfo\n\ +\n\ +In general the lower trinagular factorization is significantly faster for\n\ +sparse matrices.\n\ @seealso{cholinv, chol2inv}\n\ @end deftypefn") { octave_value_list retval; + int nargin = args.length (); + bool LLt = false; + bool vecout = false; - int nargin = args.length (); - - if (nargin != 1 || nargout > 2) + if (nargin < 1 || nargin > 3 || nargout > 3 + || (! args(0).is_sparse_type () && nargout > 2)) { print_usage (); return retval; } - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); + int n = 1; + while (n < nargin && ! error_state) + { + std::string tmp = args(n++).string_value (); + + if (! error_state ) + { + if (tmp.compare ("vector") == 0) + vecout = true; + else if (tmp.compare ("lower") == 0) + LLt = true; + else if (tmp.compare ("upper") == 0) + LLt = false; + else + error ("chol: unexpected second or third input"); + } + else + error ("chol: expecting trailing string arguments"); + } - int arg_is_empty = empty_arg ("chol", nr, nc); + if (! error_state) + { + octave_value arg = args(0); + + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + bool natural = (nargout != 3); + + int arg_is_empty = empty_arg ("chol", nr, nc); - if (arg_is_empty < 0) - return retval; - if (arg_is_empty > 0) - return octave_value (Matrix ()); + if (arg_is_empty < 0) + return retval; + if (arg_is_empty > 0) + return octave_value (Matrix ()); + + if (arg.is_sparse_type ()) + { + if (arg.is_real_type ()) + { + SparseMatrix m = arg.sparse_matrix_value (); - if (arg.is_real_type ()) - { - Matrix m = arg.matrix_value (); + if (! error_state) + { + octave_idx_type info; + SparseCHOL fact (m, info, natural); + if (nargout == 3) + if (vecout) + retval(2) = fact.perm (); + else + retval(2) = fact.Q(); - if (! error_state) - { - octave_idx_type info; - CHOL fact (m, info); - if (nargout == 2 || info == 0) + if (nargout > 1 || info == 0) + { + retval(1) = fact.P(); + if (LLt) + retval(0) = fact.L(); + else + retval(0) = fact.R(); + } + else + error ("chol: matrix not positive definite"); + } + } + else if (arg.is_complex_type ()) { - retval(1) = static_cast<double> (info); - retval(0) = fact.chol_matrix (); + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseComplexCHOL fact (m, info, natural); + + if (nargout == 3) + if (vecout) + retval(2) = fact.perm (); + else + retval(2) = fact.Q(); + + if (nargout > 1 || info == 0) + { + retval(1) = fact.P(); + if (LLt) + retval(0) = fact.L(); + else + retval(0) = fact.R(); + } + else + error ("chol: matrix not positive definite"); + } } else - error ("chol: matrix not positive definite"); + gripe_wrong_type_arg ("chol", arg); } - } - else if (arg.is_complex_type ()) - { - ComplexMatrix m = arg.complex_matrix_value (); + else + { + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); - if (! error_state) - { - octave_idx_type info; - ComplexCHOL fact (m, info); - if (nargout == 2 || info == 0) + if (! error_state) + { + octave_idx_type info; + CHOL fact (m, info); + if (nargout == 2 || info == 0) + { + retval(1) = static_cast<double> (info); + if (LLt) + retval(0) = fact.chol_matrix ().transpose (); + else + retval(0) = fact.chol_matrix (); + } + else + error ("chol: matrix not positive definite"); + } + } + else if (arg.is_complex_type ()) { - retval(1) = static_cast<double> (info); - retval(0) = fact.chol_matrix (); + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + ComplexCHOL fact (m, info); + if (nargout == 2 || info == 0) + { + retval(1) = static_cast<double> (info); + if (LLt) + retval(0) = fact.chol_matrix ().hermitian (); + else + retval(0) = fact.chol_matrix (); + } + else + error ("chol: matrix not positive definite"); + } } else - error ("chol: matrix not positive definite"); + gripe_wrong_type_arg ("chol", arg); } } - else - { - gripe_wrong_type_arg ("chol", arg); - } return retval; } @@ -141,36 +295,72 @@ retval = Matrix (); else { - if (arg.is_real_type ()) + if (arg.is_sparse_type ()) { - Matrix m = arg.matrix_value (); - - if (! error_state) + if (arg.is_real_type ()) { - octave_idx_type info; - CHOL chol (m, info); - if (info == 0) - retval = chol.inverse (); - else - error ("cholinv: matrix not positive definite"); + SparseMatrix m = arg.sparse_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } } - } - else if (arg.is_complex_type ()) - { - ComplexMatrix m = arg.complex_matrix_value (); - - if (! error_state) + else if (arg.is_complex_type ()) { - octave_idx_type info; - ComplexCHOL chol (m, info); - if (info == 0) - retval = chol.inverse (); - else - error ("cholinv: matrix not positive definite"); + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + SparseComplexCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } } + else + gripe_wrong_type_arg ("cholinv", arg); } else - gripe_wrong_type_arg ("chol", arg); + { + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + octave_idx_type info; + CHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + octave_idx_type info; + ComplexCHOL chol (m, info); + if (info == 0) + retval = chol.inverse (); + else + error ("cholinv: matrix not positive definite"); + } + } + else + gripe_wrong_type_arg ("chol", arg); + } } } else @@ -205,22 +395,44 @@ retval = Matrix (); else { - if (arg.is_real_type ()) + if (arg.is_sparse_type ()) { - Matrix r = arg.matrix_value (); + if (arg.is_real_type ()) + { + SparseMatrix r = arg.sparse_matrix_value (); - if (! error_state) - retval = chol2inv (r); - } - else if (arg.is_complex_type ()) - { - ComplexMatrix r = arg.complex_matrix_value (); + if (! error_state) + retval = chol2inv (r); + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix r = arg.sparse_complex_matrix_value (); - if (! error_state) - retval = chol2inv (r); + if (! error_state) + retval = chol2inv (r); + } + else + gripe_wrong_type_arg ("chol2inv", arg); } else - gripe_wrong_type_arg ("chol2inv", arg); + { + if (arg.is_real_type ()) + { + Matrix r = arg.matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else if (arg.is_complex_type ()) + { + ComplexMatrix r = arg.complex_matrix_value (); + + if (! error_state) + retval = chol2inv (r); + } + else + gripe_wrong_type_arg ("chol2inv", arg); + } } } else
--- a/src/DLD-FUNCTIONS/inv.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/inv.cc Fri Feb 22 15:50:51 2008 +0100 @@ -38,6 +38,12 @@ Compute the inverse of the square matrix @var{a}. Return an estimate\n\ of the reciprocal condition number if requested, otherwise warn of an\n\ ill-conditioned matrix if the reciprocal condition number is small.\n\ +\n\ +If called with a sparse matrix, then in general @var{x} will be a full\n\ +matrix, and so if possible forming the inverse of a sparse matrix should\n\ +be avoided. It is significantly more accurate and faster to do\n\ +@code{@var{y} = @var{a} \\ @var{b}}, rather than\n\ +@code{@var{y} = inv (@var{a}) * @var{b}}.\n\ @end deftypefn") { octave_value_list retval; @@ -68,63 +74,73 @@ return retval; } + octave_value result; + octave_idx_type info; + double rcond = 0.0; if (arg.is_real_type ()) { - Matrix m = arg.matrix_value (); - - if (! error_state) + if (arg.is_sparse_type ()) { - MatrixType mattyp = args(0).matrix_type (); - - octave_idx_type info; - double rcond = 0.0; - - Matrix result = m.inverse (mattyp, info, rcond, 1); - - args(0).matrix_type (mattyp); - - if (nargout > 1) - retval(1) = rcond; - - retval(0) = result; - - volatile double xrcond = rcond; - xrcond += 1.0; - if (nargout < 2 && (info == -1 || xrcond == 1.0)) - warning ("inverse: matrix singular to machine precision,\ - rcond = %g", rcond); + SparseMatrix m = arg.sparse_matrix_value (); + if (! error_state) + { + MatrixType mattyp = args(0).matrix_type (); + result = m.inverse (mattyp, info, rcond, 1); + args(0).matrix_type (mattyp); + } + } + else + { + Matrix m = arg.matrix_value (); + if (! error_state) + { + MatrixType mattyp = args(0).matrix_type (); + result = m.inverse (mattyp, info, rcond, 1); + args(0).matrix_type (mattyp); + } } } else if (arg.is_complex_type ()) { - ComplexMatrix m = arg.complex_matrix_value (); - - if (! error_state) + if (arg.is_sparse_type ()) { - MatrixType mattyp = args(0).matrix_type (); - - octave_idx_type info; - double rcond = 0.0; - - ComplexMatrix result = m.inverse (mattyp, info, rcond, 1); + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + if (! error_state) + { + MatrixType mattyp = args(0).matrix_type (); + result = m.inverse (mattyp, info, rcond, 1); + args(0).matrix_type (mattyp); + } + } + else + { + ComplexMatrix m = arg.complex_matrix_value (); + if (! error_state) + { + MatrixType mattyp = args(0).matrix_type (); + result = m.inverse (mattyp, info, rcond, 1); + args(0).matrix_type (mattyp); + } + } - args(0).matrix_type (mattyp); - if (nargout > 1) - retval(1) = rcond; - - retval(0) = result; - - volatile double xrcond = rcond; - xrcond += 1.0; - if (nargout < 2 && (info == -1 || xrcond == 1.0)) - warning ("inverse: matrix singular to machine precision,\ - rcond = %g", rcond); - } } else + gripe_wrong_type_arg ("inv", arg); + + + if (! error_state) { - gripe_wrong_type_arg ("inv", arg); + if (nargout > 1) + retval(1) = rcond; + + retval(0) = result; + + volatile double xrcond = rcond; + xrcond += 1.0; + if (nargout < 2 && (info == -1 || xrcond == 1.0)) + warning ("inverse: matrix singular to machine precision, rcond = %g", + rcond); } return retval;
--- a/src/DLD-FUNCTIONS/lu.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/lu.cc Fri Feb 22 15:50:51 2008 +0100 @@ -27,21 +27,30 @@ #include "CmplxLU.h" #include "dbleLU.h" +#include "SparseCmplxLU.h" +#include "SparsedbleLU.h" #include "defun-dld.h" #include "error.h" #include "gripes.h" #include "oct-obj.h" #include "utils.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" DEFUN_DLD (lu, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{l}, @var{u}, @var{p}] =} lu (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}] =} lu (@var{s})\n\ +@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{p}, @var{q}, @var{r}] =} lu (@var{s})\n\ +@deftypefnx {Loadable Function} {[@dots{}] =} lu (@var{s}, @var{thres})\n\ +@deftypefnx {Loadable Function} {@var{y} =} lu (@dots{})\n\ +@deftypefnx {Loadable Function} {[@dots{}] =} lu (@dots{}, 'vector')\n\ @cindex LU decomposition\n\ -Compute the LU decomposition of @var{a}, using subroutines from\n\ -@sc{Lapack}. The result is returned in a permuted form, according to\n\ -the optional return value @var{p}. For example, given the matrix\n\ -@code{a = [1, 2; 3, 4]},\n\ +Compute the LU decomposition of @var{a}. If @var{a} is full subroutines from\n\ +@sc{Lapack} are used and if @var{a} is sparse then UMFPACK is used. The\n\ +result is returned in a permuted form, according to the optional return\n\ +value @var{p}. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\ \n\ @example\n\ [l, u, p] = lu (a)\n\ @@ -68,18 +77,92 @@ @end example\n\ \n\ The matrix is not required to be square.\n\ +\n\ +Called with two or three output arguments and a spare input matrix,\n\ +then @dfn{lu} does not attempt to perform sparsity preserving column\n\ +permutations. Called with a fourth output argument, the sparsity\n\ +preserving column transformation @var{Q} is returned, such that\n\ +@code{@var{p} * @var{a} * @var{q} = @var{l} * @var{u}}.\n\ +\n\ +Called with a fifth output argument and a sparse input matrix, then\n\ +@dfn{lu} attempts to use a scaling factor @var{r} on the input matrix\n\ +such that @code{@var{p} * (@var{r} \\ @var{a}) * @var{q} = @var{l} * @var{u}}.\n\ +This typically leads to a sparser and more stable factorsation.\n\ +\n\ +An additional input argument @var{thres}, that defines the pivoting\n\ +threshold can be given. @var{thres} can be a scalar, in which case\n\ +it defines UMFPACK pivoting tolerance for both symmetric and unsymmetric\n\ +cases. If @var{thres} is a two element vector, then the first element\n\ +defines the pivoting tolerance for the unsymmetric UMFPACK pivoting\n\ +strategy and the second the symmetric strategy. By default, the values\n\ +defined by @code{spparms} are used and are by default @code{[0.1, 0.001]}.\n\ +\n\ +Given the string argument 'vector', @dfn{lu} returns the values of @var{p}\n\ +@var{q} as vector values, such that for full matrix, @code{@var{a}\n\ +(@var{p},:) = @var{l} * @var{u}}, and @code{@var{r}(@var{p},:) * @var{a}\n\ +(:, @var{q}) = @var{l} * @var{u}}.\n\ +\n\ +With two output arguments, returns the permuted forms of the upper and\n\ +lower triangular matrices, such that @code{@var{a} = @var{l} * @var{u}}.\n\ +With one output argument @var{y}, then the matrix returned by the @sc{Lapack}\n\ +routines is returned. If the input matrix is sparse then the matrix @var{l}\n\ +is embedded into @var{u} to give a return value similar to the full case.\n\ +For both full and sparse matrices, @dfn{lu} looses the permutation\n\ +information.\n\ @end deftypefn") { octave_value_list retval; + int nargin = args.length (); + bool issparse = (nargin > 0 && args(0).is_sparse_type ()); + bool scale = (nargout == 5); - int nargin = args.length (); - - if (nargin != 1 || nargout > 3) + if (nargin < 1 || (issparse && (nargin > 3 || nargout > 5)) + || (!issparse && (nargin > 2 || nargout > 3))) { print_usage (); return retval; } + bool vecout = false; + Matrix thres; + + int n = 1; + while (n < nargin && ! error_state) + { + if (args (n).is_string ()) + { + std::string tmp = args(n++).string_value (); + + if (! error_state ) + { + if (tmp.compare ("vector") == 0) + vecout = true; + else + error ("lu: unrecognized string argument"); + } + } + else + { + Matrix tmp = args(n++).matrix_value (); + + if (! error_state ) + { + if (!issparse) + error ("lu: can not define pivoting threshold for full matrices"); + else if (tmp.nelem () == 1) + { + thres.resize(1,2); + thres(0) = tmp(0); + thres(1) = tmp(0); + } + else if (tmp.nelem () == 2) + thres = tmp; + else + error ("lu: expecting 2 element vector for thres"); + } + } + } + octave_value arg = args(0); octave_idx_type nr = arg.rows (); @@ -87,18 +170,94 @@ int arg_is_empty = empty_arg ("lu", nr, nc); - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty > 0) - return octave_value_list (3, Matrix ()); + if (issparse) + { + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (5, SparseMatrix ()); + + ColumnVector Qinit; + if (nargout < 4) + { + Qinit.resize (nc); + for (octave_idx_type i = 0; i < nc; i++) + Qinit (i) = i; + } + + if (arg.is_real_type ()) + { + SparseMatrix m = arg.sparse_matrix_value (); + + switch (nargout) + { + case 0: + case 1: + case 2: + { + SparseLU fact (m, Qinit, thres, false, true); + + if (nargout < 2) + retval (0) = fact.Y (); + else + { + SparseMatrix P = fact.Pr (); + SparseMatrix L = P.transpose () * fact.L (); + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + + retval(0) = octave_value (L, + MatrixType (MatrixType::Permuted_Lower, + nr, fact.row_perm ())); + } + } + break; - if (arg.is_real_type ()) - { - Matrix m = arg.matrix_value (); + case 3: + { + SparseLU fact (m, Qinit, thres, false, true); + + if (vecout) + retval (2) = fact.Pr_vec (); + else + retval(2) = fact.Pr (); + + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + retval(0) = octave_value (fact.L (), + MatrixType (MatrixType::Lower)); + } + break; + + case 4: + default: + { + SparseLU fact (m, thres, scale); - if (! error_state) + if (scale) + retval(4) = fact.R (); + + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + } + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + retval(0) = octave_value (fact.L (), + MatrixType (MatrixType::Lower)); + } + break; + } + } + else if (arg.is_complex_type ()) { - LU fact (m); + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); switch (nargout) { @@ -106,55 +265,154 @@ case 1: case 2: { - Matrix P = fact.P (); - Matrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); - retval(0) = L; + SparseComplexLU fact (m, Qinit, thres, false, true); + + if (nargout < 2) + retval (0) = fact.Y (); + else + { + SparseMatrix P = fact.Pr (); + SparseComplexMatrix L = P.transpose () * fact.L (); + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + + retval(0) = octave_value (L, + MatrixType (MatrixType::Permuted_Lower, + nr, fact.row_perm ())); + } } break; case 3: + { + SparseComplexLU fact (m, Qinit, thres, false, true); + + if (vecout) + retval (2) = fact.Pr_vec (); + else + retval(2) = fact.Pr (); + + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + retval(0) = octave_value (fact.L (), + MatrixType (MatrixType::Lower)); + } + break; + + case 4: default: - retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); + { + SparseComplexLU fact (m, thres, scale); + + if (scale) + retval(4) = fact.R (); + + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + } + retval(1) = octave_value (fact.U (), + MatrixType (MatrixType::Upper)); + retval(0) = octave_value (fact.L (), + MatrixType (MatrixType::Lower)); + } break; } } - } - else if (arg.is_complex_type ()) - { - ComplexMatrix m = arg.complex_matrix_value (); - - if (! error_state) - { - ComplexLU fact (m); - - switch (nargout) - { - case 0: - case 1: - case 2: - { - Matrix P = fact.P (); - ComplexMatrix L = P.transpose () * fact.L (); - retval(1) = fact.U (); - retval(0) = L; - } - break; - - case 3: - default: - retval(2) = fact.P (); - retval(1) = fact.U (); - retval(0) = fact.L (); - break; - } - } + else + gripe_wrong_type_arg ("lu", arg); } else { - gripe_wrong_type_arg ("lu", arg); + if (arg_is_empty < 0) + return retval; + else if (arg_is_empty > 0) + return octave_value_list (3, Matrix ()); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + LU fact (m); + + switch (nargout) + { + case 0: + case 1: + retval(0) = fact.Y (); + break; + + case 2: + { + Matrix P = fact.P (); + Matrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + default: + { + if (vecout) + retval(2) = fact.P_vec (); + else + retval(2) = fact.P (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + ComplexLU fact (m); + + switch (nargout) + { + case 0: + case 1: + retval(0) = fact.Y (); + break; + + case 2: + { + Matrix P = fact.P (); + ComplexMatrix L = P.transpose () * fact.L (); + retval(1) = fact.U (); + retval(0) = L; + } + break; + + case 3: + default: + { + if (vecout) + retval(2) = fact.P_vec (); + else + retval(2) = fact.P (); + retval(1) = fact.U (); + retval(0) = fact.L (); + } + break; + } + } + } + else + gripe_wrong_type_arg ("lu", arg); } return retval;
--- a/src/DLD-FUNCTIONS/luinc.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/luinc.cc Fri Feb 22 15:50:51 2008 +0100 @@ -90,6 +90,9 @@ \n\ All other fields in @var{opts} are ignored. The outputs from @dfn{luinc}\n\ are the same as for @dfn{lu}.\n\ +\n\ +Given the string argument 'vector', @dfn{luinc} returns the values of @var{p}\n\ +@var{q} as vector values.\n\ @seealso{sparse, lu, cholinc}\n\ @end deftypefn") { @@ -98,15 +101,16 @@ if (nargin == 0) print_usage (); - else if (nargin != 2) + else if (nargin < 2 || nargin > 3) error ("luinc: incorrect number of arguments"); else { bool zero_level = false; bool milu = false; bool udiag = false; - bool thresh = -1; + Matrix thresh; double droptol = -1.; + bool vecout; if (args(1).is_string ()) { @@ -137,11 +141,34 @@ } if (map.contains ("thresh")) - thresh = map.contents ("thresh")(0).double_value (); + { + thresh = map.contents ("thresh")(0).matrix_value (); + + if (thresh.nelem () == 1) + { + thresh.resize(1,2); + thresh(1) = thresh(0); + } + else if (thresh.nelem () != 2) + error ("chol: expecting 2 element vector for thresh"); + } } else droptol = args(1).double_value (); + if (nargin == 3) + { + std::string tmp = args(2).string_value (); + + if (! error_state ) + { + if (tmp.compare ("vector") == 0) + vecout = true; + else + error ("luinc: unrecognized string argument"); + } + } + // FIXME Add code for zero-level factorization if (zero_level) error ("luinc: zero-level factorization not implemented"); @@ -166,7 +193,7 @@ case 1: case 2: { - SparseLU fact (sm, Qinit, thresh, true, droptol, + SparseLU fact (sm, Qinit, thresh, false, true, droptol, milu, udiag); if (! error_state) @@ -184,12 +211,15 @@ case 3: { - SparseLU fact (sm, Qinit, thresh, true, droptol, + SparseLU fact (sm, Qinit, thresh, false, true, droptol, milu, udiag); if (! error_state) { - retval(2) = fact.Pr (); + if (vecout) + retval(2) = fact.Pr_vec (); + else + retval(2) = fact.Pr (); retval(1) = octave_value (fact.U (), MatrixType (MatrixType::Upper)); retval(0) = octave_value (fact.L (), @@ -201,13 +231,21 @@ case 4: default: { - SparseLU fact (sm, Qinit, thresh, false, droptol, + SparseLU fact (sm, Qinit, thresh, false, false, droptol, milu, udiag); if (! error_state) { - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + } retval(1) = octave_value (fact.U (), MatrixType (MatrixType::Upper)); retval(0) = octave_value (fact.L (), @@ -237,7 +275,7 @@ case 1: case 2: { - SparseComplexLU fact (sm, Qinit, thresh, true, + SparseComplexLU fact (sm, Qinit, thresh, false, true, droptol, milu, udiag); @@ -256,12 +294,15 @@ case 3: { - SparseComplexLU fact (sm, Qinit, thresh, true, + SparseComplexLU fact (sm, Qinit, thresh, false, true, droptol, milu, udiag); if (! error_state) { - retval(2) = fact.Pr (); + if (vecout) + retval(2) = fact.Pr_vec (); + else + retval(2) = fact.Pr (); retval(1) = octave_value (fact.U (), MatrixType (MatrixType::Upper)); retval(0) = octave_value (fact.L (), @@ -273,13 +314,21 @@ case 4: default: { - SparseComplexLU fact (sm, Qinit, thresh, false, + SparseComplexLU fact (sm, Qinit, thresh, false, false, droptol, milu, udiag); if (! error_state) { - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); + if (vecout) + { + retval(3) = fact.Pc_vec (); + retval(2) = fact.Pr_vec (); + } + else + { + retval(3) = fact.Pc (); + retval(2) = fact.Pr (); + } retval(1) = octave_value (fact.U (), MatrixType (MatrixType::Upper)); retval(0) = octave_value (fact.L (),
--- a/src/DLD-FUNCTIONS/sparse.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/sparse.cc Fri Feb 22 15:50:51 2008 +0100 @@ -39,12 +39,6 @@ #include "ov-cx-sparse.h" #include "ov-bool-sparse.h" -static bool -is_sparse (const octave_value& arg) -{ - return (arg.is_sparse_type ()); -} - DEFUN_DLD (issparse, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {} issparse (@var{expr})\n\ @@ -57,7 +51,7 @@ return octave_value (); } else - return octave_value (is_sparse (args(0))); + return octave_value (args(0).is_sparse_type ()); } DEFUN_DLD (sparse, args, , @@ -132,7 +126,7 @@ { octave_value arg = args (0); - if (is_sparse (arg)) + if (arg.is_sparse_type ()) { if (use_complex) { @@ -387,342 +381,6 @@ return retval; } -#define SPARSE_DIM_ARG_BODY(NAME, FUNC) \ - int nargin = args.length(); \ - octave_value retval; \ - if ((nargin != 1 ) && (nargin != 2)) \ - print_usage (); \ - else { \ - int dim = (nargin == 1 ? -1 : args(1).int_value(true) - 1); \ - if (error_state) return retval; \ - if (dim < -1 || dim > 1) { \ - error (#NAME ": invalid dimension argument = %d", dim + 1); \ - return retval; \ - } \ - if (args(0).type_id () == \ - octave_sparse_matrix::static_type_id () || args(0).type_id () == \ - octave_sparse_bool_matrix::static_type_id ()) { \ - retval = args(0).sparse_matrix_value () .FUNC (dim); \ - } else if (args(0).type_id () == \ - octave_sparse_complex_matrix::static_type_id ()) { \ - retval = args(0).sparse_complex_matrix_value () .FUNC (dim); \ - } else \ - print_usage (); \ - } \ - return retval - -// PKG_ADD: dispatch ("prod", "spprod", "sparse matrix"); -// PKG_ADD: dispatch ("prod", "spprod", "sparse complex matrix"); -// PKG_ADD: dispatch ("prod", "spprod", "sparse bool matrix"); -DEFUN_DLD (spprod, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{y} =} spprod (@var{x},@var{dim})\n\ -Product of elements along dimension @var{dim}. If @var{dim} is omitted,\n\ -it defaults to 1 (column-wise products).\n\ -@seealso{spsum, spsumsq}\n\ -@end deftypefn") -{ - SPARSE_DIM_ARG_BODY (spprod, prod); -} - -// PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse matrix"); -// PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse complex matrix"); -// PKG_ADD: dispatch ("cumprod", "spcumprod", "sparse bool matrix"); -DEFUN_DLD (spcumprod, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{y} =} spcumprod (@var{x},@var{dim})\n\ -Cumulative product of elements along dimension @var{dim}. If @var{dim}\n\ -is omitted, it defaults to 1 (column-wise cumulative products).\n\ -@seealso{spcumsum}\n\ -@end deftypefn") -{ - SPARSE_DIM_ARG_BODY (spcumprod, cumprod); -} - -// PKG_ADD: dispatch ("sum", "spsum", "sparse matrix"); -// PKG_ADD: dispatch ("sum", "spsum", "sparse complex matrix"); -// PKG_ADD: dispatch ("sum", "spsum", "sparse bool matrix"); -DEFUN_DLD (spsum, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{y} =} spsum (@var{x},@var{dim})\n\ -Sum of elements along dimension @var{dim}. If @var{dim} is omitted, it\n\ -defaults to 1 (column-wise sum).\n\ -@seealso{spprod, spsumsq}\n\ -@end deftypefn") -{ - SPARSE_DIM_ARG_BODY (spsum, sum); -} - -// PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse matrix"); -// PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse complex matrix"); -// PKG_ADD: dispatch ("cumsum", "spcumsum", "sparse bool matrix"); -DEFUN_DLD (spcumsum, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{y} =} spcumsum (@var{x},@var{dim})\n\ -Cumulative sum of elements along dimension @var{dim}. If @var{dim}\n\ -is omitted, it defaults to 1 (column-wise cumulative sums).\n\ -@seealso{spcumprod}\n\ -@end deftypefn") -{ - SPARSE_DIM_ARG_BODY (spcumsum, cumsum); -} - -// PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse matrix"); -// PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse complex matrix"); -// PKG_ADD: dispatch ("sumsq", "spsumsq", "sparse bool matrix"); -DEFUN_DLD (spsumsq, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{y} =} spsumsq (@var{x},@var{dim})\n\ -Sum of squares of elements along dimension @var{dim}. If @var{dim}\n\ -is omitted, it defaults to 1 (column-wise sum of squares).\n\ -This function is equivalent to computing\n\ -@example\n\ -spsum (x .* spconj (x), dim)\n\ -@end example\n\ -but it uses less memory and avoids calling @code{spconj} if @var{x} is\n\ -real.\n\ -@seealso{spprod, spsum}\n\ -@end deftypefn") -{ - SPARSE_DIM_ARG_BODY (spsumsq, sumsq); -} - - -static octave_value -make_spdiag (const octave_value& a, const octave_value& b) -{ - octave_value retval; - - if (a.is_complex_type ()) - { - SparseComplexMatrix m = a.sparse_complex_matrix_value (); - octave_idx_type k = b.nint_value(true); - - if (error_state) - return retval; - - octave_idx_type nr = m.rows (); - octave_idx_type nc = m.columns (); - - if (nr == 0 || nc == 0) - retval = m; - else if (nr == 1 || nc == 1) - { - octave_idx_type roff = 0; - octave_idx_type coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - k = -k; - roff = k; - coff = 0; - } - - if (nr == 1) - { - octave_idx_type n = nc + k; - octave_idx_type nz = m.nzmax (); - SparseComplexMatrix r (n, n, nz); - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) - { - r.xdata (i) = m.data (i); - r.xridx (i) = j + roff; - } - r.xcidx (j+coff+1) = m.cidx(j+1); - } - for (octave_idx_type i = nc+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - else - { - octave_idx_type n = nr + k; - octave_idx_type nz = m.nzmax (); - octave_idx_type ii = 0; - octave_idx_type ir = m.ridx(0); - SparseComplexMatrix r (n, n, nz); - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type i = 0; i < nr; i++) - { - if (ir == i) - { - r.xdata (ii) = m.data (ii); - r.xridx (ii++) = ir + roff; - if (ii != nz) - ir = m.ridx (ii); - } - r.xcidx (i+coff+1) = ii; - } - for (octave_idx_type i = nr+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - } - else - { - SparseComplexMatrix r = m.diag (k); - // Don't use numel, since it can overflow for very large matrices - if (r.rows () > 0 && r.cols () > 0) - retval = r; - } - } - else if (a.is_real_type ()) - { - SparseMatrix m = a.sparse_matrix_value (); - - octave_idx_type k = b.nint_value(true); - - if (error_state) - return retval; - - octave_idx_type nr = m.rows (); - octave_idx_type nc = m.columns (); - - if (nr == 0 || nc == 0) - retval = m; - else if (nr == 1 || nc == 1) - { - octave_idx_type roff = 0; - octave_idx_type coff = 0; - if (k > 0) - { - roff = 0; - coff = k; - } - else if (k < 0) - { - k = -k; - roff = k; - coff = 0; - } - - if (nr == 1) - { - octave_idx_type n = nc + k; - octave_idx_type nz = m.nzmax (); - SparseMatrix r (n, n, nz); - - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type j = 0; j < nc; j++) - { - for (octave_idx_type i = m.cidx(j); i < m.cidx(j+1); i++) - { - r.xdata (i) = m.data (i); - r.xridx (i) = j + roff; - } - r.xcidx (j+coff+1) = m.cidx(j+1); - } - for (octave_idx_type i = nc+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - else - { - octave_idx_type n = nr + k; - octave_idx_type nz = m.nzmax (); - octave_idx_type ii = 0; - octave_idx_type ir = m.ridx(0); - SparseMatrix r (n, n, nz); - for (octave_idx_type i = 0; i < coff+1; i++) - r.xcidx (i) = 0; - for (octave_idx_type i = 0; i < nr; i++) - { - if (ir == i) - { - r.xdata (ii) = m.data (ii); - r.xridx (ii++) = ir + roff; - if (ii != nz) - ir = m.ridx (ii); - } - r.xcidx (i+coff+1) = ii; - } - for (octave_idx_type i = nr+coff+1; i < n+1; i++) - r.xcidx (i) = nz; - retval = r; - } - } - else - { - SparseMatrix r = m.diag (k); - if (r.rows () > 0 && r.cols () > 0) - retval = r; - } - } - else - gripe_wrong_type_arg ("spdiag", a); - - return retval; -} - -static octave_value -make_spdiag (const octave_value& a) -{ - octave_value retval; - octave_idx_type nr = a.rows (); - octave_idx_type nc = a.columns (); - - if (nr == 0 || nc == 0) - retval = SparseMatrix (); - else - retval = make_spdiag (a, octave_value (0.)); - - return retval; -} - -// PKG_ADD: dispatch ("diag", "spdiag", "sparse matrix"); -// PKG_ADD: dispatch ("diag", "spdiag", "sparse complex matrix"); -// PKG_ADD: dispatch ("diag", "spdiag", "sparse bool matrix"); -DEFUN_DLD (spdiag, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spdiag (@var{v}, @var{k})\n\ -Return a diagonal matrix with the sparse vector @var{v} on diagonal\n\ -@var{k}. The second argument is optional. If it is positive, the vector is\n\ -placed on the @var{k}-th super-diagonal. If it is negative, it is placed\n\ -on the @var{-k}-th sub-diagonal. The default value of @var{k} is 0, and\n\ -the vector is placed on the main diagonal. For example,\n\ -\n\ -@example\n\ -@group\n\ -spdiag ([1, 2, 3], 1)\n\ -ans =\n\ -\n\ -Compressed Column Sparse (rows=4, cols=4, nnz=3)\n\ - (1 , 2) -> 1\n\ - (2 , 3) -> 2\n\ - (3 , 4) -> 3\n\ -@end group\n\ -@end example\n\ -\n\ -@noindent\n\ -Given a matrix argument, instead of a vector, @code{spdiag} extracts the\n\ -@var{k}-th diagonal of the sparse matrix.\n\ -@seealso{diag}\n\ -@end deftypefn") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin == 1 && args(0).is_defined ()) - retval = make_spdiag (args(0)); - else if (nargin == 2 && args(0).is_defined () && args(1).is_defined ()) - retval = make_spdiag (args(0), args(1)); - else - print_usage (); - - return retval; -} - /* ;;; Local Variables: *** ;;; mode: C++ ***
--- a/src/DLD-FUNCTIONS/spchol.cc Fri Feb 22 13:47:38 2008 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,659 +0,0 @@ -/* - -Copyright (C) 2005, 2006, 2007 David Bateman -Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Andy Adler - -This file is part of Octave. - -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 3 of the License, 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 Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#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 "SparseCmplxCHOL.h" -#include "SparsedbleCHOL.h" -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" -#include "oct-spparms.h" -#include "sparse-util.h" - -static octave_value_list -sparse_chol (const octave_value_list& args, const int nargout, - const std::string& name, const bool LLt) -{ - octave_value_list retval; - int nargin = args.length (); - - if (nargin != 1 || nargout > 3) - { - print_usage (); - return retval; - } - - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - bool natural = (nargout != 3); - - int arg_is_empty = empty_arg (name.c_str(), nr, nc); - - if (arg_is_empty < 0) - return retval; - if (arg_is_empty > 0) - return octave_value (Matrix ()); - - if (arg.is_real_type ()) - { - SparseMatrix m = arg.sparse_matrix_value (); - - if (! error_state) - { - octave_idx_type info; - SparseCHOL fact (m, info, natural); - if (nargout == 3) - retval(2) = fact.Q(); - - if (nargout > 1 || info == 0) - { - retval(1) = fact.P(); - if (LLt) - retval(0) = fact.L(); - else - retval(0) = fact.R(); - } - else - error ("%s: matrix not positive definite", name.c_str()); - } - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix m = arg.sparse_complex_matrix_value (); - - if (! error_state) - { - octave_idx_type info; - SparseComplexCHOL fact (m, info, natural); - - if (nargout == 3) - retval(2) = fact.Q(); - - if (nargout > 1 || info == 0) - { - retval(1) = fact.P(); - if (LLt) - retval(0) = fact.L(); - else - retval(0) = fact.R(); - } - else - error ("%s: matrix not positive definite", name.c_str()); - } - } - else - gripe_wrong_type_arg (name.c_str(), arg); - - return retval; -} - -// PKG_ADD: dispatch ("chol", "spchol", "sparse matrix"); -// PKG_ADD: dispatch ("chol", "spchol", "sparse complex matrix"); -// PKG_ADD: dispatch ("chol", "spchol", "sparse bool matrix"); -DEFUN_DLD (spchol, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{r} =} spchol (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{r}, @var{p}] =} spchol (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{r}, @var{p}, @var{q}] =} spchol (@var{a})\n\ -@cindex Cholesky factorization\n\ -Compute the Cholesky factor, @var{r}, of the symmetric positive definite\n\ -sparse matrix @var{a}, where\n\ -@iftex\n\ -@tex\n\ -$ R^T R = A $.\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -\n\ -@example\n\ -r' * r = a.\n\ -@end example\n\ -@end ifinfo\n\ -\n\ -If called with 2 or more outputs @var{p} is the 0 when @var{r} is positive\n\ -definite and @var{p} is a positive integer otherwise.\n\ -\n\ -If called with 3 outputs then a sparsity preserving row/column permutation\n\ -is applied to @var{a} prior to the factorization. That is @var{r}\n\ -is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ -@iftex\n\ -@tex\n\ -$ R^T R = Q A Q^T$.\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -\n\ -@example\n\ -r' * r = q * a * q'.\n\ -@end example\n\ -@end ifinfo\n\ -\n\ -Note that @code{splchol} factorization is faster and uses less memory.\n\ -@seealso{spcholinv, spchol2inv, splchol}\n\ -@end deftypefn") -{ - return sparse_chol (args, nargout, "spchol", false); -} - -// PKG_ADD: dispatch ("lchol", "splchol", "sparse matrix"); -// PKG_ADD: dispatch ("lchol", "splchol", "sparse complex matrix"); -// PKG_ADD: dispatch ("lchol", "splchol", "sparse bool matrix"); -DEFUN_DLD (splchol, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{l} =} splchol (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{p}] =} splchol (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{p}, @var{q}] =} splchol (@var{a})\n\ -@cindex Cholesky factorization\n\ -Compute the Cholesky factor, @var{l}, of the symmetric positive definite\n\ -sparse matrix @var{a}, where\n\ -@iftex\n\ -@tex\n\ -$ L L^T = A $.\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -\n\ -@example\n\ -l * l' = a.\n\ -@end example\n\ -@end ifinfo\n\ -\n\ -If called with 2 or more outputs @var{p} is the 0 when @var{l} is positive\n\ -definite and @var{l} is a positive integer otherwise.\n\ -\n\ -If called with 3 outputs that a sparsity preserving row/column permutation\n\ -is applied to @var{a} prior to the factorization. That is @var{l}\n\ -is the factorization of @code{@var{a}(@var{q},@var{q})} such that\n\ -@iftex\n\ -@tex\n\ -$ L R^T = A (Q, Q)$.\n\ -@end tex\n\ -@end iftex\n\ -@ifinfo\n\ -\n\ -@example\n\ -r * r' = a (q, q).\n\ -@end example\n\ -@end ifinfo\n\ -\n\ -Note that @code{splchol} factorization is faster and uses less memory\n\ -than @code{spchol}. @code{splchol(@var{a})} is equivalent to\n\ -@code{spchol(@var{a})'}.\n\ -@seealso{spcholinv, spchol2inv, splchol}\n\ -@end deftypefn") -{ - return sparse_chol (args, nargout, "splchol", true); -} - -// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse matrix"); -// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse complex matrix"); -// PKG_ADD: dispatch ("cholinv", "spcholinv", "sparse bool matrix"); -DEFUN_DLD (spcholinv, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spcholinv (@var{a})\n\ -Use the Cholesky factorization to compute the inverse of the\n\ -sparse symmetric positive definite matrix @var{a}.\n\ -@seealso{spchol, spchol2inv}\n\ -@end deftypefn") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin == 1) - { - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - if (nr == 0 || nc == 0) - retval = Matrix (); - else - { - if (arg.is_real_type ()) - { - SparseMatrix m = arg.sparse_matrix_value (); - - if (! error_state) - { - octave_idx_type info; - SparseCHOL chol (m, info); - if (info == 0) - retval = chol.inverse (); - else - error ("spcholinv: matrix not positive definite"); - } - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix m = arg.sparse_complex_matrix_value (); - - if (! error_state) - { - octave_idx_type info; - SparseComplexCHOL chol (m, info); - if (info == 0) - retval = chol.inverse (); - else - error ("spcholinv: matrix not positive definite"); - } - } - else - gripe_wrong_type_arg ("spcholinv", arg); - } - } - else - print_usage (); - - return retval; -} - -// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse matrix"); -// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse complex matrix"); -// PKG_ADD: dispatch ("chol2inv", "spchol2inv", "sparse bool matrix"); -DEFUN_DLD (spchol2inv, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spchol2inv (@var{u})\n\ -Invert a sparse symmetric, positive definite square matrix from its\n\ -Cholesky decomposition, @var{u}. Note that @var{u} should be an\n\ -upper-triangular matrix with positive diagonal elements.\n\ -@code{chol2inv (@var{u})} provides @code{inv (@var{u}'*@var{u})} but\n\ -it is much faster than using @code{inv}.\n\ -@seealso{spchol, spcholinv}\n\ -@end deftypefn") -{ - octave_value retval; - - int nargin = args.length (); - - if (nargin == 1) - { - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - if (nr == 0 || nc == 0) - retval = Matrix (); - else - { - if (arg.is_real_type ()) - { - SparseMatrix r = arg.sparse_matrix_value (); - - if (! error_state) - retval = chol2inv (r); - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix r = arg.sparse_complex_matrix_value (); - - if (! error_state) - retval = chol2inv (r); - } - else - gripe_wrong_type_arg ("spchol2inv", arg); - } - } - else - print_usage (); - - return retval; -} - -DEFUN_DLD (symbfact, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{s}, @var{typ}, @var{mode})\n\ -\n\ -Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\ -Where\n\ -\n\ -@table @asis\n\ -@item @var{s}\n\ -@var{s} is a complex or real sparse matrix.\n\ -\n\ -@item @var{typ}\n\ -Is the type of the factorization and can be one of\n\ -\n\ -@table @code\n\ -@item sym\n\ -Factorize @var{s}. This is the default.\n\ -\n\ -@item col\n\ -Factorize @code{@var{s}' * @var{s}}.\n\ -@item row\n\ -Factorize @code{@var{s} * @var{s}'}.\n\ -@item lo\n\ -Factorize @code{@var{s}'}\n\ -@end table\n\ -\n\ -@item @var{mode}\n\ -The default is to return the Cholesky factorization for @var{r}, and if\n\ -@var{mode} is 'L', the conjugate transpose of the Cholesky factorization\n\ -is returned. The conjugate transpose version is faster and uses less\n\ -memory, but returns the same values for @var{count}, @var{h}, @var{parent}\n\ -and @var{post} outputs.\n\ -@end table\n\ -\n\ -The output variables are\n\ -\n\ -@table @asis\n\ -@item @var{count}\n\ -The row counts of the Cholesky factorization as determined by @var{typ}.\n\ -\n\ -@item @var{h}\n\ -The height of the elimination tree.\n\ -\n\ -@item @var{parent}\n\ -The elimination tree itself.\n\ -\n\ -@item @var{post}\n\ -A sparse boolean matrix whose structure is that of the Cholesky\n\ -factorization as determined by @var{typ}.\n\ -@end table\n\ -@end deftypefn") -{ - octave_value_list retval; - int nargin = args.length (); - - if (nargin < 1 || nargin > 3 || nargout > 5) - { - print_usage (); - return retval; - } - -#ifdef HAVE_CHOLMOD - - cholmod_common Common; - cholmod_common *cm = &Common; - CHOLMOD_NAME(start) (cm); - - double spu = octave_sparse_params::get_key ("spumoni"); - if (spu == 0.) - { - cm->print = -1; - cm->print_function = NULL; - } - else - { - cm->print = static_cast<int> (spu) + 2; - cm->print_function =&SparseCholPrint; - } - - cm->error_handler = &SparseCholError; - cm->complex_divide = CHOLMOD_NAME(divcomplex); - cm->hypotenuse = CHOLMOD_NAME(hypot); - - double dummy; - cholmod_sparse Astore; - cholmod_sparse *A = &Astore; - A->packed = true; - A->sorted = true; - A->nz = NULL; -#ifdef IDX_TYPE_LONG - A->itype = CHOLMOD_LONG; -#else - A->itype = CHOLMOD_INT; -#endif - A->dtype = CHOLMOD_DOUBLE; - A->stype = 1; - A->x = &dummy; - - if (args(0).is_real_type ()) - { - const SparseMatrix a = args(0).sparse_matrix_value(); - A->nrow = a.rows(); - A->ncol = a.cols(); - A->p = a.cidx(); - A->i = a.ridx(); - A->nzmax = a.nnz(); - A->xtype = CHOLMOD_REAL; - - if (a.rows() > 0 && a.cols() > 0) - A->x = a.data(); - } - else if (args(0).is_complex_type ()) - { - const SparseComplexMatrix a = args(0).sparse_complex_matrix_value(); - A->nrow = a.rows(); - A->ncol = a.cols(); - A->p = a.cidx(); - A->i = a.ridx(); - A->nzmax = a.nnz(); - A->xtype = CHOLMOD_COMPLEX; - - if (a.rows() > 0 && a.cols() > 0) - A->x = a.data(); - } - else - gripe_wrong_type_arg ("symbfact", arg(0)); - - octave_idx_type coletree = false; - octave_idx_type n = A->nrow; - - if (nargin > 1) - { - char ch; - std::string str = args(1).string_value(); - ch = tolower (str.c_str()[0]); - if (ch == 'r') - A->stype = 0; - else if (ch == 'c') - { - n = A->ncol; - coletree = true; - A->stype = 0; - } - else if (ch == 's') - A->stype = 1; - else if (ch == 's') - A->stype = -1; - else - error ("Unrecognized typ in symbolic factorization"); - } - - if (A->stype && A->nrow != A->ncol) - error ("Matrix must be square"); - - if (!error_state) - { - OCTAVE_LOCAL_BUFFER (octave_idx_type, Parent, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n); - OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n); - - cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm); - cholmod_sparse *Aup, *Alo; - - if (A->stype == 1 || coletree) - { - Aup = A ; - Alo = F ; - } - else - { - Aup = F ; - Alo = A ; - } - - CHOLMOD_NAME(etree) (Aup, Parent, cm); - - if (cm->status < CHOLMOD_OK) - { - error("matrix corrupted"); - goto symbfact_error; - } - - if (CHOLMOD_NAME(postorder) (Parent, n, NULL, Post, cm) != n) - { - error("postorder failed"); - goto symbfact_error; - } - - CHOLMOD_NAME(rowcolcounts) (Alo, NULL, 0, Parent, Post, NULL, - ColCount, First, Level, cm); - - if (cm->status < CHOLMOD_OK) - { - error("matrix corrupted"); - goto symbfact_error; - } - - if (nargout > 4) - { - cholmod_sparse *A1, *A2; - - if (A->stype == 1) - { - A1 = A; - A2 = NULL; - } - else if (A->stype == -1) - { - A1 = F; - A2 = NULL; - } - else if (coletree) - { - A1 = F; - A2 = A; - } - else - { - A1 = A; - A2 = F; - } - - // count the total number of entries in L - octave_idx_type lnz = 0 ; - for (octave_idx_type j = 0 ; j < n ; j++) - lnz += ColCount [j] ; - - - // allocate the output matrix L (pattern-only) - SparseBoolMatrix L (n, n, lnz); - - // initialize column pointers - lnz = 0; - for (octave_idx_type j = 0 ; j < n ; j++) - { - L.xcidx(j) = lnz; - lnz += ColCount [j]; - } - L.xcidx(n) = lnz; - - - /* create a copy of the column pointers */ - octave_idx_type *W = First; - for (octave_idx_type j = 0 ; j < n ; j++) - W [j] = L.xcidx(j); - - // get workspace for computing one row of L - cholmod_sparse *R = cholmod_allocate_sparse (n, 1, n, false, true, - 0, CHOLMOD_PATTERN, cm); - octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p); - octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i); - - // compute L one row at a time - for (octave_idx_type k = 0 ; k < n ; k++) - { - // get the kth row of L and store in the columns of L - CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ; - for (octave_idx_type p = 0 ; p < Rp [1] ; p++) - L.xridx (W [Ri [p]]++) = k ; - - // add the diagonal entry - L.xridx (W [k]++) = k ; - } - - // free workspace - cholmod_free_sparse (&R, cm) ; - - - // transpose L to get R, or leave as is - if (nargin < 3) - L = L.transpose (); - - // fill numerical values of L with one's - for (octave_idx_type p = 0 ; p < lnz ; p++) - L.xdata(p) = true; - - retval(4) = L; - } - - ColumnVector tmp (n); - if (nargout > 3) - { - for (octave_idx_type i = 0; i < n; i++) - tmp(i) = Post[i] + 1; - retval(3) = tmp; - } - - if (nargout > 2) - { - for (octave_idx_type i = 0; i < n; i++) - tmp(i) = Parent[i] + 1; - retval(2) = tmp; - } - - if (nargout > 1) - { - /* compute the elimination tree height */ - octave_idx_type height = 0 ; - for (int i = 0 ; i < n ; i++) - height = (height > Level[i] ? height : Level[i]); - height++ ; - retval(1) = static_cast<double> (height); - } - - for (octave_idx_type i = 0; i < n; i++) - tmp(i) = ColCount[i]; - retval(0) = tmp; - } - - symbfact_error: -#else - error ("symbfact: not available in this version of Octave"); -#endif - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/ -
--- a/src/DLD-FUNCTIONS/splu.cc Fri Feb 22 13:47:38 2008 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,506 +0,0 @@ -/* - -Copyright (C) 2004, 2005, 2006, 2007 David Bateman -Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler - -This file is part of Octave. - -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 3 of the License, 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 Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#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 "SparseCmplxLU.h" -#include "SparsedbleLU.h" -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" - -// PKG_ADD: dispatch ("lu", "splu", "sparse matrix"); -// PKG_ADD: dispatch ("lu", "splu", "sparse complex matrix"); -// PKG_ADD: dispatch ("lu", "splu", "sparse bool matrix"); -DEFUN_DLD (splu, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{l}, @var{u}] =} splu (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}] =} splu (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}, @var{Q}] =} splu (@var{a})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}, @var{Q}] =} splu (@dots{}, @var{thres})\n\ -@deftypefnx {Loadable Function} {[@var{l}, @var{u}, @var{P}] =} splu (@dots{}, @var{Q})\n\ -@cindex LU decomposition\n\ -Compute the LU decomposition of the sparse matrix @var{a}, using\n\ -subroutines from UMFPACK. The result is returned in a permuted\n\ -form, according to the optional return values @var{P} and @var{Q}.\n\ -\n\ -Called with two or three output arguments and a single input argument,\n\ -@dfn{splu} is a replacement for @dfn{lu}, and therefore the sparsity\n\ -preserving column permutations @var{Q} are not performed. Called with\n\ -a fourth output argument, the sparsity preserving column transformation\n\ -@var{Q} is returned, such that @code{@var{P} * @var{a} * @var{Q} =\n\ -@var{l} * @var{u}}.\n\ -\n\ -An additional input argument @var{thres}, that defines the pivoting\n\ -threshold can be given. Alternatively, the desired sparsity preserving\n\ -column permutations @var{Q} can be passed. Note that @var{Q} is assumed\n\ -to be fixed if there are fewer than four output arguments. Otherwise,\n\ -the updated column permutations are returned as the fourth argument.\n\ -\n\ -With two output arguments, returns the permuted forms of the upper and\n\ -lower triangular matrices, such that @code{@var{a} = @var{l} * @var{u}}.\n\ -With two or three output arguments, if a user-defined @var{Q} is given,\n\ -then @code{@var{u} * @var{Q}'} is returned. The matrix is not required to\n\ -be square.\n\ -@seealso{sparse, spinv, colamd, symamd}\n\ -@end deftypefn") -{ - octave_value_list retval; - - int nargin = args.length (); - - if (nargin < 1 || nargin > 3 || nargout > 4) - { - print_usage (); - return retval; - } - - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - int arg_is_empty = empty_arg ("splu", nr, nc); - - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty > 0) - return octave_value_list (3, SparseMatrix ()); - - ColumnVector Qinit; - bool have_Qinit = false; - double thres = -1.; - - for (int k = 1; k < nargin; k++) - { - if (args(k).is_sparse_type ()) - { - SparseMatrix tmp = args (k).sparse_matrix_value (); - - if (error_state) - { - error ("splu: Not a valid permutation/threshold"); - return retval; - } - - dim_vector dv = tmp.dims (); - - if (dv(0) == 1 && dv(1) == 1) - thres = tmp (0); - else if (dv(0) == 1 || dv(1) == 1) - { - octave_idx_type nel = tmp.numel (); - Qinit.resize (nel); - for (octave_idx_type i = 0; i < nel; i++) - Qinit (i) = tmp (i) - 1; - have_Qinit = true; - } - else - { - octave_idx_type t_nc = tmp.cols (); - - if (tmp.nzmax () != t_nc) - error ("splu: Not a valid permutation matrix"); - else - { - for (octave_idx_type i = 0; i < t_nc + 1; i++) - if (tmp.cidx(i) != i) - { - error ("splu: Not a valid permutation matrix"); - break; - } - } - - if (!error_state) - { - for (octave_idx_type i = 0; i < t_nc; i++) - if (tmp.data (i) != 1.) - { - error ("splu: Not a valid permutation matrix"); - break; - } - else - Qinit (i) = tmp.ridx (i) - 1; - } - - if (! error_state) - have_Qinit = true; - } - } - else - { - NDArray tmp = args(k).array_value (); - - if (error_state) - return retval; - - dim_vector dv = tmp.dims (); - if (dv.length () > 2) - { - error ("splu: second argument must be a vector/matrix or a scalar"); - } - else if (dv(0) == 1 && dv(1) == 1) - thres = tmp (0); - else if (dv(0) == 1 || dv(1) == 1) - { - octave_idx_type nel = tmp.numel (); - Qinit.resize (nel); - for (octave_idx_type i = 0; i < nel; i++) - Qinit (i) = tmp (i) - 1; - have_Qinit = true; - } - else - { - SparseMatrix tmp2 (tmp); - - octave_idx_type t_nc = tmp2.cols (); - - if (tmp2.nzmax () != t_nc) - error ("splu: Not a valid permutation matrix"); - else - { - for (octave_idx_type i = 0; i < t_nc + 1; i++) - if (tmp2.cidx(i) != i) - { - error ("splu: Not a valid permutation matrix"); - break; - } - } - - if (!error_state) - { - for (octave_idx_type i = 0; i < t_nc; i++) - if (tmp2.data (i) != 1.) - { - error ("splu: Not a valid permutation matrix"); - break; - } - else - Qinit (i) = tmp2.ridx (i) - 1; - } - - if (! error_state) - have_Qinit = true; - } - } - } - - if (error_state) - return retval; - - if (arg.is_real_type ()) - { - SparseMatrix m = arg.sparse_matrix_value (); - - if (nargout < 4 && ! have_Qinit) - { - octave_idx_type m_nc = m.cols (); - Qinit.resize (m_nc); - for (octave_idx_type i = 0; i < m_nc; i++) - Qinit (i) = i; - } - - if (! error_state) - { - switch (nargout) - { - case 0: - case 1: - case 2: - { - SparseLU fact (m, Qinit, thres, true); - - SparseMatrix P = fact.Pr (); - SparseMatrix L = P.transpose () * fact.L (); - if (have_Qinit) - retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), - MatrixType (MatrixType::Permuted_Upper, nc, fact.col_perm ())); - else - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - - retval(0) = octave_value (L, - MatrixType (MatrixType::Permuted_Lower, nr, fact.row_perm ())); - } - break; - - case 3: - { - SparseLU fact (m, Qinit, thres, true); - - retval(2) = fact.Pr (); - if (have_Qinit) - retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), - MatrixType (MatrixType::Permuted_Upper, nc, fact.col_perm ())); - else - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - break; - - case 4: - default: - { - if (have_Qinit) - { - SparseLU fact (m, Qinit, thres, false); - - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - else - { - SparseLU fact (m, thres); - - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - } - break; - } - } - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix m = arg.sparse_complex_matrix_value (); - - if (nargout < 4 && ! have_Qinit) - { - octave_idx_type m_nc = m.cols (); - Qinit.resize (m_nc); - for (octave_idx_type i = 0; i < m_nc; i++) - Qinit (i) = i; - } - - if (! error_state) - { - switch (nargout) - { - case 0: - case 1: - case 2: - { - SparseComplexLU fact (m, Qinit, thres, true); - - SparseMatrix P = fact.Pr (); - SparseComplexMatrix L = P.transpose () * fact.L (); - - if (have_Qinit) - retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), - MatrixType (MatrixType::Permuted_Upper, nc, fact.col_perm ())); - else - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - - retval(0) = octave_value (L, - MatrixType (MatrixType::Permuted_Lower, nr, fact.row_perm ())); - } - break; - - case 3: - { - SparseComplexLU fact (m, Qinit, thres, true); - - retval(2) = fact.Pr (); - if (have_Qinit) - retval(1) = octave_value (fact.U () * fact.Pc ().transpose (), - MatrixType (MatrixType::Permuted_Upper, nc, fact.col_perm ())); - else - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - break; - - case 4: - default: - { - if (have_Qinit) - { - SparseComplexLU fact (m, Qinit, thres, false); - - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - else - { - SparseComplexLU fact (m, thres); - - retval(3) = fact.Pc (); - retval(2) = fact.Pr (); - retval(1) = octave_value (fact.U (), - MatrixType (MatrixType::Upper)); - retval(0) = octave_value (fact.L (), - MatrixType (MatrixType::Lower)); - } - } - break; - } - } - } - else - { - gripe_wrong_type_arg ("splu", arg); - } - - return retval; -} - -// PKG_ADD: dispatch ("inv", "spinv", "sparse matrix"); -// PKG_ADD: dispatch ("inv", "spinv", "sparse complex matrix"); -// PKG_ADD: dispatch ("inv", "spinv", "sparse bool matrix"); -// PKG_ADD: dispatch ("inverse", "spinv", "sparse matrix"); -// PKG_ADD: dispatch ("inverse", "spinv", "sparse complex matrix"); -// PKG_ADD: dispatch ("inverse", "spinv", "sparse bool matrix"); -DEFUN_DLD (spinv, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{x}, @var{rcond}] = } spinv (@var{a}, @var{Q})\n\ -Compute the inverse of the sparse square matrix @var{a}. Return an estimate\n\ -of the reciprocal condition number if requested, otherwise warn of an\n\ -ill-conditioned matrix if the reciprocal condition number is small.\n\ -This function takes advantage of the sparsity of the matrix to accelerate\n\ -the calculation of the inverse.\n\ -\n\ -In general @var{x} will be a full matrix, and so if possible forming the\n\ -inverse of a sparse matrix should be avoided. It is significantly more\n\ -accurate and faster to do @code{@var{y} = @var{a} \\ @var{b}}, rather\n\ -than @code{@var{y} = spinv (@var{a}) * @var{b}}.\n\ -@end deftypefn") -{ - octave_value_list retval; - - int nargin = args.length (); - - if (nargin != 1) - { - print_usage (); - return retval; - } - - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - int arg_is_empty = empty_arg ("spinverse", nr, nc); - - if (arg_is_empty < 0) - return retval; - else if (arg_is_empty > 0) - return octave_value (Matrix ()); - - if (nr != nc) - { - gripe_square_matrix_required ("spinverse"); - return retval; - } - - if (arg.is_real_type ()) - { - - SparseMatrix m = arg.sparse_matrix_value (); - - if (! error_state) - { - MatrixType mattyp = args(0).matrix_type (); - - octave_idx_type info; - double rcond = 0.0; - SparseMatrix result = m.inverse (mattyp, info, rcond, 1); - - args(0).matrix_type (mattyp); - - if (nargout > 1) - retval(1) = rcond; - - retval(0) = result; - - volatile double xrcond = rcond; - xrcond += 1.0; - if (nargout < 2 && (info == -1 || xrcond == 1.0)) - warning ("spinverse: matrix singular to machine precision,\ - rcond = %g", rcond); - } - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix m = arg.sparse_complex_matrix_value (); - - if (! error_state) - { - MatrixType mattyp = args(0).matrix_type (); - - octave_idx_type info; - double rcond = 0.0; - - SparseComplexMatrix result = m.inverse (mattyp, info, rcond, 1); - - args(0).matrix_type (mattyp); - - if (nargout > 1) - retval(1) = rcond; - - retval(0) = result; - - volatile double xrcond = rcond; - xrcond += 1.0; - if (nargout < 2 && (info == -1 || xrcond == 1.0)) - warning ("spinverse: matrix singular to machine precision,\ - rcond = %g", rcond); - } - } - else - gripe_wrong_type_arg ("spinverse", arg); - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/DLD-FUNCTIONS/spparms.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/DLD-FUNCTIONS/spparms.cc Fri Feb 22 15:50:51 2008 +0100 @@ -72,8 +72,11 @@ use the sparsity preserving amd functions (default 1)\n\ @item piv_tol\n\ The pivot tolerance of the UMFPACK solvers (default 0.1)\n\ +@item sym_tol\n\ +The pivot tolerance of the UMFPACK symmetric solvers (default 0.001)\n\ @item bandden\n\ -?? (default 0.5)\n\ +The density of non-zero elements in a banded matrix before it is treated\n\ +by the LAPACK banded solvers (default 0.5)\n\ @item umfpack\n\ Flag whether the UMFPACK or mmd solvers are used for the LU, '\\' and\n\ '/' operations (default 1)\n\
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/symbfact.cc Fri Feb 22 15:50:51 2008 +0100 @@ -0,0 +1,365 @@ +/* + +Copyright (C) 2005, 2006, 2007 David Bateman +Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Andy Adler + +This file is part of Octave. + +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 3 of the License, 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 Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "SparseCmplxCHOL.h" +#include "SparsedbleCHOL.h" +#include "oct-spparms.h" +#include "sparse-util.h" + +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" + +DEFUN_DLD (symbfact, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {[@var{count}, @var{h}, @var{parent}, @var{post}, @var{r}] =} symbfact (@var{s}, @var{typ}, @var{mode})\n\ +\n\ +Performs a symbolic factorization analysis on the sparse matrix @var{s}.\n\ +Where\n\ +\n\ +@table @asis\n\ +@item @var{s}\n\ +@var{s} is a complex or real sparse matrix.\n\ +\n\ +@item @var{typ}\n\ +Is the type of the factorization and can be one of\n\ +\n\ +@table @code\n\ +@item sym\n\ +Factorize @var{s}. This is the default.\n\ +\n\ +@item col\n\ +Factorize @code{@var{s}' * @var{s}}.\n\ +@item row\n\ +Factorize @code{@var{s} * @var{s}'}.\n\ +@item lo\n\ +Factorize @code{@var{s}'}\n\ +@end table\n\ +\n\ +@item @var{mode}\n\ +The default is to return the Cholesky factorization for @var{r}, and if\n\ +@var{mode} is 'L', the conjugate transpose of the Cholesky factorization\n\ +is returned. The conjugate transpose version is faster and uses less\n\ +memory, but returns the same values for @var{count}, @var{h}, @var{parent}\n\ +and @var{post} outputs.\n\ +@end table\n\ +\n\ +The output variables are\n\ +\n\ +@table @asis\n\ +@item @var{count}\n\ +The row counts of the Cholesky factorization as determined by @var{typ}.\n\ +\n\ +@item @var{h}\n\ +The height of the elimination tree.\n\ +\n\ +@item @var{parent}\n\ +The elimination tree itself.\n\ +\n\ +@item @var{post}\n\ +A sparse boolean matrix whose structure is that of the Cholesky\n\ +factorization as determined by @var{typ}.\n\ +@end table\n\ +@end deftypefn") +{ + octave_value_list retval; + int nargin = args.length (); + + if (nargin < 1 || nargin > 3 || nargout > 5) + { + print_usage (); + return retval; + } + +#ifdef HAVE_CHOLMOD + + cholmod_common Common; + cholmod_common *cm = &Common; + CHOLMOD_NAME(start) (cm); + + double spu = octave_sparse_params::get_key ("spumoni"); + if (spu == 0.) + { + cm->print = -1; + cm->print_function = NULL; + } + else + { + cm->print = static_cast<int> (spu) + 2; + cm->print_function =&SparseCholPrint; + } + + cm->error_handler = &SparseCholError; + cm->complex_divide = CHOLMOD_NAME(divcomplex); + cm->hypotenuse = CHOLMOD_NAME(hypot); + + double dummy; + cholmod_sparse Astore; + cholmod_sparse *A = &Astore; + A->packed = true; + A->sorted = true; + A->nz = NULL; +#ifdef IDX_TYPE_LONG + A->itype = CHOLMOD_LONG; +#else + A->itype = CHOLMOD_INT; +#endif + A->dtype = CHOLMOD_DOUBLE; + A->stype = 1; + A->x = &dummy; + + if (args(0).is_real_type ()) + { + const SparseMatrix a = args(0).sparse_matrix_value(); + A->nrow = a.rows(); + A->ncol = a.cols(); + A->p = a.cidx(); + A->i = a.ridx(); + A->nzmax = a.nnz(); + A->xtype = CHOLMOD_REAL; + + if (a.rows() > 0 && a.cols() > 0) + A->x = a.data(); + } + else if (args(0).is_complex_type ()) + { + const SparseComplexMatrix a = args(0).sparse_complex_matrix_value(); + A->nrow = a.rows(); + A->ncol = a.cols(); + A->p = a.cidx(); + A->i = a.ridx(); + A->nzmax = a.nnz(); + A->xtype = CHOLMOD_COMPLEX; + + if (a.rows() > 0 && a.cols() > 0) + A->x = a.data(); + } + else + gripe_wrong_type_arg ("symbfact", arg(0)); + + octave_idx_type coletree = false; + octave_idx_type n = A->nrow; + + if (nargin > 1) + { + char ch; + std::string str = args(1).string_value(); + ch = tolower (str.c_str()[0]); + if (ch == 'r') + A->stype = 0; + else if (ch == 'c') + { + n = A->ncol; + coletree = true; + A->stype = 0; + } + else if (ch == 's') + A->stype = 1; + else if (ch == 's') + A->stype = -1; + else + error ("Unrecognized typ in symbolic factorization"); + } + + if (A->stype && A->nrow != A->ncol) + error ("Matrix must be square"); + + if (!error_state) + { + OCTAVE_LOCAL_BUFFER (octave_idx_type, Parent, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Post, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, ColCount, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, First, n); + OCTAVE_LOCAL_BUFFER (octave_idx_type, Level, n); + + cholmod_sparse *F = CHOLMOD_NAME(transpose) (A, 0, cm); + cholmod_sparse *Aup, *Alo; + + if (A->stype == 1 || coletree) + { + Aup = A ; + Alo = F ; + } + else + { + Aup = F ; + Alo = A ; + } + + CHOLMOD_NAME(etree) (Aup, Parent, cm); + + if (cm->status < CHOLMOD_OK) + { + error("matrix corrupted"); + goto symbfact_error; + } + + if (CHOLMOD_NAME(postorder) (Parent, n, NULL, Post, cm) != n) + { + error("postorder failed"); + goto symbfact_error; + } + + CHOLMOD_NAME(rowcolcounts) (Alo, NULL, 0, Parent, Post, NULL, + ColCount, First, Level, cm); + + if (cm->status < CHOLMOD_OK) + { + error("matrix corrupted"); + goto symbfact_error; + } + + if (nargout > 4) + { + cholmod_sparse *A1, *A2; + + if (A->stype == 1) + { + A1 = A; + A2 = NULL; + } + else if (A->stype == -1) + { + A1 = F; + A2 = NULL; + } + else if (coletree) + { + A1 = F; + A2 = A; + } + else + { + A1 = A; + A2 = F; + } + + // count the total number of entries in L + octave_idx_type lnz = 0 ; + for (octave_idx_type j = 0 ; j < n ; j++) + lnz += ColCount [j] ; + + + // allocate the output matrix L (pattern-only) + SparseBoolMatrix L (n, n, lnz); + + // initialize column pointers + lnz = 0; + for (octave_idx_type j = 0 ; j < n ; j++) + { + L.xcidx(j) = lnz; + lnz += ColCount [j]; + } + L.xcidx(n) = lnz; + + + /* create a copy of the column pointers */ + octave_idx_type *W = First; + for (octave_idx_type j = 0 ; j < n ; j++) + W [j] = L.xcidx(j); + + // get workspace for computing one row of L + cholmod_sparse *R = cholmod_allocate_sparse (n, 1, n, false, true, + 0, CHOLMOD_PATTERN, cm); + octave_idx_type *Rp = static_cast<octave_idx_type *>(R->p); + octave_idx_type *Ri = static_cast<octave_idx_type *>(R->i); + + // compute L one row at a time + for (octave_idx_type k = 0 ; k < n ; k++) + { + // get the kth row of L and store in the columns of L + CHOLMOD_NAME (row_subtree) (A1, A2, k, Parent, R, cm) ; + for (octave_idx_type p = 0 ; p < Rp [1] ; p++) + L.xridx (W [Ri [p]]++) = k ; + + // add the diagonal entry + L.xridx (W [k]++) = k ; + } + + // free workspace + cholmod_free_sparse (&R, cm) ; + + + // transpose L to get R, or leave as is + if (nargin < 3) + L = L.transpose (); + + // fill numerical values of L with one's + for (octave_idx_type p = 0 ; p < lnz ; p++) + L.xdata(p) = true; + + retval(4) = L; + } + + ColumnVector tmp (n); + if (nargout > 3) + { + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = Post[i] + 1; + retval(3) = tmp; + } + + if (nargout > 2) + { + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = Parent[i] + 1; + retval(2) = tmp; + } + + if (nargout > 1) + { + /* compute the elimination tree height */ + octave_idx_type height = 0 ; + for (int i = 0 ; i < n ; i++) + height = (height > Level[i] ? height : Level[i]); + height++ ; + retval(1) = static_cast<double> (height); + } + + for (octave_idx_type i = 0; i < n; i++) + tmp(i) = ColCount[i]; + retval(0) = tmp; + } + + symbfact_error: +#else + error ("symbfact: not available in this version of Octave"); +#endif + + return retval; +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
--- a/src/Makefile.in Fri Feb 22 13:47:38 2008 -0500 +++ b/src/Makefile.in Fri Feb 22 15:50:51 2008 +0100 @@ -70,8 +70,8 @@ givens.cc hess.cc inv.cc kron.cc lsode.cc \ lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \ - spchol.cc splu.cc spparms.cc \ - sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \ + spparms.cc sqrtm.cc svd.cc syl.cc symrcm.cc symbfact.cc \ + time.cc tsearch.cc typecast.cc \ urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \ __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \ __qp__.cc __voronoi__.cc
--- a/src/data.cc Fri Feb 22 13:47:38 2008 -0500 +++ b/src/data.cc Fri Feb 22 15:50:51 2008 +0100 @@ -455,20 +455,11 @@ if (! error_state) { - if (arg_x.is_sparse_type ()) - { - SparseMatrix x = arg_x.sparse_matrix_value (); - - if (! error_state) - retval = map_d_s (atan2, y, x); - } - else - { - NDArray x = arg_x.array_value (); - - if (! error_state) - retval = map_d_m (atan2, y, x); - } + // Even if x is sparse return a full matrix here + NDArray x = arg_x.array_value (); + + if (! error_state) + retval = map_d_m (atan2, y, x); } } else if (x_is_scalar) @@ -480,7 +471,7 @@ if (! error_state) { double x = arg_x.double_value (); - + if (! error_state) retval = map_s_d (atan2, y, x); } @@ -500,7 +491,8 @@ } else if (y_dims == x_dims) { - if (arg_y.is_sparse_type () || arg_x.is_sparse_type ()) + // Even if y is sparse return a full matrix here + if (arg_x.is_sparse_type ()) { SparseMatrix y = arg_y.sparse_matrix_value (); @@ -712,21 +704,67 @@ { \ if (dim >= -1) \ { \ - if (isnative) \ - { \ - if NATIVE_REDUCTION_1 (FCN, uint8, dim) \ - else if NATIVE_REDUCTION_1 (FCN, uint16, dim) \ - else if NATIVE_REDUCTION_1 (FCN, uint32, dim) \ - else if NATIVE_REDUCTION_1 (FCN, uint64, dim) \ - else if NATIVE_REDUCTION_1 (FCN, int8, dim) \ - else if NATIVE_REDUCTION_1 (FCN, int16, dim) \ - else if NATIVE_REDUCTION_1 (FCN, int32, dim) \ - else if NATIVE_REDUCTION_1 (FCN, int64, dim) \ - else if NATIVE_REDUCTION_1 (FCN, bool, dim) \ - else if (arg.is_char_matrix ()) \ - { \ - error (#FCN, ": invalid char type"); \ + if (arg.is_sparse_type ()) \ + { \ + if (arg.is_real_type ()) \ + { \ + SparseMatrix tmp = arg.sparse_matrix_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + else \ + { \ + SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + } \ + else \ + { \ + if (isnative) \ + { \ + if NATIVE_REDUCTION_1 (FCN, uint8, dim) \ + else if NATIVE_REDUCTION_1 (FCN, uint16, dim) \ + else if NATIVE_REDUCTION_1 (FCN, uint32, dim) \ + else if NATIVE_REDUCTION_1 (FCN, uint64, dim) \ + else if NATIVE_REDUCTION_1 (FCN, int8, dim) \ + else if NATIVE_REDUCTION_1 (FCN, int16, dim) \ + else if NATIVE_REDUCTION_1 (FCN, int32, dim) \ + else if NATIVE_REDUCTION_1 (FCN, int64, dim) \ + else if NATIVE_REDUCTION_1 (FCN, bool, dim) \ + else if (arg.is_char_matrix ()) \ + { \ + error (#FCN, ": invalid char type"); \ + } \ + else if (arg.is_complex_type ()) \ + { \ + ComplexNDArray tmp = arg.complex_array_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + else if (arg.is_real_type ()) \ + { \ + NDArray tmp = arg.array_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + else \ + { \ + gripe_wrong_type_arg (#FCN, arg); \ + return retval; \ + } \ } \ + else if (arg.is_real_type ()) \ + { \ + NDArray tmp = arg.array_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ else if (arg.is_complex_type ()) \ { \ ComplexNDArray tmp = arg.complex_array_value (); \ @@ -734,37 +772,11 @@ if (! error_state) \ retval = tmp.FCN (dim); \ } \ - else if (arg.is_real_type ()) \ - { \ - NDArray tmp = arg.array_value (); \ - \ - if (! error_state) \ - retval = tmp.FCN (dim); \ - } \ - else \ + else \ { \ gripe_wrong_type_arg (#FCN, arg); \ return retval; \ } \ - } \ - else if (arg.is_real_type ()) \ - { \ - NDArray tmp = arg.array_value (); \ - \ - if (! error_state) \ - retval = tmp.FCN (dim); \ - } \ - else if (arg.is_complex_type ()) \ - { \ - ComplexNDArray tmp = arg.complex_array_value (); \ - \ - if (! error_state) \ - retval = tmp.FCN (dim); \ - } \ - else \ - { \ - gripe_wrong_type_arg (#FCN, arg); \ - return retval; \ } \ } \ else \ @@ -795,17 +807,37 @@ { \ if (arg.is_real_type ()) \ { \ - NDArray tmp = arg.array_value (); \ + if (arg.is_sparse_type ()) \ + { \ + SparseMatrix tmp = arg.sparse_matrix_value (); \ \ - if (! error_state) \ - retval = tmp.FCN (dim); \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + else \ + { \ + NDArray tmp = arg.array_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ } \ else if (arg.is_complex_type ()) \ { \ - ComplexNDArray tmp = arg.complex_array_value (); \ + if (arg.is_sparse_type ()) \ + { \ + SparseComplexMatrix tmp = arg.sparse_complex_matrix_value (); \ \ - if (! error_state) \ - retval = tmp.FCN (dim); \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ + else \ + { \ + ComplexNDArray tmp = arg.complex_array_value (); \ + \ + if (! error_state) \ + retval = tmp.FCN (dim); \ + } \ } \ else \ { \ @@ -946,6 +978,94 @@ make_diag (const uint64NDArray& v, octave_idx_type k); #endif +template <class T> +static octave_value +make_spdiag (const T& v, octave_idx_type k) +{ + octave_value retval; + dim_vector dv = v.dims (); + octave_idx_type nr = dv (0); + octave_idx_type nc = dv (1); + + if (nr == 0 || nc == 0) + retval = T (); + else if (nr != 1 && nc != 1) + retval = v.diag (k); + else + { + octave_idx_type roff = 0; + octave_idx_type coff = 0; + if (k > 0) + { + roff = 0; + coff = k; + } + else if (k < 0) + { + roff = -k; + coff = 0; + } + + if (nr == 1) + { + octave_idx_type n = nc + std::abs (k); + octave_idx_type nz = v.nzmax (); + T r (n, n, nz); + for (octave_idx_type i = 0; i < coff+1; i++) + r.xcidx (i) = 0; + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++) + { + r.xdata (i) = v.data (i); + r.xridx (i) = j + roff; + } + r.xcidx (j+coff+1) = v.cidx(j+1); + } + for (octave_idx_type i = nc+coff+1; i < n+1; i++) + r.xcidx (i) = nz; + retval = r; + } + else + { + octave_idx_type n = nr + std::abs (k); + octave_idx_type nz = v.nzmax (); + octave_idx_type ii = 0; + octave_idx_type ir = v.ridx(0); + T r (n, n, nz); + for (octave_idx_type i = 0; i < coff+1; i++) + r.xcidx (i) = 0; + for (octave_idx_type i = 0; i < nr; i++) + { + if (ir == i) + { + r.xdata (ii) = v.data (ii); + r.xridx (ii++) = ir + roff; + if (ii != nz) + ir = v.ridx (ii); + } + r.xcidx (i+coff+1) = ii; + } + for (octave_idx_type i = nr+coff+1; i < n+1; i++) + r.xcidx (i) = nz; + retval = r; + } + } + + return retval; +} + +#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL) +static octave_value +make_spdiag (const SparseMatrix& v, octave_idx_type k); + +static octave_value +make_spdiag (const SparseComplexMatrix& v, octave_idx_type k); + +static octave_value +make_spdiag (const SparseBoolMatrix& v, octave_idx_type k); +#endif + static octave_value make_diag (const octave_value& a, octave_idx_type k) { @@ -954,17 +1074,35 @@ if (result_type == "double") { - if (a.is_real_type ()) + if (a.is_sparse_type ()) { - Matrix m = a.matrix_value (); - if (!error_state) - retval = make_diag (m, k); + if (a.is_real_type ()) + { + SparseMatrix m = a.sparse_matrix_value (); + if (!error_state) + retval = make_spdiag (m, k); + } + else + { + SparseComplexMatrix m = a.sparse_complex_matrix_value (); + if (!error_state) + retval = make_spdiag (m, k); + } } else { - ComplexMatrix m = a.complex_matrix_value (); - if (!error_state) - retval = make_diag (m, k); + if (a.is_real_type ()) + { + Matrix m = a.matrix_value (); + if (!error_state) + retval = make_diag (m, k); + } + else + { + ComplexMatrix m = a.complex_matrix_value (); + if (!error_state) + retval = make_diag (m, k); + } } } #if 0 @@ -983,9 +1121,18 @@ } else if (result_type == "logical") { - boolMatrix m = a.bool_matrix_value (); - if (!error_state) - retval = make_diag (m, k); + if (a.is_sparse_type ()) + { + SparseBoolMatrix m = a.sparse_bool_matrix_value (); + if (!error_state) + retval = make_spdiag (m, k); + } + else + { + boolMatrix m = a.bool_matrix_value (); + if (!error_state) + retval = make_diag (m, k); + } } else if (result_type == "int8") retval = make_diag (a.int8_array_value (), k);
--- a/test/ChangeLog Fri Feb 22 13:47:38 2008 -0500 +++ b/test/ChangeLog Fri Feb 22 15:50:51 2008 +0100 @@ -1,3 +1,10 @@ +2008-02-22 David Bateman <dbateman@free.fr> + + * build_sparse_tests.sh: Replaced removed sparse functions like + spdiag with their generic names. Fix lu tests for modified + syntax. Test vector and scaling or LU and chol functions. + * test_linalg.m: Change error message of failing chol/lu test. + 2008-02-19 David Bateman <dbateman@free.fr> * build_sparse_tests.sh: Replaced removed spars functions like
--- a/test/build_sparse_tests.sh Fri Feb 22 13:47:38 2008 -0500 +++ b/test/build_sparse_tests.sh Fri Feb 22 15:50:51 2008 +0100 @@ -180,11 +180,11 @@ %% of a singular matrix, which is consistent with the full matrix %% behaviour. They are therefore disabled.. %!testif HAVE_UMFPACK -%! assert(spinv(sparse([1,1;1,1+i])),sparse([1-1i,1i;1i,-1i]),10*eps); -% !error spinv( sparse( [1,1;1,1] ) ); -% !error spinv( sparse( [0,0;0,1] ) ); -% !error spinv( sparse( [0,0;0,1+i] ) ); -% !error spinv( sparse( [0,0;0,0] ) ); +%! assert(inv(sparse([1,1;1,1+i])),sparse([1-1i,1i;1i,-1i]),10*eps); +% !error inv( sparse( [1,1;1,1] ) ); +% !error inv( sparse( [0,0;0,1] ) ); +% !error inv( sparse( [0,0;0,1+i] ) ); +% !error inv( sparse( [0,0;0,0] ) ); %% error handling in constructor %!error sparse(1,[2,3],[1,2,3]); @@ -388,18 +388,18 @@ gen_matrixdiag_tests() { cat >>$TESTS <<EOF %% Matrix diagonal tests (uses af,as,bf,bs) -%!assert(spdiag(as),sparse(diag(af))) -%!assert(spdiag(bs),sparse(diag(bf))) -%!assert(spdiag(as,1),sparse(diag(af,1))) -%!assert(spdiag(bs,1),sparse(diag(bf,1))) -%!assert(spdiag(as,-1),sparse(diag(af,-1))) -%!assert(spdiag(bs,-1),sparse(diag(bf,-1))) -%!assert(spdiag(as(:)),sparse(diag(af(:)))) -%!assert(spdiag(as(:),1),sparse(diag(af(:),1))) -%!assert(spdiag(as(:),-1),sparse(diag(af(:),-1))) -%!assert(spdiag(as(:)'),sparse(diag(af(:)'))) -%!assert(spdiag(as(:)',1),sparse(diag(af(:)',1))) -%!assert(spdiag(as(:)',-1),sparse(diag(af(:)',-1))) +%!assert(diag(as),sparse(diag(af))) +%!assert(diag(bs),sparse(diag(bf))) +%!assert(diag(as,1),sparse(diag(af,1))) +%!assert(diag(bs,1),sparse(diag(bf,1))) +%!assert(diag(as,-1),sparse(diag(af,-1))) +%!assert(diag(bs,-1),sparse(diag(bf,-1))) +%!assert(diag(as(:)),sparse(diag(af(:)))) +%!assert(diag(as(:),1),sparse(diag(af(:),1))) +%!assert(diag(as(:),-1),sparse(diag(af(:),-1))) +%!assert(diag(as(:)'),sparse(diag(af(:)'))) +%!assert(diag(as(:)',1),sparse(diag(af(:)',1))) +%!assert(diag(as(:)',-1),sparse(diag(af(:)',-1))) %!assert(spdiags(as,[0,1]),[diag(af,0),diag(af,1)]) %!test [tb,tc]=spdiags(as); %! assert(spdiags(tb,tc,sparse(zeros(size(as)))),as) @@ -531,21 +531,21 @@ %!assert(!issparse(af)) %!assert(!(issparse(af)&&iscomplex(af))) %!assert(!(issparse(af)&&isreal(af))) -%!assert(spsum(as),sparse(sum(af))) -%!assert(spsum(as,1),sparse(sum(af,1))) -%!assert(spsum(as,2),sparse(sum(af,2))) -%!assert(spcumsum(as),sparse(cumsum(af))) -%!assert(spcumsum(as,1),sparse(cumsum(af,1))) -%!assert(spcumsum(as,2),sparse(cumsum(af,2))) -%!assert(spsumsq(as),sparse(sumsq(af))) -%!assert(spsumsq(as,1),sparse(sumsq(af,1))) -%!assert(spsumsq(as,2),sparse(sumsq(af,2))) -%!assert(spprod(as),sparse(prod(af))) -%!assert(spprod(as,1),sparse(prod(af,1))) -%!assert(spprod(as,2),sparse(prod(af,2))) -%!assert(spcumprod(as),sparse(cumprod(af))) -%!assert(spcumprod(as,1),sparse(cumprod(af,1))) -%!assert(spcumprod(as,2),sparse(cumprod(af,2))) +%!assert(sum(as),sparse(sum(af))) +%!assert(sum(as,1),sparse(sum(af,1))) +%!assert(sum(as,2),sparse(sum(af,2))) +%!assert(cumsum(as),sparse(cumsum(af))) +%!assert(cumsum(as,1),sparse(cumsum(af,1))) +%!assert(cumsum(as,2),sparse(cumsum(af,2))) +%!assert(sumsq(as),sparse(sumsq(af))) +%!assert(sumsq(as,1),sparse(sumsq(af,1))) +%!assert(sumsq(as,2),sparse(sumsq(af,2))) +%!assert(prod(as),sparse(prod(af))) +%!assert(prod(as,1),sparse(prod(af,1))) +%!assert(prod(as,2),sparse(prod(af,2))) +%!assert(cumprod(as),sparse(cumprod(af))) +%!assert(cumprod(as,1),sparse(cumprod(af,1))) +%!assert(cumprod(as,2),sparse(cumprod(af,2))) %!assert(min(as),sparse(min(af))) %!assert(full(min(as(:))),min(af(:))) @@ -634,20 +634,20 @@ %! assert(det(bs+speye(size(bs))),det(bf+eye(size(bf))),100*eps*abs(det(bf+eye(size(bf))))) %!testif HAVE_UMFPACK -%! [l,u]=splu(sparse([1,1;1,1])); +%! [l,u]=lu(sparse([1,1;1,1])); %! assert(l*u,[1,1;1,1],10*eps); %!testif HAVE_UMFPACK -%! [l,u]=splu(sparse([1,1;1,1+i])); +%! [l,u]=lu(sparse([1,1;1,1+i])); %! assert(l,sparse([1,2,2],[1,1,2],1),10*eps); %! assert(u,sparse([1,1,2],[1,2,2],[1,1,1i]),10*eps); %!testif HAVE_UMFPACK ;# permuted LU -%! [L,U] = splu(bs); +%! [L,U] = lu(bs); %! assert(L*U,bs,1e-10); %!testif HAVE_UMFPACK ;# simple LU + row permutations -%! [L,U,P] = splu(bs); +%! [L,U,P] = lu(bs); %! assert(P'*L*U,bs,1e-10); %! # triangularity %! [i,j,v]=find(L); @@ -656,7 +656,7 @@ %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# simple LU + row/col permutations -%! [L,U,P,Q] = splu(bs); +%! [L,U,P,Q] = lu(bs); %! assert(P'*L*U*Q',bs,1e-10); %! # triangularity %! [i,j,v]=find(L); @@ -664,18 +664,18 @@ %! [i,j,v]=find(U); %! assert(j-i>=0); -%!testif HAVE_UMFPACK ;# LU with fixed column permutation -%! [L,U,P] = splu(bs,colamd(bs)); -%! assert(P'*L*U,bs,1e-10); +%!testif HAVE_UMFPACK ;# LU with vector permutations +%! [L,U,P] = lu(bs,'vector'); +%! assert(L(P,:)*U,bs,1e-10); %! # triangularity %! [i,j,v]=find(L); %! assert(i-j>=0); -%! [i,j,v]=find(U(:,colamd(bs))); +%! [i,j,v]=find(U); %! assert(j-i>=0); -%!testif HAVE_UMFPACK ;# LU with initial column permutation -%! [L,U,P,Q] = splu(bs,colamd(bs)); -%! assert(P'*L*U*Q',bs,1e-10); +%!testif HAVE_UMFPACK ;# LU with scaling +%! [L,U,P,Q,R] = lu(bs); +%! assert(R*P'*L*U*Q',bs,1e-10); %! # triangularity %! [i,j,v]=find(L); %! assert(i-j>=0); @@ -683,7 +683,7 @@ %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# inverse -%! assert(spinv(bs)*bs,sparse(eye(rows(bs))),1e-10); +%! assert(inv(bs)*bs,sparse(eye(rows(bs))),1e-10); %!assert(bf\as',bf\af',100*eps); %!assert(bs\af',bf\af',100*eps); @@ -696,25 +696,25 @@ gen_cholesky_tests() { cat >>$TESTS <<EOF %!testif HAVE_CHOLMOD -%! assert(spchol(bs)'*spchol(bs),bs,1e-10); +%! assert(chol(bs)'*chol(bs),bs,1e-10); %!testif HAVE_CHOLMOD -%! assert(splchol(bs)*splchol(bs)',bs,1e-10); +%! assert(chol(bs,'lower')*chol(bs,'lower')',bs,1e-10); %!testif HAVE_CHOLMOD -%! assert(splchol(bs),spchol(bs)',1e-10); +%! assert(chol(bs,'lower'),chol(bs)',1e-10); %!testif HAVE_CHOLMOD ;# Return Partial Cholesky factorization -%! [RS,PS] = spchol(bs); +%! [RS,PS] = chol(bs); %! assert(RS'*RS,bs,1e-10); %! assert(PS,0); -%! [LS,PS] = splchol(bs); +%! [LS,PS] = chol(bs,'lower'); %! assert(LS*LS',bs,1e-10); %! assert(PS,0); %!testif HAVE_CHOLMOD ;# Permuted Cholesky factorization -%! [RS,PS,QS] = spchol(bs); +%! [RS,PS,QS] = chol(bs); %! assert(RS'*RS,QS*bs*QS',1e-10); %! assert(PS,0); -%! [LS,PS,QS] = splchol(bs); +%! [LS,PS,QS] = chol(bs,'lower'); %! assert(LS*LS',QS*bs*QS',1e-10); %! assert(PS,0);
--- a/test/test_linalg.m Fri Feb 22 13:47:38 2008 -0500 +++ b/test/test_linalg.m Fri Feb 22 15:50:51 2008 +0100 @@ -106,7 +106,7 @@ %!error <Invalid call to chol.*> chol (); %% test/octave.test/linalg/chol-5.m -%!error <Invalid call to chol.*> chol (1, 2); +%!error <unexpected second or third input.*> chol (1, 2); %% test/octave.test/linalg/hess-1.m %!test @@ -124,7 +124,7 @@ %!error hess ([1, 2; 3, 4; 5, 6]); %% test/octave.test/linalg/lu-1.m -%!assert(all (all (lu ([1, 2; 3, 4]) - [1/3, 1; 1, 0] < eps))); +%!assert(all (all (lu ([1, 2; 3, 4]) - [3, 4; 1/3, 2/3] < eps))); %% test/octave.test/linalg/lu-2.m %!test @@ -139,11 +139,17 @@ %! && abs (u - [3, 4; 0, 2/3]) < sqrt (eps) %! && abs (p - [0, 1; 1, 0]) < sqrt (eps))); +%!test +%! [l, u, p] = lu ([1, 2; 3, 4],'vector'); +%! assert((abs (l - [1, 0; 1/3, 1]) < sqrt (eps) +%! && abs (u - [3, 4; 0, 2/3]) < sqrt (eps) +%! && abs (p - [2;1]) < sqrt (eps))); + %% test/octave.test/linalg/lu-4.m %!error <Invalid call to lu.*> lu (); %% test/octave.test/linalg/lu-5.m -%!error <Invalid call to lu.*> lu ([1, 2; 3, 4], 2); +%!error lu ([1, 2; 3, 4], 2); %% test/octave.test/linalg/lu-6.m %!test