Mercurial > octave-nkf
changeset 5610:9761b7d24e9e
[project @ 2006-02-09 09:12:02 by dbateman]
line wrap: on
line diff
--- a/ChangeLog Thu Feb 09 02:30:00 2006 +0000 +++ b/ChangeLog Thu Feb 09 09:12:03 2006 +0000 @@ -1,3 +1,9 @@ +2006-02-09 David Bateman <dbateman@free.fr> + + * configure.in: Fix for probe of colamd, cccolamd and metis. New test + for the presence of cxsparse. + Makeconf.in: Include CXSPARSE_LIBS. + 2006-01-19 John W. Eaton <jwe@octave.org> * configure.in: Use $WITH_PCRE instead of $HAVE_PCRE in shell test.
--- a/Makeconf.in Thu Feb 09 02:30:00 2006 +0000 +++ b/Makeconf.in Thu Feb 09 09:12:03 2006 +0000 @@ -200,6 +200,7 @@ COLAMD_LIBS = @COLAMD_LIBS@ CCOLAMD_LIBS = @CCOLAMD_LIBS@ CHOLMOD_LIBS = @CHOLMOD_LIBS@ +CXSPARSE_LIBS = @CXSPARSE_LIBS@ LIBS = @LIBS@ USE_64_BIT_IDX_T = @USE_64_BIT_IDX_T@
--- a/PROJECTS Thu Feb 09 02:30:00 2006 +0000 +++ b/PROJECTS Thu Feb 09 09:12:03 2006 +0000 @@ -107,14 +107,11 @@ Sparse Matrices: --------------- - * QR factorization functions, also for use in lssolve functions. Write - svds function based on this. Write sprank function based on svds. + * Improve QR factorization functions, using idea based on CSPARSE + cs_dmsol.m - * Once dmperm is implemented, perhaps use the technique to detect permuted - triangular matrices with permutations in both rows and cols. - - * Implement fourth argument to the sprand and sprandn that the leading - brand implements. + * Implement fourth argument to the sprand and sprandn, and addition + arguments to sprandsym that the leading brand implements. * Sparse logical indexing in idx_vector class so that something like "a=sprandn(1e6,1e6,1e-6); a(a<1) = 0" won't cause a memory overflow. @@ -122,13 +119,13 @@ * Write the rest of the sparse docs * The algo in TOMS 582 is perfect for symrcm function. However, this is - under the ACM license and can't be used in a GPL program. + under the ACM license and can't be used in a GPL program. An alternative is that PETSC is GPL compatiable and has a symrcm implemented from the original SPARSPAK. Its not clear that this is legal to me as I have found no clarification of the original license of SPARSPAK. As PETSC has had this code for over 10 years without - problem, we can perhaps assume that there is no issues. Maybe need + problem, we can perhaps assume that there are no issues. Maybe need to contact PETSC people or the SPARSPAK people at uni of waterloo to check issues. @@ -155,18 +152,10 @@ way to treat this but a complete rewrite of the sparse assignment functions... - * Port the sparse testing code into the DejaGNU testing code. Add tests - for spkron, matrix_type, luinc.. - * Other missing Functions - eigs - - dmperm Use Tim Davis' BTF code when available - symmmd Superseded by symamd - colmmd Superseded by colamd - - sprandsym - - etreeplot - - treeplot - - gplot - treelayout - cholinc - condest
--- a/configure.in Thu Feb 09 02:30:00 2006 +0000 +++ b/configure.in Thu Feb 09 09:12:03 2006 +0000 @@ -29,7 +29,7 @@ EXTERN_CXXFLAGS="$CXXFLAGS" AC_INIT -AC_REVISION($Revision: 1.498 $) +AC_REVISION($Revision: 1.499 $) AC_PREREQ(2.57) AC_CONFIG_SRCDIR([src/octave.cc]) AC_CONFIG_HEADER(config.h) @@ -808,7 +808,7 @@ if test "$with_colamd" = yes; then warn_colamd="COLAMD not found. This will result in some lack of functionality for sparse matrices." with_colamd=no - AC_CHECK_HEADERS([ufsparse/colamd.h umfpack/colamd.h colamd.h], [ + AC_CHECK_HEADERS([ufsparse/colamd.h colamd/colamd.h colamd.h], [ AC_CHECK_LIB(colamd, colamd, [COLAMD_LIBS="-lcolamd"; with_colamd=yes]) if test "$with_colamd" = yes; then AC_DEFINE(HAVE_COLAMD, 1, [Define if the COLAMD library is used.]) @@ -831,7 +831,7 @@ if test "$with_ccolamd" = yes; then warn_ccolamd="CCOLAMD not found. This will result in some lack of functionality for sparse matrices." with_ccolamd=no - AC_CHECK_HEADERS([ufsparse/ccolamd.h umfpack/ccolamd.h ccolamd.h], [ + AC_CHECK_HEADERS([ufsparse/ccolamd.h ccolamd/ccolamd.h ccolamd.h], [ AC_CHECK_LIB(ccolamd, ccolamd, [CCOLAMD_LIBS="-lccolamd"; with_ccolamd=yes]) if test "$with_ccolamd" = yes; then AC_DEFINE(HAVE_CCOLAMD, 1, [Define if the CCOLAMD library is used.]) @@ -855,8 +855,8 @@ test "$with_ccolamd" = yes && test "$with_amd" = yes; then warn_cholmod="CHOLMOD not found. This will result in some lack of functionality for sparse matrices." with_cholmod=no - AC_CHECK_HEADERS([ufsparse/cholmod.h umfpack/cholmod.h cholmod.h], [ - AC_CHECK_HEADERS([metis/metis.h ufsparse/metis.h umfpack/metis.h metis.h], [ + AC_CHECK_HEADERS([ufsparse/cholmod.h cholmod/cholmod.h cholmod.h], [ + AC_CHECK_HEADERS([ufsparse/metis.h metis/metis.h metis.h], [ AC_CHECK_LIB(metis, METIS_NodeND, with_metis=yes, with_metis=no) break], with_metis=no) @@ -865,14 +865,14 @@ AC_DEFINE(HAVE_METIS, 1, [Define if the METIS library is used.]) AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod -lmetis"; with_cholmod=yes], [ - AC_CHECK_LIB(cholmod_start, cholmod, + AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod -cblas -lmetis"; with_cholmod=yes], [], $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis)], $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS -lmetis) else AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod"; with_cholmod=yes], [ - AC_CHECK_LIB(cholmod_start, cholmod, [CHOLMOD_LIBS="-lcholmod -cblas"; + AC_CHECK_LIB(cholmod, cholmod_start, [CHOLMOD_LIBS="-lcholmod -cblas"; with_cholmod=yes], [], $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS)], $AMD_LIBS $COLAMD_LIBS $CCOLAMD_LIBS $BLAS_LIBS $FLIBS) @@ -889,6 +889,29 @@ AC_MSG_WARN($warn_cholmod) fi +CXSPARSE_LIBS= +AC_SUBST(CXSPARSE_LIBS) + +AC_ARG_WITH(cxsparse, + [ --without-cxsparse don't use CXSparse, disable some sparse functionality], + with_cxsparse=$withval, with_cxsparse=yes) + +if test "$with_cxsparse" = yes; then + warn_cxsparse="CXSparse not found. This will result in some lack of functionality for sparse matrices." + with_cxsparse=no + AC_CHECK_HEADERS([ufsparse/cxs.h cxsparse/cxs.h cxs.h], [ + AC_CHECK_LIB(cxspack, cs_sqr_di, [CXSPARSE_LIBS="-lcxspack"; with_cxsparse=yes]) + if test "$with_cxsparse" = yes; then + AC_DEFINE(HAVE_CXSPARSE, 1, [Define if the CXSparse library is used.]) + warn_cxsparse= + fi + break]) +fi + +if test -n "$warn_cxsparse"; then + AC_MSG_WARN($warn_cxsparse) +fi + ### Handle shared library options. ### Enable creation of static libraries. @@ -1810,6 +1833,7 @@ COLAMD libraries: $COLAMD_LIBS CCOLAMD libraries: $CCOLAMD_LIBS CHOLMOD libraries: $CHOLMOD_LIBS + CXSPARSE libraries: $CXSPARSE_LIBS HDF5 libraries: $HDF5_LIBS MPI libraries: $MPI_LIBS LIBS: $LIBS
--- a/liboctave/CSparse.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/CSparse.cc Thu Feb 09 09:12:03 2006 +0000 @@ -43,6 +43,7 @@ #include "oct-sparse.h" #include "sparse-util.h" #include "SparseCmplxCHOL.h" +#include "SparseCmplxQR.h" #include "oct-sort.h" @@ -610,7 +611,7 @@ SparseComplexMatrix SparseComplexMatrix::dinverse (SparseType &mattyp, octave_idx_type& info, - double& rcond, const bool force, + double& rcond, const bool, const bool calccond) const { SparseComplexMatrix retval; @@ -664,7 +665,7 @@ SparseComplexMatrix SparseComplexMatrix::tinverse (SparseType &mattyp, octave_idx_type& info, - double& rcond, const bool force, + double& rcond, const bool, const bool calccond) const { SparseComplexMatrix retval; @@ -902,7 +903,7 @@ SparseComplexMatrix SparseComplexMatrix::inverse (SparseType& mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const + double& rcond, int, int calc_cond) const { int typ = mattype.type (false); SparseComplexMatrix ret; @@ -979,7 +980,7 @@ } ComplexDET -SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int calc_cond) const +SparseComplexMatrix::determinant (octave_idx_type& err, double& rcond, int) const { ComplexDET retval; #ifdef HAVE_UMFPACK @@ -4936,7 +4937,6 @@ Symbolic, &Numeric, control, info) ; UMFPACK_ZNAME (free_symbolic) (&Symbolic) ; -#ifdef HAVE_LSSOLVE rcond = Info (UMFPACK_RCOND); volatile double rcond_plus_one = rcond + 1.0; @@ -4955,9 +4955,7 @@ rcond); } - else -#endif - if (status < 0) + else if (status < 0) { (*current_liboctave_error_handler) ("SparseComplexMatrix::solve numeric factorization failed"); @@ -5117,9 +5115,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_dense *X; @@ -5216,25 +5212,6 @@ } } -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); @@ -5396,9 +5373,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_sparse *X; @@ -5530,25 +5505,6 @@ retval.maybe_compress (); -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); @@ -5700,9 +5656,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_dense *X; @@ -5777,25 +5731,6 @@ } } -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseComplexMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_ZNAME (report_info) (control, info); UMFPACK_ZNAME (free_numeric) (&Numeric); @@ -5957,9 +5892,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_sparse *X; @@ -6070,7 +6003,6 @@ retval.maybe_compress (); -#ifndef HAVE_LSSOLVE rcond = Info (UMFPACK_RCOND); volatile double rcond_plus_one = rcond + 1.0; @@ -6087,7 +6019,6 @@ rcond); } -#endif UMFPACK_ZNAME (report_info) (control, info); @@ -6580,12 +6511,9 @@ } ComplexMatrix -SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return ComplexMatrix (); +SparseComplexMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type&) const +{ + return qrsolve (*this, b, info); } SparseComplexMatrix @@ -6605,12 +6533,9 @@ SparseComplexMatrix SparseComplexMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return SparseComplexMatrix (); + octave_idx_type&) const +{ + return qrsolve (*this, b, info); } ComplexMatrix @@ -6630,12 +6555,9 @@ ComplexMatrix SparseComplexMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return ComplexMatrix (); + octave_idx_type&) const +{ + return qrsolve (*this, b, info); } SparseComplexMatrix @@ -6655,12 +6577,9 @@ SparseComplexMatrix SparseComplexMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return SparseComplexMatrix (); + octave_idx_type&) const +{ + return qrsolve (*this, b, info); } ComplexColumnVector @@ -6681,10 +6600,8 @@ ComplexColumnVector SparseComplexMatrix::lssolve (const ColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return ComplexColumnVector (); + Matrix tmp (b); + return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0)); } ComplexColumnVector @@ -6706,10 +6623,8 @@ SparseComplexMatrix::lssolve (const ComplexColumnVector& b, octave_idx_type& info, octave_idx_type& rank) const { - info = -1; - (*current_liboctave_error_handler) - ("SparseComplexMatrix::lssolve not implemented yet"); - return ComplexColumnVector (); + ComplexMatrix tmp (b); + return lssolve (tmp, info, rank).column (static_cast<octave_idx_type> (0)); } // unary operations
--- a/liboctave/ChangeLog Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/ChangeLog Thu Feb 09 09:12:03 2006 +0000 @@ -1,3 +1,20 @@ +2006-02-09 David Bateman <dbateman@free.fr> + + * SparseQR.cc: new file for real sparse QR class. + * SparseQR.h: declaration. + * SparseCmplxQR.cc: new file for complex sparse QR class. + * SparseCmplxQR.h: declaration. + * dSparse.cc (dinverse,tinverse,inverse): Remove unused input args. + (factorize, fsolve): Enable code code lssolve. + (lssolve): disable unused args, write based in above sparse QR class. + * CSparse.cc (dinverse,tinverse,inverse): Remove unused input args. + (factorize, fsolve): Enable code code lssolve. + (lssolve): disable unused args, write based in above sparse QR class. + * oct-sparse.h: fix location of colamd, ccolamd and metis headers. + Include CXSparse headers. + * Makefile.in (MATRIX_INC): Include SparseQR.h and SparseCmplxQR.h. + (MATRIX_SRC): Include SparseQR.cc and SparseCmplxQR.cc. + 2006-02-08 John W. Eaton <jwe@octave.org> * Array-util.h (calc_permutated_idx): Delete.
--- a/liboctave/Makefile.in Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/Makefile.in Thu Feb 09 09:12:03 2006 +0000 @@ -37,7 +37,7 @@ boolSparse.h CSparse.h dSparse.h MSparse-defs.h MSparse.h \ Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \ sparse-base-chol.h SparseCmplxCHOL.h SparsedbleCHOL.h \ - Sparse-op-defs.h SparseType.h \ + SparseCmplxQR.h SparseQR.h Sparse-op-defs.h SparseType.h \ int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \ int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \ intNDArray.h @@ -96,7 +96,8 @@ dbleDET.cc dbleHESS.cc dbleLU.cc dbleQR.cc dbleQRP.cc \ dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \ MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \ - SparseCmplxCHOL.cc SparsedbleCHOL.cc SparseType.cc \ + SparseCmplxCHOL.cc SparsedbleCHOL.cc \ + SparseCmplxQR.cc SparseQR.cc SparseType.cc \ int8NDArray.cc uint8NDArray.cc int16NDArray.cc uint16NDArray.cc \ int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/SparseCmplxQR.cc Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,676 @@ +/* + +Copyright (C) 2005 David Bateman + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif +#include <vector> + +#include "lo-error.h" +#include "SparseCmplxQR.h" + +SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep +(const SparseComplexMatrix& a, int order) +{ +#ifdef HAVE_CXSPARSE + // cast away const on A, with full knowledge that CSparse won't touch it + CXSPARSE_ZNAME (cs) A; + A.nzmax = a.nnz (); + A.m = a.rows (); + A.n = a.cols (); + nrows = A.m; + // Cast away const on A, with full knowledge that CSparse won't touch it + // Prevents the methods below making a copy of the data. + A.p = const_cast<octave_idx_type *>(a.cidx ()); + A.i = const_cast<octave_idx_type *>(a.ridx ()); + A.x = const_cast<double _Complex *>(reinterpret_cast<const double _Complex *> + (a.data ())); + A.nz = -1; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + S = CXSPARSE_ZNAME (cs_sqr) (&A, order, 1); + N = CXSPARSE_ZNAME (cs_qr) (&A, S); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + if (!N) + (*current_liboctave_error_handler) + ("SparseComplexQR: sparse matrix QR factorization filled"); + count = 1; +#else + (*current_liboctave_error_handler) + ("SparseComplexQR: sparse matrix QR factorization not implemented"); +#endif +} + +SparseComplexQR::SparseComplexQR_rep::~SparseComplexQR_rep (void) +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_ZNAME (cs_sfree) (S); + CXSPARSE_ZNAME (cs_nfree) (N); +#endif +} + +SparseComplexMatrix +SparseComplexQR::SparseComplexQR_rep::V (void) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from V and sort + // XXX FIXME XXX Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_dropzeros) (N->L); + CXSPARSE_ZNAME (cs) *D = CXSPARSE_ZNAME (cs_transpose) (N->L, 1); + CXSPARSE_ZNAME (cs_spfree) (N->L); + N->L = CXSPARSE_ZNAME (cs_transpose) (D, 1); + CXSPARSE_ZNAME (cs_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->L->n; + octave_idx_type nz = N->L->nzmax; + SparseComplexMatrix ret (N->L->m, nc, nz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->L->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->L->i[j]; + ret.xdata (j) = reinterpret_cast<Complex *>(N->L->x)[j]; + } + return ret; +#else + return SparseComplexMatrix (); +#endif +} + +ColumnVector +SparseComplexQR::SparseComplexQR_rep::Pinv (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret(N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem(i) = S->Pinv[i]; + return ret; +#else + return ColumnVector (); +#endif +} + +ColumnVector +SparseComplexQR::SparseComplexQR_rep::P (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret(N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem(S->Pinv[i]) = i; + return ret; +#else + return ColumnVector (); +#endif +} + +SparseComplexMatrix +SparseComplexQR::SparseComplexQR_rep::R (const bool econ) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from R and sort + // XXX FIXME XXX Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_dropzeros) (N->U); + CXSPARSE_ZNAME (cs) *D = CXSPARSE_ZNAME (cs_transpose) (N->U, 1); + CXSPARSE_ZNAME (cs_spfree) (N->U); + N->U = CXSPARSE_ZNAME (cs_transpose) (D, 1); + CXSPARSE_ZNAME (cs_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->U->n; + octave_idx_type nz = N->U->nzmax; + SparseComplexMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->U->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->U->i[j]; + ret.xdata (j) = reinterpret_cast<Complex *>(N->U->x)[j]; + } + return ret; +#else + return SparseComplexMatrix (); +#endif +} + +ComplexMatrix +SparseComplexQR::SparseComplexQR_rep::C (const ComplexMatrix &b) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type b_nr = b.rows(); + octave_idx_type b_nc = b.cols(); + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + const double _Complex *bvec = + reinterpret_cast<const double _Complex *>(b.fortran_vec()); + ComplexMatrix ret(b_nr,b_nc); + Complex *vec = ret.fortran_vec(); + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + else + { + OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) + { + OCTAVE_QUIT; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_ipvec) (b_nr, S->Pinv, bvec + idx, + reinterpret_cast<double _Complex *>(buf)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type i = 0; i < nm; i++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) + (N->L, i, N->B[i], reinterpret_cast<double _Complex *>(buf)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < b_nr; i++) + vec[i+idx] = buf[i]; + } + } + return ret; +#else + return ComplexMatrix (); +#endif +} + +ComplexMatrix +qrsolve(const SparseComplexMatrix&a, const Matrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + ComplexMatrix x; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseComplexQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + double _Complex *vec = reinterpret_cast<double _Complex *> + (x.fortran_vec()); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, q.S()->m2); + OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_ipvec) + (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_ZNAME (cs_ipvec) (nc, q.S()->Q, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + else + { + SparseComplexMatrix at = a.hermitian(); + SparseComplexQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + double _Complex *vec = reinterpret_cast<double _Complex *> + (x.fortran_vec()); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, + nc > q.S()->m2 ? nc : q.S()->m2); + OCTAVE_LOCAL_BUFFER (Complex, Xx, b_nr); + OCTAVE_LOCAL_BUFFER (Complex, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) + (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); + CXSPARSE_ZNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + CXSPARSE_ZNAME (cs_happly) + (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) (nc, q.S()->Pinv, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + + return x; +#else + return ComplexMatrix (); +#endif +} + +SparseComplexMatrix +qrsolve(const SparseComplexMatrix&a, const SparseMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + SparseComplexMatrix x; + volatile octave_idx_type ii, x_nz; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseComplexQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_ipvec) + (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_ZNAME (cs_ipvec) (nc, q.S()->Q, buf, + reinterpret_cast<double _Complex *>(Xx)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + else + { + SparseComplexMatrix at = a.hermitian(); + SparseComplexQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, + nc > q.S()->m2 ? nc : q.S()->m2); + OCTAVE_LOCAL_BUFFER (Complex, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) + (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); + CXSPARSE_ZNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) + (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) (nc, q.S()->Pinv, buf, + reinterpret_cast<double _Complex *>(Xx)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + + x.maybe_compress (); + return x; +#else + return SparseComplexMatrix (); +#endif +} + +ComplexMatrix +qrsolve(const SparseComplexMatrix&a, const ComplexMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + const double _Complex *bvec = + reinterpret_cast<const double _Complex *>(b.fortran_vec()); + ComplexMatrix x; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseComplexQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + double _Complex *vec = reinterpret_cast<double _Complex *> + (x.fortran_vec()); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_ZNAME (cs_ipvec) (nc, q.S()->Q, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + else + { + SparseComplexMatrix at = a.hermitian(); + SparseComplexQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + double _Complex *vec = reinterpret_cast<double _Complex *> + (x.fortran_vec()); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, + nc > q.S()->m2 ? nc : q.S()->m2); + OCTAVE_LOCAL_BUFFER (Complex, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) (nr, q.S()->Q, bvec + bidx, buf); + CXSPARSE_ZNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) + (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) (nc, q.S()->Pinv, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + + return x; +#else + return ComplexMatrix (); +#endif +} + +SparseComplexMatrix +qrsolve(const SparseComplexMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + SparseComplexMatrix x; + volatile octave_idx_type ii, x_nz; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseComplexQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_ipvec) + (nr, q.S()->Pinv, reinterpret_cast<double _Complex *>(Xx), buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_ZNAME (cs_ipvec) (nc, q.S()->Q, buf, + reinterpret_cast<double _Complex *>(Xx)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + else + { + SparseComplexMatrix at = a.hermitian(); + SparseComplexQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (Complex, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double _Complex, buf, + nc > q.S()->m2 ? nc : q.S()->m2); + OCTAVE_LOCAL_BUFFER (Complex, B, nr); + for (octave_idx_type i = 0; i < nr; i++) + B[i] = conj (reinterpret_cast<Complex *>(q.N()->B) [i]); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) + (nr, q.S()->Q, reinterpret_cast<double _Complex *>(Xx), buf); + CXSPARSE_ZNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_happly) + (q.N()->L, j, reinterpret_cast<double _Complex *>(B)[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_ZNAME (cs_pvec) (nc, q.S()->Pinv, buf, + reinterpret_cast<double _Complex *>(Xx)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + + x.maybe_compress (); + return x; +#else + return SparseComplexMatrix (); +#endif +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/ +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/SparseCmplxQR.h Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,146 @@ +/* + +Copyright (C) 2005 David Bateman + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (sparse_cmplx_QR_h) +#define sparse_cmplx_QR_h 1 + +#include <iostream> + +#include "dMatrix.h" +#include "CMatrix.h" +#include "dSparse.h" +#include "CSparse.h" +#include "oct-sparse.h" + +#ifdef IDX_TYPE_LONG +#define CXSPARSE_ZNAME(name) name ## _cl +#else +#define CXSPARSE_ZNAME(name) name ## _ci +#endif + +class +SparseComplexQR +{ +protected: + class SparseComplexQR_rep + { + public: + SparseComplexQR_rep (const SparseComplexMatrix& a, int order); + + ~SparseComplexQR_rep (void); +#ifdef HAVE_CXSPARSE + bool ok (void) const { return (N && S); } +#else + bool ok (void) const { return false; } +#endif + SparseComplexMatrix V (void) const; + + ColumnVector Pinv (void) const; + + ColumnVector P (void) const; + + SparseComplexMatrix R (const bool econ) const; + + ComplexMatrix C (const ComplexMatrix &b) const; + + int count; + + octave_idx_type nrows; +#ifdef HAVE_CXSPARSE + CXSPARSE_ZNAME (css) *S; + + CXSPARSE_ZNAME (csn) *N; +#endif + }; +private: + SparseComplexQR_rep *rep; + +public: + SparseComplexQR (void) : + rep (new SparseComplexQR_rep (SparseComplexMatrix(), -1)) { } + + SparseComplexQR (const SparseComplexMatrix& a, int order = -1) : + rep (new SparseComplexQR_rep (a, order)) { } + + SparseComplexQR (const SparseComplexQR& a) : rep (a.rep) { rep->count++; } + + ~SparseComplexQR (void) + { + if (--rep->count <= 0) + delete rep; + } + + SparseComplexQR& operator = (const SparseComplexQR& a) + { + if (this != &a) + { + if (--rep->count <= 0) + delete rep; + + rep = a.rep; + rep->count++; + } + return *this; + } + + bool ok (void) const { return rep->ok(); } + + SparseComplexMatrix V (void) const { return rep->V(); } + + ColumnVector Pinv (void) const { return rep->P(); } + + ColumnVector P (void) const { return rep->P(); } + + SparseComplexMatrix R (const bool econ = false) const + { return rep->R(econ); } + + ComplexMatrix C (const ComplexMatrix &b) const { return rep->C(b); } + + friend ComplexMatrix qrsolve (const SparseComplexMatrix &a, const Matrix &b, + octave_idx_type &info); + + friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, + const SparseMatrix &b, + octave_idx_type &info); + + friend ComplexMatrix qrsolve (const SparseComplexMatrix &a, + const ComplexMatrix &b, + octave_idx_type &info); + + friend SparseComplexMatrix qrsolve (const SparseComplexMatrix &a, + const SparseComplexMatrix &b, + octave_idx_type &info); + +protected: +#ifdef HAVE_CXSPARSE + CXSPARSE_ZNAME (css) * S (void) { return rep->S; } + + CXSPARSE_ZNAME (csn) * N (void) { return rep->N; } +#endif +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/SparseQR.cc Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,715 @@ +/* + +Copyright (C) 2005 David Bateman + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif +#include <vector> + +#include "lo-error.h" +#include "SparseQR.h" + +SparseQR::SparseQR_rep::SparseQR_rep (const SparseMatrix& a, int order) +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME (cs) A; + A.nzmax = a.nzmax (); + A.m = a.rows (); + A.n = a.cols (); + nrows = A.m; + // Cast away const on A, with full knowledge that CSparse won't touch it + // Prevents the methods below making a copy of the data. + A.p = const_cast<octave_idx_type *>(a.cidx ()); + A.i = const_cast<octave_idx_type *>(a.ridx ()); + A.x = const_cast<double *>(a.data ()); + A.nz = -1; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + S = CXSPARSE_DNAME (cs_sqr) (&A, order, 1); + N = CXSPARSE_DNAME (cs_qr) (&A, S); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + if (!N) + (*current_liboctave_error_handler) + ("SparseQR: sparse matrix QR factorization filled"); + count = 1; +#else + (*current_liboctave_error_handler) + ("SparseQR: sparse matrix QR factorization not implemented"); +#endif +} + +SparseQR::SparseQR_rep::~SparseQR_rep (void) +{ +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME (cs_sfree) (S); + CXSPARSE_DNAME (cs_nfree) (N); +#endif +} + +SparseMatrix +SparseQR::SparseQR_rep::V (void) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from V and sort + // XXX FIXME XXX Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_dropzeros) (N->L); + CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->L, 1); + CXSPARSE_DNAME (cs_spfree) (N->L); + N->L = CXSPARSE_DNAME (cs_transpose) (D, 1); + CXSPARSE_DNAME (cs_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->L->n; + octave_idx_type nz = N->L->nzmax; + SparseMatrix ret (N->L->m, nc, nz); + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->L->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->L->i[j]; + ret.xdata (j) = N->L->x[j]; + } + return ret; +#else + return SparseMatrix (); +#endif +} + +ColumnVector +SparseQR::SparseQR_rep::Pinv (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret(N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem(i) = S->Pinv[i]; + return ret; +#else + return ColumnVector (); +#endif +} + +ColumnVector +SparseQR::SparseQR_rep::P (void) const +{ +#ifdef HAVE_CXSPARSE + ColumnVector ret(N->L->m); + for (octave_idx_type i = 0; i < N->L->m; i++) + ret.xelem(S->Pinv[i]) = i; + return ret; +#else + return ColumnVector (); +#endif +} + +SparseMatrix +SparseQR::SparseQR_rep::R (const bool econ) const +{ +#ifdef HAVE_CXSPARSE + // Drop zeros from R and sort + // XXX FIXME XXX Is the double transpose to sort necessary? + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_dropzeros) (N->U); + CXSPARSE_DNAME (cs) *D = CXSPARSE_DNAME (cs_transpose) (N->U, 1); + CXSPARSE_DNAME (cs_spfree) (N->U); + N->U = CXSPARSE_DNAME (cs_transpose) (D, 1); + CXSPARSE_DNAME (cs_spfree) (D); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + octave_idx_type nc = N->U->n; + octave_idx_type nz = N->U->nzmax; + + SparseMatrix ret ((econ ? (nc > nrows ? nrows : nc) : nrows), nc, nz); + + for (octave_idx_type j = 0; j < nc+1; j++) + ret.xcidx (j) = N->U->p[j]; + for (octave_idx_type j = 0; j < nz; j++) + { + ret.xridx (j) = N->U->i[j]; + ret.xdata (j) = N->U->x[j]; + } + return ret; +#else + return SparseMatrix (); +#endif +} + +Matrix +SparseQR::SparseQR_rep::C (const Matrix &b) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type b_nr = b.rows(); + octave_idx_type b_nc = b.cols(); + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + const double *bvec = b.fortran_vec(); + Matrix ret(b_nr,b_nc); + double *vec = ret.fortran_vec(); + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + else + { + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < b_nc; j++, idx+=b_nr) + { + OCTAVE_QUIT; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + // cast away const on bvec, with full knowledge that CSparse + // won't touch it + CXSPARSE_DNAME (cs_ipvec) (b_nr, S->Pinv, bvec + idx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (volatile octave_idx_type i = 0; i < nm; i++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (N->L, i, N->B[i], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < b_nr; i++) + vec[i+idx] = buf[i]; + } + } + return ret; +#else + return Matrix (); +#endif +} + +Matrix +qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + const double *bvec = b.fortran_vec(); + Matrix x; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return Matrix(); + } + x.resize(nc, b_nc); + double *vec = x.fortran_vec(); + OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + // cast away const on bvec, with full knowledge that CSparse + // won't touch it + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, bvec + bidx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + else + { + SparseMatrix at = a.hermitian(); + SparseQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return Matrix(); + } + x.resize(nc, b_nc); + double *vec = x.fortran_vec(); + OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0, bidx = 0; i < b_nc; + i++, idx+=nc, bidx+=b_nr) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + // cast away const on bvec, with full knowledge that CSparse + // won't touch it + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, bvec + bidx, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, vec + idx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + } + + return x; +#else + return Matrix (); +#endif +} + +SparseMatrix +qrsolve(const SparseMatrix&a, const SparseMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nr = b.rows(); + octave_idx_type b_nc = b.cols(); + SparseMatrix x; + volatile octave_idx_type ii, x_nz; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return SparseMatrix(); + } + x = SparseMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + double tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + else + { + SparseMatrix at = a.hermitian(); + SparseQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return SparseMatrix(); + } + x = SparseMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + Xx[j] = b.xelem(j,i); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + double tmp = Xx[j]; + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + + x.maybe_compress (); + return x; +#else + return SparseMatrix (); +#endif +} + +ComplexMatrix +qrsolve(const SparseMatrix&a, const ComplexMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nc = b.cols(); + octave_idx_type b_nr = b.rows(); + ComplexMatrix x; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + Complex *vec = x.fortran_vec(); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); + + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = 0; j < nc; j++) + vec[j+idx] = Complex (Xx[j], Xz[j]); + } + } + else + { + SparseMatrix at = a.hermitian(); + SparseQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return ComplexMatrix(); + } + x.resize(nc, b_nc); + Complex *vec = x.fortran_vec(); + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (octave_idx_type j = 0; j < nc; j++) + vec[j+idx] = Complex (Xx[j], Xz[j]); + } + } + + return x; +#else + return ComplexMatrix (); +#endif +} + +SparseComplexMatrix +qrsolve(const SparseMatrix&a, const SparseComplexMatrix &b, octave_idx_type &info) +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nr = a.rows(); + octave_idx_type nc = a.cols(); + octave_idx_type b_nr = b.rows(); + octave_idx_type b_nc = b.cols(); + SparseComplexMatrix x; + volatile octave_idx_type ii, x_nz; + info = 0; + + if (nr < 1 || nc < 1 || nr != b_nr) + (*current_liboctave_error_handler) + ("matrix dimension mismatch in solution of minimum norm problem"); + else if (nr >= nc) + { + SparseQR q (a, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xx, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xx); + CXSPARSE_DNAME (cs_ipvec) (nr, q.S()->Pinv, Xz, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_usolve) (q.N()->U, buf); + CXSPARSE_DNAME (cs_ipvec) (nc, q.S()->Q, buf, Xz); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Complex (Xx[j], Xz[j]); + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + else + { + SparseMatrix at = a.hermitian(); + SparseQR q (at, 2); + if (! q.ok ()) + { + info = -1; + return SparseComplexMatrix(); + } + x = SparseComplexMatrix (nc, b_nc, b.nzmax()); + x.xcidx(0) = 0; + x_nz = b.nzmax(); + ii = 0; + OCTAVE_LOCAL_BUFFER (double, Xx, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, Xz, (b_nr > nc ? b_nr : nc)); + OCTAVE_LOCAL_BUFFER (double, buf, nc > q.S()->m2 ? nc : q.S()->m2); + for (volatile octave_idx_type i = 0, idx = 0; i < b_nc; i++, idx+=nc) + { + OCTAVE_QUIT; + for (octave_idx_type j = 0; j < b_nr; j++) + { + Complex c = b.xelem (j,i); + Xx[j] = std::real (c); + Xz[j] = std::imag (c); + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xx, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xx); + CXSPARSE_DNAME (cs_pvec) (nr, q.S()->Q, Xz, buf); + CXSPARSE_DNAME (cs_utsolve) (q.N()->U, buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + for (volatile octave_idx_type j = nr-1; j >= 0; j--) + { + OCTAVE_QUIT; + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_happly) (q.N()->L, j, q.N()->B[j], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + CXSPARSE_DNAME (cs_pvec) (nc, q.S()->Pinv, buf, Xz); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + + for (octave_idx_type j = 0; j < nc; j++) + { + Complex tmp = Complex (Xx[j], Xz[j]); + if (tmp != 0.0) + { + if (ii == x_nz) + { + // Resize the sparse matrix + octave_idx_type sz = x_nz * (b_nc - i) / b_nc; + sz = (sz > 10 ? sz : 10) + x_nz; + x.change_capacity (sz); + x_nz = sz; + } + x.xdata(ii) = tmp; + x.xridx(ii++) = j; + } + } + x.xcidx(i+1) = ii; + } + } + + x.maybe_compress (); + return x; +#else + return SparseComplexMatrix (); +#endif +} + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/liboctave/SparseQR.h Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,142 @@ +/* + +Copyright (C) 2005 David Bateman + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#if !defined (sparse_QR_h) +#define sparse_QR_h 1 + +#include <iostream> + +#include "dMatrix.h" +#include "CMatrix.h" +#include "dSparse.h" +#include "CSparse.h" +#include "oct-sparse.h" + +#ifdef IDX_TYPE_LONG +#define CXSPARSE_DNAME(name) name ## _dl +#else +#define CXSPARSE_DNAME(name) name ## _di +#endif + +class +SparseQR +{ +protected: + class SparseQR_rep + { + public: + SparseQR_rep (const SparseMatrix& a, int order); + + ~SparseQR_rep (void); +#ifdef HAVE_CXSPARSE + bool ok (void) const { return (N && S); } +#else + bool ok (void) const { return false; } +#endif + SparseMatrix V (void) const; + + ColumnVector Pinv (void) const; + + ColumnVector P (void) const; + + SparseMatrix R (const bool econ) const; + + Matrix C (const Matrix &b) const; + + int count; + + octave_idx_type nrows; +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME (css) *S; + + CXSPARSE_DNAME (csn) *N; +#endif + }; +private: + SparseQR_rep *rep; + +public: + SparseQR (void) : rep (new SparseQR_rep (SparseMatrix(), -1)) { } + + SparseQR (const SparseMatrix& a, int order = -1) : + rep (new SparseQR_rep (a, order)) { } + + SparseQR (const SparseQR& a) : rep (a.rep) { rep->count++; } + + ~SparseQR (void) + { + if (--rep->count <= 0) + delete rep; + } + + SparseQR& operator = (const SparseQR& a) + { + if (this != &a) + { + if (--rep->count <= 0) + delete rep; + + rep = a.rep; + rep->count++; + } + return *this; + } + + bool ok (void) const { return rep->ok(); } + + SparseMatrix V (void) const { return rep->V(); } + + ColumnVector Pinv (void) const { return rep->P(); } + + ColumnVector P (void) const { return rep->P(); } + + SparseMatrix R (const bool econ = false) const { return rep->R(econ); } + + Matrix C (const Matrix &b) const { return rep->C(b); } + + friend Matrix qrsolve (const SparseMatrix &a, const Matrix &b, + octave_idx_type &info); + + friend SparseMatrix qrsolve (const SparseMatrix &a, const SparseMatrix &b, + octave_idx_type &info); + + friend ComplexMatrix qrsolve (const SparseMatrix &a, const ComplexMatrix &b, + octave_idx_type &info); + + friend SparseComplexMatrix qrsolve (const SparseMatrix &a, + const SparseComplexMatrix &b, + octave_idx_type &info); + +protected: +#ifdef HAVE_CXSPARSE + CXSPARSE_DNAME (css) * S (void) { return rep->S; } + + CXSPARSE_DNAME (csn) * N (void) { return rep->N; } +#endif +}; + +#endif + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/liboctave/SparseType.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/SparseType.cc Thu Feb 09 09:12:03 2006 +0000 @@ -806,7 +806,7 @@ ("Hermitian/Symmetric Tridiagonal Sparse Matrix"); else if (typ == SparseType::Rectangular) (*current_liboctave_warning_handler) - ("Rectangular Sparse Matrix"); + ("Rectangular/Singular Sparse Matrix"); else if (typ == SparseType::Full) (*current_liboctave_warning_handler) ("Full Sparse Matrix");
--- a/liboctave/dSparse.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/dSparse.cc Thu Feb 09 09:12:03 2006 +0000 @@ -44,6 +44,7 @@ #include "oct-sparse.h" #include "sparse-util.h" #include "SparsedbleCHOL.h" +#include "SparseQR.h" #include "oct-sort.h" @@ -697,7 +698,7 @@ SparseMatrix SparseMatrix::dinverse (SparseType &mattyp, octave_idx_type& info, - double& rcond, const bool force, + double& rcond, const bool, const bool calccond) const { SparseMatrix retval; @@ -751,7 +752,7 @@ SparseMatrix SparseMatrix::tinverse (SparseType &mattyp, octave_idx_type& info, - double& rcond, const bool force, + double& rcond, const bool, const bool calccond) const { SparseMatrix retval; @@ -989,7 +990,7 @@ SparseMatrix SparseMatrix::inverse (SparseType &mattype, octave_idx_type& info, - double& rcond, int force, int calc_cond) const + double& rcond, int, int calc_cond) const { int typ = mattype.type (false); SparseMatrix ret; @@ -5154,7 +5155,6 @@ &Numeric, control, info) ; UMFPACK_DNAME (free_symbolic) (&Symbolic) ; -#ifdef HAVE_LSSOLVE rcond = Info (UMFPACK_RCOND); volatile double rcond_plus_one = rcond + 1.0; @@ -5173,9 +5173,7 @@ rcond); } - else -#endif - if (status < 0) + else if (status < 0) { (*current_liboctave_error_handler) ("SparseMatrix::solve numeric factorization failed"); @@ -5337,9 +5335,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_dense *X; @@ -5410,25 +5406,6 @@ } } -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_DNAME (report_info) (control, info); UMFPACK_DNAME (free_numeric) (&Numeric); @@ -5589,9 +5566,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_sparse *X; @@ -5699,25 +5674,6 @@ retval.maybe_compress (); -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_DNAME (report_info) (control, info); UMFPACK_DNAME (free_numeric) (&Numeric); @@ -5868,9 +5824,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_dense *X; @@ -5960,25 +5914,6 @@ retval (i, j) = Complex (Xx[i], Xz[i]); } -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_DNAME (report_info) (control, info); UMFPACK_DNAME (free_numeric) (&Numeric); @@ -6140,9 +6075,7 @@ ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", rcond); -#ifdef HAVE_LSSOLVE return retval; -#endif } cholmod_sparse *X; @@ -6261,25 +6194,6 @@ retval.maybe_compress (); -#ifndef HAVE_LSSOLVE - rcond = Info (UMFPACK_RCOND); - volatile double rcond_plus_one = rcond + 1.0; - - if (status == UMFPACK_WARNING_singular_matrix || - rcond_plus_one == 1.0 || xisnan (rcond)) - { - err = -2; - - if (sing_handler) - sing_handler (rcond); - else - (*current_liboctave_error_handler) - ("SparseMatrix::solve matrix singular to machine precision, rcond = %g", - rcond); - - } -#endif - UMFPACK_DNAME (report_info) (control, info); UMFPACK_DNAME (free_numeric) (&Numeric); @@ -6761,12 +6675,9 @@ } Matrix -SparseMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseMatrix::lssolve not implemented yet"); - return Matrix (); +SparseMatrix::lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type&) const +{ + return qrsolve (*this, b, info); } SparseMatrix @@ -6785,12 +6696,9 @@ } SparseMatrix -SparseMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info, octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseMatrix::lssolve not implemented yet"); - return SparseMatrix (); +SparseMatrix::lssolve (const SparseMatrix& b, octave_idx_type& info, octave_idx_type&) const +{ + return qrsolve (*this, b, info); } ComplexMatrix @@ -6809,12 +6717,9 @@ } ComplexMatrix -SparseMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseMatrix::lssolve not implemented yet"); - return ComplexMatrix (); +SparseMatrix::lssolve (const ComplexMatrix& b, octave_idx_type& info, octave_idx_type&) const +{ + return qrsolve (*this, b, info); } SparseComplexMatrix @@ -6834,12 +6739,9 @@ SparseComplexMatrix SparseMatrix::lssolve (const SparseComplexMatrix& b, octave_idx_type& info, - octave_idx_type& rank) const -{ - info = -1; - (*current_liboctave_error_handler) - ("SparseMatrix::lssolve not implemented yet"); - return SparseComplexMatrix (); + octave_idx_type&) const +{ + return qrsolve (*this, b, info); } ColumnVector
--- a/liboctave/oct-sparse.h Thu Feb 09 02:30:00 2006 +0000 +++ b/liboctave/oct-sparse.h Thu Feb 09 09:12:03 2006 +0000 @@ -42,26 +42,24 @@ #if defined (HAVE_UFSPARSE_COLAMD_H) #include <ufsparse/colamd.h> -#elif defined (HAVE_UMFPACK_COLAMD_H) -#include <umfpack/colamd.h> +#elif defined (HAVE_COLAMD_COLAMD_H) +#include <colamd/colamd.h> #elif defined (HAVE_COLAMD_H) #include <colamd.h> #endif #if defined (HAVE_UFSPARSE_CCOLAMD_H) #include <ufsparse/ccolamd.h> -#elif defined (HAVE_UMFPACK_CCOLAMD_H) -#include <umfpack/ccolamd.h> +#elif defined (HAVE_CCOLAMD_CCOLAMD_H) +#include <ccolamd/ccolamd.h> #elif defined (HAVE_CCOLAMD_H) #include <ccolamd.h> #endif -#if defined (HAVE_METIS_METIS_H) +#if defined (HAVE_UFSPARSE_METIS_H) +#include <ufsparse/metis.h> +#elif defined (HAVE_METIS_METIS_H) #include <metis/metis.h> -#elif defined (HAVE_UFSPARSE_METIS_H) -#include <ufsparse/metis.h> -#elif defined (HAVE_UMFPACK_METIS_H) -#include <umfpack/metis.h> #elif defined (HAVE_METIS_H) #include <metis.h> #endif @@ -69,11 +67,19 @@ #if defined (HAVE_UFSPARSE_CHOLMOD_H) #include <ufsparse/cholmod.h> #elif defined (HAVE_UMFPACK_CHOLMOD_H) -#include <umfpack/cholmod.h> +#include <cholmod/cholmod.h> #elif defined (HAVE_CHOLMOD_H) #include <cholmod.h> #endif +#if defined (HAVE_UFSPARSE_CXS_H) +#include <ufsparse/cxs.h> +#elif defined (HAVE_CXSPARSE_CXS_H) +#include <cxsparse/cxs.h> +#elif defined (HAVE_CXS_H) +#include <cxs.h> +#endif + #if (defined (HAVE_UFSPARSE_CHOLMOD_H) \ || defined (HAVE_UMFPACK_CHOLMOD_H) \ || defined (HAVE_CHOLMOD_H))
--- a/scripts/ChangeLog Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/ChangeLog Thu Feb 09 09:12:03 2006 +0000 @@ -1,3 +1,15 @@ +2006-02-09 David Bateman <dbateman@free.fr> + + * general/triu.m: Minimum change to allow sparse matrix. More needed + for arbitrary user type. + * general/tril.m: ditto. + * sparse/sprand.m : Doc fix. + * sparse/sprandn.m: Ditto. + * sparse/sprandsym.m: New function. + * audio/setaudio.m, general/cart2pol.m, general/cart2sph.m, + general/pol2cart.m, general/sph2cart.m, signal/freqz_plot.m: + Update for syntax error for latest texinfo.tex file. + 2006-02-02 John W. Eaton <jwe@octave.org> * plot/grid.m: Append ";\n" to "set grid" command.
--- a/scripts/audio/setaudio.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/audio/setaudio.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,7 +18,7 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} setaudio ([@var{w_type} [, @var{value}]]) +## @deftypefn {Function File} {} setaudio ([@var{w_type} [, @var{value}]]) ## Execute the shell command @samp{mixer [@var{w_type} [, @var{value}]]} ## @end deftypefn
--- a/scripts/general/cart2pol.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/cart2pol.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,8 +18,8 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} [@var{theta}, @var{r}] = cart2pol (@var{x}, @var{y}) -## @deftypefnx {Function File} {} [@var{theta}, @var{r}, @var{z}] = cart2pol (@var{x}, @var{y}, @var{z}) +## @deftypefn {Function File} {[@var{theta}, @var{r}] =} cart2pol (@var{x}, @var{y}) +## @deftypefnx {Function File} {[@var{theta}, @var{r}, @var{z}] =} cart2pol (@var{x}, @var{y}, @var{z}) ## Transform cartesian to polar or cylindrical coordinates. ## @var{x}, @var{y} (and @var{z}) must be of same shape. ## @var{theta} describes the angle relative to the x - axis.
--- a/scripts/general/cart2sph.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/cart2sph.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,7 +18,7 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} [@var{theta}, @var{phi}, @var{r}] = cart2sph (@var{x}, @var{y}, @var{z}) +## @deftypefn {Function File} {[@var{theta}, @var{phi}, @var{r}] =} cart2sph (@var{x}, @var{y}, @var{z}) ## Transform cartesian to spherical coordinates. ## @var{x}, @var{y} and @var{z} must be of same shape. ## @var{theta} describes the angle relative to the x - axis.
--- a/scripts/general/pol2cart.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/pol2cart.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,8 +18,8 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} [@var{x}, @var{y}] = pol2cart (@var{theta}, @var{r}) -## @deftypefnx {Function File} {} [@var{x}, @var{y}, @var{z}] = pol2cart (@var{theta}, @var{r}, @var{z}) +## @deftypefn {Function File} {[@var{x}, @var{y}] =} pol2cart (@var{theta}, @var{r}) +## @deftypefnx {Function File} {[@var{x}, @var{y}, @var{z}] =} pol2cart (@var{theta}, @var{r}, @var{z}) ## Transform polar or cylindrical to cartesian coordinates. ## @var{theta}, @var{r} (and @var{z}) must be of same shape. ## @var{theta} describes the angle relative to the x - axis.
--- a/scripts/general/sph2cart.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/sph2cart.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,7 +18,7 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} {} [@var{x}, @var{y}, @var{z}] = sph2cart (@var{theta}, @var{phi}, @var{r}) +## @deftypefn {Function File} {[@var{x}, @var{y}, @var{z}] =} sph2cart (@var{theta}, @var{phi}, @var{r}) ## Transform spherical to cartesian coordinates. ## @var{x}, @var{y} and @var{z} must be of same shape. ## @var{theta} describes the angle relative to the x-axis.
--- a/scripts/general/tril.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/tril.m Thu Feb 09 09:12:03 2006 +0000 @@ -69,7 +69,11 @@ if (nargin > 0) [nr, nc] = size (x); - retval = zeros (nr, nc, class (x)); + if (isa (x, "sparse")) + retval = sparse (nr, nc); + else + retval = zeros (nr, nc, class (x)); + endif endif if (nargin == 1)
--- a/scripts/general/triu.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/general/triu.m Thu Feb 09 09:12:03 2006 +0000 @@ -28,7 +28,11 @@ if (nargin > 0) [nr, nc] = size (x); - retval = zeros (nr, nc, class (x)); + if (isa (x, "sparse")) + retval = sparse (nr, nc); + else + retval = zeros (nr, nc, class (x)); + endif endif if (nargin == 1)
--- a/scripts/signal/freqz_plot.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/signal/freqz_plot.m Thu Feb 09 09:12:03 2006 +0000 @@ -18,7 +18,7 @@ ## 02110-1301, USA. ## -*- texinfo -*- -## @deftypefn {Function File} freqz_plot (@var{w}, @var{h}) +## @deftypefn {Function File} {} freqz_plot (@var{w}, @var{h}) ## Plot the pass band, stop band and phase response of @var{h}. ## @end deftypefn
--- a/scripts/sparse/sprand.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/sparse/sprand.m Thu Feb 09 09:12:03 2006 +0000 @@ -7,8 +7,8 @@ ## @deftypefnx {Function File} {} sprand (@var{s}) ## Generate a random sparse matrix. The size of the matrix will be ## @var{m} by @var{n}, with a density of values given by @var{d}. -## @var{d} should be between 0 and 1. Values will be normally -## distributed with mean of zero and variance 1. +## @var{d} should be between 0 and 1. Values will be uniformly +## distributed between 0 and 1. ## ## Note: sometimes the actual density may be a bit smaller than @var{d}. ## This is unlikely to happen for large really sparse matrices.
--- a/scripts/sparse/sprandn.m Thu Feb 09 02:30:00 2006 +0000 +++ b/scripts/sparse/sprandn.m Thu Feb 09 09:12:03 2006 +0000 @@ -3,8 +3,8 @@ ## This program is free software and is in the public domain ## -*- texinfo -*- -## @deftypefn {Function File} {} sprand (@var{m}, @var{n}, @var{d}) -## @deftypefnx {Function File} {} sprand (@var{s}) +## @deftypefn {Function File} {} sprandn (@var{m}, @var{n}, @var{d}) +## @deftypefnx {Function File} {} sprandn (@var{s}) ## Generate a random sparse matrix. The size of the matrix will be ## @var{m} by @var{n}, with a density of values given by @var{d}. ## @var{d} should be between 0 and 1. Values will be normally
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/sparse/sprandsym.m Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,75 @@ +## Copyright (C) 2004 David Bateman & Andy Adler +## +## This program is free software; you can redistribute it and/or modify +## it under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 2 of the License, or +## (at your option) any later version. +## +## This program is distributed in the hope that it will be useful, +## but WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with this program; if not, write to the Free Software +## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA +## 02110-1301 USA + +## -*- texinfo -*- +## @deftypefn {Function File} {} sprandsym (@var{n}, @var{d}) +## @deftypefnx {Function File} {} sprandsym (@var{s}) +## Generate a symmetric random sparse matrix. The size of the matrix will be +## @var{n} by @var{n}, with a density of values given by @var{d}. +## @var{d} should be between 0 and 1. Values will be normally +## distributed with mean of zero and variance 1. +## +## Note: sometimes the actual density may be a bit smaller than @var{d}. +## This is unlikely to happen for large really sparse matrices. +## +## If called with a single matrix argument, a random sparse matrix is +## generated wherever the matrix @var{S} is non-zero in its lower +## triangular part. +## @end deftypefn +## @seealso{sprand, sprandn} + +function S = sprandsym(n,d) + if nargin == 1 + [i,j,v,nr,nc] = spfind(tril(n)); + S = sparse(i,j,randn(size(v)),nr,nc); + S = S + tril(S,-1)'; + elseif nargin == 2 + m1 = floor(n/2); + n1 = m1 + 1; + mn1 = m1*n1; + k1 = round(d*mn1); + idx1=unique(fix(rand(min(k1*1.01,k1+10),1)*mn1))+1; + # idx contains random numbers in [1,mn] + # generate 1% or 10 more random values than necessary + # in order to reduce the probability that there are less than k + # distinct values; + # maybe a better strategy could be used + # but I don't think it's worth the price + k1 = min(length(idx1),k1); # actual number of entries in S + j1 = floor((idx1(1:k1)-1)/m1); + i1 = idx1(1:k1) - j1*m1; + + n2 = ceil(n/2); + nn2 = n2*n2; + k2 = round(d*nn2); + idx2=unique(fix(rand(min(k2*1.01,k1+10),1)*nn2))+1; + k2 = min(length(idx2),k2); + j2 = floor((idx2(1:k2)-1)/n2); + i2 = idx2(1:k2) - j2*n2; + + if isempty(i1) && isempty(i2) + S = sparse(n,n); + else + S1 = sparse(i1,j1+1,randn(k1,1),m1,n1); + S = [tril(S1), sparse(m1,m1); ... + sparse(i2,j2+1,randn(k2,1),n2,n2), triu(S1,1)']; + S = S + tril(S,-1)'; + endif + else + usage("sprandsym(n,density) OR sprandsym(S)"); + endif +endfunction
--- a/src/ChangeLog Thu Feb 09 02:30:00 2006 +0000 +++ b/src/ChangeLog Thu Feb 09 09:12:03 2006 +0000 @@ -1,3 +1,15 @@ +2006-02-09 David Bateman <dbateman@free.fr> + + * DLD-FUNCTIONS/spqr.cc: New file for sparse QR and dmperm based on + CSparse. + * DLD-FUNCTIONS/matrix_type.cc (Fmatrix_type): dintinguish between + rectangular and singular matrices. Add tests. + * DLD-FUNCTIONS/luinc.cc: Add tests. + * DLD-FUNCTIONS/spkron.cc: Ditto. + * Makefile.in (DLD_XSRC): Add spqr.cc. + (OCT_LINK_DEPS): Add CSSPARSE_LIBS. + * sparse-xdiv.h: Remove conditio of lssolve. + 2006-02-08 John W. Eaton <jwe@octave.org> * parse.y (frob_function): Clear ID_NAME from top_level symbol
--- a/src/DLD-FUNCTIONS/luinc.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/src/DLD-FUNCTIONS/luinc.cc Thu Feb 09 09:12:03 2006 +0000 @@ -279,6 +279,28 @@ } /* + +%!test +%! a=sparse([1,2,0,0;0,1,2,0;1e-14,0,3,0;0,0,0,1]); +%! [l,u]=luinc(a,1e-10); +%! assert(l*u, sparse([1,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]),1e-10); +%! opts.droptol=1e-10; +%! [l,u]=luinc(a,opts); +%! assert(l*u, sparse([1,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]),1e-10); + +%!test +%! a=sparse([1i,2,0,0;0,1,2,0;1e-14,0,3,0;0,0,0,1]); +%! [l,u]=luinc(a,1e-10); +%! assert(l*u, sparse([1i,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]),1e-10); +%! opts.droptol=1e-10; +%! [l,u]=luinc(a,opts); +%! assert(l*u, sparse([1i,2,0,0;0,1,2,0;0,0,3,0;0,0,0,1]),1e-10); + +%!error splu(sparse([1i,2,0,0;0,1,2,0;1e-14,0,3,0;0,0,0,1])); + +*/ + +/* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: ***
--- a/src/DLD-FUNCTIONS/matrix_type.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/src/DLD-FUNCTIONS/matrix_type.cc Thu Feb 09 09:12:03 2006 +0000 @@ -157,6 +157,13 @@ retval = octave_value ("Tridiagonal Positive Definite"); else if (typ == SparseType::Hermitian) retval = octave_value ("Positive Definite"); + else if (typ == SparseType::Rectangular) + { + if (args(0).rows() == args(0).columns()) + retval = octave_value ("Singular"); + else + retval = octave_value ("Rectangular"); + } else if (typ == SparseType::Full) retval = octave_value ("Full"); else @@ -284,6 +291,66 @@ } /* + +%!assert(matrix_type(speye(10,10)),"Diagonal"); +%!assert(matrix_type(speye(10,10)([2:10,1],:)),"Permuted Diagonal"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[1;sparse(9,1);1]]),"Upper"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[1;sparse(9,1);1]](:,[2,1,3:11])),"Permuted Upper"); +%!assert(matrix_type([speye(10,10),sparse(10,1);1,sparse(1,9),1]),"Lower"); +%!assert(matrix_type([speye(10,10),sparse(10,1);1,sparse(1,9),1]([2,1,3:11],:)),"Permuted Lower"); +%!test +%! bnd=spparms("bandden"); +%! spparms("bandden",0.5); +%! a = spdiags(randn(10,3),[-1,0,1],10,10); +%! assert(matrix_type(a),"Tridiagonal"); +%! assert(matrix_type(abs(a')+abs(a)),"Tridiagonal Positive Definite"); +%! spparms("bandden",bnd); +%!test +%! bnd=spparms("bandden"); +%! spparms("bandden",0.5); +%! a = spdiags(randn(10,4),[-2:1],10,10); +%! assert(matrix_type(a),"Banded"); +%! assert(matrix_type(a'*a),"Banded Positive Definite"); +%! spparms("bandden",bnd); +%!test +%! a=[speye(10,10),[sparse(9,1);1];-1,sparse(1,9),1]; +%! assert(matrix_type(a),"Full"); +%! assert(matrix_type(a'*a),"Positive Definite"); +%!assert(matrix_type(speye(10,11)),"Rectangular"); + +%!assert(matrix_type(1i*speye(10,10)),"Diagonal"); +%!assert(matrix_type(1i*speye(10,10)([2:10,1],:)),"Permuted Diagonal"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[1i;sparse(9,1);1]]),"Upper"); +%!assert(matrix_type([[speye(10,10);sparse(1,10)],[1i;sparse(9,1);1]](:,[2,1,3:11])),"Permuted Upper"); +%!assert(matrix_type([speye(10,10),sparse(10,1);1i,sparse(1,9),1]),"Lower"); +%!assert(matrix_type([speye(10,10),sparse(10,1);1i,sparse(1,9),1]([2,1,3:11],:)),"Permuted Lower"); +%!test +%! bnd=spparms("bandden"); +%! spparms("bandden",0.5); +%! assert(matrix_type(spdiags(1i*randn(10,3),[-1,0,1],10,10)),"Tridiagonal"); +%! a = 1i*randn(9,1);a=[[a;0],ones(10,1),[0;-a]]; +%! assert(matrix_type(spdiags(a,[-1,0,1],10,10)),"Tridiagonal Positive Definite"); +%! spparms("bandden",bnd); +%!test +%! bnd=spparms("bandden"); +%! spparms("bandden",0.5); +%! assert(matrix_type(spdiags(1i*randn(10,4),[-2:1],10,10)),"Banded"); +%! a = 1i*randn(9,2);a=[[a;[0,0]],ones(10,1),[[0;-a(:,2)],[0;0;-a(1:8,1)]]]; +%! assert(matrix_type(spdiags(a,[-2:2],10,10)),"Banded Positive Definite"); +%! spparms("bandden",bnd); +%!test +%! a=[speye(10,10),[sparse(9,1);1i];-1,sparse(1,9),1]; +%! assert(matrix_type(a),"Full"); +%! assert(matrix_type(a'*a),"Positive Definite"); +%!assert(matrix_type(speye(10,11)),"Rectangular"); + +%!test +%! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular"); +%! assert(matrix_type(a),"Singular"); + +*/ + +/* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: ***
--- a/src/DLD-FUNCTIONS/spkron.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/src/DLD-FUNCTIONS/spkron.cc Thu Feb 09 09:12:03 2006 +0000 @@ -138,3 +138,16 @@ return retval; } + +/* + +%!assert(spkron(spdiag([1,2,3]),spdiag([1,2,3])),sparse(kron(diag([1,2,3]),diag([1,2,3])))) +%!assert(spkron(spdiag([1i,2,3]),spdiag([1i,2,3])),sparse(kron(diag([1i,2,3]),diag([1i,2,3])))) + +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/spqr.cc Thu Feb 09 09:12:03 2006 +0000 @@ -0,0 +1,339 @@ +/* + +Copyright (C) 2005 David Bateman +Copyright (C) 1998-2005 Andy Adler + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 2, or (at your option) any +later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with this program; see the file COPYING. If not, write to the +Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, +Boston, MA 02110-1301, USA. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#if HAVE_CXSPARSE +#include <cxsparse/cxs.h> +#endif + +#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" +#include "SparseQR.h" +#include "SparseCmplxQR.h" + +#ifdef IDX_TYPE_LONG +#define CSSPARSE_NAME(name) name ## _dl +#else +#define CSSPARSE_NAME(name) name ## _di +#endif + +// PKG_ADD: dispatch ("qr", "spqr", "sparse matrix"); +// PKG_ADD: dispatch ("qr", "spqr", "sparse complex matrix"); +// PKG_ADD: dispatch ("qr", "spqr", "sparse bool matrix"); +DEFUN_DLD (spqr, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{r} =} spqr (@var{a})\n\ +@deftypefnx {Loadable Function} {@var{r} =} spqr (@var{a},0)\n\ +@deftypefnx {Loadable Function} {[@var{c}, @var{r}] =} spqr (@var{a},@var{b})\n\ +@deftypefnx {Loadable Function} {[@var{c}, @var{r}] =} spqr (@var{a},@var{b},0)\n\ +@cindex QR factorization\n\ +Compute the sparse QR factorization of @var{a}, using @sc{CSparse}.\n\ +As the matrix @var{Q} is in general a full matrix, this function returns\n\ +the @var{Q}-less factorization @var{r} of @var{a}, such that\n\ +@code{@var{r} = chol (@var{a}' * @var{a})}.\n\ +\n\ +If the final argument is the scalar @code{0} and the number of rows is\n\ +larger than the number of columns, then an economy factorization is\n\ +returned. That is @var{r} will have only @code{size (@var{a},1)} rows.\n\ +\n\ +If an additional matrix @var{b} is supplied, then @code{spqr} returns\n\ +@var{c}, where @code{@var{c} = @var{q}' * @var{b}}. This allows the\n\ +least squares approximation of @code{@var{a} \\ @var{b}} to be calculated\n\ +as\n\ +\n\ +@example\n\ +[@var{c},@var{r}] = spqr (@var{a},@var{b})\n\ +@var{x} = @var{r} \\ @var{c}\n\ +@end example\n\ +\n\ +@end deftypefn\n\ +@seealso{spchol, qr}") +{ + int nargin = args.length (); + octave_value_list retval; + bool economy = false; + bool is_cmplx = false; + bool have_b = false; + + if (nargin < 1 || nargin > 3) + print_usage ("spqr"); + else + { + if (args(0).is_complex_type ()) + is_cmplx = true; + if (nargin > 1) + { + have_b = true; + if (args(nargin-1).is_scalar_type ()) + { + int val = args(nargin-1).int_value (); + if (val == 0) + { + economy = true; + have_b = (nargin > 2); + } + } + if (have_b && args(1).is_complex_type ()) + is_cmplx = true; + } + + if (!error_state) + { + if (have_b && nargout < 2) + error ("spqr: incorrect number of output arguments"); + else if (is_cmplx) + { + SparseComplexQR q (args(0).sparse_complex_matrix_value ()); + if (!error_state) + { + if (have_b) + { + retval(1) = q.R (economy); + retval(0) = q.C (args(1).complex_matrix_value ()); + } + else + retval(0) = q.R (economy); + } + } + else + { + SparseQR q (args(0).sparse_matrix_value ()); + if (!error_state) + { + if (have_b) + { + retval(1) = q.R (economy); + retval(0) = q.C (args(1).matrix_value ()); + } + else + retval(0) = q.R (economy); + } + } + } + } + return retval; +} + +/* + +The deactivated tests below can't be tested till rectangular back-subs is +implemented for sparse matrices. + +%!test +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! r = spqr(a); +%! assert(r'*r,a'*a,1e-10) + +%!test +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! q = symamd(a); +%! a = a(q,q); +%! r = spqr(a); +%! assert(r'*r,a'*a,1e-10) + +%!test +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! [c,r] = spqr(a,ones(n,1)); +%! assert (r\c,full(a)\ones(n,1),10e-10) + +%!test +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! b = randn(n,2); +%! [c,r] = spqr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%% Test under-determined systems!! +%!#test +%! n = 20; d= 0.2; +%! a = sprandn(n,n+1,d)+speye(n,n+1); +%! b = randn(n,2); +%! [c,r] = spqr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%!test +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! r = spqr(a); +%! assert(r'*r,a'*a,1e-10) + +%!test +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! q = symamd(a); +%! a = a(q,q); +%! r = spqr(a); +%! assert(r'*r,a'*a,1e-10) + +%!test +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! [c,r] = spqr(a,ones(n,1)); +%! assert (r\c,full(a)\ones(n,1),10e-10) + +%!test +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! b = randn(n,2); +%! [c,r] = spqr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%% Test under-determined systems!! +%!#test +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n+1,d)+speye(n,n+1); +%! b = randn(n,2); +%! [c,r] = spqr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%!error spqr(sprandn(10,10,0.2),ones(10,1)); + +*/ + +static RowVector +put_int (octave_idx_type *p, octave_idx_type n) +{ + RowVector ret (n); + for (octave_idx_type i = 0; i < n; i++) + ret.xelem(i) = p[i] + 1; + return ret; +} + +DEFUN_DLD (dmperm, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{p} =} dmperm (@var{s})\n\ +@deftypefnx {Loadable Function} {[@var{p}. @var{q}. @var{r}, @var{s}] =} dmperm (@var{s})\n\ +\n\ +@cindex Dulmage-Mendelsohn decomposition\n\ +Perform a Deulmage-Mendelsohn permutation on the sparse matrix @var{s}.\n\ +With a single output argument @dfn{dmperm} performs the row permutations\n\ +@var{p} such that @code{@var{s} (@var{p},:)} has no zero elements on the\n\ +diagonal.\n\ +\n\ +Called with two or more output arguments, returns the row and column\n\ +permutations, such that @code{@var{s} (@var{p}, @var{q})} is in block\n\ +triangular form. The values of @var{r} and @var{s} define the boundaries\n\ +of the blocks. If @var{s} is square then @code{@var{r} == @var{s}}.\n\ +\n\ +The method used is described in: A. Pothen & C.-J. Fan. Computing the block\n\ +triangular form of a sparse matrix. ACM Trans. Math. Software,\n\ +16(4):303-324, 1990.\n\ +@end deftypefn\n\ +@seealso{colamd,ccolamd}") +{ + int nargin = args.length(); + octave_value_list retval; + +#if HAVE_CXSPARSE + if (nargin != 1) + { + print_usage ("dmperm"); + return retval; + } + + octave_value arg = args(0); + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + SparseMatrix m; + SparseComplexMatrix cm; + CSSPARSE_NAME (cs) csm; + csm.m = nr; + csm.n = nc; + csm.x = NULL; + csm.nz = -1; + + if (arg.is_real_type ()) + { + m = arg.sparse_matrix_value (); + csm.nzmax = m.nnz(); + csm.p = m.xcidx (); + csm.i = m.xridx (); + } + else + { + cm = arg.sparse_complex_matrix_value (); + csm.nzmax = cm.nnz(); + csm.p = cm.xcidx (); + csm.i = cm.xridx (); + } + + if (!error_state) + { + if (nargout <= 1) + { + octave_idx_type *jmatch = CSSPARSE_NAME (cs_maxtrans) (&csm); + retval(0) = put_int (jmatch + nr, nc); + CSSPARSE_NAME (cs_free) (jmatch); + } + else + { + CSSPARSE_NAME (csd) *dm = CSSPARSE_NAME(cs_dmperm) (&csm); + //retval(5) = put_int (dm->rr, 5); + //retval(4) = put_int (dm->cc, 5); + retval(3) = put_int (dm->S, dm->nb+1); + retval(2) = put_int (dm->R, dm->nb+1); + retval(1) = put_int (dm->Q, nc); + retval(0) = put_int (dm->P, nr); + CSSPARSE_NAME (cs_dfree) (dm); + } + } +#else + error ("dmperm: not available in this version of Octave"); +#endif + + return retval; +} + +/* + +%!test +%! n=20; +%! a=speye(n,n);a=a(randperm(n),:); +%! assert(a(dmperm(a),:),speye(n)) + +%!test +%! n=20; +%! d=0.2; +%! a=tril(sprandn(n,n,d),-1)+speye(n,n); +%! a=a(randperm(n),randperm(n)); +%! [p,q,r,s]=dmperm(a); +%! assert(tril(a(p,q),-1),sparse(n,n)) + +*/ + +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/Makefile.in Thu Feb 09 02:30:00 2006 +0000 +++ b/src/Makefile.in Thu Feb 09 09:12:03 2006 +0000 @@ -51,7 +51,7 @@ lpsolve.cc lsode.cc lu.cc luinc.cc matrix_type.cc minmax.cc \ pinv.cc qr.cc quad.cc qz.cc rand.cc regexp.cc schur.cc sort.cc \ sparse.cc spchol.cc spdet.cc spkron.cc splu.cc spparms.cc \ - sqrtm.cc svd.cc syl.cc time.cc \ + spqr.cc sqrtm.cc svd.cc syl.cc time.cc \ __gnuplot_raw__.l __glpk__.cc __qp__.cc DLD_SRC := $(addprefix DLD-FUNCTIONS/, $(DLD_XSRC)) @@ -252,7 +252,8 @@ OCT_LINK_DEPS = \ -L../libcruft $(LIBCRUFT) -L../liboctave $(LIBOCTAVE) \ -L. $(LIBOCTINTERP) $(CHOLMOD_LIBS) $(UMFPACK_LIBS) $(AMD_LIBS) \ - $(COLAMD_LIBS) $(CCOLAMD_LIBS) $(BLAS_LIBS) $(FFTW_LIBS) $(LIBS) $(FLIBS) + $(COLAMD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ + $(FFTW_LIBS) $(LIBS) $(FLIBS) DISTFILES = Makefile.in ChangeLog mkdefs mkops mkgendoc \ DOCSTRINGS mkbuiltins mk-errno-list mk-pkg-add \ @@ -313,7 +314,7 @@ $(OCTAVE_LFLAGS) \ $(OCTAVE_LIBS) \ $(LEXLIB) $(UMFPACK_LIBS) $(AMD_LIBS) $(COLAMD_LIBS) \ - $(CHOLMOD_LIBS) $(CCOLAMD_LIBS) $(BLAS_LIBS) \ + $(CHOLMOD_LIBS) $(CCOLAMD_LIBS) $(CXSPARSE_LIBS) $(BLAS_LIBS) \ $(FFTW_LIBS) $(LIBS) $(FLIBS) stmp-pic: pic
--- a/src/sparse-xdiv.cc Thu Feb 09 02:30:00 2006 +0000 +++ b/src/sparse-xdiv.cc Thu Feb 09 09:12:03 2006 +0000 @@ -39,12 +39,7 @@ static inline bool result_ok (octave_idx_type info) { -#ifdef HAVE_LSSOLVE return (info != -2 && info != -1); -#else - // If the matrix is singular, who cares as we don't have QR based solver yet - return (info != -1); -#endif } static void
--- a/test/ChangeLog Thu Feb 09 02:30:00 2006 +0000 +++ b/test/ChangeLog Thu Feb 09 09:12:03 2006 +0000 @@ -1,3 +1,7 @@ +2006-02-09 David Bateman <dbateman@free.fr> + + * test/build_sparse_tests.sh: Add tests for sparse QR solvers. + 2006-01-21 David Bateman <dbateman@free.fr> * build_sparsetest.sh: Add new un-ordered indexing, assignment and
--- a/test/build_sparse_tests.sh Thu Feb 09 02:30:00 2006 +0000 +++ b/test/build_sparse_tests.sh Thu Feb 09 09:12:03 2006 +0000 @@ -914,6 +914,26 @@ %!assert(sparse(tcs\xs),sparse(tcf\xf),1e-10); EOF + +cat >>$TESTS <<EOF +%% QR solver tests + +%!function f(a, sz, feps) +%! b = randn(sz); x = a \b; +%! assert (a * x, b, feps); +%! b = randn(sz)+1i*randn(sz); x = a \ b; +%! assert (a * x, b, feps); +%! b = sprandn(sz(1),sz(2),0.2); x = a \b; +%! assert (sparse(a * x), b, feps); +%! b = sprandn(sz(1),sz(2),0.2)+1i*sprandn(sz(1),sz(2),0.2); x = a \b; +%! assert (sparse(a * x), b, feps); +%!test +%! a = alpha*sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10); +%! ## Test this by forcing matrix_type +%! a = alpha*sprandn(10,10,0.2)+speye(10,10); matrix_type(a, "Singular"); +%! f(a,[10,2],1e-10); + +EOF }