Mercurial > octave-antonio
diff src/DLD-FUNCTIONS/spqr.cc @ 5610:9761b7d24e9e
[project @ 2006-02-09 09:12:02 by dbateman]
author | dbateman |
---|---|
date | Thu, 09 Feb 2006 09:12:03 +0000 |
parents | |
children | 44c592934077 |
line wrap: on
line diff
--- /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: *** +*/