Mercurial > octave-nkf
changeset 7505:f5005d9510f4
Remove dispatched sparse functions and treat in the generic versions of the functions
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Wed, 20 Feb 2008 15:52:11 -0500 |
| parents | ddcf233d765b |
| children | 798b0a00e80c |
| files | doc/ChangeLog doc/interpreter/sparse.txi liboctave/ChangeLog liboctave/SparseCmplxQR.cc liboctave/SparseCmplxQR.h liboctave/SparseQR.cc liboctave/SparseQR.h scripts/ChangeLog scripts/sparse/colperm.m scripts/sparse/nonzeros.m scripts/sparse/spdiags.m scripts/sparse/spfun.m scripts/sparse/spones.m scripts/sparse/sprand.m scripts/sparse/sprandn.m scripts/sparse/sprandsym.m scripts/sparse/spy.m src/ChangeLog src/DLD-FUNCTIONS/det.cc src/DLD-FUNCTIONS/dmperm.cc src/DLD-FUNCTIONS/find.cc src/DLD-FUNCTIONS/minmax.cc src/DLD-FUNCTIONS/qr.cc src/DLD-FUNCTIONS/sparse.cc src/DLD-FUNCTIONS/spdet.cc src/DLD-FUNCTIONS/spfind.cc src/DLD-FUNCTIONS/spqr.cc src/Makefile.in src/data.cc test/ChangeLog test/build_sparse_tests.sh |
| diffstat | 31 files changed, 1455 insertions(+), 1424 deletions(-) [+] |
line wrap: on
line diff
--- a/doc/ChangeLog Wed Feb 20 14:56:22 2008 -0500 +++ b/doc/ChangeLog Wed Feb 20 15:52:11 2008 -0500 @@ -1,3 +1,8 @@ +2008-02-20 David Bateman <dbateman@free.fr> + + * interpreter/sparse.txi: Remove references to spmin, spmax, + spatan2, spfind, spqr and spdet. + 2008-02-19 Carlo de Falco <kingcrimson@tiscali.it> * interpreter/package.txi: Improve INDEX file documentation.
--- a/doc/interpreter/sparse.txi Wed Feb 20 14:56:22 2008 -0500 +++ b/doc/interpreter/sparse.txi Wed Feb 20 15:52:11 2008 -0500 @@ -213,8 +213,6 @@ creates a sparse (@var{n}+1)-by-(@var{n}+1) sparse matrix with a single diagonal defined. -@DOCSTRING(spatan2) - @DOCSTRING(spcumprod) @DOCSTRING(spcumsum) @@ -310,8 +308,6 @@ @DOCSTRING(spconvert) -@DOCSTRING(spfind) - The above problem can be avoided in oct-files. However, the construction of a sparse matrix from an oct-file is more complex than can be discussed in this brief introduction, and you are referred to chapter @@ -472,7 +468,7 @@ @dfn{sprandn}, @dfn{sprandsym} @item Sparse matrix conversion: - @dfn{full}, @dfn{sparse}, @dfn{spconvert}, @dfn{spfind} + @dfn{full}, @dfn{sparse}, @dfn{spconvert} @item Manipulate sparse matrices @dfn{issparse}, @dfn{nnz}, @dfn{nonzeros}, @dfn{nzmax}, @@ -489,8 +485,8 @@ @item Linear algebra: @dfn{matrix_type}, @dfn{spchol}, @dfn{cpcholinv}, - @dfn{spchol2inv}, @dfn{spdet}, @dfn{spinv}, - @dfn{splchol}, @dfn{splu}, @dfn{spqr}, @dfn{normest}, @dfn{condest}, + @dfn{spchol2inv}, @dfn{spinv}, + @dfn{splchol}, @dfn{splu}, @dfn{normest}, @dfn{condest}, @dfn{sprank} @c @dfn{spaugment} @c @dfn{eigs}, @dfn{svds} but these are in octave-forge for now @@ -503,8 +499,7 @@ @item Miscellaneous: @dfn{spparms}, @dfn{symbfact}, @dfn{spstats}, @dfn{spprod}, @dfn{spcumsum}, @dfn{spsum}, - @dfn{spsumsq}, @dfn{spmin}, @dfn{spmax}, @dfn{spatan2}, - @dfn{spdiag} + @dfn{spsumsq}, @dfn{spdiag} @end table In addition all of the standard Octave mapper functions (ie. basic @@ -844,8 +839,6 @@ @DOCSTRING(spchol2inv) -@DOCSTRING(spdet) - @DOCSTRING(spinv) @DOCSTRING(splchol) @@ -854,8 +847,6 @@ @DOCSTRING(spparms) -@DOCSTRING(spqr) - @DOCSTRING(sprank) @DOCSTRING(symbfact)
--- a/liboctave/ChangeLog Wed Feb 20 14:56:22 2008 -0500 +++ b/liboctave/ChangeLog Wed Feb 20 15:52:11 2008 -0500 @@ -1,3 +1,14 @@ +2008-02-20 David Bateman <dbateman@free.fr> + + * SparseComplexQR.cc (ComplexMatrix + SparseComplexQR::SparseComplexQR_rep::Q + (void) const): New method. + * SparseComplexQR.h (ComplexMatrix + SparseComplexQR::SparseComplexQR_rep::Q + (void) const): Declare it. + * SparseQR.cc (Matrix SparseQR::SparseQR_rep::Q (void) const): ditto. + * SparseQR.h (Matrix SparseQR::SparseQR_rep::Q (void) const): ditto. + 2008-02-20 John W. Eaton <jwe@octave.org> * boolNDArray.h (boolNDArray (const Array2<bool>&)): Delete.
--- a/liboctave/SparseCmplxQR.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/liboctave/SparseCmplxQR.cc Wed Feb 20 15:52:11 2008 -0500 @@ -39,10 +39,12 @@ cs_complex_t *buf = reinterpret_cast<cs_complex_t *> (buf ## tmp); #define OCTAVE_C99_ZERO (0. + 0.iF) +#define OCTAVE_C99_ONE (1. + 0.iF) #else #define OCTAVE_C99_COMPLEX(buf, n) \ OCTAVE_LOCAL_BUFFER (cs_complex_t, buf, (n)); #define OCTAVE_C99_ZERO cs_complex_t(0., 0.); +#define OCTAVE_C99_ONE cs_complex_t(1., 0.); #endif SparseComplexQR::SparseComplexQR_rep::SparseComplexQR_rep @@ -232,6 +234,57 @@ } ComplexMatrix +SparseComplexQR::SparseComplexQR_rep::Q (void) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + ComplexMatrix ret(nr, nr); + Complex *vec = ret.fortran_vec(); + if (nr < 0 || nc < 0) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + else if (nr == 0 || nc == 0) + ret = ComplexMatrix (nc, nr, Complex (0.0, 0.0)); + else + { + OCTAVE_C99_COMPLEX (bvec, nr); + for (octave_idx_type i = 0; i < nr; i++) + bvec[i] = OCTAVE_C99_ZERO; + OCTAVE_LOCAL_BUFFER (Complex, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) + { + OCTAVE_QUIT; + bvec[j] = OCTAVE_C99_ONE; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; +#if defined(CS_VER) && (CS_VER >= 2) + CXSPARSE_ZNAME (_ipvec) + (S->pinv, bvec, reinterpret_cast<cs_complex_t *>(buf), nr); +#else + CXSPARSE_ZNAME (_ipvec) + (nr, S->Pinv, bvec, reinterpret_cast<cs_complex_t *>(buf)); +#endif + 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 (_happly) + (N->L, i, N->B[i], reinterpret_cast<cs_complex_t *>(buf)); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < nr; i++) + vec[i+idx] = buf[i]; + bvec[j] = OCTAVE_C99_ZERO; + } + } + return ret.hermitian (); +#else + return ComplexMatrix (); +#endif +} + +ComplexMatrix qrsolve(const SparseComplexMatrix&a, const Matrix &b, octave_idx_type &info) { info = -1;
--- a/liboctave/SparseCmplxQR.h Wed Feb 20 14:56:22 2008 -0500 +++ b/liboctave/SparseCmplxQR.h Wed Feb 20 15:52:11 2008 -0500 @@ -63,6 +63,8 @@ ComplexMatrix C (const ComplexMatrix &b) const; + ComplexMatrix Q (void) const; + int count; octave_idx_type nrows; @@ -116,6 +118,8 @@ ComplexMatrix C (const ComplexMatrix &b) const { return rep->C(b); } + ComplexMatrix Q (void) const { return rep->Q(); } + friend ComplexMatrix qrsolve (const SparseComplexMatrix &a, const Matrix &b, octave_idx_type &info);
--- a/liboctave/SparseQR.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/liboctave/SparseQR.cc Wed Feb 20 15:52:11 2008 -0500 @@ -215,6 +215,57 @@ } Matrix +SparseQR::SparseQR_rep::Q (void) const +{ +#ifdef HAVE_CXSPARSE + octave_idx_type nc = N->L->n; + octave_idx_type nr = nrows; + Matrix ret (nr, nr); + double *vec = ret.fortran_vec(); + if (nr < 0 || nc < 0) + (*current_liboctave_error_handler) ("matrix dimension mismatch"); + else if (nr == 0 || nc == 0) + ret = Matrix (nc, nr, 0.0); + else + { + OCTAVE_LOCAL_BUFFER (double, bvec, nr + 1); + for (octave_idx_type i = 0; i < nr; i++) + bvec[i] = 0.; + OCTAVE_LOCAL_BUFFER (double, buf, S->m2); + for (volatile octave_idx_type j = 0, idx = 0; j < nr; j++, idx+=nr) + { + OCTAVE_QUIT; + bvec[j] = 1.0; + for (octave_idx_type i = nr; i < S->m2; i++) + buf[i] = 0.; + volatile octave_idx_type nm = (nr < nc ? nr : nc); + BEGIN_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; +#if defined(CS_VER) && (CS_VER >= 2) + CXSPARSE_DNAME (_ipvec) (S->pinv, bvec, buf, nr); +#else + CXSPARSE_DNAME (_ipvec) (nr, S->Pinv, bvec, buf); +#endif + 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 (_happly) (N->L, i, N->B[i], buf); + END_INTERRUPT_IMMEDIATELY_IN_FOREIGN_CODE; + } + for (octave_idx_type i = 0; i < nr; i++) + vec[i+idx] = buf[i]; + bvec[j] = 0.0; + } + } + return ret.transpose (); +#else + return Matrix (); +#endif +} + +Matrix qrsolve(const SparseMatrix&a, const Matrix &b, octave_idx_type& info) { info = -1;
--- a/liboctave/SparseQR.h Wed Feb 20 14:56:22 2008 -0500 +++ b/liboctave/SparseQR.h Wed Feb 20 15:52:11 2008 -0500 @@ -63,6 +63,8 @@ Matrix C (const Matrix &b) const; + Matrix Q (void) const; + int count; octave_idx_type nrows; @@ -114,6 +116,8 @@ Matrix C (const Matrix &b) const { return rep->C(b); } + Matrix Q (void) const { return rep->Q(); } + friend Matrix qrsolve (const SparseMatrix &a, const Matrix &b, octave_idx_type &info);
--- a/scripts/ChangeLog Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/ChangeLog Wed Feb 20 15:52:11 2008 -0500 @@ -2,6 +2,13 @@ * strings/strcat.m: Detect cellstr args. +2008-02-20 David Bateman <dbateman@free.fr> + + * sparse/colperm.m, sparse/nonzero.m, sparse/spdiags.m, + sparse/spfun.m, sparse/spones.m, sparse/sprand.m, + sparse/sprandn.m, sparse/sprandsym.m, sparse/spy.m: Use generic + version of find rather than spfind. + 2008-02-19 Ben Abbott <bpabbott@mac.com> * miscellaneous/edit.m: New option EDITINPLACE. Prefer file list
--- a/scripts/sparse/colperm.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/colperm.m Wed Feb 20 15:52:11 2008 -0500 @@ -31,7 +31,7 @@ print_usage (); endif - [i, j] = spfind (s); + [i, j] = find (s); idx = find (diff ([j; Inf]) != 0); [dummy, p] = sort (idx - [0; idx(1:(end-1))]); endfunction
--- a/scripts/sparse/nonzeros.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/nonzeros.m Wed Feb 20 15:52:11 2008 -0500 @@ -27,11 +27,7 @@ print_usage (); endif - if (issparse (s)) - [i, j, t] = spfind (s); - else - [i, j, t] = find (s); - endif + [i, j, t] = find (s); endfunction %!assert(nonzeros([1,2;3,0]),[1;3;2])
--- a/scripts/sparse/spdiags.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/spdiags.m Wed Feb 20 15:52:11 2008 -0500 @@ -56,7 +56,9 @@ if (nargin == 1 || nargin == 2) ## extract nonzero diagonals of v into A,c - [i, j, v, nr, nc] = spfind (v); + [i, j, v] = find (v); + [nr, nc] = size (v); + if (nargin == 1) ## c contains the active diagonals c = unique (j-i);
--- a/scripts/sparse/spfun.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/spfun.m Wed Feb 20 15:52:11 2008 -0500 @@ -30,12 +30,8 @@ print_usage (); endif - if (issparse (s)) - [i,j,v,m,n] = spfind (s); - else - [i, j, v] = find (s); - [m, n] = size (s); - endif + [i, j, v] = find (s); + [m, n] = size (s); if (isa (f, "function_handle") || isa (f, "inline function")) t = sparse (i, j, f(v), m, n);
--- a/scripts/sparse/spones.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/spones.m Wed Feb 20 15:52:11 2008 -0500 @@ -28,12 +28,8 @@ print_usage (); endif - if (issparse (s)) - [i, j, v, m, n] = spfind (s); - else - [i, j, v] = find (s); - [m, n] = size (s); - endif + [i, j, v] = find (s); + [m, n] = size (s); s = sparse (i, j, 1, m, n);
--- a/scripts/sparse/sprand.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/sprand.m Wed Feb 20 15:52:11 2008 -0500 @@ -47,7 +47,8 @@ function S = sprand (m, n, d) if (nargin == 1) - [i, j, v, nr, nc] = spfind (m); + [i, j, v] = find (m); + [nr, nc] = size (m); S = sparse (i, j, rand (size (v)), nr, nc); elseif (nargin == 3) mn = n*m;
--- a/scripts/sparse/sprandn.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/sprandn.m Wed Feb 20 15:52:11 2008 -0500 @@ -39,7 +39,8 @@ function S = sprandn (m, n, d) if (nargin == 1) - [i, j, v, nr, nc] = spfind (m); + [i, j, v] = find (m); + [nr, nc] = size (m); S = sparse (i, j, randn (size (v)), nr, nc); elseif (nargin == 3) mn = m*n;
--- a/scripts/sparse/sprandsym.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/sprandsym.m Wed Feb 20 15:52:11 2008 -0500 @@ -35,7 +35,8 @@ function S = sprandsym (n, d) if (nargin == 1) - [i, j, v, nr, nc] = spfind (tril (n)); + [i, j, v] = find (tril (n)); + [nr, nc] = size (n); S = sparse (i, j, randn (size (v)), nr, nc); S = S + tril (S, -1)'; elseif (nargin == 2)
--- a/scripts/sparse/spy.m Wed Feb 20 14:56:22 2008 -0500 +++ b/scripts/sparse/spy.m Wed Feb 20 15:52:11 2008 -0500 @@ -50,12 +50,8 @@ endif endfor - if (issparse (S)) - [i, j, s, m, n] = spfind (S); - else - [i, j, s] = find (S); - [m, n] = size (S); - endif + [i, j, s] = find (S); + [m, n] = size (S); if (isnan (markersize)) plot (j, i, LineSpec);
--- a/src/ChangeLog Wed Feb 20 14:56:22 2008 -0500 +++ b/src/ChangeLog Wed Feb 20 15:52:11 2008 -0500 @@ -4,6 +4,21 @@ 2008-02-20 David Bateman <dbateman@free.fr> + * DLD-FUNCTIONS/det.cc, DLD-FUNCTIONS/find.cc, + * DLD-FUNCTIONS/minmax.cc, DLD-FUNCTIONS/qr.cc: + Treat sparse matrices. + + * DLD-FUNCTIONS/sparse.cc (Fspmin, Fspmax, Fatan2): Remove functions. + + * DLD-FUNCTIONS/dmperm.cc: Rename from spqr.cc. + (Fspqr): Delete function. + + * DLD-FUNCTIONS/spqr.cc, DLD-FUNCTIONS/spdet.cc, + DLD-FUNCTIONS/spfind.cc: Remove. + + * Makefile.in (DLD_XSRC): Add dmperm.cc to the list. + Delete det.cc, find.cc, minmax.cc, and qr.cc from the list. + * Makefile.in (OV_SRC): Remove ov-mapper.cc. (OV_INCLUDES): Remove ov-mapper.h. (DEFUN_PATTERN): No longer accept DEFUN_MAPPER as valid.
--- a/src/DLD-FUNCTIONS/det.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/DLD-FUNCTIONS/det.cc Wed Feb 20 15:52:11 2008 -0500 @@ -37,8 +37,9 @@ DEFUN_DLD (det, args, , "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{d}, @var{rcond}] = } det (@var{a})\n\ -Compute the determinant of @var{a} using @sc{Lapack}. Return an estimate\n\ -of the reciprocal condition number if requested.\n\ +Compute the determinant of @var{a} using @sc{Lapack} for full and UMFPACK\n\ +for sparse matrices. Return an estimate of the reciprocal condition number\n\ +if requested.\n\ @end deftypefn") { octave_value_list retval; @@ -53,8 +54,8 @@ octave_value arg = args(0); - int nr = arg.rows (); - int nc = arg.columns (); + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); if (nr == 0 && nc == 0) { @@ -76,39 +77,67 @@ if (arg.is_real_type ()) { - Matrix m = arg.matrix_value (); - - if (! error_state) + octave_idx_type info; + double rcond = 0.0; + // Always compute rcond, so we can detect numerically + // singular matrices. + if (arg.is_sparse_type ()) { - // Always compute rcond, so we can detect numerically - // singular matrices. - - octave_idx_type info; - double rcond = 0.0; - DET det = m.determinant (info, rcond); - retval(1) = rcond; - volatile double xrcond = rcond; - xrcond += 1.0; - retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ()); + SparseMatrix m = arg.sparse_matrix_value (); + if (! error_state) + { + DET det = m.determinant (info, rcond); + retval(1) = rcond; + volatile double xrcond = rcond; + xrcond += 1.0; + retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ()); + } + } + else + { + Matrix m = arg.matrix_value (); + if (! error_state) + { + DET det = m.determinant (info, rcond); + retval(1) = rcond; + volatile double xrcond = rcond; + xrcond += 1.0; + retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ()); + } } } else if (arg.is_complex_type ()) { - ComplexMatrix m = arg.complex_matrix_value (); - - if (! error_state) + octave_idx_type info; + double rcond = 0.0; + // Always compute rcond, so we can detect numerically + // singular matrices. + if (arg.is_sparse_type ()) { - // Always compute rcond, so we can detect numerically - // singular matrices. + SparseComplexMatrix m = arg.sparse_complex_matrix_value (); + if (! error_state) + { + ComplexDET det = m.determinant (info, rcond); + retval(1) = rcond; + volatile double xrcond = rcond; + xrcond += 1.0; + retval(0) = ((info == -1 || xrcond == 1.0) + ? Complex (0.0) : det.value ()); + } + } + else + { + ComplexMatrix m = arg.complex_matrix_value (); + if (! error_state) + { + ComplexDET det = m.determinant (info, rcond); + retval(1) = rcond; + volatile double xrcond = rcond; + xrcond += 1.0; + retval(0) = ((info == -1 || xrcond == 1.0) + ? Complex (0.0) : det.value ()); - octave_idx_type info; - double rcond = 0.0; - ComplexDET det = m.determinant (info, rcond); - retval(1) = rcond; - volatile double xrcond = rcond; - xrcond += 1.0; - retval(0) = ((info == -1 || xrcond == 1.0) - ? Complex (0.0) : det.value ()); + } } } else
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/DLD-FUNCTIONS/dmperm.cc Wed Feb 20 15:52:11 2008 -0500 @@ -0,0 +1,236 @@ +/* + +Copyright (C) 2005, 2006, 2007 David Bateman +Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Andy Adler + +This file is part of Octave. + +Octave is free software; you can redistribute it and/or modify it +under the terms of the GNU General Public License as published by the +Free Software Foundation; either version 3 of the License, or (at your +option) any later version. + +Octave is distributed in the hope that it will be useful, but WITHOUT +ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or +FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License +for more details. + +You should have received a copy of the GNU General Public License +along with Octave; see the file COPYING. If not, see +<http://www.gnu.org/licenses/>. + +*/ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include "defun-dld.h" +#include "error.h" +#include "gripes.h" +#include "oct-obj.h" +#include "utils.h" + +#include "oct-sparse.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" +#include "SparseQR.h" +#include "SparseCmplxQR.h" + +#ifdef IDX_TYPE_LONG +#define CXSPARSE_NAME(name) cs_dl ## name +#else +#define CXSPARSE_NAME(name) cs_di ## name +#endif + +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; +} + +#if HAVE_CXSPARSE +static octave_value_list +dmperm_internal (bool rank, const octave_value arg, int nargout) +{ + octave_value_list retval; + octave_idx_type nr = arg.rows (); + octave_idx_type nc = arg.columns (); + SparseMatrix m; + SparseComplexMatrix cm; + CXSPARSE_NAME () 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 || rank) + { +#if defined(CS_VER) && (CS_VER >= 2) + octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm, 0); +#else + octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm); +#endif + if (rank) + { + octave_idx_type r = 0; + for (octave_idx_type i = 0; i < nc; i++) + if (jmatch[nr+i] >= 0) + r++; + retval(0) = static_cast<double>(r); + } + else + retval(0) = put_int (jmatch + nr, nc); + CXSPARSE_NAME (_free) (jmatch); + } + else + { +#if defined(CS_VER) && (CS_VER >= 2) + CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm, 0); +#else + CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm); +#endif + + //retval(5) = put_int (dm->rr, 5); + //retval(4) = put_int (dm->cc, 5); +#if defined(CS_VER) && (CS_VER >= 2) + 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); +#else + 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); +#endif + CXSPARSE_NAME (_dfree) (dm); + } + } + return retval; +} +#endif + +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 Dulmage-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\ +@seealso{colamd, ccolamd}\n\ +@end deftypefn") +{ + int nargin = args.length(); + octave_value_list retval; + + if (nargin != 1) + { + print_usage (); + return retval; + } + +#if HAVE_CXSPARSE + retval = dmperm_internal (false, args(0), nargout); +#else + error ("dmperm: not available in this version of Octave"); +#endif + + return retval; +} + +/* + +%!testif HAVE_CXSPARSE +%! n=20; +%! a=speye(n,n);a=a(randperm(n),:); +%! assert(a(dmperm(a),:),speye(n)) + +%!testif HAVE_CXSPARSE +%! 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)) + +*/ + +DEFUN_DLD (sprank, args, nargout, + "-*- texinfo -*-\n\ +@deftypefn {Loadable Function} {@var{p} =} sprank (@var{s})\n\ +\n\ +@cindex Structural Rank\n\ +Calculates the structural rank of a sparse matrix @var{s}. Note that\n\ +only the structure of the matrix is used in this calculation based on\n\ +a Dulmage-Mendelsohn to block triangular form. As such the numerical\n\ +rank of the matrix @var{s} is bounded by @code{sprank (@var{s}) >=\n\ +rank (@var{s})}. Ignoring floating point errors @code{sprank (@var{s}) ==\n\ +rank (@var{s})}.\n\ +@seealso{dmperm}\n\ +@end deftypefn") +{ + int nargin = args.length(); + octave_value_list retval; + + if (nargin != 1) + { + print_usage (); + return retval; + } + +#if HAVE_CXSPARSE + retval = dmperm_internal (true, args(0), nargout); +#else + error ("sprank: not available in this version of Octave"); +#endif + + return retval; +} + +/* + +%!error(sprank(1,2)); +%!testif HAVE_CXSPARSE +%! assert(sprank(speye(20)), 20) +%!testif HAVE_CXSPARSE +%! assert(sprank([1,0,2,0;2,0,4,0]),2) + +*/ +/* +;;; Local Variables: *** +;;; mode: C++ *** +;;; End: *** +*/
--- a/src/DLD-FUNCTIONS/find.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/DLD-FUNCTIONS/find.cc Wed Feb 20 15:52:11 2008 -0500 @@ -38,8 +38,8 @@ template <typename T> octave_value_list -find_nonzero_elem_idx (const T& nda, int nargout, octave_idx_type n_to_find, - int direction) +find_nonzero_elem_idx (const Array<T>& nda, int nargout, + octave_idx_type n_to_find, int direction) { octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ()); @@ -104,7 +104,7 @@ Matrix i_idx (result_nr, result_nc); Matrix j_idx (result_nr, result_nc); - T val (dim_vector (result_nr, result_nc)); + ArrayN<T> val (dim_vector (result_nr, result_nc)); if (count > 0) { @@ -172,10 +172,159 @@ return retval; } -template octave_value_list find_nonzero_elem_idx (const NDArray&, int, +template octave_value_list find_nonzero_elem_idx (const Array<double>&, int, + octave_idx_type, int); + +template octave_value_list find_nonzero_elem_idx (const Array<Complex>&, int, octave_idx_type, int); -template octave_value_list find_nonzero_elem_idx (const ComplexNDArray&, int, +template <typename T> +octave_value_list +find_nonzero_elem_idx (const Sparse<T>& v, int nargout, + octave_idx_type n_to_find, int direction) +{ + octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ()); + + + octave_idx_type nc = v.cols(); + octave_idx_type nr = v.rows(); + octave_idx_type nz = v.nnz(); + + // Search in the default range. + octave_idx_type start_nc = -1; + octave_idx_type end_nc = -1; + octave_idx_type count; + + // Search for the range to search + if (n_to_find < 0) + { + start_nc = 0; + end_nc = nc; + n_to_find = nz; + count = nz; + } + else if (direction > 0) + { + for (octave_idx_type j = 0; j < nc; j++) + { + OCTAVE_QUIT; + if (v.cidx(j) == 0 && v.cidx(j+1) != 0) + start_nc = j; + if (v.cidx(j+1) >= n_to_find) + { + end_nc = j + 1; + break; + } + } + } + else + { + for (octave_idx_type j = nc; j > 0; j--) + { + OCTAVE_QUIT; + if (v.cidx(j) == nz && v.cidx(j-1) != nz) + end_nc = j; + if (nz - v.cidx(j-1) >= n_to_find) + { + start_nc = j - 1; + break; + } + } + } + + count = (n_to_find > v.cidx(end_nc) - v.cidx(start_nc) ? + v.cidx(end_nc) - v.cidx(start_nc) : n_to_find); + + // If the original argument was a row vector, force a row vector of + // the overall indices to be returned. But see below for scalar + // case... + + octave_idx_type result_nr = count; + octave_idx_type result_nc = 1; + + bool scalar_arg = false; + + if (v.rows () == 1) + { + result_nr = 1; + result_nc = count; + + scalar_arg = (v.columns () == 1); + } + + Matrix idx (result_nr, result_nc); + + Matrix i_idx (result_nr, result_nc); + Matrix j_idx (result_nr, result_nc); + + ArrayN<T> val (dim_vector (result_nr, result_nc)); + + if (count > 0) + { + // Search for elements to return. Only search the region where + // there are elements to be found using the count that we want + // to find. + for (octave_idx_type j = start_nc, cx = 0; j < end_nc; j++) + for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++ ) + { + OCTAVE_QUIT; + if (direction < 0 && i < nz - count) + continue; + i_idx (cx) = static_cast<double> (v.ridx(i) + 1); + j_idx (cx) = static_cast<double> (j + 1); + idx (cx) = j * nr + v.ridx(i) + 1; + val (cx) = v.data(i); + cx++; + if (cx == count) + break; + } + } + else if (scalar_arg) + { + idx.resize (0, 0); + + i_idx.resize (0, 0); + j_idx.resize (0, 0); + + val.resize (dim_vector (0, 0)); + } + + switch (nargout) + { + case 0: + case 1: + retval(0) = idx; + break; + + case 5: + retval(4) = nc; + // Fall through + + case 4: + retval(3) = nr; + // Fall through + + case 3: + retval(2) = val; + // Fall through! + + case 2: + retval(1) = j_idx; + retval(0) = i_idx; + break; + + default: + panic_impossible (); + break; + } + + return retval; +} + +template octave_value_list find_nonzero_elem_idx (const Sparse<double>&, int, + octave_idx_type, int); + +template octave_value_list find_nonzero_elem_idx (const Sparse<Complex>&, int, octave_idx_type, int); DEFUN_DLD (find, args, nargout, @@ -224,6 +373,19 @@ If three inputs are given, @var{direction} should be one of \"first\" or\n\ \"last\" indicating that it should start counting found elements from the\n\ first or last element.\n\ +\n\ +Note that this function is particularly useful for sparse matrices, as\n\ +it extracts the non-zero elements as vectors, which can then be used to\n\ +create the original matrix. For example,\n\ +\n\ +@example\n\ +@group\n\ +sz = size(a);\n\ +[i, j, v] = find (a);\n\ +b = sparse(i, j, v, sz(1), sz(2));\n\ +@end group\n\ +@end example\n\ +@seealso{sparse}\n\ @end deftypefn") { octave_value_list retval; @@ -273,30 +435,57 @@ octave_value arg = args(0); - if (arg.is_real_type ()) - { - NDArray nda = arg.array_value (); - - if (! error_state) - retval = find_nonzero_elem_idx (nda, nargout, n_to_find, direction); - } - else if (arg.is_complex_type ()) + if (arg.is_sparse_type ()) { - ComplexNDArray cnda = arg.complex_array_value (); + if (arg.is_real_type ()) + { + SparseMatrix v = arg.sparse_matrix_value (); - if (! error_state) - retval = find_nonzero_elem_idx (cnda, nargout, n_to_find, direction); - } - else if (arg.is_string ()) - { - charNDArray cnda = arg.char_array_value (); + if (! error_state) + retval = find_nonzero_elem_idx (v, nargout, + n_to_find, direction); + } + else if (arg.is_complex_type ()) + { + SparseComplexMatrix v = arg.sparse_complex_matrix_value (); - if (! error_state) - retval = find_nonzero_elem_idx (cnda, nargout, n_to_find, direction); + if (! error_state) + retval = find_nonzero_elem_idx (v, nargout, + n_to_find, direction); + } + else + gripe_wrong_type_arg ("find", arg); } else { - gripe_wrong_type_arg ("find", arg); + if (arg.is_real_type ()) + { + NDArray nda = arg.array_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (nda, nargout, + n_to_find, direction); + } + else if (arg.is_complex_type ()) + { + ComplexNDArray cnda = arg.complex_array_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (cnda, nargout, + n_to_find, direction); + } + else if (arg.is_string ()) + { + charNDArray cnda = arg.char_array_value (); + + if (! error_state) + retval = find_nonzero_elem_idx (cnda, nargout, + n_to_find, direction); + } + else + { + gripe_wrong_type_arg ("find", arg); + } } return retval;
--- a/src/DLD-FUNCTIONS/minmax.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/DLD-FUNCTIONS/minmax.cc Wed Feb 20 15:52:11 2008 -0500 @@ -38,6 +38,8 @@ #include "oct-obj.h" #include "ov-cx-mat.h" +#include "ov-re-sparse.h" +#include "ov-cx-sparse.h" #define MINMAX_DOUBLE_BODY(FCN) \ { \ @@ -310,6 +312,155 @@ } \ } +#define MINMAX_SPARSE_BODY(FCN) \ +{ \ + bool single_arg = (nargin == 1) || arg2.is_empty(); \ + \ + if (single_arg && (nargout == 1 || nargout == 0)) \ + { \ + if (arg1.type_id () == octave_sparse_matrix::static_type_id ()) \ + retval(0) = arg1.sparse_matrix_value () .FCN (dim); \ + else if (arg1.type_id () == \ + octave_sparse_complex_matrix::static_type_id ()) \ + retval(0) = arg1.sparse_complex_matrix_value () .FCN (dim); \ + else \ + gripe_wrong_type_arg (#FCN, arg1); \ + } \ + else if (single_arg && nargout == 2) \ + { \ + Array2<octave_idx_type> index; \ + \ + if (arg1.type_id () == octave_sparse_matrix::static_type_id ()) \ + retval(0) = arg1.sparse_matrix_value () .FCN (index, dim); \ + else if (arg1.type_id () == \ + octave_sparse_complex_matrix::static_type_id ()) \ + retval(0) = arg1.sparse_complex_matrix_value () .FCN (index, dim); \ + else \ + gripe_wrong_type_arg (#FCN, arg1); \ + \ + octave_idx_type len = index.numel (); \ + \ + if (len > 0) \ + { \ + double nan_val = lo_ieee_nan_value (); \ + \ + NDArray idx (index.dims ()); \ + \ + for (octave_idx_type i = 0; i < len; i++) \ + { \ + OCTAVE_QUIT; \ + octave_idx_type tmp = index.elem (i) + 1; \ + idx.elem (i) = (tmp <= 0) \ + ? nan_val : static_cast<double> (tmp); \ + } \ + \ + retval(1) = idx; \ + } \ + else \ + retval(1) = NDArray (); \ + } \ + else \ + { \ + int arg1_is_scalar = arg1.is_scalar_type (); \ + int arg2_is_scalar = arg2.is_scalar_type (); \ + \ + int arg1_is_complex = arg1.is_complex_type (); \ + int arg2_is_complex = arg2.is_complex_type (); \ + \ + if (arg1_is_scalar) \ + { \ + if (arg1_is_complex || arg2_is_complex) \ + { \ + Complex c1 = arg1.complex_value (); \ + \ + SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseComplexMatrix result = FCN (c1, m2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + else \ + { \ + double d1 = arg1.double_value (); \ + SparseMatrix m2 = arg2.sparse_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseMatrix result = FCN (d1, m2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + } \ + else if (arg2_is_scalar) \ + { \ + if (arg1_is_complex || arg2_is_complex) \ + { \ + SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \ + \ + if (! error_state) \ + { \ + Complex c2 = arg2.complex_value (); \ + SparseComplexMatrix result = FCN (m1, c2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + else \ + { \ + SparseMatrix m1 = arg1.sparse_matrix_value (); \ + \ + if (! error_state) \ + { \ + double d2 = arg2.double_value (); \ + SparseMatrix result = FCN (m1, d2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + } \ + else \ + { \ + if (arg1_is_complex || arg2_is_complex) \ + { \ + SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseComplexMatrix result = FCN (m1, m2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + } \ + else \ + { \ + SparseMatrix m1 = arg1.sparse_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseMatrix m2 = arg2.sparse_matrix_value (); \ + \ + if (! error_state) \ + { \ + SparseMatrix result = FCN (m1, m2); \ + if (! error_state) \ + retval(0) = result; \ + } \ + } \ + } \ + } \ + } \ +} + + #define MINMAX_BODY(FCN) \ \ octave_value_list retval; \ @@ -388,6 +539,8 @@ else if (arg1.is_int64_type ()) \ MINMAX_INT_BODY (FCN, int64) \ } \ + else if (arg1.is_sparse_type ()) \ + MINMAX_SPARSE_BODY (FCN) \ else \ MINMAX_DOUBLE_BODY (FCN) \ \
--- a/src/DLD-FUNCTIONS/qr.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/DLD-FUNCTIONS/qr.cc Wed Feb 20 15:52:11 2008 -0500 @@ -28,6 +28,9 @@ #include "CmplxQRP.h" #include "dbleQR.h" #include "dbleQRP.h" +#include "SparseQR.h" +#include "SparseCmplxQR.h" + #include "defun-dld.h" #include "error.h" @@ -55,6 +58,7 @@ DEFUN_DLD (qr, args, nargout, "-*- texinfo -*-\n\ @deftypefn {Loadable Function} {[@var{q}, @var{r}, @var{p}] =} qr (@var{a})\n\ +@deftypefnx {Loadable Function} {[@var{q}, @var{r}, @var{p}] =} qr (@var{a}, '0')\n\ @cindex QR factorization\n\ Compute the QR factorization of @var{a}, using standard @sc{Lapack}\n\ subroutines. For example, given the matrix @code{a = [1, 2; 3, 4]},\n\ @@ -114,10 +118,14 @@ upper triangular.\n\ @end ifinfo\n\ \n\ -The permuted QR factorization @code{[@var{q}, @var{r}, @var{p}] =\n\ -qr (@var{a})} forms the QR factorization such that the diagonal\n\ -entries of @code{r} are decreasing in magnitude order. For example,\n\ -given the matrix @code{a = [1, 2; 3, 4]},\n\ +If given a second argument of '0', @code{qr} returns an economy-sized\n\ +QR factorization, omitting zero rows of @var{R} and the corresponding\n\ +columns of @var{Q}.\n\ +\n\ +If the matrix @var{a} is full, the permuted QR factorization\n\ +@code{[@var{q}, @var{r}, @var{p}] = qr (@var{a})} forms the QR factorization\n\ +such that the diagonal entries of @code{r} are decreasing in magnitude\n\ +order. For example,given the matrix @code{a = [1, 2; 3, 4]},\n\ \n\ @example\n\ [q, r, p] = qr(a)\n\ @@ -147,16 +155,31 @@ factorization allows the construction of an orthogonal basis of\n\ @code{span (a)}.\n\ \n\ -If given a second argument, @code{qr} returns an economy-sized\n\ -QR factorization, omitting zero rows of @var{R} and\n\ -the corresponding columns of @var{Q}.\n\ +If the matrix @var{a} is sparse, then compute the sparse QR factorization\n\ +of @var{a}, using @sc{CSparse}. As the matrix @var{Q} is in general a full\n\ +matrix, this function returns the @var{Q}-less factorization @var{r} of\n\ +@var{a}, such that @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{qr} 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\ @end deftypefn") { octave_value_list retval; int nargin = args.length (); - if (nargin < 1 || nargin > 2 || nargout > 3) + if (nargin < 1 || nargin > (args(0).is_sparse_type() ? 3 : 2) || nargout > 3) { print_usage (); return retval; @@ -171,88 +194,243 @@ else if (arg_is_empty > 0) return octave_value_list (3, Matrix ()); - QR::type type = (nargout == 0 || nargout == 1) ? QR::raw - : (nargin == 2 ? QR::economy : QR::std); - - if (arg.is_real_type ()) + if (arg.is_sparse_type ()) { - Matrix m = arg.matrix_value (); - - if (! error_state) - { - switch (nargout) - { - case 0: - case 1: - { - QR fact (m, type); - retval(0) = fact.R (); - } - break; - - case 2: - { - QR fact (m, type); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - break; + bool economy = false; + bool is_cmplx = false; + int have_b = 0; - default: - { - QRP fact (m, type); - retval(2) = fact.P (); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - break; + if (arg.is_complex_type ()) + is_cmplx = true; + if (nargin > 1) + { + have_b = 1; + if (args(nargin-1).is_scalar_type ()) + { + int val = args(nargin-1).int_value (); + if (val == 0) + { + economy = true; + have_b = (nargin > 2 ? 2 : 0); + } } + if (have_b > 0 && args(have_b).is_complex_type ()) + is_cmplx = true; } - } - else if (arg.is_complex_type ()) - { - ComplexMatrix m = arg.complex_matrix_value (); - - if (! error_state) + + if (!error_state) { - switch (nargout) + if (have_b && nargout < 2) + error ("qr: incorrect number of output arguments"); + else if (is_cmplx) { - case 0: - case 1: - { - ComplexQR fact (m, type); - retval(0) = fact.R (); - } - break; - - case 2: - { - ComplexQR fact (m, type); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - break; - - default: - { - ComplexQRP fact (m, type); - retval(2) = fact.P (); - retval(1) = fact.R (); - retval(0) = fact.Q (); - } - break; + SparseComplexQR q (arg.sparse_complex_matrix_value ()); + if (!error_state) + { + if (have_b > 0) + { + retval(1) = q.R (economy); + retval(0) = q.C (args(have_b).complex_matrix_value ()); + if (arg.rows() < arg.columns()) + warning ("qr: non minimum norm solution for under-determined problem"); + } + else if (nargout > 1) + { + retval(1) = q.R (economy); + retval(0) = q.Q (); + } + else + retval(0) = q.R (economy); + } + } + else + { + SparseQR q (arg.sparse_matrix_value ()); + if (!error_state) + { + if (have_b > 0) + { + retval(1) = q.R (economy); + retval(0) = q.C (args(have_b).matrix_value ()); + if (args(0).rows() < args(0).columns()) + warning ("qr: non minimum norm solution for under-determined problem"); + } + else if (nargout > 1) + { + retval(1) = q.R (economy); + retval(0) = q.Q (); + } + else + retval(0) = q.R (economy); + } } } } else { - gripe_wrong_type_arg ("qr", arg); + QR::type type = (nargout == 0 || nargout == 1) ? QR::raw + : (nargin == 2 ? QR::economy : QR::std); + + if (arg.is_real_type ()) + { + Matrix m = arg.matrix_value (); + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + { + QR fact (m, type); + retval(0) = fact.R (); + } + break; + + case 2: + { + QR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + break; + + default: + { + QRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + break; + } + } + } + else if (arg.is_complex_type ()) + { + ComplexMatrix m = arg.complex_matrix_value (); + + if (! error_state) + { + switch (nargout) + { + case 0: + case 1: + { + ComplexQR fact (m, type); + retval(0) = fact.R (); + } + break; + + case 2: + { + ComplexQR fact (m, type); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + break; + + default: + { + ComplexQRP fact (m, type); + retval(2) = fact.P (); + retval(1) = fact.R (); + retval(0) = fact.Q (); + } + break; + } + } + } + else + { + gripe_wrong_type_arg ("qr", arg); + } } return retval; } /* + +The deactivated tests below can't be tested till rectangular back-subs is +implemented for sparse matrices. + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! r = qr(a); +%! assert(r'*r,a'*a,1e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! q = symamd(a); +%! a = a(q,q); +%! r = qr(a); +%! assert(r'*r,a'*a,1e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! [c,r] = qr(a,ones(n,1)); +%! assert (r\c,full(a)\ones(n,1),10e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = sprandn(n,n,d)+speye(n,n); +%! b = randn(n,2); +%! [c,r] = qr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%% Test under-determined systems!! +%!#testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = sprandn(n,n+1,d)+speye(n,n+1); +%! b = randn(n,2); +%! [c,r] = qr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! r = qr(a); +%! assert(r'*r,a'*a,1e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! q = symamd(a); +%! a = a(q,q); +%! r = qr(a); +%! assert(r'*r,a'*a,1e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! [c,r] = qr(a,ones(n,1)); +%! assert (r\c,full(a)\ones(n,1),10e-10) + +%!testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n,d)+speye(n,n); +%! b = randn(n,2); +%! [c,r] = qr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%% Test under-determined systems!! +%!#testif HAVE_CXSPARSE +%! n = 20; d= 0.2; +%! a = 1i*sprandn(n,n+1,d)+speye(n,n+1); +%! b = randn(n,2); +%! [c,r] = qr(a,b); +%! assert (r\c,full(a)\b,10e-10) + +%!error qr(sprandn(10,10,0.2),ones(10,1)); + +*/ + + +/* ;;; Local Variables: *** ;;; mode: C++ *** ;;; End: ***
--- a/src/DLD-FUNCTIONS/sparse.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/DLD-FUNCTIONS/sparse.cc Wed Feb 20 15:52:11 2008 -0500 @@ -487,360 +487,6 @@ SPARSE_DIM_ARG_BODY (spsumsq, sumsq); } -#define MINMAX_BODY(FCN) \ - \ - octave_value_list retval; \ - \ - int nargin = args.length (); \ - \ - if (nargin < 1 || nargin > 3 || nargout > 2) \ - { \ - print_usage (); \ - return retval; \ - } \ - \ - octave_value arg1; \ - octave_value arg2; \ - octave_value arg3; \ - \ - switch (nargin) \ - { \ - case 3: \ - arg3 = args(2); \ - \ - case 2: \ - arg2 = args(1); \ - \ - case 1: \ - arg1 = args(0); \ - break; \ - \ - default: \ - panic_impossible (); \ - break; \ - } \ - \ - int dim; \ - dim_vector dv = arg1.dims (); \ - if (error_state) \ - { \ - gripe_wrong_type_arg (#FCN, arg1); \ - return retval; \ - } \ - \ - if (nargin == 3) \ - { \ - dim = arg3.nint_value () - 1; \ - if (dim < 0 || dim >= dv.length ()) \ - { \ - error ("%s: invalid dimension", #FCN); \ - return retval; \ - } \ - } \ - else \ - { \ - dim = 0; \ - while ((dim < dv.length ()) && (dv (dim) <= 1)) \ - dim++; \ - if (dim == dv.length ()) \ - dim = 0; \ - } \ - \ - bool single_arg = (nargin == 1) || arg2.is_empty(); \ - \ - if (single_arg && (nargout == 1 || nargout == 0)) \ - { \ - if (arg1.type_id () == octave_sparse_matrix::static_type_id ()) \ - retval(0) = arg1.sparse_matrix_value () .FCN (dim); \ - else if (arg1.type_id () == \ - octave_sparse_complex_matrix::static_type_id ()) \ - retval(0) = arg1.sparse_complex_matrix_value () .FCN (dim); \ - else \ - gripe_wrong_type_arg (#FCN, arg1); \ - } \ - else if (single_arg && nargout == 2) \ - { \ - Array2<octave_idx_type> index; \ - \ - if (arg1.type_id () == octave_sparse_matrix::static_type_id ()) \ - retval(0) = arg1.sparse_matrix_value () .FCN (index, dim); \ - else if (arg1.type_id () == \ - octave_sparse_complex_matrix::static_type_id ()) \ - retval(0) = arg1.sparse_complex_matrix_value () .FCN (index, dim); \ - else \ - gripe_wrong_type_arg (#FCN, arg1); \ - \ - octave_idx_type len = index.numel (); \ - \ - if (len > 0) \ - { \ - double nan_val = lo_ieee_nan_value (); \ - \ - NDArray idx (index.dims ()); \ - \ - for (octave_idx_type i = 0; i < len; i++) \ - { \ - OCTAVE_QUIT; \ - octave_idx_type tmp = index.elem (i) + 1; \ - idx.elem (i) = (tmp <= 0) \ - ? nan_val : static_cast<double> (tmp); \ - } \ - \ - retval(1) = idx; \ - } \ - else \ - retval(1) = NDArray (); \ - } \ - else \ - { \ - int arg1_is_scalar = arg1.is_scalar_type (); \ - int arg2_is_scalar = arg2.is_scalar_type (); \ - \ - int arg1_is_complex = arg1.is_complex_type (); \ - int arg2_is_complex = arg2.is_complex_type (); \ - \ - if (arg1_is_scalar) \ - { \ - if (arg1_is_complex || arg2_is_complex) \ - { \ - Complex c1 = arg1.complex_value (); \ - \ - SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseComplexMatrix result = FCN (c1, m2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - else \ - { \ - double d1 = arg1.double_value (); \ - SparseMatrix m2 = arg2.sparse_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseMatrix result = FCN (d1, m2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - } \ - else if (arg2_is_scalar) \ - { \ - if (arg1_is_complex || arg2_is_complex) \ - { \ - SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \ - \ - if (! error_state) \ - { \ - Complex c2 = arg2.complex_value (); \ - SparseComplexMatrix result = FCN (m1, c2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - else \ - { \ - SparseMatrix m1 = arg1.sparse_matrix_value (); \ - \ - if (! error_state) \ - { \ - double d2 = arg2.double_value (); \ - SparseMatrix result = FCN (m1, d2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - } \ - else \ - { \ - if (arg1_is_complex || arg2_is_complex) \ - { \ - SparseComplexMatrix m1 = arg1.sparse_complex_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseComplexMatrix m2 = arg2.sparse_complex_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseComplexMatrix result = FCN (m1, m2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - } \ - else \ - { \ - SparseMatrix m1 = arg1.sparse_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseMatrix m2 = arg2.sparse_matrix_value (); \ - \ - if (! error_state) \ - { \ - SparseMatrix result = FCN (m1, m2); \ - if (! error_state) \ - retval(0) = result; \ - } \ - } \ - } \ - } \ - } \ - \ - return retval - -// PKG_ADD: dispatch ("min", "spmin", "sparse matrix"); -// PKG_ADD: dispatch ("min", "spmin", "sparse complex matrix"); -// PKG_ADD: dispatch ("min", "spmin", "sparse bool matrix"); -DEFUN_DLD (spmin, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Mapping Function} {} spmin (@var{x}, @var{y}, @var{dim})\n\ -@deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} spmin (@var{x})\n\ -@cindex Utility Functions\n\ -For a vector argument, return the minimum value. For a matrix\n\ -argument, return the minimum value from each column, as a row\n\ -vector, or over the dimension @var{dim} if defined. For two matrices\n\ -(or a matrix and scalar), return the pair-wise minimum.\n\ -Thus,\n\ -\n\ -@example\n\ -min (min (@var{x}))\n\ -@end example\n\ -\n\ -@noindent\n\ -returns the smallest element of @var{x}, and\n\ -\n\ -@example\n\ -@group\n\ -min (2:5, pi)\n\ - @result{} 2.0000 3.0000 3.1416 3.1416\n\ -@end group\n\ -@end example\n\ -@noindent\n\ -compares each element of the range @code{2:5} with @code{pi}, and\n\ -returns a row vector of the minimum values.\n\ -\n\ -For complex arguments, the magnitude of the elements are used for\n\ -comparison.\n\ -\n\ -If called with one input and two output arguments,\n\ -@code{min} also returns the first index of the\n\ -minimum value(s). Thus,\n\ -\n\ -@example\n\ -@group\n\ -[x, ix] = min ([1, 3, 0, 2, 5])\n\ - @result{} x = 0\n\ - ix = 3\n\ -@end group\n\ -@end example\n\ -@end deftypefn") -{ - MINMAX_BODY (min); -} - -// PKG_ADD: dispatch ("max", "spmax", "sparse matrix"); -// PKG_ADD: dispatch ("max", "spmax", "sparse complex matrix"); -// PKG_ADD: dispatch ("max", "spmax", "sparse bool matrix"); -DEFUN_DLD (spmax, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Mapping Function} {} spmax (@var{x}, @var{y}, @var{dim})\n\ -@deftypefnx {Mapping Function} {[@var{w}, @var{iw}] =} spmax (@var{x})\n\ -@cindex Utility Functions\n\ -For a vector argument, return the maximum value. For a matrix\n\ -argument, return the maximum value from each column, as a row\n\ -vector, or over the dimension @var{dim} if defined. For two matrices\n\ -(or a matrix and scalar), return the pair-wise maximum.\n\ -Thus,\n\ -\n\ -@example\n\ -max (max (@var{x}))\n\ -@end example\n\ -\n\ -@noindent\n\ -returns the largest element of @var{x}, and\n\ -\n\ -@example\n\ -@group\n\ -max (2:5, pi)\n\ - @result{} 3.1416 3.1416 4.0000 5.0000\n\ -@end group\n\ -@end example\n\ -@noindent\n\ -compares each element of the range @code{2:5} with @code{pi}, and\n\ -returns a row vector of the maximum values.\n\ -\n\ -For complex arguments, the magnitude of the elements are used for\n\ -comparison.\n\ -\n\ -If called with one input and two output arguments,\n\ -@code{max} also returns the first index of the\n\ -maximum value(s). Thus,\n\ -\n\ -@example\n\ -@group\n\ -[x, ix] = max ([1, 3, 5, 2, 5])\n\ - @result{} x = 5\n\ - ix = 3\n\ -@end group\n\ -@end example\n\ -@end deftypefn") -{ - MINMAX_BODY (max); -} - -// PKG_ADD: dispatch ("atan2", "spatan2", "sparse matrix"); -// PKG_ADD: dispatch ("atan2", "spatan2", "sparse complex matrix"); -// PKG_ADD: dispatch ("atan2", "spatan2", "sparse bool matrix"); -DEFUN_DLD (spatan2, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spatan2 (@var{y}, @var{x})\n\ -Compute atan (Y / X) for corresponding sparse matrix elements of Y and X.\n\ -The result is in range -pi to pi.\n\ -@end deftypefn") -{ - octave_value retval; - int nargin = args.length (); - if (nargin == 2) { - SparseMatrix a, b; - double da, db; - bool is_double_a = false; - bool is_double_b = false; - - if (args(0).is_scalar_type ()) - { - is_double_a = true; - da = args(0).double_value(); - } - else - a = args(0).sparse_matrix_value (); - - if (args(1).is_scalar_type ()) - { - is_double_b = true; - db = args(1).double_value(); - } - else - b = args(1).sparse_matrix_value (); - - if (is_double_a && is_double_b) - retval = Matrix (1, 1, atan2(da, db)); - else if (is_double_a) - retval = atan2 (da, b); - else if (is_double_b) - retval = atan2 (a, db); - else - retval = atan2 (a, b); - - } else - print_usage (); - - return retval; -} static octave_value make_spdiag (const octave_value& a, const octave_value& b)
--- a/src/DLD-FUNCTIONS/spdet.cc Wed Feb 20 14:56:22 2008 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,132 +0,0 @@ -/* - -Copyright (C) 2004, 2005, 2006, 2007 David Bateman -Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004 Andy Adler - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "defun-dld.h" -#include "error.h" -#include "gripes.h" -#include "oct-obj.h" -#include "utils.h" - -#include "dbleDET.h" -#include "CmplxDET.h" - -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" - -// PKG_ADD: dispatch ("det", "spdet", "sparse matrix"); -// PKG_ADD: dispatch ("det", "spdet", "sparse complex matrix"); -// PKG_ADD: dispatch ("det", "spdet", "sparse bool matrix"); -DEFUN_DLD (spdet, args, , - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {[@var{d}, @var{rcond}] = } spdet (@var{a})\n\ -Compute the determinant of sparse matrix @var{a} using UMFPACK. Return\n\ -an estimate of the reciprocal condition number if requested.\n\ -@end deftypefn") -{ - octave_value_list retval; - - int nargin = args.length (); - - if (nargin != 1) - { - print_usage (); - return retval; - } - - octave_value arg = args(0); - - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - - if (nr == 0 && nc == 0) - { - retval(0) = 1.0; - return retval; - } - - int arg_is_empty = empty_arg ("spdet", nr, nc); - if (arg_is_empty < 0) - return retval; - if (arg_is_empty > 0) - return octave_value (Matrix (1, 1, 1.0)); - - if (nr != nc) - { - gripe_square_matrix_required ("spdet"); - return retval; - } - - if (arg.is_real_type ()) - { - SparseMatrix m = args(0).sparse_matrix_value (); - - if (! error_state) - { - // Always compute rcond, so we can detect numerically - // singular matrices. - - octave_idx_type info; - double rcond = 0.0; - DET det = m.determinant (info, rcond); - retval(1) = rcond; - volatile double xrcond = rcond; - xrcond += 1.0; - retval(0) = ((info == -1 || xrcond == 1.0) ? 0.0 : det.value ()); - } - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix m = args(0).sparse_complex_matrix_value (); - - if (! error_state) - { - // Always compute rcond, so we can detect numerically - // singular matrices. - - octave_idx_type info; - double rcond = 0.0; - ComplexDET det = m.determinant (info, rcond); - retval(1) = rcond; - volatile double xrcond = rcond; - xrcond += 1.0; - retval(0) = ((info == -1 || xrcond == 1.0) - ? Complex (0.0) : det.value ()); - } - } - else - { - gripe_wrong_type_arg ("spdet", arg); - } - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/DLD-FUNCTIONS/spfind.cc Wed Feb 20 14:56:22 2008 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,289 +0,0 @@ -/* - -Copyright (C) 2006, 2007 David Bateman - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <cstdlib> -#include <string> - -#include "variables.h" -#include "utils.h" -#include "pager.h" -#include "defun-dld.h" -#include "gripes.h" -#include "quit.h" - -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" - -template <typename T, typename M> -octave_value_list -sparse_find_non_zero_elem_idx (const T& v, M& val, int nargout, - octave_idx_type n_to_find, int direction) -{ - octave_value_list retval ((nargout == 0 ? 1 : nargout), Matrix ()); - - - octave_idx_type nc = v.cols(); - octave_idx_type nr = v.rows(); - octave_idx_type nz = v.nnz(); - - // Search in the default range. - octave_idx_type start_nc = -1; - octave_idx_type end_nc = -1; - octave_idx_type count; - - // Search for the range to search - if (n_to_find < 0) - { - start_nc = 0; - end_nc = nc; - n_to_find = nz; - count = nz; - } - else if (direction > 0) - { - for (octave_idx_type j = 0; j < nc; j++) - { - OCTAVE_QUIT; - if (v.cidx(j) == 0 && v.cidx(j+1) != 0) - start_nc = j; - if (v.cidx(j+1) >= n_to_find) - { - end_nc = j + 1; - break; - } - } - } - else - { - for (octave_idx_type j = nc; j > 0; j--) - { - OCTAVE_QUIT; - if (v.cidx(j) == nz && v.cidx(j-1) != nz) - end_nc = j; - if (nz - v.cidx(j-1) >= n_to_find) - { - start_nc = j - 1; - break; - } - } - } - - count = (n_to_find > v.cidx(end_nc) - v.cidx(start_nc) ? - v.cidx(end_nc) - v.cidx(start_nc) : n_to_find); - - // If the original argument was a row vector, force a row vector of - // the overall indices to be returned. But see below for scalar - // case... - - octave_idx_type result_nr = count; - octave_idx_type result_nc = 1; - - bool scalar_arg = false; - - if (v.rows () == 1) - { - result_nr = 1; - result_nc = count; - - scalar_arg = (v.columns () == 1); - } - - Matrix idx (result_nr, result_nc); - - Matrix i_idx (result_nr, result_nc); - Matrix j_idx (result_nr, result_nc); - - val = M(result_nr, result_nc); - - if (count > 0) - { - // Search for elements to return. Only search the region where - // there are elements to be found using the count that we want - // to find. - for (octave_idx_type j = start_nc, cx = 0; j < end_nc; j++) - for (octave_idx_type i = v.cidx(j); i < v.cidx(j+1); i++ ) - { - OCTAVE_QUIT; - if (direction < 0 && i < nz - count) - continue; - i_idx (cx) = static_cast<double> (v.ridx(i) + 1); - j_idx (cx) = static_cast<double> (j + 1); - idx (cx) = j * nr + v.ridx(i) + 1; - val (cx) = v.data(i); - cx++; - if (cx == count) - break; - } - } - else if (scalar_arg) - { - idx.resize (0, 0); - - i_idx.resize (0, 0); - j_idx.resize (0, 0); - - val.resize (0, 0); - } - - switch (nargout) - { - case 0: - case 1: - retval(0) = idx; - break; - - case 5: - retval(4) = nc; - // Fall through - - case 4: - retval(3) = nr; - // Fall through - - case 3: - retval(2) = val; - // Fall through! - - case 2: - retval(1) = j_idx; - retval(0) = i_idx; - break; - - default: - panic_impossible (); - break; - } - - return retval; -} - -template octave_value_list sparse_find_non_zero_elem_idx - (const SparseMatrix&, Matrix&, int, octave_idx_type, int); - -template octave_value_list sparse_find_non_zero_elem_idx - (const SparseComplexMatrix&, ComplexMatrix&, int, octave_idx_type, int); - - -// PKG_ADD: dispatch ("find", "spfind", "sparse matrix"); -// PKG_ADD: dispatch ("find", "spfind", "sparse complex matrix"); -// PKG_ADD: dispatch ("find", "spfind", "sparse bool matrix"); -DEFUN_DLD (spfind, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {} spfind (@var{x})\n\ -@deftypefnx {Loadable Function} {} spfind (@var{x}, @var{n})\n\ -@deftypefnx {Loadable Function} {} spfind (@var{x}, @var{n}, @var{direction})\n\ -@deftypefnx {Loadable Function} {[@var{i}, @var{j}, @var{v}} spfind (@dots{})\n\ -\n\ -A sparse version of the @code{find} function. Please see the @code{find}\n\ -for details of its use.\n\ -\n\ -Note that this function is particularly useful for sparse matrices, as\n\ -it extracts the non-zero elements as vectors, which can then be used to\n\ -create the original matrix. For example,\n\ -\n\ -@example\n\ -@group\n\ -sz = size(a);\n\ -[i, j, v] = spfind (a);\n\ -b = sparse(i, j, v, sz(1), sz(2));\n\ -@end group\n\ -@end example\n\ -@seealso{sparse}\n\ -@end deftypefn") -{ - octave_value_list retval; - int nargin = args.length (); - - if (nargin > 3 || nargin < 1) - { - print_usage (); - return retval; - } - - // Setup the default options. - octave_idx_type n_to_find = -1; - if (nargin > 1) - { - n_to_find = args(1).int_value (); - if (error_state) - { - error ("find: expecting second argument to be an integer"); - return retval; - } - } - - // Direction to do the searching (1 == forward, -1 == reverse). - int direction = 1; - if (nargin > 2) - { - direction = 0; - - std::string s_arg = args(2).string_value (); - - if (! error_state) - { - if (s_arg == "first") - direction = 1; - else if (s_arg == "last") - direction = -1; - } - - if (direction == 0) - { - error ("find: expecting third argument to be \"first\" or \"last\""); - return retval; - } - } - - octave_value arg = args(0); - - if (arg.is_real_type ()) - { - SparseMatrix v = arg.sparse_matrix_value (); - Matrix val; - - if (! error_state) - retval = sparse_find_non_zero_elem_idx (v, val, nargout, n_to_find, direction); - } - else if (arg.is_complex_type ()) - { - SparseComplexMatrix v = arg.sparse_complex_matrix_value (); - ComplexMatrix val; - - if (! error_state) - retval = sparse_find_non_zero_elem_idx (v, val, nargout, n_to_find, direction); - } - else - gripe_wrong_type_arg ("spfind", arg); - - return retval; -} - -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/DLD-FUNCTIONS/spqr.cc Wed Feb 20 14:56:22 2008 -0500 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,415 +0,0 @@ -/* - -Copyright (C) 2005, 2006, 2007 David Bateman -Copyright (C) 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005 Andy Adler - -This file is part of Octave. - -Octave is free software; you can redistribute it and/or modify it -under the terms of the GNU General Public License as published by the -Free Software Foundation; either version 3 of the License, or (at your -option) any later version. - -Octave is distributed in the hope that it will be useful, but WITHOUT -ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or -FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License -for more details. - -You should have received a copy of the GNU General Public License -along with Octave; see the file COPYING. If not, see -<http://www.gnu.org/licenses/>. - -*/ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include "defun-dld.h" -#include "error.h" -#include "gripes.h" -#include "oct-obj.h" -#include "utils.h" - -#include "oct-sparse.h" -#include "ov-re-sparse.h" -#include "ov-cx-sparse.h" -#include "SparseQR.h" -#include "SparseCmplxQR.h" - -#ifdef IDX_TYPE_LONG -#define CXSPARSE_NAME(name) cs_dl ## name -#else -#define CXSPARSE_NAME(name) cs_di ## name -#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\ -@seealso{spchol, qr}\n\ -@end deftypefn") -{ - 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 (); - 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 ()); - if (args(0).rows() < args(0).columns()) - warning ("spqr: non minimum norm solution for under-determined problem"); - } - 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 ()); - if (args(0).rows() < args(0).columns()) - warning ("spqr: non minimum norm solution for under-determined problem"); - } - 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. - -%!testif HAVE_CXSPARSE -%! n = 20; d= 0.2; -%! a = sprandn(n,n,d)+speye(n,n); -%! r = spqr(a); -%! assert(r'*r,a'*a,1e-10) - -%!testif HAVE_CXSPARSE -%! 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) - -%!testif HAVE_CXSPARSE -%! 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) - -%!testif HAVE_CXSPARSE -%! 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!! -%!#testif HAVE_CXSPARSE -%! 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) - -%!testif HAVE_CXSPARSE -%! n = 20; d= 0.2; -%! a = 1i*sprandn(n,n,d)+speye(n,n); -%! r = spqr(a); -%! assert(r'*r,a'*a,1e-10) - -%!testif HAVE_CXSPARSE -%! 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) - -%!testif HAVE_CXSPARSE -%! 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) - -%!testif HAVE_CXSPARSE -%! 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!! -%!#testif HAVE_CXSPARSE -%! 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; -} - -#if HAVE_CXSPARSE -static octave_value_list -dmperm_internal (bool rank, const octave_value arg, int nargout) -{ - octave_value_list retval; - octave_idx_type nr = arg.rows (); - octave_idx_type nc = arg.columns (); - SparseMatrix m; - SparseComplexMatrix cm; - CXSPARSE_NAME () 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 || rank) - { -#if defined(CS_VER) && (CS_VER >= 2) - octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm, 0); -#else - octave_idx_type *jmatch = CXSPARSE_NAME (_maxtrans) (&csm); -#endif - if (rank) - { - octave_idx_type r = 0; - for (octave_idx_type i = 0; i < nc; i++) - if (jmatch[nr+i] >= 0) - r++; - retval(0) = static_cast<double>(r); - } - else - retval(0) = put_int (jmatch + nr, nc); - CXSPARSE_NAME (_free) (jmatch); - } - else - { -#if defined(CS_VER) && (CS_VER >= 2) - CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm, 0); -#else - CXSPARSE_NAME (d) *dm = CXSPARSE_NAME(_dmperm) (&csm); -#endif - - //retval(5) = put_int (dm->rr, 5); - //retval(4) = put_int (dm->cc, 5); -#if defined(CS_VER) && (CS_VER >= 2) - 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); -#else - 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); -#endif - CXSPARSE_NAME (_dfree) (dm); - } - } - return retval; -} -#endif - -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 Dulmage-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\ -@seealso{colamd, ccolamd}\n\ -@end deftypefn") -{ - int nargin = args.length(); - octave_value_list retval; - - if (nargin != 1) - { - print_usage (); - return retval; - } - -#if HAVE_CXSPARSE - retval = dmperm_internal (false, args(0), nargout); -#else - error ("dmperm: not available in this version of Octave"); -#endif - - return retval; -} - -/* - -%!testif HAVE_CXSPARSE -%! n=20; -%! a=speye(n,n);a=a(randperm(n),:); -%! assert(a(dmperm(a),:),speye(n)) - -%!testif HAVE_CXSPARSE -%! 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)) - -*/ - -DEFUN_DLD (sprank, args, nargout, - "-*- texinfo -*-\n\ -@deftypefn {Loadable Function} {@var{p} =} sprank (@var{s})\n\ -\n\ -@cindex Structural Rank\n\ -Calculates the structural rank of a sparse matrix @var{s}. Note that\n\ -only the structure of the matrix is used in this calculation based on\n\ -a Dulmage-Mendelsohn to block triangular form. As such the numerical\n\ -rank of the matrix @var{s} is bounded by @code{sprank (@var{s}) >=\n\ -rank (@var{s})}. Ignoring floating point errors @code{sprank (@var{s}) ==\n\ -rank (@var{s})}.\n\ -@seealso{dmperm}\n\ -@end deftypefn") -{ - int nargin = args.length(); - octave_value_list retval; - - if (nargin != 1) - { - print_usage (); - return retval; - } - -#if HAVE_CXSPARSE - retval = dmperm_internal (true, args(0), nargout); -#else - error ("sprank: not available in this version of Octave"); -#endif - - return retval; -} - -/* - -%!error(sprank(1,2)); -%!testif HAVE_CXSPARSE -%! assert(sprank(speye(20)), 20) -%!testif HAVE_CXSPARSE -%! assert(sprank([1,0,2,0;2,0,4,0]),2) - -*/ -/* -;;; Local Variables: *** -;;; mode: C++ *** -;;; End: *** -*/
--- a/src/Makefile.in Wed Feb 20 14:56:22 2008 -0500 +++ b/src/Makefile.in Wed Feb 20 15:52:11 2008 -0500 @@ -64,13 +64,13 @@ DLD_XSRC := balance.cc besselj.cc betainc.cc bsxfun.cc cellfun.cc chol.cc \ ccolamd.cc colamd.cc colloc.cc conv2.cc convhulln.cc daspk.cc \ - dasrt.cc dassl.cc det.cc dispatch.cc eig.cc expm.cc \ + dasrt.cc dassl.cc det.cc dispatch.cc dmperm.cc eig.cc expm.cc \ fft.cc fft2.cc fftn.cc fftw.cc filter.cc find.cc fsolve.cc \ gammainc.cc gcd.cc getgrent.cc getpwent.cc getrusage.cc \ givens.cc hess.cc inv.cc kron.cc lsode.cc \ lu.cc luinc.cc matrix_type.cc md5sum.cc minmax.cc pinv.cc qr.cc \ quad.cc qz.cc rand.cc regexp.cc schur.cc sparse.cc \ - spchol.cc spdet.cc spfind.cc splu.cc spparms.cc spqr.cc \ + spchol.cc splu.cc spparms.cc \ sqrtm.cc svd.cc syl.cc symrcm.cc time.cc tsearch.cc typecast.cc \ urlwrite.cc __contourc__.cc __delaunayn__.cc __dsearchn__.cc \ __glpk__.cc __lin_interpn__.cc __pchip_deriv__.cc \
--- a/src/data.cc Wed Feb 20 14:56:22 2008 -0500 +++ b/src/data.cc Wed Feb 20 15:52:11 2008 -0500 @@ -189,6 +189,227 @@ return retval; } +static SparseMatrix +map_d_s (d_dd_fcn f, double x, const SparseMatrix& y) +{ + octave_idx_type nr = y.rows (); + octave_idx_type nc = y.columns (); + SparseMatrix retval; + double f_zero = f (x, 0.); + + if (f_zero != 0) + { + retval = SparseMatrix (nr, nc, f_zero); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++) + { + OCTAVE_QUIT; + retval.data (y.ridx(i) + j * nr) = f (x, y.data (i)); + } + + retval.maybe_compress (true); + } + else + { + octave_idx_type nz = y.nnz (); + retval = SparseMatrix (nr, nc, nz); + octave_idx_type ii = 0; + retval.cidx (ii) = 0; + + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = y.cidx (j); i < y.cidx (j+1); i++) + { + OCTAVE_QUIT; + double val = f (x, y.data (i)); + + if (val != 0.0) + { + retval.data (ii) = val; + retval.ridx (ii++) = y.ridx (i); + } + } + retval.cidx (j + 1) = ii; + } + + retval.maybe_compress (false); + } + + return retval; +} + +static SparseMatrix +map_s_d (d_dd_fcn f, const SparseMatrix& x, double y) +{ + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.columns (); + SparseMatrix retval; + double f_zero = f (0., y); + + if (f_zero != 0) + { + retval = SparseMatrix (nr, nc, f_zero); + + for (octave_idx_type j = 0; j < nc; j++) + for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++) + { + OCTAVE_QUIT; + retval.data (x.ridx(i) + j * nr) = f (x.data (i), y); + } + + retval.maybe_compress (true); + } + else + { + octave_idx_type nz = x.nnz (); + retval = SparseMatrix (nr, nc, nz); + octave_idx_type ii = 0; + retval.cidx (ii) = 0; + + for (octave_idx_type j = 0; j < nc; j++) + { + for (octave_idx_type i = x.cidx (j); i < x.cidx (j+1); i++) + { + OCTAVE_QUIT; + double val = f (x.data (i), y); + + if (val != 0.0) + { + retval.data (ii) = val; + retval.ridx (ii++) = x.ridx (i); + } + } + retval.cidx (j + 1) = ii; + } + + retval.maybe_compress (false); + } + + return retval; +} + +static SparseMatrix +map_s_s (d_dd_fcn f, const SparseMatrix& x, const SparseMatrix& y) +{ + octave_idx_type nr = x.rows (); + octave_idx_type nc = x.columns (); + + octave_idx_type y_nr = y.rows (); + octave_idx_type y_nc = y.columns (); + + assert (nr == y_nr && nc == y_nc); + + SparseMatrix retval; + double f_zero = f (0., 0.); + + if (f_zero != 0) + { + retval = SparseMatrix (nr, nc, f_zero); + octave_idx_type k1 = 0, k2 = 0; + + for (octave_idx_type j = 0; j < nc; j++) + { + while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1)) + { + OCTAVE_QUIT; + if (k1 >= x.cidx(j+1)) + { + retval.data (y.ridx(k2) + j * nr) = f (0.0, y.data (k2)); + k2++; + } + else if (k2 >= y.cidx(j+1)) + { + retval.data (x.ridx(k1) + j * nr) = f (x.data (k1), 0.0); + k1++; + } + else + { + octave_idx_type rx = x.ridx(k1); + octave_idx_type ry = y.ridx(k2); + + if (rx < ry) + { + retval.data (rx + j * nr) = f (x.data (k1), 0.0); + k1++; + } + else if (rx > ry) + { + retval.data (ry + j * nr) = f (0.0, y.data (k2)); + k2++; + } + else + { + retval.data (ry + j * nr) = f (x.data (k1), y.data (k2)); + k1++; + k2++; + } + } + } + } + + retval.maybe_compress (true); + } + else + { + octave_idx_type nz = y.nnz (); + retval = SparseMatrix (nr, nc, nz); + octave_idx_type ii = 0; + retval.cidx (ii) = 0; + octave_idx_type k1 = 0, k2 = 0; + + for (octave_idx_type j = 0; j < nc; j++) + { + while (k1 < x.cidx(j+1) && k2 < y.cidx(j+1)) + { + OCTAVE_QUIT; + double val; + octave_idx_type r; + if (k1 >= x.cidx(j+1)) + { + r = y.ridx (k2); + val = f (0.0, y.data (k2++)); + } + else if (k2 >= y.cidx(j+1)) + { + r = x.ridx (k1); + val = f (x.data (k1++), 0.0); + } + else + { + octave_idx_type rx = x.ridx(k1); + octave_idx_type ry = y.ridx(k2); + + if (rx < ry) + { + r = x.ridx (k1); + val = f (x.data (k1++), 0.0); + } + else if (rx > ry) + { + r = y.ridx (k2); + val = f (0.0, y.data (k2++)); + } + else + { + r = y.ridx (k2); + val = f (x.data (k1++), y.data (k2++)); + } + } + if (val != 0.0) + { + retval.data (ii) = val; + retval.ridx (ii++) = r; + } + } + } + + retval.maybe_compress (false); + } + + return retval; +} + DEFUN (atan2, args, , "-*- texinfo -*-\n\ @deftypefn {Mapping Function} {} atan2 (@var{y}, @var{x})\n\ @@ -244,34 +465,74 @@ if (! error_state) { - Matrix x = arg_x.matrix_value (); - - if (! error_state) - retval = map_d_m (atan2, y, x); + if (arg_x.is_sparse_type ()) + { + SparseMatrix x = arg_x.sparse_matrix_value (); + + if (! error_state) + retval = map_d_s (atan2, y, x); + } + else + { + Matrix x = arg_x.matrix_value (); + + if (! error_state) + retval = map_d_m (atan2, y, x); + } } } else if (x_is_scalar) { - Matrix y = arg_y.matrix_value (); - - if (! error_state) + if (arg_y.is_sparse_type ()) { - double x = arg_x.double_value (); + SparseMatrix y = arg_y.sparse_matrix_value (); if (! error_state) - retval = map_m_d (atan2, y, x); + { + double x = arg_x.double_value (); + + if (! error_state) + retval = map_s_d (atan2, y, x); + } + } + else + { + Matrix y = arg_y.matrix_value (); + + if (! error_state) + { + double x = arg_x.double_value (); + + if (! error_state) + retval = map_m_d (atan2, y, x); + } } } else if (y_nr == x_nr && y_nc == x_nc) { - Matrix y = arg_y.matrix_value (); - - if (! error_state) + if (arg_y.is_sparse_type () || arg_x.is_sparse_type ()) { - Matrix x = arg_x.matrix_value (); + SparseMatrix y = arg_y.sparse_matrix_value (); if (! error_state) - retval = map_m_m (atan2, y, x); + { + SparseMatrix x = arg_x.sparse_matrix_value (); + + if (! error_state) + retval = map_s_s (atan2, y, x); + } + } + else + { + Matrix y = arg_y.matrix_value (); + + if (! error_state) + { + Matrix x = arg_x.matrix_value (); + + if (! error_state) + retval = map_m_m (atan2, y, x); + } } } else @@ -289,7 +550,7 @@ @deftypefn {Mapping Function} {} fmod (@var{x}, @var{y})\n\ Compute the floating point remainder of dividing @var{x} by @var{y}\n\ using the C library function @code{fmod}. The result has the same\n\ -sign as @var{x}. If @var{y} is zero, the result implementation-defined.\n\ +sign as @var{x}. If @var{y} is zero, the result is implementation-defined.\n\ @end deftypefn") { octave_value retval; @@ -336,34 +597,74 @@ if (! error_state) { - Matrix x = arg_x.matrix_value (); - - if (! error_state) - retval = map_m_d (fmod, x, y); + if (arg_x.is_sparse_type ()) + { + SparseMatrix x = arg_x.sparse_matrix_value (); + + if (! error_state) + retval = map_s_d (fmod, x, y); + } + else + { + Matrix x = arg_x.matrix_value (); + + if (! error_state) + retval = map_m_d (fmod, x, y); + } } } else if (x_is_scalar) { - Matrix y = arg_y.matrix_value (); - - if (! error_state) + if (arg_y.is_sparse_type ()) { - double x = arg_x.double_value (); + SparseMatrix y = arg_y.sparse_matrix_value (); if (! error_state) - retval = map_d_m (fmod, x, y); + { + double x = arg_x.double_value (); + + if (! error_state) + retval = map_d_s (fmod, x, y); + } + } + else + { + Matrix y = arg_y.matrix_value (); + + if (! error_state) + { + double x = arg_x.double_value (); + + if (! error_state) + retval = map_d_m (fmod, x, y); + } } } else if (y_nr == x_nr && y_nc == x_nc) { - Matrix y = arg_y.matrix_value (); - - if (! error_state) + if (arg_y.is_sparse_type () || arg_x.is_sparse_type ()) { - Matrix x = arg_x.matrix_value (); + SparseMatrix y = arg_y.sparse_matrix_value (); if (! error_state) - retval = map_m_m (fmod, x, y); + { + SparseMatrix x = arg_x.sparse_matrix_value (); + + if (! error_state) + retval = map_s_s (fmod, x, y); + } + } + else + { + Matrix y = arg_y.matrix_value (); + + if (! error_state) + { + Matrix x = arg_x.matrix_value (); + + if (! error_state) + retval = map_m_m (fmod, x, y); + } } } else
--- a/test/ChangeLog Wed Feb 20 14:56:22 2008 -0500 +++ b/test/ChangeLog Wed Feb 20 15:52:11 2008 -0500 @@ -1,3 +1,8 @@ +2008-02-19 David Bateman <dbateman@free.fr> + + * build_sparse_tests.sh: Replaced removed spars functions like + spmin, with their generic names. + 2008-01-22 John W. Eaton <jwe@octave.org> * test_poly.m, test_set.m, test_stats.m: Delete files with no tests.
--- a/test/build_sparse_tests.sh Wed Feb 20 14:56:22 2008 -0500 +++ b/test/build_sparse_tests.sh Wed Feb 20 15:52:11 2008 -0500 @@ -547,26 +547,26 @@ %!assert(spcumprod(as,1),sparse(cumprod(af,1))) %!assert(spcumprod(as,2),sparse(cumprod(af,2))) -%!assert(spmin(as),sparse(min(af))) -%!assert(full(spmin(as(:))),min(af(:))) -%!assert(spmin(as,[],1),sparse(min(af,[],1))) -%!assert(spmin(as,[],2),sparse(min(af,[],2))) -%!assert(spmin(as,[],1),sparse(min(af,[],1))) -%!assert(spmin(as,0),sparse(min(af,0))) -%!assert(spmin(as,bs),sparse(min(af,bf))) -%!assert(spmax(as),sparse(max(af))) -%!assert(full(spmax(as(:))),max(af(:))) -%!assert(spmax(as,[],1),sparse(max(af,[],1))) -%!assert(spmax(as,[],2),sparse(max(af,[],2))) -%!assert(spmax(as,[],1),sparse(max(af,[],1))) -%!assert(spmax(as,0),sparse(max(af,0))) -%!assert(spmax(as,bs),sparse(max(af,bf))) +%!assert(min(as),sparse(min(af))) +%!assert(full(min(as(:))),min(af(:))) +%!assert(min(as,[],1),sparse(min(af,[],1))) +%!assert(min(as,[],2),sparse(min(af,[],2))) +%!assert(min(as,[],1),sparse(min(af,[],1))) +%!assert(min(as,0),sparse(min(af,0))) +%!assert(min(as,bs),sparse(min(af,bf))) +%!assert(max(as),sparse(max(af))) +%!assert(full(max(as(:))),max(af(:))) +%!assert(max(as,[],1),sparse(max(af,[],1))) +%!assert(max(as,[],2),sparse(max(af,[],2))) +%!assert(max(as,[],1),sparse(max(af,[],1))) +%!assert(max(as,0),sparse(max(af,0))) +%!assert(max(as,bs),sparse(max(af,bf))) %!assert(as==as) %!assert(as==af) %!assert(af==as) %!test -%! [ii,jj,vv,nr,nc] = spfind(as); +%! [ii,jj,vv,nr,nc] = find(as); %! assert(af,full(sparse(ii,jj,vv,nr,nc))); %!assert(nnz(as),sum(af(:)!=0)) %!assert(nnz(as),nnz(af)) @@ -581,21 +581,21 @@ %!error [i,j]=size(af);as(i-1,j+1); %!error [i,j]=size(af);as(i+1,j-1); %!test -%! [Is,Js,Vs] = spfind(as); +%! [Is,Js,Vs] = find(as); %! [If,Jf,Vf] = find(af); %! assert(Is,If); %! assert(Js,Jf); %! assert(Vs,Vf); %!error as(0,1); %!error as(1,0); -%!assert(spfind(as),find(af)) +%!assert(find(as),find(af)) %!test -%! [i,j,v] = spfind(as); +%! [i,j,v] = find(as); %! [m,n] = size(as); %! x = sparse(i,j,v,m,n); %! assert(x,as); %!test -%! [i,j,v,m,n] = spfind(as); +%! [i,j,v,m,n] = find(as); %! x = sparse(i,j,v,m,n); %! assert(x,as); %!assert(issparse(horzcat(as,as))); @@ -631,7 +631,7 @@ cat >>$TESTS <<EOF %!testif HAVE_UMFPACK -%! assert(spdet(bs+speye(size(bs))),det(bf+eye(size(bf))),100*eps*abs(det(bf+eye(size(bf))))) +%! assert(det(bs+speye(size(bs))),det(bf+eye(size(bf))),100*eps*abs(det(bf+eye(size(bf))))) %!testif HAVE_UMFPACK %! [l,u]=splu(sparse([1,1;1,1])); @@ -650,36 +650,36 @@ %! [L,U,P] = splu(bs); %! assert(P'*L*U,bs,1e-10); %! # triangularity -%! [i,j,v]=spfind(L); +%! [i,j,v]=find(L); %! assert(i-j>=0); -%! [i,j,v]=spfind(U); +%! [i,j,v]=find(U); %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# simple LU + row/col permutations %! [L,U,P,Q] = splu(bs); %! assert(P'*L*U*Q',bs,1e-10); %! # triangularity -%! [i,j,v]=spfind(L); +%! [i,j,v]=find(L); %! assert(i-j>=0); -%! [i,j,v]=spfind(U); +%! [i,j,v]=find(U); %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# LU with fixed column permutation %! [L,U,P] = splu(bs,colamd(bs)); %! assert(P'*L*U,bs,1e-10); %! # triangularity -%! [i,j,v]=spfind(L); +%! [i,j,v]=find(L); %! assert(i-j>=0); -%! [i,j,v]=spfind(U(:,colamd(bs))); +%! [i,j,v]=find(U(:,colamd(bs))); %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# LU with initial column permutation %! [L,U,P,Q] = splu(bs,colamd(bs)); %! assert(P'*L*U*Q',bs,1e-10); %! # triangularity -%! [i,j,v]=spfind(L); +%! [i,j,v]=find(L); %! assert(i-j>=0); -%! [i,j,v]=spfind(U); +%! [i,j,v]=find(U); %! assert(j-i>=0); %!testif HAVE_UMFPACK ;# inverse @@ -1020,20 +1020,20 @@ %!test %! us = alpha*[speye(11,9),[1;sparse(8,1);1;0]]; %!testif HAVE_CXSPARSE -%! [c,r] = spqr (us, xf); +%! [c,r] = qr (us, xf); %! assert(us\xf,r\c,100*eps) %!testif HAVE_CXSPARSE -%! [c,r] = spqr (us, xs); +%! [c,r] = qr (us, xs); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(us\xs,r\c,100*eps) %!test %! pus = us(:,[1:8,10,9]); %!testif HAVE_CXSPARSE -%! [c,r] = spqr (pus, xf); +%! [c,r] = qr (pus, xf); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(pus\xf,r\c,100*eps) %!testif HAVE_CXSPARSE -%! [c,r] = spqr (pus, xs); +%! [c,r] = qr (pus, xs); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(pus\xs,r\c,100*eps) %!test @@ -1051,20 +1051,20 @@ %! xf = beta * ones(12,2); %! xs = speye(12,12); %!testif HAVE_CXSPARSE -%! [c,r] = spqr (ls, xf); +%! [c,r] = qr (ls, xf); %! assert(ls\xf,r\c,100*eps) %!testif HAVE_CXSPARSE -%! [c,r] = spqr (ls, xs); +%! [c,r] = qr (ls, xs); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(ls\xs,r\c,100*eps) %!testif HAVE_CXSPARSE %! pls = ls(:,[1:8,10,9]); %!testif HAVE_CXSPARSE -%! [c,r] = spqr (pls, xf); +%! [c,r] = qr (pls, xf); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(pls\xf,r\c,100*eps) %!testif HAVE_CXSPARSE -%! [c,r] = spqr (pls, xs); +%! [c,r] = qr (pls, xs); %! r = matrix_type(r,"Singular"); ## Force Matrix Type %! assert(pls\xs,r\c,100*eps)