[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: 2.9.5 and 2.1.73
From: |
David Bateman |
Subject: |
Re: 2.9.5 and 2.1.73 |
Date: |
Tue, 14 Mar 2006 18:12:32 +0100 |
User-agent: |
Mozilla Thunderbird 1.0.6-7.5.20060mdk (X11/20050322) |
With the mentioned attachment this time..
D.
David Bateman wrote:
>John W. Eaton wrote:
>
>
>
>>I'd like to make a new set of snapshots soon (this week, if
>>possible). Are there any outstanding bugs that you think must be
>>fixed before I do that?
>>
>>Thanks,
>>
>>jwe
>>
>>
>>
>>
>>
>>
>Yes. I have quite a few changes, in particular to the sparse rectangular
>solver code I recently introduced. In fact in is partially disabled in
>the current CVS. I have a couple of issues (ie seg-faults) with the code
>at the moment but will send it hopefully before the end of the week.
>
>There are also three features I have patches for I'd like to add for 2.9.5
>
>* Disable the calculation of the condition number for the diagonal and
>triangular sparse solvers. Given the manner in which these solvers work,
>its generally not possible to pass these singular matrices. You have to
>force the matrix type to be able to. The calculation of the condition
>number can often take orders of magnitude longer than the calculation of
>the solution for these two cases, and so disabling the makes the speed
>comparable (or slightly faster) than matlab.
>* Special case sparse permutations like A = B(p,q) where p and q are
>strict permutations (ie the elements of p and q are unique). This gives
>several orders of magnitude is speedup of the above indexing operation
>* Sparse Dulmange-Mendelsohn solver. I'd like to have this in 2.9.5 so
>that my Octave 2006 paper can refer to 2.9.5 as the version the
>benchmarks are created against. However I have a question on how to use
>C++ template functions in octave. My problem is I need to write
>"dmsolve<SparseMatrix,SparseMatrix,SparseMatrix>(...)" rather than
>"dmsolve(....)" to call the specific version of a template function. I
>can find no instances of the use of such calls to specific template
>functions in octave and so I'm not sure of the best way to include the
>code. I attach my current version of the code as a separate function
>dmsol, so if you can comment on the best way to include this I'd
>appreciate it..
>
>Cheers
>David
>
>
>
>
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
/*
Copyright (C) 2006 David Bateman
Octave is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2, or (at your option) any
later version.
Octave is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with this program; see the file COPYING. If not, write to the
Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
Boston, MA 02110-1301, USA.
*/
//#ifdef HAVE_CONFIG_H
#include <config.h>
//#endif
#include "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 "MArray2.h"
#include "MSparse.h"
#include "SparseQR.h"
#include "SparseCmplxQR.h"
#include "SparseType.h"
#include "oct-sort.h"
template <class T>
static MSparse<T>
dmsolve_extract (const MSparse<T> &A, const octave_idx_type *Pinv,
const octave_idx_type *Q, octave_idx_type rst,
octave_idx_type rend, octave_idx_type cst,
octave_idx_type cend, octave_idx_type maxnz = -1,
bool lazy = false)
{
octave_idx_type nz = (rend - rst) * (cend - cst);
maxnz = (maxnz < 0 ? A.nnz () : maxnz);
MSparse<T> B (rend - rst, cend - cst, (nz < maxnz ? nz : maxnz));
// Some sparse functions can support lazy indexing (where elements
// in the row are in no particular order), even though octave in
// general can't. For those functions that can using it is a big
// win here in terms of speed.
if (lazy)
{
nz = 0;
for (octave_idx_type j = cst ; j < cend ; j++)
{
octave_idx_type qq = (Q ? Q [j] : j);
B.xcidx (j - cst) = nz;
for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
{
OCTAVE_QUIT;
octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
if (r >= rst && r < rend)
{
B.xdata (nz) = A.data (p);
B.xridx (nz++) = r - rst ;
}
}
}
B.xcidx (cend - cst) = nz ;
}
else
{
OCTAVE_LOCAL_BUFFER (T, X, rend - rst);
octave_sort<octave_idx_type> sort;
octave_idx_type *ri = B.xridx();
nz = 0;
for (octave_idx_type j = cst ; j < cend ; j++)
{
octave_idx_type qq = (Q ? Q [j] : j);
B.xcidx (j - cst) = nz;
for (octave_idx_type p = A.cidx(qq) ; p < A.cidx (qq+1) ; p++)
{
OCTAVE_QUIT;
octave_idx_type r = (Pinv ? Pinv [A.ridx (p)] : A.ridx (p));
if (r >= rst && r < rend)
{
X [r-rst] = A.data (p);
B.xridx (nz++) = r - rst ;
}
}
sort.sort (ri + B.xcidx (j - cst), nz - B.xcidx (j - cst));
for (octave_idx_type p = B.cidx (j - cst); p < nz; p++)
B.xdata (p) = X [B.xridx (p)];
}
B.xcidx (cend - cst) = nz ;
}
return B;
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static MSparse<double>
dmsolve_extract (const MSparse<double> &A, const octave_idx_type *Pinv,
const octave_idx_type *Q, octave_idx_type rst,
octave_idx_type rend, octave_idx_type cst,
octave_idx_type cend, octave_idx_type maxnz,
bool lazy);
static MSparse<Complex>
dmsolve_extract (const MSparse<Complex> &A, const octave_idx_type *Pinv,
const octave_idx_type *Q, octave_idx_type rst,
octave_idx_type rend, octave_idx_type cst,
octave_idx_type cend, octave_idx_type maxnz,
bool lazy);
#endif
template <class T>
static MArray2<T>
dmsolve_extract (const MArray2<T> &m, const octave_idx_type *,
const octave_idx_type *, octave_idx_type r1,
octave_idx_type r2, octave_idx_type c1,
octave_idx_type c2)
{
r2 -= 1;
c2 -= 1;
if (r1 > r2) { octave_idx_type tmp = r1; r1 = r2; r2 = tmp; }
if (c1 > c2) { octave_idx_type tmp = c1; c1 = c2; c2 = tmp; }
octave_idx_type new_r = r2 - r1 + 1;
octave_idx_type new_c = c2 - c1 + 1;
MArray2<T> result (new_r, new_c);
for (octave_idx_type j = 0; j < new_c; j++)
for (octave_idx_type i = 0; i < new_r; i++)
result.xelem (i, j) = m.elem (r1+i, c1+j);
return result;
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static MArray2<double>
dmsolve_extract (const MArray2<double> &m, const octave_idx_type *,
const octave_idx_type *, octave_idx_type r1,
octave_idx_type r2, octave_idx_type c1,
octave_idx_type c2)
static MArray2<Complex>
dmsolve_extract (const MArray2<Complex> &m, const octave_idx_type *,
const octave_idx_type *, octave_idx_type r1,
octave_idx_type r2, octave_idx_type c1,
octave_idx_type c2)
#endif
template <class T>
static void
dmsolve_insert (MArray2<T> &a, const MArray2<T> &b, const octave_idx_type *Q,
octave_idx_type r, octave_idx_type c)
{
T *ax = a.fortran_vec();
const T *bx = b.fortran_vec();
octave_idx_type anr = a.rows();
octave_idx_type nr = b.rows();
octave_idx_type nc = b.cols();
for (octave_idx_type j = 0; j < nc; j++)
{
octave_idx_type aoff = (c + j) * anr;
octave_idx_type boff = j * nr;
for (octave_idx_type i = 0; i < nr; i++)
{
OCTAVE_QUIT;
ax [Q [r + i] + aoff] = bx [i + boff];
}
}
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static void
dmsolve_insert (MArray2<double> &a, const MArray2<double> &b,
const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
static void
dmsolve_insert (MArray2<Complex> &a, const MArray2<Complex> &b,
const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
#endif
template <class T>
static void
dmsolve_insert (MSparse<T> &a, const MSparse<T> &b, const octave_idx_type *Q,
octave_idx_type r, octave_idx_type c)
{
octave_idx_type b_rows = b.rows ();
octave_idx_type b_cols = b.cols ();
octave_idx_type nr = a.rows ();
octave_idx_type nc = a.cols ();
OCTAVE_LOCAL_BUFFER (octave_idx_type, Qinv, nr);
for (octave_idx_type i = 0; i < nr; i++)
Qinv [Q [i]] = i;
// First count the number of elements in the final array
octave_idx_type nel = a.xcidx(c) + b.nnz ();
if (c + b_cols < nc)
nel += a.xcidx(nc) - a.xcidx(c + b_cols);
for (octave_idx_type i = c; i < c + b_cols; i++)
for (octave_idx_type j = a.xcidx(i); j < a.xcidx(i+1); j++)
if (Qinv [a.xridx(j)] < r || Qinv [a.xridx(j)] >= r + b_rows)
nel++;
OCTAVE_LOCAL_BUFFER (T, X, nr);
octave_sort<octave_idx_type> sort;
MSparse<T> tmp (a);
a = MSparse<T> (nr, nc, nel);
octave_idx_type *ri = a.xridx();
for (octave_idx_type i = 0; i < tmp.cidx(c); i++)
{
a.xdata(i) = tmp.xdata(i);
a.xridx(i) = tmp.xridx(i);
}
for (octave_idx_type i = 0; i < c + 1; i++)
a.xcidx(i) = tmp.xcidx(i);
octave_idx_type ii = a.xcidx(c);
for (octave_idx_type i = c; i < c + b_cols; i++)
{
OCTAVE_QUIT;
for (octave_idx_type j = tmp.xcidx(i); j < tmp.xcidx(i+1); j++)
if (Qinv [tmp.xridx(j)] < r || Qinv [tmp.xridx(j)] >= r + b_rows)
{
X [tmp.xridx(j)] = tmp.xdata(j);
a.xridx(ii++) = tmp.xridx(j);
}
OCTAVE_QUIT;
for (octave_idx_type j = b.cidx(i-c); j < b.cidx(i-c+1); j++)
{
X [Q [r + b.ridx(j)]] = b.data(j);
a.xridx(ii++) = Q [r + b.ridx(j)];
}
sort.sort (ri + a.xcidx (i), ii - a.xcidx (i));
for (octave_idx_type p = a.xcidx (i); p < ii; p++)
a.xdata (p) = X [a.xridx (p)];
a.xcidx(i+1) = ii;
}
for (octave_idx_type i = c + b_cols; i < nc; i++)
{
for (octave_idx_type j = tmp.xcidx(i); j < tmp.cidx(i+1); j++)
{
a.xdata(ii) = tmp.xdata(j);
a.xridx(ii++) = tmp.xridx(j);
}
a.xcidx(i+1) = ii;
}
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static void
dmsolve_insert (MSparse<double> &a, const SparseMatrix &b,
const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
static void
dmsolve_insert (MSparse<Complex> &a, const MSparse<Complex> &b,
const octave_idx_type *Q, octave_idx_type r, octave_idx_type c);
#endif
template <class T, class RT>
static void
dmsolve_permute (MArray2<RT> &a, const MArray2<T>& b, const octave_idx_type *p)
{
octave_idx_type b_nr = b.rows ();
octave_idx_type b_nc = b.cols ();
const T *Bx = b.fortran_vec();
a.resize(b_nr, b_nc);
RT *Btx = a.fortran_vec();
for (octave_idx_type j = 0; j < b_nc; j++)
{
octave_idx_type off = j * b_nr;
for (octave_idx_type i = 0; i < b_nr; i++)
{
OCTAVE_QUIT;
Btx [p [i] + off] = Bx [ i + off];
}
}
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static void
dmsolve_permute (MArray2<double> &a, const MArray2<double>& b,
const octave_idx_type *p);
static void
dmsolve_permute (MArray2<Complex> &a, const MArray2<double>& b,
const octave_idx_type *p);
static void
dmsolve_permute (MArray2<Complex> &a, const MArray2<Complex>& b,
const octave_idx_type *p);
#endif
template <class T, class RT>
static void
dmsolve_permute (MSparse<RT> &a, const MSparse<T>& b, const octave_idx_type *p)
{
octave_idx_type b_nr = b.rows ();
octave_idx_type b_nc = b.cols ();
octave_idx_type b_nz = b.nnz ();
octave_idx_type nz = 0;
a = MSparse<RT> (b_nr, b_nc, b_nz);
octave_sort<octave_idx_type> sort;
octave_idx_type *ri = a.xridx();
OCTAVE_LOCAL_BUFFER (RT, X, b_nr);
a.xcidx(0) = 0;
for (octave_idx_type j = 0; j < b_nc; j++)
{
for (octave_idx_type i = b.cidx(j); i < b.cidx(j+1); i++)
{
OCTAVE_QUIT;
octave_idx_type r = p [b.ridx (i)];
X [r] = b.data (i);
a.xridx(nz++) = p [b.ridx (i)];
}
sort.sort (ri + a.xcidx (j), nz - a.xcidx (j));
for (octave_idx_type i = a.cidx (j); i < nz; i++)
{
OCTAVE_QUIT;
a.xdata (i) = X [a.xridx (i)];
}
a.xcidx(j+1) = nz;
}
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
static void
dmsolve_permute (MSparse<double> &a, const MSparse<double>& b,
const octave_idx_type *p);
static void
dmsolve_permute (MSparse<Complex> &a, const MSparse<double>& b,
const octave_idx_type *p);
static void
dmsolve_permute (MSparse<Complex> &a, const MSparse<Complex>& b,
const octave_idx_type *p);
#endif
static void
solve_singularity_warning (double rcond)
{
// Dummy singularity handler so that LU solver doesn't flag
// an error for numerically rank defficient matrices
}
template <class RT, class ST, class T>
RT
dmsolve (const ST &a, const T &b, octave_idx_type &info)
{
octave_idx_type nr = a.rows ();
octave_idx_type nc = a.cols ();
octave_idx_type nz = a.nnz ();
octave_idx_type b_nr = b.rows ();
octave_idx_type b_nc = b.cols ();
RT retval;
if (nr < 1 || nc < 1 || nr != b_nr)
(*current_liboctave_error_handler)
("matrix dimension mismatch in solution of minimum norm problem");
else
{
octave_idx_type nnz_remaining = a.nnz ();
CXSPARSE_DNAME () csm;
csm.m = nr;
csm.n = nc;
csm.x = NULL;
csm.nz = -1;
csm.nzmax = a.nnz ();
// Cast away const on A, with full knowledge that CSparse won't touch it.
// Prevents the methods below making a copy of the data.
csm.p = const_cast<octave_idx_type *>(a.cidx ());
csm.i = const_cast<octave_idx_type *>(a.ridx ());
CXSPARSE_DNAME (d) *dm = CXSPARSE_DNAME(_dmperm) (&csm);
octave_idx_type *p = dm->P;
octave_idx_type *q = dm->Q;
OCTAVE_LOCAL_BUFFER (octave_idx_type, pinv, nr);
for (octave_idx_type i = 0; i < nr; i++)
pinv [p [i]] = i;
RT btmp;
dmsolve_permute (btmp, b, pinv);
SparseType mtyp (SparseType::Full);
info = 0;
double rcond = 0.0;
retval.resize (nc, b_nc);
// Leading over-determined block
if (dm->rr [2] < nr && dm->cc [3] < nc)
{
ST m = dmsolve_extract (a, pinv, q, dm->rr [2], nr, dm->cc [3], nc,
nnz_remaining, true);
nnz_remaining -= m.nnz();
RT mtmp =
qrsolve (m, dmsolve_extract (btmp, NULL, NULL, dm->rr[2], b_nr, 0,
b_nc), info);
dmsolve_insert (retval, mtmp, q, dm->cc [3], 0);
if (dm->rr [2] > 0 && !info && !error_state)
{
m = dmsolve_extract (a, pinv, q, 0, dm->rr [2],
dm->cc [3], nc, nnz_remaining, true);
nnz_remaining -= m.nnz();
RT ctmp = dmsolve_extract (btmp, NULL, NULL, 0,
dm->rr[2], 0, b_nc);
btmp.insert (ctmp - m * mtmp, 0, 0);
}
}
// Structurally non-singular blocks
// XXX FIXME XXX Should use fine Dulmange-Mendelsohn decomposition here.
if (dm->rr [1] < dm->rr [2] && dm->cc [2] < dm->cc [3] &&
!info && !error_state)
{
ST m = dmsolve_extract (a, pinv, q, dm->rr [1], dm->rr [2],
dm->cc [2], dm->cc [3], nnz_remaining, false);
nnz_remaining -= m.nnz();
RT btmp2 = dmsolve_extract (btmp, NULL, NULL, dm->rr [1], dm->rr [2],
0, b_nc);
RT mtmp = m.solve (mtyp, btmp2, info, rcond,
solve_singularity_warning);
if (info != 0)
{
info = 0;
mtmp = qrsolve (m, btmp2, info);
}
dmsolve_insert (retval, mtmp, q, dm->cc [2], 0);
if (dm->rr [1] > 0 && !info && !error_state)
{
m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], dm->cc [2],
dm->cc [3], nnz_remaining, true);
nnz_remaining -= m.nnz();
RT ctmp = dmsolve_extract (btmp, NULL, NULL, 0,
dm->rr[1], 0, b_nc);
btmp.insert (ctmp - m * mtmp, 0, 0);
}
}
// Trailing under-determined block
if (dm->rr [1] > 0 && dm->cc [2] > 0 && !info && !error_state)
{
ST m = dmsolve_extract (a, pinv, q, 0, dm->rr [1], 0,
dm->cc [2], nnz_remaining, true);
RT mtmp =
qrsolve (m, dmsolve_extract(btmp, NULL, NULL, 0, dm->rr [1] , 0,
b_nc), info);
dmsolve_insert (retval, mtmp, q, 0, 0);
}
CXSPARSE_DNAME (_dfree) (dm);
}
return retval;
}
#if !defined (CXX_NEW_FRIEND_TEMPLATE_DECL)
extern Matrix
dmsolve (const SparseMatrix &a, const Matrix &b,
octave_idx_type &info);
extern ComplexMatrix
dmsolve (const SparseMatrix &a, const ComplexMatrix &b,
octave_idx_type &info);
extern ComplexMatrix
dmsolve (const SparseComplexMatrix &a, const Matrix &b,
octave_idx_type &info);
extern ComplexMatrix
dmsolve (const SparseComplexMatrix &a, const ComplexMatrix &b,
octave_idx_type &info);
extern SparseMatrix
dmsolve (const SparseMatrix &a, const SparseMatrix &b,
octave_idx_type &info);
extern SparseComplexMatrix
dmsolve (const SparseMatrix &a, const SparseComplexMatrix &b,
octave_idx_type &info);
extern SparseComplexMatrix
dmsolve (const SparseComplexMatrix &a, const SparseMatrix &b,
octave_idx_type &info);
extern SparseComplexMatrix
dmsolve (const SparseComplexMatrix &a, const SparseComplexMatrix &b,
octave_idx_type &info);
#endif
DEFUN_DLD (dmsol, args, , "x = dmsol (A, b)")
{
octave_value retval;
int nargin = args.length();
octave_idx_type info;
if (nargin !=2)
{
print_usage ("dmsol");
return retval;
}
if (args(0).is_real_type ())
{
const SparseMatrix m = args(0).sparse_matrix_value ();
if (args(1).is_real_type ())
{
if (args(1).class_name () == "sparse")
{
const SparseMatrix b = args(1).sparse_matrix_value();
if (!error_state)
retval = dmsolve <SparseMatrix, SparseMatrix,
SparseMatrix> (m, b, info);
}
else
{
const Matrix b = args(1).matrix_value();
if (!error_state)
retval = dmsolve <Matrix, SparseMatrix, Matrix>(m, b, info);
}
}
else
{
if (args(1).class_name () == "sparse")
{
const SparseComplexMatrix b =
args(1).sparse_complex_matrix_value();
if (!error_state)
retval = dmsolve <SparseComplexMatrix, SparseMatrix,
SparseComplexMatrix> (m, b, info);
}
else
{
const ComplexMatrix b = args(1).complex_matrix_value();
if (!error_state)
retval = dmsolve <ComplexMatrix, SparseMatrix,
ComplexMatrix> (m, b, info);
}
}
}
else
{
const SparseComplexMatrix m = args(0).sparse_complex_matrix_value ();
if (args(1).is_real_type ())
{
if (args(1).class_name () == "sparse")
{
const SparseMatrix b = args(1).sparse_matrix_value();
if (!error_state)
retval = dmsolve <SparseComplexMatrix, SparseComplexMatrix,
SparseMatrix> (m, b, info);
}
else
{
const Matrix b = args(1).matrix_value();
if (!error_state)
retval = dmsolve <ComplexMatrix, SparseComplexMatrix,
Matrix> (m, b, info);
}
}
else
{
if (args(1).class_name () == "sparse")
{
const SparseComplexMatrix b =
args(1).sparse_complex_matrix_value();
if (!error_state)
retval = dmsolve <SparseComplexMatrix, SparseComplexMatrix,
SparseComplexMatrix> (m, b, info);
}
else
{
const ComplexMatrix b = args(1).complex_matrix_value();
if (!error_state)
retval = dmsolve <ComplexMatrix, SparseComplexMatrix,
ComplexMatrix> (m, b, info);
}
}
}
return retval;
}
/*
%!function f(a, sz, feps)
%! b = randn(sz); x = dmsol(a,b);
%! assert (a * x, b, feps);
%! b = randn(sz)+1i*randn(sz); x = dmsol(a,b);
%! assert (a * x, b, feps);
%! b = sprandn(sz(1),sz(2),0.2); x = dmsol(a,b);
%! assert (sparse(a * x), b, feps);
%! b = sprandn(sz(1),sz(2),0.2)+1i*sprandn(sz(1),sz(2),0.2); x = dmsol(a,b);
%! assert (sparse(a * x), b, feps);
%!test
%! a = sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10);
%! ## Test this by forcing matrix_type
%! a = sprandn(10,10,0.2)+speye(10,10); f(a,[10,2],1e-10);
%!test
%! a = 1i*sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10);
%! ## Test this by forcing matrix_type
%! a = 1i*sprandn(10,10,0.2)+speye(10,10); f(a,[10,2],1e-10);
*/
DEFUN_DLD (qrsol, args, nargout, "x = qrsol (A, b)")
{
octave_value retval;
int nargin = args.length();
octave_idx_type info;
if (nargin !=2)
{
print_usage ("qrsol");
return retval;
}
if (args(0).is_real_type ())
{
SparseMatrix m = args(0).sparse_matrix_value ();
if (args(1).is_real_type ())
{
if (args(1).class_name () == "sparse")
retval = qrsolve (m, args(1).sparse_matrix_value(), info);
else
retval = qrsolve (m, args(1).matrix_value(), info);
}
else
{
if (args(1).class_name () == "sparse")
retval = qrsolve (m, args(1).sparse_complex_matrix_value(), info);
else
retval = qrsolve (m, args(1).complex_matrix_value(), info);
}
}
else
{
SparseComplexMatrix m = args(0).sparse_complex_matrix_value ();
if (args(1).is_real_type ())
{
if (args(1).class_name () == "sparse")
retval = qrsolve (m, args(1).sparse_matrix_value(), info);
else
retval = qrsolve (m, args(1).matrix_value(), info);
}
else
{
if (args(1).class_name () == "sparse")
retval = qrsolve (m, args(1).sparse_complex_matrix_value(), info);
else
retval = qrsolve (m, args(1).complex_matrix_value(), info);
}
}
return retval;
}
/*
%!function f(a, sz, feps)
%! b = randn(sz); x = qrsol(a,b);
%! assert (a * x, b, feps);
%! b = randn(sz)+1i*randn(sz); x = qrsol(a,b);
%! assert (a * x, b, feps);
%! b = sprandn(sz(1),sz(2),0.2); x = qrsol(a,b);
%! assert (sparse(a * x), b, feps);
%! b = sprandn(sz(1),sz(2),0.2)+1i*sprandn(sz(1),sz(2),0.2); x = qrsol(a,b);
%! assert (sparse(a * x), b, feps);
%!test
%! a = sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10);
%! ## Test this by forcing matrix_type
%! a = sprandn(10,10,0.2)+speye(10,10); f(a,[10,2],1e-10);
%!test
%! a = 1i*sprandn(10,11,0.2)+speye(10,11); f(a,[10,2],1e-10);
%! ## Test this by forcing matrix_type
%! a = 1i*sprandn(10,10,0.2)+speye(10,10); f(a,[10,2],1e-10);
*/
/*
;;; Local Variables: ***
;;; mode: C++ ***
;;; End: ***
*/
Re: 2.9.5 and 2.1.73, David Bateman, 2006/03/14
- Re: 2.9.5 and 2.1.73,
David Bateman <=
Message not available
Re: 2.9.5 and 2.1.73, Colin Ingram, 2006/03/23
Re: 2.9.5 and 2.1.73, David Bateman, 2006/03/14
Re: 2.9.5 and 2.1.73, David Bateman, 2006/03/14
Re: 2.9.5 and 2.1.73, David Bateman, 2006/03/14