[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
matrix type caching, triangualar and Cholesky solvers
|
From: |
David Bateman |
|
Subject: |
matrix type caching, triangualar and Cholesky solvers |
|
Date: |
Thu, 27 Apr 2006 00:06:50 +0200 |
|
User-agent: |
Mozilla Thunderbird 1.0.6-7.5.20060mdk (X11/20050322) |
Here is a patch that allows whether a full matrix is upper, lower or
positive definite to be probed, cached and the method of solving linear
equations with this matrix with the "/" or "\" operators automatically
chosen based on this information. This seems to be about 40% for
positive definite matrices and about 100 times for triangular matrices
using atlas on my machine. I tried to test all possible cases of the use
of this code before submitting it, but as there is no testing code
whatsoever for the "\" and "/" full operators yet, it wasn't easy to
propose appropriately adapted test code in the absence existing test code.
This should supersede the chol.oct file from octave-forge and the
ov-tri-mat type as well.. John do you want it committed? Does someone
else want to try it out?
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
Index: liboctave/CMatrix.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/CMatrix.cc,v
retrieving revision 1.116
diff -c -r1.116 CMatrix.cc
*** liboctave/CMatrix.cc 24 Apr 2006 19:13:07 -0000 1.116
--- liboctave/CMatrix.cc 26 Apr 2006 21:33:56 -0000
***************
*** 112,117 ****
--- 112,154 ----
const octave_idx_type&, double*, double&,
octave_idx_type&,
Complex*, const octave_idx_type&, double*,
octave_idx_type&);
+ F77_RET_T
+ F77_FUNC (zpotrf, ZPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ Complex*, const octave_idx_type&,
+ octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (zpocon, ZPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ Complex*, const octave_idx_type&, const double&,
+ double&, Complex*, double*,
+ octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (zpotrs, ZPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const octave_idx_type&, const Complex*,
+ const octave_idx_type&, Complex*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (ztrcon, ZTRCON) (F77_CONST_CHAR_ARG_DECL,
F77_CONST_CHAR_ARG_DECL,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const Complex*, const octave_idx_type&, double&,
+ Complex*, double*, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (ztrtrs, ZTRTRS) (F77_CONST_CHAR_ARG_DECL,
F77_CONST_CHAR_ARG_DECL,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const octave_idx_type&, const Complex*,
+ const octave_idx_type&, Complex*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
// Note that the original complex fft routines were not written for
// double complex arguments. They have been modified by adding an
// implicit double precision (a-h,o-z) statement at the beginning of
***************
*** 1491,1645 ****
}
ComplexMatrix
! ComplexMatrix::solve (const Matrix& b) const
{
octave_idx_type info;
double rcond;
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
{
double rcond;
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond)
const
{
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler) const
{
ComplexMatrix tmp (b);
! return solve (tmp, info, rcond, sing_handler);
}
ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b) const
{
octave_idx_type info;
double rcond;
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
{
double rcond;
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double&
rcond) const
{
! return solve (b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double&
rcond,
! solve_singularity_handler sing_handler) const
{
ComplexMatrix retval;
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
! if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch in solution of linear equations");
! else
{
! info = 0;
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
! ComplexMatrix atmp = *this;
! Complex *tmp_data = atmp.fortran_vec ();
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (2 * nc);
! double *prz = rz.fortran_vec ();
! // Calculate the norm of the matrix, for later use.
! double anorm =
atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
! F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision");
! }
! else
! {
! // Now calculate the condition number for non-singular matrix.
! char job = '1';
! F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, prz, info
! F77_CHAR_ARG_LEN (1)));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgecon");
! if (info != 0)
! info = -2;
! volatile double rcond_plus_one = rcond + 1.0;
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! else
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
! octave_idx_type b_nc = b.cols ();
! job = 'N';
! F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! pipvt, result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgetrs");
! }
! }
! }
! }
!
! return retval;
}
ComplexColumnVector
--- 1528,2165 ----
}
ComplexMatrix
! ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler,
! bool calc_cond) const
! {
! ComplexMatrix retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch solution of linear equations");
! else
! {
! volatile int typ = mattype.type ();
!
! if (typ == MatrixType::Permuted_Upper ||
! typ == MatrixType::Upper)
! {
! octave_idx_type b_nc = b.cols ();
! rcond = 1.;
! info = 0;
!
! if (typ == MatrixType::Permuted_Upper)
! {
! (*current_liboctave_error_handler)
! ("Permuted triangular matrix not implemented");
! }
! else
! {
! const Complex *tmp_data = fortran_vec ();
!
! if (calc_cond)
! {
! char norm = '1';
! char uplo = 'U';
! char dia = 'N';
!
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (nc);
! double *prz = rz.fortran_vec ();
!
! F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
! F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, tmp_data, nr, rcond,
! pz, prz, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in ztrcon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
!
! char uplo = 'U';
! char trans = 'N';
! char dia = 'N';
!
! F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&trans, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, b_nc, tmp_data, nr,
! result, nr, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrtrs");
! }
! }
! }
! else
! (*current_liboctave_error_handler) ("incorrect matrix type");
! }
!
! return retval;
! }
!
! ComplexMatrix
! ComplexMatrix::ltsolve (MatrixType &mattype, const ComplexMatrix& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler,
! bool calc_cond) const
! {
! ComplexMatrix retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch solution of linear equations");
! else
! {
! volatile int typ = mattype.type ();
!
! if (typ == MatrixType::Permuted_Lower ||
! typ == MatrixType::Lower)
! {
! octave_idx_type b_nc = b.cols ();
! rcond = 1.;
! info = 0;
!
! if (typ == MatrixType::Permuted_Lower)
! {
! (*current_liboctave_error_handler)
! ("Permuted triangular matrix not implemented");
! }
! else
! {
! const Complex *tmp_data = fortran_vec ();
!
! if (calc_cond)
! {
! char norm = '1';
! char uplo = 'L';
! char dia = 'N';
!
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (nc);
! double *prz = rz.fortran_vec ();
!
! F77_XFCN (ztrcon, ZTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
! F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, tmp_data, nr, rcond,
! pz, prz, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in ztrcon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
!
! char uplo = 'L';
! char trans = 'N';
! char dia = 'N';
!
! F77_XFCN (ztrtrs, ZTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&trans, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, b_nc, tmp_data, nr,
! result, nr, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrtrs");
! }
! }
! }
! else
! (*current_liboctave_error_handler) ("incorrect matrix type");
! }
!
! return retval;
! }
!
! ComplexMatrix
! ComplexMatrix::fsolve (MatrixType &mattype, const ComplexMatrix& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler,
! bool calc_cond) const
! {
! ComplexMatrix retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != nc || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch in solution of linear equations");
! else
! {
! volatile int typ = mattype.type ();
!
! // Calculate the norm of the matrix, for later use.
! double anorm = -1.;
!
! if (typ == MatrixType::Hermitian)
! {
! info = 0;
! char job = 'L';
! ComplexMatrix atmp = *this;
! Complex *tmp_data = atmp.fortran_vec ();
! anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! F77_XFCN (zpotrf, ZPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
! tmp_data, nr, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zpotrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
!
! mattype.mark_as_unsymmetric ();
! typ = MatrixType::Full;
! }
! else
! {
! if (calc_cond)
! {
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (nc);
! double *prz = rz.fortran_vec ();
!
! F77_XFCN (zpocon, ZPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, tmp_data, nr, anorm,
! rcond, pz, prz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zpocon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond =
%g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
!
! octave_idx_type b_nc = b.cols ();
!
! F77_XFCN (zpotrs, ZPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zpotrs");
! }
! else
! {
! mattype.mark_as_unsymmetric ();
! typ = MatrixType::Full;
! }
! }
! }
! }
!
! if (typ == MatrixType::Full)
! {
! info = 0;
!
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
!
! ComplexMatrix atmp = *this;
! Complex *tmp_data = atmp.fortran_vec ();
!
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (2 * nc);
! double *prz = rz.fortran_vec ();
!
! // Calculate the norm of the matrix, for later use.
! if (anorm < 0.)
! anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgetrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision");
!
! mattype.mark_as_rectangular ();
! }
! else
! {
! if (calc_cond)
! {
! // Now calculate the condition number for
! // non-singular matrix.
! char job = '1';
! F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, prz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgecon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond =
%g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
!
! octave_idx_type b_nc = b.cols ();
!
! char job = 'N';
! F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! pipvt, result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgetrs");
! }
! else
! mattype.mark_as_rectangular ();
! }
! }
! }
! }
!
! return retval;
! }
!
! ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const Matrix& b) const
{
octave_idx_type info;
double rcond;
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const Matrix& b,
! octave_idx_type& info) const
{
double rcond;
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
! double& rcond) const
{
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type&
info,
! double& rcond, solve_singularity_handler sing_handler,
! bool singular_fallback) const
{
ComplexMatrix tmp (b);
! return solve (typ, tmp, info, rcond, sing_handler, singular_fallback);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b) const
{
octave_idx_type info;
double rcond;
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
! octave_idx_type& info) const
{
double rcond;
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &typ, const ComplexMatrix& b,
! octave_idx_type& info, double& rcond) const
{
! return solve (typ, b, info, rcond, 0);
}
ComplexMatrix
! ComplexMatrix::solve (MatrixType &mattype, const ComplexMatrix& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler,
! bool singular_fallback) const
{
ComplexMatrix retval;
+ int typ = mattype.type ();
! if (typ == MatrixType::Unknown)
! typ = mattype.type (*this);
! // Only calculate the condition number for LU/Cholesky
! if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
! retval = utsolve (mattype, b, info, rcond, sing_handler, false);
! else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
! retval = ltsolve (mattype, b, info, rcond, sing_handler, false);
! else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
! retval = fsolve (mattype, b, info, rcond, sing_handler, true);
! else if (typ != MatrixType::Rectangular)
{
! (*current_liboctave_error_handler) ("unknown matrix type");
! return ComplexMatrix ();
! }
! // Rectangular or one of the above solvers flags a singular matrix
! if (singular_fallback && mattype.type () == MatrixType::Rectangular)
! {
! octave_idx_type rank;
! retval = lssolve (b, info, rank);
! }
! return retval;
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b) const
! {
! octave_idx_type info;
! double rcond;
! return solve (typ, ComplexColumnVector (b), info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
! octave_idx_type& info) const
! {
! double rcond;
! return solve (typ, ComplexColumnVector (b), info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
! octave_idx_type& info, double& rcond) const
! {
! return solve (typ, ComplexColumnVector (b), info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ColumnVector& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler) const
! {
! return solve (typ, ComplexColumnVector (b), info, rcond, sing_handler);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
! {
! octave_idx_type info;
! double rcond;
! return solve (typ, b, info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
! octave_idx_type& info) const
! {
! double rcond;
! return solve (typ, b, info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
! octave_idx_type& info, double& rcond) const
! {
! return solve (typ, b, info, rcond, 0);
! }
! ComplexColumnVector
! ComplexMatrix::solve (MatrixType &typ, const ComplexColumnVector& b,
! octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler) const
! {
! ComplexMatrix tmp (b);
! return solve (typ, tmp, info, rcond,
sing_handler).column(static_cast<octave_idx_type> (0));
! }
! ComplexMatrix
! ComplexMatrix::solve (const Matrix& b) const
! {
! octave_idx_type info;
! double rcond;
! return solve (b, info, rcond, 0);
! }
! ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info) const
! {
! double rcond;
! return solve (b, info, rcond, 0);
! }
!
! ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond)
const
! {
! return solve (b, info, rcond, 0);
! }
! ComplexMatrix
! ComplexMatrix::solve (const Matrix& b, octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler) const
! {
! ComplexMatrix tmp (b);
! return solve (tmp, info, rcond, sing_handler);
! }
! ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b) const
! {
! octave_idx_type info;
! double rcond;
! return solve (b, info, rcond, 0);
! }
! ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info) const
! {
! double rcond;
! return solve (b, info, rcond, 0);
! }
!
! ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double&
rcond) const
! {
! return solve (b, info, rcond, 0);
! }
!
! ComplexMatrix
! ComplexMatrix::solve (const ComplexMatrix& b, octave_idx_type& info, double&
rcond,
! solve_singularity_handler sing_handler) const
! {
! MatrixType mattype (*this);
! return solve (b, info, rcond, sing_handler);
}
ComplexColumnVector
***************
*** 1658,1670 ****
}
ComplexColumnVector
! ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double&
rcond) const
{
return solve (ComplexColumnVector (b), info, rcond, 0);
}
ComplexColumnVector
! ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info, double&
rcond,
solve_singularity_handler sing_handler) const
{
return solve (ComplexColumnVector (b), info, rcond, sing_handler);
--- 2178,2192 ----
}
ComplexColumnVector
! ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
! double& rcond) const
{
return solve (ComplexColumnVector (b), info, rcond, 0);
}
ComplexColumnVector
! ComplexMatrix::solve (const ColumnVector& b, octave_idx_type& info,
! double& rcond,
solve_singularity_handler sing_handler) const
{
return solve (ComplexColumnVector (b), info, rcond, sing_handler);
***************
*** 1697,1796 ****
double& rcond,
solve_singularity_handler sing_handler) const
{
! ComplexColumnVector retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch in solution of linear equations");
! else
! {
! info = 0;
!
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
!
! ComplexMatrix atmp = *this;
! Complex *tmp_data = atmp.fortran_vec ();
!
! Array<Complex> z (2 * nc);
! Complex *pz = z.fortran_vec ();
! Array<double> rz (2 * nc);
! double *prz = rz.fortran_vec ();
!
! // Calculate the norm of the matrix, for later use.
! double anorm =
atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! F77_XFCN (zgetrf, ZGETRF, (nr, nr, tmp_data, nr, pipvt, info));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler) ("unrecoverable error in zgetrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! else
! {
! // Now calculate the condition number for non-singular matrix.
! char job = '1';
! F77_XFCN (zgecon, ZGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, prz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgecon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! else
! {
! retval = b;
! Complex *result = retval.fortran_vec ();
!
! job = 'N';
! F77_XFCN (zgetrs, ZGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, 1, tmp_data, nr, pipvt,
! result, b.length(), info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in zgetrs");
!
! }
! }
! }
! }
! return retval;
}
ComplexMatrix
--- 2219,2226 ----
double& rcond,
solve_singularity_handler sing_handler) const
{
! MatrixType mattype (*this);
! return solve (mattype, b, info, rcond, sing_handler);
}
ComplexMatrix
Index: liboctave/CMatrix.h
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/CMatrix.h,v
retrieving revision 1.56
diff -c -r1.56 CMatrix.h
*** liboctave/CMatrix.h 27 Mar 2006 22:26:19 -0000 1.56
--- liboctave/CMatrix.h 26 Apr 2006 21:33:56 -0000
***************
*** 26,31 ****
--- 26,32 ----
#include "MArray2.h"
#include "MDiagArray2.h"
+ #include "MatrixType.h"
#include "mx-defs.h"
#include "mx-op-defs.h"
***************
*** 150,155 ****
--- 151,216 ----
ComplexDET determinant (octave_idx_type& info) const;
ComplexDET determinant (octave_idx_type& info, double& rcond, int calc_cond
= 1) const;
+ private:
+ // Upper triangular matrix solvers
+ ComplexMatrix utsolve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ // Lower triangular matrix solvers
+ ComplexMatrix ltsolve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ // Full matrix solvers (umfpack/cholesky)
+ ComplexMatrix fsolve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ public:
+ // Generic interface to solver with no probing of type
+ ComplexMatrix solve (MatrixType &typ, const Matrix& b) const;
+ ComplexMatrix solve (MatrixType &typ, const Matrix& b,
+ octave_idx_type& info) const;
+ ComplexMatrix solve (MatrixType &typ, const Matrix& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexMatrix solve (MatrixType &typ, const Matrix& b, octave_idx_type&
info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool singular_fallback = true) const;
+
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool singular_fallback = true) const;
+
+ ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b) const;
+ ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info) const;
+ ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler) const;
+
+ ComplexColumnVector solve (MatrixType &typ,
+ const ComplexColumnVector& b) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler) const;
+
+ // Generic interface to solver with probing of type
ComplexMatrix solve (const Matrix& b) const;
ComplexMatrix solve (const Matrix& b, octave_idx_type& info) const;
ComplexMatrix solve (const Matrix& b, octave_idx_type& info, double& rcond)
const;
Index: liboctave/Makefile.in
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/Makefile.in,v
retrieving revision 1.215
diff -c -r1.215 Makefile.in
*** liboctave/Makefile.in 17 Apr 2006 05:05:16 -0000 1.215
--- liboctave/Makefile.in 26 Apr 2006 21:33:57 -0000
***************
*** 39,45 ****
Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \
sparse-base-chol.h sparse-dmsolve.cc SparseCmplxCHOL.h \
SparsedbleCHOL.h SparseCmplxQR.h SparseQR.h Sparse-op-defs.h \
! SparseType.h \
int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \
int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \
intNDArray.h
--- 39,45 ----
Sparse.h sparse-base-lu.h SparseCmplxLU.h SparsedbleLU.h \
sparse-base-chol.h sparse-dmsolve.cc SparseCmplxCHOL.h \
SparsedbleCHOL.h SparseCmplxQR.h SparseQR.h Sparse-op-defs.h \
! SparseType.h MatrixType.h \
int8NDArray.h uint8NDArray.h int16NDArray.h uint16NDArray.h \
int32NDArray.h uint32NDArray.h int64NDArray.h uint64NDArray.h \
intNDArray.h
***************
*** 101,107 ****
dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \
MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \
SparseCmplxCHOL.cc SparsedbleCHOL.cc \
! SparseCmplxQR.cc SparseQR.cc SparseType.cc \
int8NDArray.cc uint8NDArray.cc int16NDArray.cc uint16NDArray.cc \
int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc
--- 101,107 ----
dbleSCHUR.cc dbleSVD.cc boolSparse.cc CSparse.cc dSparse.cc \
MSparse.cc Sparse.cc SparseCmplxLU.cc SparsedbleLU.cc \
SparseCmplxCHOL.cc SparsedbleCHOL.cc \
! SparseCmplxQR.cc SparseQR.cc SparseType.cc MatrixType.cc \
int8NDArray.cc uint8NDArray.cc int16NDArray.cc uint16NDArray.cc \
int32NDArray.cc uint32NDArray.cc int64NDArray.cc uint64NDArray.cc
Index: liboctave/dMatrix.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/dMatrix.cc,v
retrieving revision 1.122
diff -c -r1.122 dMatrix.cc
*** liboctave/dMatrix.cc 24 Apr 2006 19:13:07 -0000 1.122
--- liboctave/dMatrix.cc 26 Apr 2006 21:33:58 -0000
***************
*** 108,113 ****
--- 108,148 ----
const octave_idx_type&, double*, double&,
octave_idx_type&,
double*, const octave_idx_type&, octave_idx_type&);
+ F77_RET_T
+ F77_FUNC (dpotrf, DPOTRF) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ double *, const octave_idx_type&,
+ octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dpocon, DPOCON) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ double*, const octave_idx_type&, const double&,
+ double&, double*, octave_idx_type*,
+ octave_idx_type& F77_CHAR_ARG_LEN_DECL);
+ F77_RET_T
+ F77_FUNC (dpotrs, DPOTRS) (F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const octave_idx_type&, const double*,
+ const octave_idx_type&, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL);
+
+ F77_RET_T
+ F77_FUNC (dtrcon, DTRCON) (F77_CONST_CHAR_ARG_DECL,
F77_CONST_CHAR_ARG_DECL,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const double*, const octave_idx_type&, double&,
+ double*, octave_idx_type*, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+ F77_RET_T
+ F77_FUNC (dtrtrs, DTRTRS) (F77_CONST_CHAR_ARG_DECL,
F77_CONST_CHAR_ARG_DECL,
+ F77_CONST_CHAR_ARG_DECL, const octave_idx_type&,
+ const octave_idx_type&, const double*,
+ const octave_idx_type&, double*,
+ const octave_idx_type&, octave_idx_type&
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL
+ F77_CHAR_ARG_LEN_DECL);
+
// Note that the original complex fft routines were not written for
// double complex arguments. They have been modified by adding an
// implicit double precision (a-h,o-z) statement at the beginning of
***************
*** 1154,1182 ****
}
Matrix
! Matrix::solve (const Matrix& b) const
{
! octave_idx_type info;
! double rcond;
! return solve (b, info, rcond, 0);
! }
! Matrix
! Matrix::solve (const Matrix& b, octave_idx_type& info) const
! {
! double rcond;
! return solve (b, info, rcond, 0);
}
Matrix
! Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const
{
! return solve (b, info, rcond, 0);
}
Matrix
! Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond,
! solve_singularity_handler sing_handler) const
{
Matrix retval;
--- 1189,1409 ----
}
Matrix
! Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
! double& rcond, solve_singularity_handler sing_handler,
! bool calc_cond) const
{
! Matrix retval;
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch solution of linear equations");
! else
! {
! volatile int typ = mattype.type ();
!
! if (typ == MatrixType::Permuted_Upper ||
! typ == MatrixType::Upper)
! {
! octave_idx_type b_nc = b.cols ();
! rcond = 1.;
! info = 0;
!
! if (typ == MatrixType::Permuted_Upper)
! {
! (*current_liboctave_error_handler)
! ("Permuted triangular matrix not implemented");
! }
! else
! {
! const double *tmp_data = fortran_vec ();
!
! if (calc_cond)
! {
! char norm = '1';
! char uplo = 'U';
! char dia = 'N';
!
! Array<double> z (3 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
!
! F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
! F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, tmp_data, nr, rcond,
! pz, piz, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrcon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! double *result = retval.fortran_vec ();
!
! char uplo = 'U';
! char trans = 'N';
! char dia = 'N';
!
! F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&trans, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, b_nc, tmp_data, nr,
! result, nr, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrtrs");
! }
! }
! }
! else
! (*current_liboctave_error_handler) ("incorrect matrix type");
! }
!
! return retval;
}
Matrix
! Matrix::ltsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
! double& rcond, solve_singularity_handler sing_handler,
! bool calc_cond) const
{
! Matrix retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != b.rows ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch solution of linear equations");
! else
! {
! volatile int typ = mattype.type ();
!
! if (typ == MatrixType::Permuted_Lower ||
! typ == MatrixType::Lower)
! {
! octave_idx_type b_nc = b.cols ();
! rcond = 1.;
! info = 0;
!
! if (typ == MatrixType::Permuted_Lower)
! {
! (*current_liboctave_error_handler)
! ("Permuted triangular matrix not implemented");
! }
! else
! {
! const double *tmp_data = fortran_vec ();
!
! if (calc_cond)
! {
! char norm = '1';
! char uplo = 'L';
! char dia = 'N';
!
! Array<double> z (3 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
!
! F77_XFCN (dtrcon, DTRCON, (F77_CONST_CHAR_ARG2 (&norm, 1),
! F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, tmp_data, nr, rcond,
! pz, piz, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrcon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! double *result = retval.fortran_vec ();
!
! char uplo = 'L';
! char trans = 'N';
! char dia = 'N';
!
! F77_XFCN (dtrtrs, DTRTRS, (F77_CONST_CHAR_ARG2 (&uplo, 1),
! F77_CONST_CHAR_ARG2 (&trans, 1),
! F77_CONST_CHAR_ARG2 (&dia, 1),
! nr, b_nc, tmp_data, nr,
! result, nr, info
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dtrtrs");
! }
! }
! }
! else
! (*current_liboctave_error_handler) ("incorrect matrix type");
! }
!
! return retval;
}
Matrix
! Matrix::fsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
! double& rcond, solve_singularity_handler sing_handler,
! bool calc_cond) const
{
Matrix retval;
***************
*** 1188,1247 ****
("matrix dimension mismatch solution of linear equations");
else
{
! info = 0;
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
! Matrix atmp = *this;
! double *tmp_data = atmp.fortran_vec ();
! Array<double> z (4 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
! // Calculate the norm of the matrix, for later use.
! double anorm =
atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
! F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
! else
{
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision");
! }
! else
! {
! // Now calculate the condition number for non-singular matrix.
! char job = '1';
! F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, piz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgecon");
!
! if (info != 0)
! info = -2;
! volatile double rcond_plus_one = rcond + 1.0;
! if (rcond_plus_one == 1.0 || xisnan (rcond))
{
info = -2;
--- 1415,1539 ----
("matrix dimension mismatch solution of linear equations");
else
{
! volatile int typ = mattype.type ();
!
! // Calculate the norm of the matrix, for later use.
! double anorm = -1.;
! if (typ == MatrixType::Hermitian)
! {
! info = 0;
! char job = 'L';
! Matrix atmp = *this;
! double *tmp_data = atmp.fortran_vec ();
! anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! F77_XFCN (dpotrf, DPOTRF, (F77_CONST_CHAR_ARG2 (&job, 1), nr,
! tmp_data, nr, info
! F77_CHAR_ARG_LEN (1)));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dpotrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
! {
! info = -2;
! mattype.mark_as_unsymmetric ();
! typ = MatrixType::Full;
! }
! else
! {
! if (calc_cond)
! {
! Array<double> z (3 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
!
! F77_XFCN (dpocon, DPOCON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, tmp_data, nr, anorm,
! rcond, pz, piz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dpocon");
!
! if (info != 0)
! info = -2;
! volatile double rcond_plus_one = rcond + 1.0;
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond =
%g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! double *result = retval.fortran_vec ();
!
! octave_idx_type b_nc = b.cols ();
!
! F77_XFCN (dpotrs, DPOTRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dpotrs");
! }
! else
! {
! mattype.mark_as_unsymmetric ();
! typ = MatrixType::Full;
! }
! }
! }
! }
!
! if (typ == MatrixType::Full)
{
! info = 0;
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
! Matrix atmp = *this;
! double *tmp_data = atmp.fortran_vec ();
! if(anorm < 0.)
! anorm = atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! Array<double> z (4 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
! F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgetrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info != 0)
{
info = -2;
***************
*** 1249,1281 ****
sing_handler (rcond);
else
(*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
}
! else
{
! retval = b;
! double *result = retval.fortran_vec ();
! octave_idx_type b_nc = b.cols ();
! job = 'N';
! F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! pipvt, result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgetrs");
}
}
}
}
return retval;
}
ComplexMatrix
Matrix::solve (const ComplexMatrix& b) const
{
--- 1541,1785 ----
sing_handler (rcond);
else
(*current_liboctave_error_handler)
! ("matrix singular to machine precision");
!
! mattype.mark_as_rectangular ();
}
! else
{
! if (calc_cond)
! {
! // Now calculate the condition number for
! // non-singular matrix.
! char job = '1';
! F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, piz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgecon");
!
! if (info != 0)
! info = -2;
! volatile double rcond_plus_one = rcond + 1.0;
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond =
%g",
! rcond);
! }
! }
!
! if (info == 0)
! {
! retval = b;
! double *result = retval.fortran_vec ();
!
! octave_idx_type b_nc = b.cols ();
!
! char job = 'N';
! F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, b_nc, tmp_data, nr,
! pipvt, result, b.rows(), info
! F77_CHAR_ARG_LEN (1)));
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgetrs");
! }
! else
! mattype.mark_as_rectangular ();
}
}
}
+ else if (typ != MatrixType::Hermitian)
+ (*current_liboctave_error_handler) ("incorrect matrix type");
}
return retval;
}
+ Matrix
+ Matrix::solve (MatrixType &typ, const Matrix& b) const
+ {
+ octave_idx_type info;
+ double rcond;
+ return solve (typ, b, info, rcond, 0);
+ }
+
+ Matrix
+ Matrix::solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond) const
+ {
+ return solve (typ, b, info, rcond, 0);
+ }
+
+ Matrix
+ Matrix::solve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool singular_fallback) const
+ {
+ Matrix retval;
+ int typ = mattype.type ();
+
+ if (typ == MatrixType::Unknown)
+ typ = mattype.type (*this);
+
+ // Only calculate the condition number for LU/Cholesky
+ if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+ retval = utsolve (mattype, b, info, rcond, sing_handler, false);
+ else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+ retval = ltsolve (mattype, b, info, rcond, sing_handler, false);
+ else if (typ == MatrixType::Full || typ == MatrixType::Hermitian)
+ retval = fsolve (mattype, b, info, rcond, sing_handler, true);
+ else if (typ != MatrixType::Rectangular)
+ {
+ (*current_liboctave_error_handler) ("unknown matrix type");
+ return Matrix ();
+ }
+
+ // Rectangular or one of the above solvers flags a singular matrix
+ if (singular_fallback && mattype.type () == MatrixType::Rectangular)
+ {
+ octave_idx_type rank;
+ retval = lssolve (b, info, rank);
+ }
+
+ return retval;
+ }
+
+ ComplexMatrix
+ Matrix::solve (MatrixType &typ, const ComplexMatrix& b) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b);
+ }
+
+ ComplexMatrix
+ Matrix::solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b, info);
+ }
+
+ ComplexMatrix
+ Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
+ double& rcond) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b, info, rcond);
+ }
+
+ ComplexMatrix
+ Matrix::solve (MatrixType &typ, const ComplexMatrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool singular_fallback) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b, info, rcond, sing_handler, singular_fallback);
+ }
+
+ ColumnVector
+ Matrix::solve (MatrixType &typ, const ColumnVector& b) const
+ {
+ octave_idx_type info; double rcond;
+ return solve (typ, b, info, rcond);
+ }
+
+ ColumnVector
+ Matrix::solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info) const
+ {
+ double rcond;
+ return solve (typ, b, info, rcond);
+ }
+
+ ColumnVector
+ Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info,
+ double& rcond) const
+ {
+ return solve (typ, b, info, rcond, 0);
+ }
+
+ ColumnVector
+ Matrix::solve (MatrixType &typ, const ColumnVector& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler) const
+ {
+ Matrix tmp (b);
+ return solve (typ, tmp, info, rcond,
sing_handler).column(static_cast<octave_idx_type> (0));
+ }
+
+ ComplexColumnVector
+ Matrix::solve (MatrixType &typ, const ComplexColumnVector& b) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b);
+ }
+
+ ComplexColumnVector
+ Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b, info);
+ }
+
+ ComplexColumnVector
+ Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve (typ, b, info, rcond);
+ }
+
+ ComplexColumnVector
+ Matrix::solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler) const
+ {
+ ComplexMatrix tmp (*this);
+ return tmp.solve(typ, b, info, rcond, sing_handler);
+ }
+
+ Matrix
+ Matrix::solve (const Matrix& b) const
+ {
+ octave_idx_type info;
+ double rcond;
+ return solve (b, info, rcond, 0);
+ }
+
+ Matrix
+ Matrix::solve (const Matrix& b, octave_idx_type& info) const
+ {
+ double rcond;
+ return solve (b, info, rcond, 0);
+ }
+
+ Matrix
+ Matrix::solve (const Matrix& b, octave_idx_type& info, double& rcond) const
+ {
+ return solve (b, info, rcond, 0);
+ }
+
+ Matrix
+ Matrix::solve (const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler) const
+ {
+ MatrixType mattype (*this);
+ return solve (mattype, b, info, rcond, sing_handler);
+ }
+
ComplexMatrix
Matrix::solve (const ComplexMatrix& b) const
{
***************
*** 1329,1428 ****
Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
solve_singularity_handler sing_handler) const
{
! ColumnVector retval;
!
! octave_idx_type nr = rows ();
! octave_idx_type nc = cols ();
!
! if (nr == 0 || nc == 0 || nr != nc || nr != b.length ())
! (*current_liboctave_error_handler)
! ("matrix dimension mismatch solution of linear equations");
! else
! {
! info = 0;
!
! Array<octave_idx_type> ipvt (nr);
! octave_idx_type *pipvt = ipvt.fortran_vec ();
!
! Matrix atmp = *this;
! double *tmp_data = atmp.fortran_vec ();
!
! Array<double> z (4 * nc);
! double *pz = z.fortran_vec ();
! Array<octave_idx_type> iz (nc);
! octave_idx_type *piz = iz.fortran_vec ();
!
! // Calculate the norm of the matrix, for later use.
! double anorm =
atmp.abs().sum().row(static_cast<octave_idx_type>(0)).max();
!
! F77_XFCN (dgetrf, DGETRF, (nr, nr, tmp_data, nr, pipvt, info));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler) ("unrecoverable error in dgetrf");
! else
! {
! // Throw-away extra info LAPACK gives so as to not change output.
! rcond = 0.0;
! if (info > 0)
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision");
!
! }
! else
! {
! // Now calculate the condition number for non-singular matrix.
! char job = '1';
! F77_XFCN (dgecon, DGECON, (F77_CONST_CHAR_ARG2 (&job, 1),
! nc, tmp_data, nr, anorm,
! rcond, pz, piz, info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgecon");
!
! if (info != 0)
! info = -2;
!
! volatile double rcond_plus_one = rcond + 1.0;
!
! if (rcond_plus_one == 1.0 || xisnan (rcond))
! {
! info = -2;
!
! if (sing_handler)
! sing_handler (rcond);
! else
! (*current_liboctave_error_handler)
! ("matrix singular to machine precision, rcond = %g",
! rcond);
! }
! else
! {
! retval = b;
! double *result = retval.fortran_vec ();
!
! job = 'N';
! F77_XFCN (dgetrs, DGETRS, (F77_CONST_CHAR_ARG2 (&job, 1),
! nr, 1, tmp_data, nr, pipvt,
! result, b.length(), info
! F77_CHAR_ARG_LEN (1)));
!
! if (f77_exception_encountered)
! (*current_liboctave_error_handler)
! ("unrecoverable error in dgetrs");
! }
! }
! }
! }
!
! return retval;
}
ComplexColumnVector
--- 1833,1840 ----
Matrix::solve (const ColumnVector& b, octave_idx_type& info, double& rcond,
solve_singularity_handler sing_handler) const
{
! MatrixType mattype (*this);
! return solve (mattype, b, info, rcond, sing_handler);
}
ComplexColumnVector
Index: liboctave/dMatrix.h
===================================================================
RCS file: /usr/local/cvsroot/octave/liboctave/dMatrix.h,v
retrieving revision 1.62
diff -c -r1.62 dMatrix.h
*** liboctave/dMatrix.h 27 Mar 2006 22:26:20 -0000 1.62
--- liboctave/dMatrix.h 26 Apr 2006 21:33:58 -0000
***************
*** 26,31 ****
--- 26,32 ----
#include "MArray2.h"
#include "MDiagArray2.h"
+ #include "MatrixType.h"
#include "mx-defs.h"
#include "mx-op-defs.h"
***************
*** 122,127 ****
--- 123,184 ----
DET determinant (octave_idx_type& info) const;
DET determinant (octave_idx_type& info, double& rcond, int calc_cond = 1)
const;
+ private:
+ // Upper triangular matrix solvers
+ Matrix utsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ // Lower triangular matrix solvers
+ Matrix ltsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ // Full matrix solvers (lu/cholesky)
+ Matrix fsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
+ public:
+ // Generic interface to solver with no probing of type
+ Matrix solve (MatrixType &typ, const Matrix& b) const;
+ Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info)
const;
+ Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond) const;
+ Matrix solve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool singular_fallback = true) const;
+
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexMatrix solve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool singular_fallback = true) const;
+
+ ColumnVector solve (MatrixType &typ, const ColumnVector& b) const;
+ ColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info) const;
+ ColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info, double& rcond) const;
+ ColumnVector solve (MatrixType &typ, const ColumnVector& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler) const;
+
+ ComplexColumnVector solve (MatrixType &typ,
+ const ComplexColumnVector& b) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond) const;
+ ComplexColumnVector solve (MatrixType &typ, const ComplexColumnVector& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler) const;
+
+ // Generic interface to solver with probing of type
Matrix solve (const Matrix& b) const;
Matrix solve (const Matrix& b, octave_idx_type& info) const;
Matrix solve (const Matrix& b, octave_idx_type& info, double& rcond) const;
***************
*** 148,153 ****
--- 205,211 ----
double& rcond,
solve_singularity_handler sing_handler) const;
+ // Singular solvers
Matrix lssolve (const Matrix& b) const;
Matrix lssolve (const Matrix& b, octave_idx_type& info) const;
Matrix lssolve (const Matrix& b, octave_idx_type& info, octave_idx_type&
rank) const;
Index: src/ov-base-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-base-mat.h,v
retrieving revision 1.37
diff -c -r1.37 ov-base-mat.h
*** src/ov-base-mat.h 13 Apr 2006 13:04:32 -0000 1.37
--- src/ov-base-mat.h 26 Apr 2006 21:33:59 -0000
***************
*** 37,42 ****
--- 37,44 ----
#include "ov-base.h"
#include "ov-typeinfo.h"
+ #include "MatrixType.h"
+
class Octave_map;
class tree_walker;
***************
*** 50,66 ****
public:
octave_base_matrix (void)
! : octave_base_value () { }
octave_base_matrix (const MT& m)
! : octave_base_value (), matrix (m)
{
if (matrix.ndims () == 0)
matrix.resize (dim_vector (0, 0));
}
octave_base_matrix (const octave_base_matrix& m)
! : octave_base_value (), matrix (m.matrix) { }
~octave_base_matrix (void) { }
--- 52,75 ----
public:
octave_base_matrix (void)
! : octave_base_value (), typ (MatrixType ()) { }
octave_base_matrix (const MT& m)
! : octave_base_value (), matrix (m), typ (MatrixType ())
! {
! if (matrix.ndims () == 0)
! matrix.resize (dim_vector (0, 0));
! }
!
! octave_base_matrix (const MT& m, const MatrixType& t)
! : octave_base_value (), matrix (m), typ (t)
{
if (matrix.ndims () == 0)
matrix.resize (dim_vector (0, 0));
}
octave_base_matrix (const octave_base_matrix& m)
! : octave_base_value (), matrix (m.matrix), typ (m.typ) { }
~octave_base_matrix (void) { }
***************
*** 107,112 ****
--- 116,125 ----
octave_value all (int dim = 0) const { return matrix.all (dim); }
octave_value any (int dim = 0) const { return matrix.any (dim); }
+ MatrixType matrix_type (void) const { return typ; }
+ MatrixType matrix_type (const MatrixType& _typ) const
+ { MatrixType ret = typ; typ = _typ; return ret; }
+
bool is_matrix_type (void) const { return true; }
bool is_numeric_type (void) const { return true; }
***************
*** 126,131 ****
--- 139,146 ----
protected:
MT matrix;
+
+ mutable MatrixType typ;
};
#endif
Index: src/ov-bool-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-bool-mat.h,v
retrieving revision 1.42
diff -c -r1.42 ov-bool-mat.h
*** src/ov-bool-mat.h 13 Apr 2006 13:04:33 -0000 1.42
--- src/ov-bool-mat.h 26 Apr 2006 21:33:59 -0000
***************
*** 38,43 ****
--- 38,45 ----
#include "ov-base-mat.h"
#include "ov-typeinfo.h"
+ #include "MatrixType.h"
+
class Octave_map;
class octave_value_list;
***************
*** 59,64 ****
--- 61,69 ----
octave_bool_matrix (const boolMatrix& bm)
: octave_base_matrix<boolNDArray> (bm) { }
+ octave_bool_matrix (const boolMatrix& bm, const MatrixType& t)
+ : octave_base_matrix<boolNDArray> (bm, t) { }
+
octave_bool_matrix (const Array2<bool>& a)
: octave_base_matrix<boolNDArray> (a) { }
Index: src/ov-cx-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-cx-mat.h,v
retrieving revision 1.39
diff -c -r1.39 ov-cx-mat.h
*** src/ov-cx-mat.h 13 Apr 2006 13:04:33 -0000 1.39
--- src/ov-cx-mat.h 26 Apr 2006 21:33:59 -0000
***************
*** 39,44 ****
--- 39,46 ----
#include "ov-base-mat.h"
#include "ov-typeinfo.h"
+ #include "MatrixType.h"
+
class Octave_map;
class octave_value_list;
***************
*** 60,65 ****
--- 62,70 ----
octave_complex_matrix (const ComplexMatrix& m)
: octave_base_matrix<ComplexNDArray> (m) { }
+ octave_complex_matrix (const ComplexMatrix& m, const MatrixType& t)
+ : octave_base_matrix<ComplexNDArray> (m, t) { }
+
octave_complex_matrix (const ArrayN<Complex>& m)
: octave_base_matrix<ComplexNDArray> (ComplexNDArray (m)) { }
Index: src/ov-re-mat.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov-re-mat.h,v
retrieving revision 1.50
diff -c -r1.50 ov-re-mat.h
*** src/ov-re-mat.h 13 Apr 2006 13:04:33 -0000 1.50
--- src/ov-re-mat.h 26 Apr 2006 21:33:59 -0000
***************
*** 40,45 ****
--- 40,47 ----
#include "ov-base-mat.h"
#include "ov-typeinfo.h"
+ #include "MatrixType.h"
+
class Octave_map;
class octave_value_list;
***************
*** 58,63 ****
--- 60,68 ----
octave_matrix (const Matrix& m)
: octave_base_matrix<NDArray> (m) { }
+ octave_matrix (const Matrix& m, const MatrixType& t)
+ : octave_base_matrix<NDArray> (m, t) { }
+
octave_matrix (const NDArray& nda)
: octave_base_matrix<NDArray> (nda) { }
Index: src/ov.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov.cc,v
retrieving revision 1.140
diff -c -r1.140 ov.cc
*** src/ov.cc 24 Apr 2006 19:13:10 -0000 1.140
--- src/ov.cc 26 Apr 2006 21:33:59 -0000
***************
*** 383,390 ****
{
}
! octave_value::octave_value (const Matrix& m)
! : rep (new octave_matrix (m))
{
maybe_mutate ();
}
--- 383,390 ----
{
}
! octave_value::octave_value (const Matrix& m, const MatrixType& t)
! : rep (new octave_matrix (m, t))
{
maybe_mutate ();
}
***************
*** 425,432 ****
maybe_mutate ();
}
! octave_value::octave_value (const ComplexMatrix& m)
! : rep (new octave_complex_matrix (m))
{
maybe_mutate ();
}
--- 425,432 ----
maybe_mutate ();
}
! octave_value::octave_value (const ComplexMatrix& m, const MatrixType& t)
! : rep (new octave_complex_matrix (m, t))
{
maybe_mutate ();
}
***************
*** 466,473 ****
{
}
! octave_value::octave_value (const boolMatrix& bm)
! : rep (new octave_bool_matrix (bm))
{
maybe_mutate ();
}
--- 466,473 ----
{
}
! octave_value::octave_value (const boolMatrix& bm, const MatrixType& t)
! : rep (new octave_bool_matrix (bm, t))
{
maybe_mutate ();
}
Index: src/ov.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/ov.h,v
retrieving revision 1.138
diff -c -r1.138 ov.h
*** src/ov.h 24 Apr 2006 19:13:10 -0000 1.138
--- src/ov.h 26 Apr 2006 21:34:00 -0000
***************
*** 41,46 ****
--- 41,47 ----
#include "oct-time.h"
#include "str-vec.h"
#include "SparseType.h"
+ #include "MatrixType.h"
class Cell;
class streamoff_array;
***************
*** 163,183 ****
octave_value (double d);
octave_value (const ArrayN<octave_value>& a, bool is_cs_list = false);
octave_value (const Cell& c, bool is_cs_list = false);
! octave_value (const Matrix& m);
octave_value (const NDArray& nda);
octave_value (const ArrayN<double>& m);
octave_value (const DiagMatrix& d);
octave_value (const RowVector& v);
octave_value (const ColumnVector& v);
octave_value (const Complex& C);
! octave_value (const ComplexMatrix& m);
octave_value (const ComplexNDArray& cnda);
octave_value (const ArrayN<Complex>& m);
octave_value (const ComplexDiagMatrix& d);
octave_value (const ComplexRowVector& v);
octave_value (const ComplexColumnVector& v);
octave_value (bool b);
! octave_value (const boolMatrix& bm);
octave_value (const boolNDArray& bnda);
octave_value (char c, char type = '"');
octave_value (const char *s, char type = '"');
--- 164,184 ----
octave_value (double d);
octave_value (const ArrayN<octave_value>& a, bool is_cs_list = false);
octave_value (const Cell& c, bool is_cs_list = false);
! octave_value (const Matrix& m, const MatrixType& t = MatrixType());
octave_value (const NDArray& nda);
octave_value (const ArrayN<double>& m);
octave_value (const DiagMatrix& d);
octave_value (const RowVector& v);
octave_value (const ColumnVector& v);
octave_value (const Complex& C);
! octave_value (const ComplexMatrix& m, const MatrixType& t = MatrixType());
octave_value (const ComplexNDArray& cnda);
octave_value (const ArrayN<Complex>& m);
octave_value (const ComplexDiagMatrix& d);
octave_value (const ComplexRowVector& v);
octave_value (const ComplexColumnVector& v);
octave_value (bool b);
! octave_value (const boolMatrix& bm, const MatrixType& t = MatrixType());
octave_value (const boolNDArray& bnda);
octave_value (char c, char type = '"');
octave_value (const char *s, char type = '"');
Index: src/xdiv.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/xdiv.cc,v
retrieving revision 1.31
diff -c -r1.31 xdiv.cc
*** src/xdiv.cc 26 Apr 2005 19:24:35 -0000 1.31
--- src/xdiv.cc 26 Apr 2006 21:34:00 -0000
***************
*** 118,230 ****
// -*- 1 -*-
Matrix
! xdiv (const Matrix& a, const Matrix& b)
{
if (! mx_div_conform (a, b))
return Matrix ();
Matrix atmp = a.transpose ();
Matrix btmp = b.transpose ();
octave_idx_type info;
! if (btmp.rows () == btmp.columns ())
! {
! double rcond = 0.0;
!
! Matrix result
! = btmp.solve (atmp, info, rcond, solve_singularity_warning);
! if (result_ok (info))
! return Matrix (result.transpose ());
! }
!
! octave_idx_type rank;
! Matrix result = btmp.lssolve (atmp, info, rank);
return result.transpose ();
}
// -*- 2 -*-
ComplexMatrix
! xdiv (const Matrix& a, const ComplexMatrix& b)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
Matrix atmp = a.transpose ();
ComplexMatrix btmp = b.hermitian ();
octave_idx_type info;
! if (btmp.rows () == btmp.columns ())
! {
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (atmp, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result.hermitian ();
! }
!
! octave_idx_type rank;
! ComplexMatrix result = btmp.lssolve (atmp, info, rank);
return result.hermitian ();
}
// -*- 3 -*-
ComplexMatrix
! xdiv (const ComplexMatrix& a, const Matrix& b)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
ComplexMatrix atmp = a.hermitian ();
Matrix btmp = b.transpose ();
octave_idx_type info;
! if (btmp.rows () == btmp.columns ())
! {
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (atmp, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result.hermitian ();
! }
!
! octave_idx_type rank;
! ComplexMatrix result = btmp.lssolve (atmp, info, rank);
return result.hermitian ();
}
// -*- 4 -*-
ComplexMatrix
! xdiv (const ComplexMatrix& a, const ComplexMatrix& b)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
ComplexMatrix atmp = a.hermitian ();
ComplexMatrix btmp = b.hermitian ();
octave_idx_type info;
! if (btmp.rows () == btmp.columns ())
! {
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (atmp, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result.hermitian ();
! }
!
! octave_idx_type rank;
! ComplexMatrix result = btmp.lssolve (atmp, info, rank);
return result.hermitian ();
}
--- 118,202 ----
// -*- 1 -*-
Matrix
! xdiv (const Matrix& a, const Matrix& b, MatrixType &typ)
{
if (! mx_div_conform (a, b))
return Matrix ();
Matrix atmp = a.transpose ();
Matrix btmp = b.transpose ();
+ MatrixType btyp = typ.transpose ();
octave_idx_type info;
! double rcond = 0.0;
! Matrix result
! = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
+ typ = btyp.transpose ();
return result.transpose ();
}
// -*- 2 -*-
ComplexMatrix
! xdiv (const Matrix& a, const ComplexMatrix& b, MatrixType &typ)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
Matrix atmp = a.transpose ();
ComplexMatrix btmp = b.hermitian ();
+ MatrixType btyp = typ.transpose ();
octave_idx_type info;
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
+ typ = btyp.transpose ();
return result.hermitian ();
}
// -*- 3 -*-
ComplexMatrix
! xdiv (const ComplexMatrix& a, const Matrix& b, MatrixType &typ)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
ComplexMatrix atmp = a.hermitian ();
Matrix btmp = b.transpose ();
+ MatrixType btyp = typ.transpose ();
octave_idx_type info;
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
+ typ = btyp.transpose ();
return result.hermitian ();
}
// -*- 4 -*-
ComplexMatrix
! xdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType &typ)
{
if (! mx_div_conform (a, b))
return ComplexMatrix ();
ComplexMatrix atmp = a.hermitian ();
ComplexMatrix btmp = b.hermitian ();
+ MatrixType btyp = typ.transpose ();
octave_idx_type info;
! double rcond = 0.0;
! ComplexMatrix result
! = btmp.solve (btyp, atmp, info, rcond, solve_singularity_warning);
+ typ = btyp.transpose ();
return result.hermitian ();
}
***************
*** 385,478 ****
// -*- 1 -*-
Matrix
! xleftdiv (const Matrix& a, const Matrix& b)
{
if (! mx_leftdiv_conform (a, b))
return Matrix ();
octave_idx_type info;
! if (a.rows () == a.columns ())
! {
! double rcond = 0.0;
!
! Matrix result
! = a.solve (b, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result;
! }
!
! octave_idx_type rank;
! return a.lssolve (b, info, rank);
}
// -*- 2 -*-
ComplexMatrix
! xleftdiv (const Matrix& a, const ComplexMatrix& b)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! if (a.rows () == a.columns ())
! {
! double rcond = 0.0;
!
! ComplexMatrix result
! = a.solve (b, info, rcond, solve_singularity_warning);
! if (result_ok (info))
! return result;
! }
!
! octave_idx_type rank;
! return a.lssolve (b, info, rank);
}
// -*- 3 -*-
ComplexMatrix
! xleftdiv (const ComplexMatrix& a, const Matrix& b)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! if (a.rows () == a.columns ())
! {
! double rcond = 0.0;
!
! ComplexMatrix result
! = a.solve (b, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result;
! }
!
! octave_idx_type rank;
! return a.lssolve (b, info, rank);
}
// -*- 4 -*-
ComplexMatrix
! xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! if (a.rows () == a.columns ())
! {
! double rcond = 0.0;
!
! ComplexMatrix result
! = a.solve (b, info, rcond, solve_singularity_warning);
!
! if (result_ok (info))
! return result;
! }
!
! octave_idx_type rank;
! return a.lssolve (b, info, rank);
}
/*
--- 357,407 ----
// -*- 1 -*-
Matrix
! xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ)
{
if (! mx_leftdiv_conform (a, b))
return Matrix ();
octave_idx_type info;
! double rcond = 0.0;
! return a.solve (typ, b, info, rcond, solve_singularity_warning);
}
// -*- 2 -*-
ComplexMatrix
! xleftdiv (const Matrix& a, const ComplexMatrix& b, MatrixType &typ)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! double rcond = 0.0;
! return a.solve (typ, b, info, rcond, solve_singularity_warning);
}
// -*- 3 -*-
ComplexMatrix
! xleftdiv (const ComplexMatrix& a, const Matrix& b, MatrixType &typ)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! double rcond = 0.0;
! return a.solve (typ, b, info, rcond, solve_singularity_warning);
}
// -*- 4 -*-
ComplexMatrix
! xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b, MatrixType &typ)
{
if (! mx_leftdiv_conform (a, b))
return ComplexMatrix ();
octave_idx_type info;
! double rcond = 0.0;
! return a.solve (typ, b, info, rcond, solve_singularity_warning);
}
/*
Index: src/xdiv.h
===================================================================
RCS file: /usr/local/cvsroot/octave/src/xdiv.h,v
retrieving revision 1.13
diff -c -r1.13 xdiv.h
*** src/xdiv.h 26 Apr 2005 19:24:35 -0000 1.13
--- src/xdiv.h 26 Apr 2006 21:34:00 -0000
***************
*** 25,30 ****
--- 25,31 ----
#define octave_xdiv_h 1
#include "oct-cmplx.h"
+ #include "MatrixType.h"
class Matrix;
class ComplexMatrix;
***************
*** 32,41 ****
class NDArray;
class ComplexNDArray;
! extern Matrix xdiv (const Matrix& a, const Matrix& b);
! extern ComplexMatrix xdiv (const Matrix& a, const ComplexMatrix& b);
! extern ComplexMatrix xdiv (const ComplexMatrix& a, const Matrix& b);
! extern ComplexMatrix xdiv (const ComplexMatrix& a, const ComplexMatrix& b);
extern Matrix x_el_div (double a, const Matrix& b);
extern ComplexMatrix x_el_div (double a, const ComplexMatrix& b);
--- 33,45 ----
class NDArray;
class ComplexNDArray;
! extern Matrix xdiv (const Matrix& a, const Matrix& b, MatrixType &typ);
! extern ComplexMatrix xdiv (const Matrix& a, const ComplexMatrix& b,
! MatrixType &typ);
! extern ComplexMatrix xdiv (const ComplexMatrix& a, const Matrix& b,
! MatrixType &typ);
! extern ComplexMatrix xdiv (const ComplexMatrix& a, const ComplexMatrix& b,
! MatrixType &typ);
extern Matrix x_el_div (double a, const Matrix& b);
extern ComplexMatrix x_el_div (double a, const ComplexMatrix& b);
***************
*** 47,56 ****
extern ComplexNDArray x_el_div (const Complex a, const NDArray& b);
extern ComplexNDArray x_el_div (const Complex a, const ComplexNDArray& b);
! extern Matrix xleftdiv (const Matrix& a, const Matrix& b);
! extern ComplexMatrix xleftdiv (const Matrix& a, const ComplexMatrix& b);
! extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const Matrix& b);
! extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const ComplexMatrix&
b);
#endif
--- 51,63 ----
extern ComplexNDArray x_el_div (const Complex a, const NDArray& b);
extern ComplexNDArray x_el_div (const Complex a, const ComplexNDArray& b);
! extern Matrix xleftdiv (const Matrix& a, const Matrix& b, MatrixType &typ);
! extern ComplexMatrix xleftdiv (const Matrix& a, const ComplexMatrix& b,
! MatrixType &typ);
! extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const Matrix& b,
! MatrixType &typ);
! extern ComplexMatrix xleftdiv (const ComplexMatrix& a, const ComplexMatrix& b,
! MatrixType &typ);
#endif
Index: src/DLD-FUNCTIONS/matrix_type.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/DLD-FUNCTIONS/matrix_type.cc,v
retrieving revision 1.13
diff -c -r1.13 matrix_type.cc
*** src/DLD-FUNCTIONS/matrix_type.cc 24 Apr 2006 19:13:11 -0000 1.13
--- src/DLD-FUNCTIONS/matrix_type.cc 26 Apr 2006 21:34:00 -0000
***************
*** 30,38 ****
--- 30,41 ----
#include "ov.h"
#include "defun-dld.h"
#include "error.h"
+ #include "ov-re-mat.h"
+ #include "ov-cx-mat.h"
#include "ov-re-sparse.h"
#include "ov-cx-sparse.h"
#include "SparseType.h"
+ #include "MatrixType.h"
DEFUN_DLD (matrix_type, args, ,
"-*- texinfo -*-\n\
***************
*** 300,306 ****
}
}
else
! error ("matrix_type: Only sparse matrices treated at the moment");
}
return retval;
--- 303,456 ----
}
}
else
! {
! if (nargin == 1)
! {
! MatrixType mattyp;
!
! if (args(0).type_name () == "complex matrix" )
! {
! const octave_complex_matrix& rep
! = dynamic_cast<const octave_complex_matrix&>
(args(0).get_rep ());
!
! mattyp = rep.matrix_type ();
!
! if (mattyp.is_unknown ())
! {
! ComplexMatrix m = args(0).complex_matrix_value ();
! if (!error_state)
! {
! mattyp = MatrixType (m);
! rep.matrix_type (mattyp);
! }
! }
! }
! else
! {
! const octave_matrix& rep
! = dynamic_cast<const octave_matrix&> (args(0).get_rep ());
!
! mattyp = rep.matrix_type ();
!
! if (mattyp.is_unknown ())
! {
! Matrix m = args(0).matrix_value ();
! if (!error_state)
! {
! mattyp = MatrixType (m);
! rep.matrix_type (mattyp);
! }
! }
! }
!
! int typ = mattyp.type ();
!
! if (typ == MatrixType::Upper)
! retval = octave_value ("Upper");
! else if (typ == MatrixType::Permuted_Upper)
! retval = octave_value ("Permuted Upper");
! else if (typ == MatrixType::Lower)
! retval = octave_value ("Lower");
! else if (typ == MatrixType::Permuted_Lower)
! retval = octave_value ("Permuted Lower");
! else if (typ == MatrixType::Hermitian)
! retval = octave_value ("Positive Definite");
! else if (typ == MatrixType::Rectangular)
! {
! if (args(0).rows() == args(0).columns())
! retval = octave_value ("Singular");
! else
! retval = octave_value ("Rectangular");
! }
! else if (typ == MatrixType::Full)
! retval = octave_value ("Full");
! else
! // This should never happen!!!
! retval = octave_value ("Unknown");
! }
! else
! {
! // Ok, we're changing the matrix type
! std::string str_typ = args(1).string_value ();
!
! // FIXME -- why do I have to explicitly call the constructor?
! MatrixType mattyp = MatrixType ();
!
! if (error_state)
! error ("Matrix type must be a string");
! else
! {
! // Use STL function to convert to lower case
! std::transform (str_typ.begin (), str_typ.end (),
! str_typ.begin (), tolower);
!
! if (str_typ == "upper")
! mattyp.mark_as_upper_triangular ();
! else if (str_typ == "lower")
! mattyp.mark_as_lower_triangular ();
! else if (str_typ == "positive definite")
! {
! mattyp.mark_as_full ();
! mattyp.mark_as_symmetric ();
! }
! else if (str_typ == "singular")
! mattyp.mark_as_rectangular ();
! else if (str_typ == "full")
! mattyp.mark_as_full ();
! else if (str_typ == "unknown")
! mattyp.invalidate_type ();
! else
! error ("matrix_type: Unknown matrix type %s",
str_typ.c_str());
!
! if (! error_state)
! {
! if (nargin == 3 && (str_typ == "upper"
! || str_typ == "lower"))
! {
! const ColumnVector perm =
! ColumnVector (args (2).vector_value ());
!
! if (error_state)
! error ("matrix_type: Invalid permutation vector");
! else
! {
! octave_idx_type len = perm.length ();
! dim_vector dv = args(0).dims ();
!
! if (len != dv(0))
! error ("matrix_type: Invalid permutation
vector");
! else
! {
! OCTAVE_LOCAL_BUFFER (octave_idx_type, p, len);
!
! for (octave_idx_type i = 0; i < len; i++)
! p[i] = static_cast<octave_idx_type> (perm
(i)) - 1;
!
! if (str_typ == "upper")
! mattyp.mark_as_permuted (len, p);
! else
! mattyp.mark_as_permuted (len, p);
! }
! }
! }
! else if (nargin != 2)
! error ("matrix_type: Invalid number of arguments");
!
! if (! error_state)
! {
! // Set the matrix type
! if (args(0).type_name () == "complex matrix" )
! retval =
! octave_value (args(0).complex_matrix_value (),
! mattyp);
! else
! retval = octave_value (args(0).matrix_value (),
! mattyp);
! }
! }
! }
! }
! }
}
return retval;
***************
*** 393,398 ****
--- 543,568 ----
%! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular");
%! assert(matrix_type(a),"Singular");
+ %!assert(matrix_type(triu(ones(10,10))),"Upper");
+ %!assert(matrix_type(triu(ones(10,10),-1)),"Full");
+ %!assert(matrix_type(tril(ones(10,10))),"Lower");
+ %!assert(matrix_type(tril(ones(10,10),1)),"Full");
+ %!assert(matrix_type(10*eye(10,10) + ones(10,10)), "Positive Definite");
+ %!assert(matrix_type(ones(11,10)),"Rectangular")
+ %!test
+ %! a = matrix_type(ones(10,10),"Singular");
+ %! assert(matrix_type(a),"Singular");
+
+ %!assert(matrix_type(triu(1i*ones(10,10))),"Upper");
+ %!assert(matrix_type(triu(1i*ones(10,10),-1)),"Full");
+ %!assert(matrix_type(tril(1i*ones(10,10))),"Lower");
+ %!assert(matrix_type(tril(1i*ones(10,10),1)),"Full");
+ %!assert(matrix_type(10*eye(10,10) + 1i*triu(ones(10,10),1)
-1i*tril(ones(10,10),-1)), "Positive Definite");
+ %!assert(matrix_type(ones(11,10)),"Rectangular")
+ %!test
+ %! a = matrix_type(ones(10,10),"Singular");
+ %! assert(matrix_type(a),"Singular");
+
*/
/*
Index: src/OPERATORS/op-cm-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-cm.cc,v
retrieving revision 1.15
diff -c -r1.15 op-cm-cm.cc
*** src/OPERATORS/op-cm-cm.cc 26 Apr 2005 19:24:35 -0000 1.15
--- src/OPERATORS/op-cm-cm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 75,81 ****
DEFNDBINOP_OP (sub, complex_matrix, complex_matrix, complex_array,
complex_array, -)
DEFBINOP_OP (mul, complex_matrix, complex_matrix, *)
! DEFBINOP_FN (div, complex_matrix, complex_matrix, xdiv)
DEFBINOPX (pow, complex_matrix, complex_matrix)
{
--- 75,92 ----
DEFNDBINOP_OP (sub, complex_matrix, complex_matrix, complex_array,
complex_array, -)
DEFBINOP_OP (mul, complex_matrix, complex_matrix, *)
!
! DEFBINOP (div, complex_matrix, complex_matrix)
! {
! CAST_BINOP_ARGS (const octave_complex_matrix&, const
octave_complex_matrix&);
! MatrixType typ = v2.matrix_type ();
!
! ComplexMatrix ret = xdiv (v1.complex_matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
! }
DEFBINOPX (pow, complex_matrix, complex_matrix)
{
***************
*** 83,89 ****
return octave_value ();
}
! DEFBINOP_FN (ldiv, complex_matrix, complex_matrix, xleftdiv)
DEFNDBINOP_FN (lt, complex_matrix, complex_matrix, complex_array,
complex_array, mx_el_lt)
DEFNDBINOP_FN (le, complex_matrix, complex_matrix, complex_array,
complex_array, mx_el_le)
--- 94,110 ----
return octave_value ();
}
! DEFBINOP (ldiv, complex_matrix, complex_matrix)
! {
! CAST_BINOP_ARGS (const octave_complex_matrix&, const
octave_complex_matrix&);
! MatrixType typ = v1.matrix_type ();
!
! ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
! }
DEFNDBINOP_FN (lt, complex_matrix, complex_matrix, complex_array,
complex_array, mx_el_lt)
DEFNDBINOP_FN (le, complex_matrix, complex_matrix, complex_array,
complex_array, mx_el_le)
Index: src/OPERATORS/op-cm-cs.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-cs.cc,v
retrieving revision 1.13
diff -c -r1.13 op-cm-cs.cc
*** src/OPERATORS/op-cm-cs.cc 26 Apr 2005 19:24:35 -0000 1.13
--- src/OPERATORS/op-cm-cs.cc 26 Apr 2006 21:34:01 -0000
***************
*** 61,68 ****
ComplexMatrix m1 = v1.complex_matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
! return octave_value (xleftdiv (m1, m2));
}
DEFNDBINOP_FN (lt, complex_matrix, complex, complex_array, complex, mx_el_lt)
--- 61,71 ----
ComplexMatrix m1 = v1.complex_matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (m1, m2, typ);
! v1.matrix_type (typ);
! return ret;
}
DEFNDBINOP_FN (lt, complex_matrix, complex, complex_array, complex, mx_el_lt)
Index: src/OPERATORS/op-cm-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-m.cc,v
retrieving revision 1.14
diff -c -r1.14 op-cm-m.cc
*** src/OPERATORS/op-cm-m.cc 26 Apr 2005 19:24:35 -0000 1.14
--- src/OPERATORS/op-cm-m.cc 26 Apr 2006 21:34:01 -0000
***************
*** 46,52 ****
DEFNDBINOP_OP (sub, complex_matrix, matrix, complex_array, array, -)
DEFBINOP_OP (mul, complex_matrix, matrix, *)
! DEFBINOP_FN (div, complex_matrix, matrix, xdiv)
DEFBINOPX (pow, complex_matrix, matrix)
{
--- 46,64 ----
DEFNDBINOP_OP (sub, complex_matrix, matrix, complex_array, array, -)
DEFBINOP_OP (mul, complex_matrix, matrix, *)
!
! DEFBINOP (div, complex_matrix, matrix)
! {
! CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_matrix&);
! MatrixType typ = v2.matrix_type ();
!
! ComplexMatrix ret = xdiv (v1.complex_matrix_value (),
! v2.matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
! }
!
DEFBINOPX (pow, complex_matrix, matrix)
{
***************
*** 54,60 ****
return octave_value ();
}
! DEFBINOP_FN (ldiv, complex_matrix, matrix, xleftdiv)
DEFNDBINOP_FN (lt, complex_matrix, matrix, complex_array, array, mx_el_lt)
DEFNDBINOP_FN (le, complex_matrix, matrix, complex_array, array, mx_el_le)
--- 66,82 ----
return octave_value ();
}
! DEFBINOP (ldiv, complex_matrix, matrix)
! {
! CAST_BINOP_ARGS (const octave_complex_matrix&, const octave_matrix&);
! MatrixType typ = v1.matrix_type ();
!
! ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (),
! v2.matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
! }
DEFNDBINOP_FN (lt, complex_matrix, matrix, complex_array, array, mx_el_lt)
DEFNDBINOP_FN (le, complex_matrix, matrix, complex_array, array, mx_el_le)
Index: src/OPERATORS/op-cm-s.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-s.cc,v
retrieving revision 1.14
diff -c -r1.14 op-cm-s.cc
*** src/OPERATORS/op-cm-s.cc 26 Apr 2005 19:24:35 -0000 1.14
--- src/OPERATORS/op-cm-s.cc 26 Apr 2006 21:34:01 -0000
***************
*** 65,72 ****
ComplexMatrix m1 = v1.complex_matrix_value ();
Matrix m2 = v2.matrix_value ();
! return octave_value (xleftdiv (m1, m2));
}
DEFNDBINOP_FN (lt, complex_matrix, scalar, complex_array, scalar, mx_el_lt)
--- 65,76 ----
ComplexMatrix m1 = v1.complex_matrix_value ();
Matrix m2 = v2.matrix_value ();
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (m1, m2, typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFNDBINOP_FN (lt, complex_matrix, scalar, complex_array, scalar, mx_el_lt)
Index: src/OPERATORS/op-cm-scm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-scm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-cm-scm.cc
*** src/OPERATORS/op-cm-scm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-cm-scm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 69,77 ****
{
CAST_BINOP_ARGS (const octave_complex_matrix&,
const octave_sparse_complex_matrix&);
! return xleftdiv (v1.complex_matrix_value (),
! v2.complex_matrix_value ());
}
DEFBINOP_FN (lt, complex_matrix, sparse_complex_matrix, mx_el_lt)
--- 69,81 ----
{
CAST_BINOP_ARGS (const octave_complex_matrix&,
const octave_sparse_complex_matrix&);
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (lt, complex_matrix, sparse_complex_matrix, mx_el_lt)
Index: src/OPERATORS/op-cm-sm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cm-sm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-cm-sm.cc
*** src/OPERATORS/op-cm-sm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-cm-sm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 68,75 ****
{
CAST_BINOP_ARGS (const octave_complex_matrix&,
const octave_sparse_matrix&);
! return xleftdiv (v1.complex_matrix_value (), v2.matrix_value ());
}
DEFBINOP_FN (lt, complex_matrix, sparse_matrix, mx_el_lt)
--- 68,80 ----
{
CAST_BINOP_ARGS (const octave_complex_matrix&,
const octave_sparse_matrix&);
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (v1.complex_matrix_value (),
! v2.matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (lt, complex_matrix, sparse_matrix, mx_el_lt)
Index: src/OPERATORS/op-cs-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cs-cm.cc,v
retrieving revision 1.13
diff -c -r1.13 op-cs-cm.cc
*** src/OPERATORS/op-cs-cm.cc 26 Apr 2005 19:24:35 -0000 1.13
--- src/OPERATORS/op-cs-cm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 47,54 ****
ComplexMatrix m1 = v1.complex_matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
! return octave_value (xdiv (m1, m2));
}
DEFBINOP_FN (pow, complex, complex_matrix, xpow)
--- 47,58 ----
ComplexMatrix m1 = v1.complex_matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (m1, m2, typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (pow, complex, complex_matrix, xpow)
Index: src/OPERATORS/op-cs-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-cs-m.cc,v
retrieving revision 1.14
diff -c -r1.14 op-cs-m.cc
*** src/OPERATORS/op-cs-m.cc 26 Apr 2005 19:24:35 -0000 1.14
--- src/OPERATORS/op-cs-m.cc 26 Apr 2006 21:34:01 -0000
***************
*** 53,60 ****
ComplexMatrix m1 = v1.complex_matrix_value ();
Matrix m2 = v2.matrix_value ();
! return octave_value (xdiv (m1, m2));
}
DEFBINOP_FN (pow, complex, matrix, xpow)
--- 53,64 ----
ComplexMatrix m1 = v1.complex_matrix_value ();
Matrix m2 = v2.matrix_value ();
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (m1, m2, typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (pow, complex, matrix, xpow)
Index: src/OPERATORS/op-m-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-cm.cc,v
retrieving revision 1.14
diff -c -r1.14 op-m-cm.cc
*** src/OPERATORS/op-m-cm.cc 26 Apr 2005 19:24:35 -0000 1.14
--- src/OPERATORS/op-m-cm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 46,52 ****
DEFNDBINOP_OP (sub, matrix, complex_matrix, array, complex_array, -)
DEFBINOP_OP (mul, matrix, complex_matrix, *)
! DEFBINOP_FN (div, matrix, complex_matrix, xdiv)
DEFBINOPX (pow, matrix, complex_matrix)
{
--- 46,63 ----
DEFNDBINOP_OP (sub, matrix, complex_matrix, array, complex_array, -)
DEFBINOP_OP (mul, matrix, complex_matrix, *)
!
! DEFBINOP (div, matrix, complex_matrix)
! {
! CAST_BINOP_ARGS (const octave_matrix&, const octave_complex_matrix&);
! MatrixType typ = v2.matrix_type ();
!
! ComplexMatrix ret = xdiv (v1.matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
! }
DEFBINOPX (pow, matrix, complex_matrix)
{
***************
*** 54,60 ****
return octave_value ();
}
! DEFBINOP_FN (ldiv, matrix, complex_matrix, xleftdiv)
DEFNDBINOP_FN (lt, matrix, complex_matrix, array, complex_array, mx_el_lt)
DEFNDBINOP_FN (le, matrix, complex_matrix, array, complex_array, mx_el_le)
--- 65,81 ----
return octave_value ();
}
! DEFBINOP (ldiv, matrix, complex_matrix)
! {
! CAST_BINOP_ARGS (const octave_matrix&, const octave_complex_matrix&);
! MatrixType typ = v1.matrix_type ();
!
! ComplexMatrix ret = xleftdiv (v1.matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
! }
DEFNDBINOP_FN (lt, matrix, complex_matrix, array, complex_array, mx_el_lt)
DEFNDBINOP_FN (le, matrix, complex_matrix, array, complex_array, mx_el_le)
Index: src/OPERATORS/op-m-cs.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-cs.cc,v
retrieving revision 1.14
diff -c -r1.14 op-m-cs.cc
*** src/OPERATORS/op-m-cs.cc 26 Apr 2005 19:24:35 -0000 1.14
--- src/OPERATORS/op-m-cs.cc 26 Apr 2006 21:34:01 -0000
***************
*** 67,74 ****
Matrix m1 = v1.matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
! return octave_value (xleftdiv (m1, m2));
}
DEFNDBINOP_FN (lt, matrix, complex, array, complex, mx_el_lt)
--- 67,78 ----
Matrix m1 = v1.matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (m1, m2, typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFNDBINOP_FN (lt, matrix, complex, array, complex, mx_el_lt)
Index: src/OPERATORS/op-m-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-m.cc,v
retrieving revision 1.19
diff -c -r1.19 op-m-m.cc
*** src/OPERATORS/op-m-m.cc 26 Apr 2005 19:24:35 -0000 1.19
--- src/OPERATORS/op-m-m.cc 26 Apr 2006 21:34:01 -0000
***************
*** 62,68 ****
DEFNDBINOP_OP (sub, matrix, matrix, array, array, -)
DEFBINOP_OP (mul, matrix, matrix, *)
! DEFBINOP_FN (div, matrix, matrix, xdiv)
DEFBINOPX (pow, matrix, matrix)
{
--- 62,78 ----
DEFNDBINOP_OP (sub, matrix, matrix, array, array, -)
DEFBINOP_OP (mul, matrix, matrix, *)
!
! DEFBINOP (div, matrix, matrix)
! {
! CAST_BINOP_ARGS (const octave_matrix&, const octave_matrix&);
! MatrixType typ = v2.matrix_type ();
!
! Matrix ret = xdiv (v1.matrix_value (), v2.matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
! }
DEFBINOPX (pow, matrix, matrix)
{
***************
*** 70,76 ****
return octave_value ();
}
! DEFBINOP_FN (ldiv, matrix, matrix, xleftdiv)
DEFNDBINOP_FN (lt, matrix, matrix, array, array, mx_el_lt)
DEFNDBINOP_FN (le, matrix, matrix, array, array, mx_el_le)
--- 80,95 ----
return octave_value ();
}
! DEFBINOP (ldiv, matrix, matrix)
! {
! CAST_BINOP_ARGS (const octave_matrix&, const octave_matrix&);
! MatrixType typ = v1.matrix_type ();
!
! Matrix ret = xleftdiv (v1.matrix_value (), v2.matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
! }
DEFNDBINOP_FN (lt, matrix, matrix, array, array, mx_el_lt)
DEFNDBINOP_FN (le, matrix, matrix, array, array, mx_el_le)
Index: src/OPERATORS/op-m-s.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-s.cc,v
retrieving revision 1.12
diff -c -r1.12 op-m-s.cc
*** src/OPERATORS/op-m-s.cc 26 Apr 2005 19:24:35 -0000 1.12
--- src/OPERATORS/op-m-s.cc 26 Apr 2006 21:34:01 -0000
***************
*** 61,68 ****
Matrix m1 = v1.matrix_value ();
Matrix m2 = v2.matrix_value ();
! return octave_value (xleftdiv (m1, m2));
}
DEFNDBINOP_FN (lt, matrix, scalar, array, scalar, mx_el_lt)
--- 61,72 ----
Matrix m1 = v1.matrix_value ();
Matrix m2 = v2.matrix_value ();
+ MatrixType typ = v1.matrix_type ();
! Matrix ret = xleftdiv (m1, m2, typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFNDBINOP_FN (lt, matrix, scalar, array, scalar, mx_el_lt)
Index: src/OPERATORS/op-m-scm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-scm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-m-scm.cc
*** src/OPERATORS/op-m-scm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-m-scm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 69,76 ****
{
CAST_BINOP_ARGS (const octave_matrix&,
const octave_sparse_complex_matrix&);
! return xleftdiv (v1.matrix_value (), v2.complex_matrix_value ());
}
DEFBINOP_FN (lt, matrix, sparse_complex_matrix, mx_el_lt)
--- 69,81 ----
{
CAST_BINOP_ARGS (const octave_matrix&,
const octave_sparse_complex_matrix&);
+ MatrixType typ = v1.matrix_type ();
! ComplexMatrix ret = xleftdiv (v1.matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (lt, matrix, sparse_complex_matrix, mx_el_lt)
Index: src/OPERATORS/op-m-sm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-m-sm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-m-sm.cc
*** src/OPERATORS/op-m-sm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-m-sm.cc 26 Apr 2006 21:34:01 -0000
***************
*** 65,72 ****
DEFBINOP (ldiv, matrix, sparse_matrix)
{
CAST_BINOP_ARGS (const octave_matrix&, const octave_sparse_matrix&);
! return xleftdiv (v1.matrix_value (), v2.matrix_value ());
}
DEFBINOP_FN (lt, matrix, sparse_matrix, mx_el_lt)
--- 65,76 ----
DEFBINOP (ldiv, matrix, sparse_matrix)
{
CAST_BINOP_ARGS (const octave_matrix&, const octave_sparse_matrix&);
+ MatrixType typ = v1.matrix_type ();
! Matrix ret = xleftdiv (v1.matrix_value (), v2.matrix_value (), typ);
!
! v1.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (lt, matrix, sparse_matrix, mx_el_lt)
Index: src/OPERATORS/op-s-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-s-cm.cc,v
retrieving revision 1.15
diff -c -r1.15 op-s-cm.cc
*** src/OPERATORS/op-s-cm.cc 26 Apr 2005 19:24:35 -0000 1.15
--- src/OPERATORS/op-s-cm.cc 26 Apr 2006 21:34:02 -0000
***************
*** 53,60 ****
Matrix m1 = v1.matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
! return octave_value (xdiv (m1, m2));
}
DEFBINOP_FN (pow, scalar, complex_matrix, xpow)
--- 53,64 ----
Matrix m1 = v1.matrix_value ();
ComplexMatrix m2 = v2.complex_matrix_value ();
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (m1, m2, typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (pow, scalar, complex_matrix, xpow)
Index: src/OPERATORS/op-s-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-s-m.cc,v
retrieving revision 1.12
diff -c -r1.12 op-s-m.cc
*** src/OPERATORS/op-s-m.cc 26 Apr 2005 19:24:35 -0000 1.12
--- src/OPERATORS/op-s-m.cc 26 Apr 2006 21:34:02 -0000
***************
*** 47,54 ****
Matrix m1 = v1.matrix_value ();
Matrix m2 = v2.matrix_value ();
! return octave_value (xdiv (m1, m2));
}
DEFBINOP_FN (pow, scalar, matrix, xpow)
--- 47,58 ----
Matrix m1 = v1.matrix_value ();
Matrix m2 = v2.matrix_value ();
+ MatrixType typ = v2.matrix_type ();
! Matrix ret = xdiv (m1, m2, typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOP_FN (pow, scalar, matrix, xpow)
Index: src/OPERATORS/op-scm-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-scm-cm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-scm-cm.cc
*** src/OPERATORS/op-scm-cm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-scm-cm.cc 26 Apr 2006 21:34:02 -0000
***************
*** 49,56 ****
{
CAST_BINOP_ARGS (const octave_sparse_complex_matrix&,
const octave_complex_matrix&);
! return xdiv (v1.complex_matrix_value (), v2.complex_matrix_value ());
}
DEFBINOPX (pow, sparse_complex_matrix, complex_matrix)
--- 49,61 ----
{
CAST_BINOP_ARGS (const octave_sparse_complex_matrix&,
const octave_complex_matrix&);
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (v1.complex_matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOPX (pow, sparse_complex_matrix, complex_matrix)
Index: src/OPERATORS/op-scm-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-scm-m.cc,v
retrieving revision 1.6
diff -c -r1.6 op-scm-m.cc
*** src/OPERATORS/op-scm-m.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-scm-m.cc 26 Apr 2006 21:34:02 -0000
***************
*** 50,57 ****
{
CAST_BINOP_ARGS (const octave_sparse_complex_matrix&,
const octave_matrix&);
! return xdiv (v1.complex_matrix_value (), v2.matrix_value ());
}
DEFBINOPX (pow, sparse_complex_matrix, matrix)
--- 50,62 ----
{
CAST_BINOP_ARGS (const octave_sparse_complex_matrix&,
const octave_matrix&);
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (v1.complex_matrix_value (),
! v2.matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOPX (pow, sparse_complex_matrix, matrix)
Index: src/OPERATORS/op-sm-cm.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-sm-cm.cc,v
retrieving revision 1.6
diff -c -r1.6 op-sm-cm.cc
*** src/OPERATORS/op-sm-cm.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-sm-cm.cc 26 Apr 2006 21:34:02 -0000
***************
*** 49,56 ****
{
CAST_BINOP_ARGS (const octave_sparse_matrix&,
const octave_complex_matrix&);
! return xdiv (v1.matrix_value (), v2.complex_matrix_value ());
}
DEFBINOPX (pow, sparse_matrix, complex_matrix)
--- 49,61 ----
{
CAST_BINOP_ARGS (const octave_sparse_matrix&,
const octave_complex_matrix&);
+ MatrixType typ = v2.matrix_type ();
! ComplexMatrix ret = xdiv (v1.matrix_value (),
! v2.complex_matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOPX (pow, sparse_matrix, complex_matrix)
Index: src/OPERATORS/op-sm-m.cc
===================================================================
RCS file: /usr/local/cvsroot/octave/src/OPERATORS/op-sm-m.cc,v
retrieving revision 1.6
diff -c -r1.6 op-sm-m.cc
*** src/OPERATORS/op-sm-m.cc 14 Apr 2006 04:01:40 -0000 1.6
--- src/OPERATORS/op-sm-m.cc 26 Apr 2006 21:34:02 -0000
***************
*** 48,55 ****
DEFBINOP (div, sparse_matrix, matrix)
{
CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_matrix&);
! return xdiv (v1.matrix_value (), v2.matrix_value ());
}
DEFBINOPX (pow, sparse_matrix, matrix)
--- 48,59 ----
DEFBINOP (div, sparse_matrix, matrix)
{
CAST_BINOP_ARGS (const octave_sparse_matrix&, const octave_matrix&);
+ MatrixType typ = v2.matrix_type ();
! Matrix ret = xdiv (v1.matrix_value (), v2.matrix_value (), typ);
!
! v2.matrix_type (typ);
! return ret;
}
DEFBINOPX (pow, sparse_matrix, matrix)
*** ./libcruft/lapack/dpocon.f.orig 2006-04-26 15:09:53.884432128 +0200
--- ./libcruft/lapack/dpocon.f 2006-04-26 15:08:44.049048720 +0200
***************
*** 0 ****
--- 1,173 ----
+ SUBROUTINE DPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, IWORK,
+ $ INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, N
+ DOUBLE PRECISION ANORM, RCOND
+ * ..
+ * .. Array Arguments ..
+ INTEGER IWORK( * )
+ DOUBLE PRECISION A( LDA, * ), WORK( * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * DPOCON estimates the reciprocal of the condition number (in the
+ * 1-norm) of a real symmetric positive definite matrix using the
+ * Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF.
+ *
+ * An estimate is obtained for norm(inv(A)), and the reciprocal of the
+ * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': Upper triangle of A is stored;
+ * = 'L': Lower triangle of A is stored.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * A (input) DOUBLE PRECISION array, dimension (LDA,N)
+ * The triangular factor U or L from the Cholesky factorization
+ * A = U**T*U or A = L*L**T, as computed by DPOTRF.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * ANORM (input) DOUBLE PRECISION
+ * The 1-norm (or infinity-norm) of the symmetric matrix A.
+ *
+ * RCOND (output) DOUBLE PRECISION
+ * The reciprocal of the condition number of the matrix A,
+ * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
+ * estimate of the 1-norm of inv(A) computed in this routine.
+ *
+ * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
+ *
+ * IWORK (workspace) INTEGER array, dimension (N)
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ONE, ZERO
+ PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL UPPER
+ CHARACTER NORMIN
+ INTEGER IX, KASE
+ DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL LSAME, IDAMAX, DLAMCH
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL DLACON, DLATRS, DRSCL, XERBLA
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC ABS, MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -4
+ ELSE IF( ANORM.LT.ZERO ) THEN
+ INFO = -5
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DPOCON', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ RCOND = ZERO
+ IF( N.EQ.0 ) THEN
+ RCOND = ONE
+ RETURN
+ ELSE IF( ANORM.EQ.ZERO ) THEN
+ RETURN
+ END IF
+ *
+ SMLNUM = DLAMCH( 'Safe minimum' )
+ *
+ * Estimate the 1-norm of inv(A).
+ *
+ KASE = 0
+ NORMIN = 'N'
+ 10 CONTINUE
+ CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE )
+ IF( KASE.NE.0 ) THEN
+ IF( UPPER ) THEN
+ *
+ * Multiply by inv(U').
+ *
+ CALL DLATRS( 'Upper', 'Transpose', 'Non-unit', NORMIN, N, A,
+ $ LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
+ NORMIN = 'Y'
+ *
+ * Multiply by inv(U).
+ *
+ CALL DLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
+ $ A, LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
+ ELSE
+ *
+ * Multiply by inv(L).
+ *
+ CALL DLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
+ $ A, LDA, WORK, SCALEL, WORK( 2*N+1 ), INFO )
+ NORMIN = 'Y'
+ *
+ * Multiply by inv(L').
+ *
+ CALL DLATRS( 'Lower', 'Transpose', 'Non-unit', NORMIN, N, A,
+ $ LDA, WORK, SCALEU, WORK( 2*N+1 ), INFO )
+ END IF
+ *
+ * Multiply by 1/SCALE if doing so will not cause overflow.
+ *
+ SCALE = SCALEL*SCALEU
+ IF( SCALE.NE.ONE ) THEN
+ IX = IDAMAX( N, WORK, 1 )
+ IF( SCALE.LT.ABS( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
+ $ GO TO 20
+ CALL DRSCL( N, SCALE, WORK, 1 )
+ END IF
+ GO TO 10
+ END IF
+ *
+ * Compute the estimate of the reciprocal condition number.
+ *
+ IF( AINVNM.NE.ZERO )
+ $ RCOND = ( ONE / AINVNM ) / ANORM
+ *
+ 20 CONTINUE
+ RETURN
+ *
+ * End of DPOCON
+ *
+ END
*** ./libcruft/lapack/dpotrs.f.orig 2006-04-26 15:09:53.884432128 +0200
--- ./libcruft/lapack/dpotrs.f 2006-04-26 15:08:44.053048112 +0200
***************
*** 0 ****
--- 1,133 ----
+ SUBROUTINE DPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ * ..
+ * .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * DPOTRS solves a system of linear equations A*X = B with a symmetric
+ * positive definite matrix A using the Cholesky factorization
+ * A = U**T*U or A = L*L**T computed by DPOTRF.
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': Upper triangle of A is stored;
+ * = 'L': Lower triangle of A is stored.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * NRHS (input) INTEGER
+ * The number of right hand sides, i.e., the number of columns
+ * of the matrix B. NRHS >= 0.
+ *
+ * A (input) DOUBLE PRECISION array, dimension (LDA,N)
+ * The triangular factor U or L from the Cholesky factorization
+ * A = U**T*U or A = L*L**T, as computed by DPOTRF.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
+ * On entry, the right hand side matrix B.
+ * On exit, the solution matrix X.
+ *
+ * LDB (input) INTEGER
+ * The leading dimension of the array B. LDB >= max(1,N).
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ONE
+ PARAMETER ( ONE = 1.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL UPPER
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL DTRSM, XERBLA
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DPOTRS', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 .OR. NRHS.EQ.0 )
+ $ RETURN
+ *
+ IF( UPPER ) THEN
+ *
+ * Solve A*X = B where A = U'*U.
+ *
+ * Solve U'*X = B, overwriting B with X.
+ *
+ CALL DTRSM( 'Left', 'Upper', 'Transpose', 'Non-unit', N, NRHS,
+ $ ONE, A, LDA, B, LDB )
+ *
+ * Solve U*X = B, overwriting B with X.
+ *
+ CALL DTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
+ $ NRHS, ONE, A, LDA, B, LDB )
+ ELSE
+ *
+ * Solve A*X = B where A = L*L'.
+ *
+ * Solve L*X = B, overwriting B with X.
+ *
+ CALL DTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
+ $ NRHS, ONE, A, LDA, B, LDB )
+ *
+ * Solve L'*X = B, overwriting B with X.
+ *
+ CALL DTRSM( 'Left', 'Lower', 'Transpose', 'Non-unit', N, NRHS,
+ $ ONE, A, LDA, B, LDB )
+ END IF
+ *
+ RETURN
+ *
+ * End of DPOTRS
+ *
+ END
*** ./libcruft/lapack/dtrtrs.f.orig 2006-04-26 15:09:53.913427720 +0200
--- ./libcruft/lapack/dtrtrs.f 2006-04-26 15:08:44.059047200 +0200
***************
*** 0 ****
--- 1,148 ----
+ SUBROUTINE DTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB,
+ $ INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER DIAG, TRANS, UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ * ..
+ * .. Array Arguments ..
+ DOUBLE PRECISION A( LDA, * ), B( LDB, * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * DTRTRS solves a triangular system of the form
+ *
+ * A * X = B or A**T * X = B,
+ *
+ * where A is a triangular matrix of order N, and B is an N-by-NRHS
+ * matrix. A check is made to verify that A is nonsingular.
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': A is upper triangular;
+ * = 'L': A is lower triangular.
+ *
+ * TRANS (input) CHARACTER*1
+ * Specifies the form of the system of equations:
+ * = 'N': A * X = B (No transpose)
+ * = 'T': A**T * X = B (Transpose)
+ * = 'C': A**H * X = B (Conjugate transpose = Transpose)
+ *
+ * DIAG (input) CHARACTER*1
+ * = 'N': A is non-unit triangular;
+ * = 'U': A is unit triangular.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * NRHS (input) INTEGER
+ * The number of right hand sides, i.e., the number of columns
+ * of the matrix B. NRHS >= 0.
+ *
+ * A (input) DOUBLE PRECISION array, dimension (LDA,N)
+ * The triangular matrix A. If UPLO = 'U', the leading N-by-N
+ * upper triangular part of the array A contains the upper
+ * triangular matrix, and the strictly lower triangular part of
+ * A is not referenced. If UPLO = 'L', the leading N-by-N lower
+ * triangular part of the array A contains the lower triangular
+ * matrix, and the strictly upper triangular part of A is not
+ * referenced. If DIAG = 'U', the diagonal elements of A are
+ * also not referenced and are assumed to be 1.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * B (input/output) DOUBLE PRECISION array, dimension (LDB,NRHS)
+ * On entry, the right hand side matrix B.
+ * On exit, if INFO = 0, the solution matrix X.
+ *
+ * LDB (input) INTEGER
+ * The leading dimension of the array B. LDB >= max(1,N).
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ * > 0: if INFO = i, the i-th diagonal element of A is zero,
+ * indicating that the matrix is singular and the solutions
+ * X have not been computed.
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ZERO, ONE
+ PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL NOUNIT
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL DTRSM, XERBLA
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ NOUNIT = LSAME( DIAG, 'N' )
+ IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
+ $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DTRTRS', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 )
+ $ RETURN
+ *
+ * Check for singularity.
+ *
+ IF( NOUNIT ) THEN
+ DO 10 INFO = 1, N
+ IF( A( INFO, INFO ).EQ.ZERO )
+ $ RETURN
+ 10 CONTINUE
+ END IF
+ INFO = 0
+ *
+ * Solve A * x = b or A' * x = b.
+ *
+ CALL DTRSM( 'Left', UPLO, TRANS, DIAG, N, NRHS, ONE, A, LDA, B,
+ $ LDB )
+ *
+ RETURN
+ *
+ * End of DTRTRS
+ *
+ END
*** ./libcruft/lapack/dtrcon.f.orig 2006-04-26 15:09:53.914427568 +0200
--- ./libcruft/lapack/dtrcon.f 2006-04-26 15:08:44.063046592 +0200
***************
*** 0 ****
--- 1,193 ----
+ SUBROUTINE DTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK,
+ $ IWORK, INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER DIAG, NORM, UPLO
+ INTEGER INFO, LDA, N
+ DOUBLE PRECISION RCOND
+ * ..
+ * .. Array Arguments ..
+ INTEGER IWORK( * )
+ DOUBLE PRECISION A( LDA, * ), WORK( * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * DTRCON estimates the reciprocal of the condition number of a
+ * triangular matrix A, in either the 1-norm or the infinity-norm.
+ *
+ * The norm of A is computed and an estimate is obtained for
+ * norm(inv(A)), then the reciprocal of the condition number is
+ * computed as
+ * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
+ *
+ * Arguments
+ * =========
+ *
+ * NORM (input) CHARACTER*1
+ * Specifies whether the 1-norm condition number or the
+ * infinity-norm condition number is required:
+ * = '1' or 'O': 1-norm;
+ * = 'I': Infinity-norm.
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': A is upper triangular;
+ * = 'L': A is lower triangular.
+ *
+ * DIAG (input) CHARACTER*1
+ * = 'N': A is non-unit triangular;
+ * = 'U': A is unit triangular.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * A (input) DOUBLE PRECISION array, dimension (LDA,N)
+ * The triangular matrix A. If UPLO = 'U', the leading N-by-N
+ * upper triangular part of the array A contains the upper
+ * triangular matrix, and the strictly lower triangular part of
+ * A is not referenced. If UPLO = 'L', the leading N-by-N lower
+ * triangular part of the array A contains the lower triangular
+ * matrix, and the strictly upper triangular part of A is not
+ * referenced. If DIAG = 'U', the diagonal elements of A are
+ * also not referenced and are assumed to be 1.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * RCOND (output) DOUBLE PRECISION
+ * The reciprocal of the condition number of the matrix A,
+ * computed as RCOND = 1/(norm(A) * norm(inv(A))).
+ *
+ * WORK (workspace) DOUBLE PRECISION array, dimension (3*N)
+ *
+ * IWORK (workspace) INTEGER array, dimension (N)
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ONE, ZERO
+ PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL NOUNIT, ONENRM, UPPER
+ CHARACTER NORMIN
+ INTEGER IX, KASE, KASE1
+ DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IDAMAX
+ DOUBLE PRECISION DLAMCH, DLANTR
+ EXTERNAL LSAME, IDAMAX, DLAMCH, DLANTR
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL DLACON, DLATRS, DRSCL, XERBLA
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
+ NOUNIT = LSAME( DIAG, 'N' )
+ *
+ IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'DTRCON', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 ) THEN
+ RCOND = ONE
+ RETURN
+ END IF
+ *
+ RCOND = ZERO
+ SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) )
+ *
+ * Compute the norm of the triangular matrix A.
+ *
+ ANORM = DLANTR( NORM, UPLO, DIAG, N, N, A, LDA, WORK )
+ *
+ * Continue only if ANORM > 0.
+ *
+ IF( ANORM.GT.ZERO ) THEN
+ *
+ * Estimate the norm of the inverse of A.
+ *
+ AINVNM = ZERO
+ NORMIN = 'N'
+ IF( ONENRM ) THEN
+ KASE1 = 1
+ ELSE
+ KASE1 = 2
+ END IF
+ KASE = 0
+ 10 CONTINUE
+ CALL DLACON( N, WORK( N+1 ), WORK, IWORK, AINVNM, KASE )
+ IF( KASE.NE.0 ) THEN
+ IF( KASE.EQ.KASE1 ) THEN
+ *
+ * Multiply by inv(A).
+ *
+ CALL DLATRS( UPLO, 'No transpose', DIAG, NORMIN, N, A,
+ $ LDA, WORK, SCALE, WORK( 2*N+1 ), INFO )
+ ELSE
+ *
+ * Multiply by inv(A').
+ *
+ CALL DLATRS( UPLO, 'Transpose', DIAG, NORMIN, N, A, LDA,
+ $ WORK, SCALE, WORK( 2*N+1 ), INFO )
+ END IF
+ NORMIN = 'Y'
+ *
+ * Multiply by 1/SCALE if doing so will not cause overflow.
+ *
+ IF( SCALE.NE.ONE ) THEN
+ IX = IDAMAX( N, WORK, 1 )
+ XNORM = ABS( WORK( IX ) )
+ IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
+ $ GO TO 20
+ CALL DRSCL( N, SCALE, WORK, 1 )
+ END IF
+ GO TO 10
+ END IF
+ *
+ * Compute the estimate of the reciprocal condition number.
+ *
+ IF( AINVNM.NE.ZERO )
+ $ RCOND = ( ONE / ANORM ) / AINVNM
+ END IF
+ *
+ 20 CONTINUE
+ RETURN
+ *
+ * End of DTRCON
+ *
+ END
*** ./libcruft/lapack/zpocon.f.orig 2006-04-26 15:09:53.914427568 +0200
--- ./libcruft/lapack/zpocon.f 2006-04-26 15:08:44.067045984 +0200
***************
*** 0 ****
--- 1,180 ----
+ SUBROUTINE ZPOCON( UPLO, N, A, LDA, ANORM, RCOND, WORK, RWORK,
+ $ INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, N
+ DOUBLE PRECISION ANORM, RCOND
+ * ..
+ * .. Array Arguments ..
+ DOUBLE PRECISION RWORK( * )
+ COMPLEX*16 A( LDA, * ), WORK( * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * ZPOCON estimates the reciprocal of the condition number (in the
+ * 1-norm) of a complex Hermitian positive definite matrix using the
+ * Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF.
+ *
+ * An estimate is obtained for norm(inv(A)), and the reciprocal of the
+ * condition number is computed as RCOND = 1 / (ANORM * norm(inv(A))).
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': Upper triangle of A is stored;
+ * = 'L': Lower triangle of A is stored.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * A (input) COMPLEX*16 array, dimension (LDA,N)
+ * The triangular factor U or L from the Cholesky factorization
+ * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * ANORM (input) DOUBLE PRECISION
+ * The 1-norm (or infinity-norm) of the Hermitian matrix A.
+ *
+ * RCOND (output) DOUBLE PRECISION
+ * The reciprocal of the condition number of the matrix A,
+ * computed as RCOND = 1/(ANORM * AINVNM), where AINVNM is an
+ * estimate of the 1-norm of inv(A) computed in this routine.
+ *
+ * WORK (workspace) COMPLEX*16 array, dimension (2*N)
+ *
+ * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ONE, ZERO
+ PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL UPPER
+ CHARACTER NORMIN
+ INTEGER IX, KASE
+ DOUBLE PRECISION AINVNM, SCALE, SCALEL, SCALEU, SMLNUM
+ COMPLEX*16 ZDUM
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IZAMAX
+ DOUBLE PRECISION DLAMCH
+ EXTERNAL LSAME, IZAMAX, DLAMCH
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL XERBLA, ZDRSCL, ZLACON, ZLATRS
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, DIMAG, MAX
+ * ..
+ * .. Statement Functions ..
+ DOUBLE PRECISION CABS1
+ * ..
+ * .. Statement Function definitions ..
+ CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -4
+ ELSE IF( ANORM.LT.ZERO ) THEN
+ INFO = -5
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZPOCON', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ RCOND = ZERO
+ IF( N.EQ.0 ) THEN
+ RCOND = ONE
+ RETURN
+ ELSE IF( ANORM.EQ.ZERO ) THEN
+ RETURN
+ END IF
+ *
+ SMLNUM = DLAMCH( 'Safe minimum' )
+ *
+ * Estimate the 1-norm of inv(A).
+ *
+ KASE = 0
+ NORMIN = 'N'
+ 10 CONTINUE
+ CALL ZLACON( N, WORK( N+1 ), WORK, AINVNM, KASE )
+ IF( KASE.NE.0 ) THEN
+ IF( UPPER ) THEN
+ *
+ * Multiply by inv(U').
+ *
+ CALL ZLATRS( 'Upper', 'Conjugate transpose', 'Non-unit',
+ $ NORMIN, N, A, LDA, WORK, SCALEL, RWORK, INFO )
+ NORMIN = 'Y'
+ *
+ * Multiply by inv(U).
+ *
+ CALL ZLATRS( 'Upper', 'No transpose', 'Non-unit', NORMIN, N,
+ $ A, LDA, WORK, SCALEU, RWORK, INFO )
+ ELSE
+ *
+ * Multiply by inv(L).
+ *
+ CALL ZLATRS( 'Lower', 'No transpose', 'Non-unit', NORMIN, N,
+ $ A, LDA, WORK, SCALEL, RWORK, INFO )
+ NORMIN = 'Y'
+ *
+ * Multiply by inv(L').
+ *
+ CALL ZLATRS( 'Lower', 'Conjugate transpose', 'Non-unit',
+ $ NORMIN, N, A, LDA, WORK, SCALEU, RWORK, INFO )
+ END IF
+ *
+ * Multiply by 1/SCALE if doing so will not cause overflow.
+ *
+ SCALE = SCALEL*SCALEU
+ IF( SCALE.NE.ONE ) THEN
+ IX = IZAMAX( N, WORK, 1 )
+ IF( SCALE.LT.CABS1( WORK( IX ) )*SMLNUM .OR. SCALE.EQ.ZERO )
+ $ GO TO 20
+ CALL ZDRSCL( N, SCALE, WORK, 1 )
+ END IF
+ GO TO 10
+ END IF
+ *
+ * Compute the estimate of the reciprocal condition number.
+ *
+ IF( AINVNM.NE.ZERO )
+ $ RCOND = ( ONE / AINVNM ) / ANORM
+ *
+ 20 CONTINUE
+ RETURN
+ *
+ * End of ZPOCON
+ *
+ END
*** ./libcruft/lapack/zpotrs.f.orig 2006-04-26 15:09:53.914427568 +0200
--- ./libcruft/lapack/zpotrs.f 2006-04-26 15:08:44.068045832 +0200
***************
*** 0 ****
--- 1,133 ----
+ SUBROUTINE ZPOTRS( UPLO, N, NRHS, A, LDA, B, LDB, INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * September 30, 1994
+ *
+ * .. Scalar Arguments ..
+ CHARACTER UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ * ..
+ * .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * ZPOTRS solves a system of linear equations A*X = B with a Hermitian
+ * positive definite matrix A using the Cholesky factorization
+ * A = U**H*U or A = L*L**H computed by ZPOTRF.
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': Upper triangle of A is stored;
+ * = 'L': Lower triangle of A is stored.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * NRHS (input) INTEGER
+ * The number of right hand sides, i.e., the number of columns
+ * of the matrix B. NRHS >= 0.
+ *
+ * A (input) COMPLEX*16 array, dimension (LDA,N)
+ * The triangular factor U or L from the Cholesky factorization
+ * A = U**H*U or A = L*L**H, as computed by ZPOTRF.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
+ * On entry, the right hand side matrix B.
+ * On exit, the solution matrix X.
+ *
+ * LDB (input) INTEGER
+ * The leading dimension of the array B. LDB >= max(1,N).
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ COMPLEX*16 ONE
+ PARAMETER ( ONE = ( 1.0D+0, 0.0D+0 ) )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL UPPER
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL XERBLA, ZTRSM
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -2
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -3
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -5
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZPOTRS', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 .OR. NRHS.EQ.0 )
+ $ RETURN
+ *
+ IF( UPPER ) THEN
+ *
+ * Solve A*X = B where A = U'*U.
+ *
+ * Solve U'*X = B, overwriting B with X.
+ *
+ CALL ZTRSM( 'Left', 'Upper', 'Conjugate transpose', 'Non-unit',
+ $ N, NRHS, ONE, A, LDA, B, LDB )
+ *
+ * Solve U*X = B, overwriting B with X.
+ *
+ CALL ZTRSM( 'Left', 'Upper', 'No transpose', 'Non-unit', N,
+ $ NRHS, ONE, A, LDA, B, LDB )
+ ELSE
+ *
+ * Solve A*X = B where A = L*L'.
+ *
+ * Solve L*X = B, overwriting B with X.
+ *
+ CALL ZTRSM( 'Left', 'Lower', 'No transpose', 'Non-unit', N,
+ $ NRHS, ONE, A, LDA, B, LDB )
+ *
+ * Solve L'*X = B, overwriting B with X.
+ *
+ CALL ZTRSM( 'Left', 'Lower', 'Conjugate transpose', 'Non-unit',
+ $ N, NRHS, ONE, A, LDA, B, LDB )
+ END IF
+ *
+ RETURN
+ *
+ * End of ZPOTRS
+ *
+ END
*** ./libcruft/lapack/ztrtrs.f.orig 2006-04-26 15:09:53.914427568 +0200
--- ./libcruft/lapack/ztrtrs.f 2006-04-26 15:08:44.075044768 +0200
***************
*** 0 ****
--- 1,149 ----
+ SUBROUTINE ZTRTRS( UPLO, TRANS, DIAG, N, NRHS, A, LDA, B, LDB,
+ $ INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * September 30, 1994
+ *
+ * .. Scalar Arguments ..
+ CHARACTER DIAG, TRANS, UPLO
+ INTEGER INFO, LDA, LDB, N, NRHS
+ * ..
+ * .. Array Arguments ..
+ COMPLEX*16 A( LDA, * ), B( LDB, * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * ZTRTRS solves a triangular system of the form
+ *
+ * A * X = B, A**T * X = B, or A**H * X = B,
+ *
+ * where A is a triangular matrix of order N, and B is an N-by-NRHS
+ * matrix. A check is made to verify that A is nonsingular.
+ *
+ * Arguments
+ * =========
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': A is upper triangular;
+ * = 'L': A is lower triangular.
+ *
+ * TRANS (input) CHARACTER*1
+ * Specifies the form of the system of equations:
+ * = 'N': A * X = B (No transpose)
+ * = 'T': A**T * X = B (Transpose)
+ * = 'C': A**H * X = B (Conjugate transpose)
+ *
+ * DIAG (input) CHARACTER*1
+ * = 'N': A is non-unit triangular;
+ * = 'U': A is unit triangular.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * NRHS (input) INTEGER
+ * The number of right hand sides, i.e., the number of columns
+ * of the matrix B. NRHS >= 0.
+ *
+ * A (input) COMPLEX*16 array, dimension (LDA,N)
+ * The triangular matrix A. If UPLO = 'U', the leading N-by-N
+ * upper triangular part of the array A contains the upper
+ * triangular matrix, and the strictly lower triangular part of
+ * A is not referenced. If UPLO = 'L', the leading N-by-N lower
+ * triangular part of the array A contains the lower triangular
+ * matrix, and the strictly upper triangular part of A is not
+ * referenced. If DIAG = 'U', the diagonal elements of A are
+ * also not referenced and are assumed to be 1.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
+ * On entry, the right hand side matrix B.
+ * On exit, if INFO = 0, the solution matrix X.
+ *
+ * LDB (input) INTEGER
+ * The leading dimension of the array B. LDB >= max(1,N).
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ * > 0: if INFO = i, the i-th diagonal element of A is zero,
+ * indicating that the matrix is singular and the solutions
+ * X have not been computed.
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ COMPLEX*16 ZERO, ONE
+ PARAMETER ( ZERO = ( 0.0D+0, 0.0D+0 ),
+ $ ONE = ( 1.0D+0, 0.0D+0 ) )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL NOUNIT
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ EXTERNAL LSAME
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL XERBLA, ZTRSM
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC MAX
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ NOUNIT = LSAME( DIAG, 'N' )
+ IF( .NOT.LSAME( UPLO, 'U' ) .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.LSAME( TRANS, 'N' ) .AND. .NOT.
+ $ LSAME( TRANS, 'T' ) .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( NRHS.LT.0 ) THEN
+ INFO = -5
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -7
+ ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
+ INFO = -9
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZTRTRS', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 )
+ $ RETURN
+ *
+ * Check for singularity.
+ *
+ IF( NOUNIT ) THEN
+ DO 10 INFO = 1, N
+ IF( A( INFO, INFO ).EQ.ZERO )
+ $ RETURN
+ 10 CONTINUE
+ END IF
+ INFO = 0
+ *
+ * Solve A * x = b, A**T * x = b, or A**H * x = b.
+ *
+ CALL ZTRSM( 'Left', UPLO, TRANS, DIAG, N, NRHS, ONE, A, LDA, B,
+ $ LDB )
+ *
+ RETURN
+ *
+ * End of ZTRTRS
+ *
+ END
*** ./libcruft/lapack/ztrcon.f.orig 2006-04-26 15:09:53.914427568 +0200
--- ./libcruft/lapack/ztrcon.f 2006-04-26 15:08:44.083043552 +0200
***************
*** 0 ****
--- 1,200 ----
+ SUBROUTINE ZTRCON( NORM, UPLO, DIAG, N, A, LDA, RCOND, WORK,
+ $ RWORK, INFO )
+ *
+ * -- LAPACK routine (version 3.0) --
+ * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
+ * Courant Institute, Argonne National Lab, and Rice University
+ * March 31, 1993
+ *
+ * .. Scalar Arguments ..
+ CHARACTER DIAG, NORM, UPLO
+ INTEGER INFO, LDA, N
+ DOUBLE PRECISION RCOND
+ * ..
+ * .. Array Arguments ..
+ DOUBLE PRECISION RWORK( * )
+ COMPLEX*16 A( LDA, * ), WORK( * )
+ * ..
+ *
+ * Purpose
+ * =======
+ *
+ * ZTRCON estimates the reciprocal of the condition number of a
+ * triangular matrix A, in either the 1-norm or the infinity-norm.
+ *
+ * The norm of A is computed and an estimate is obtained for
+ * norm(inv(A)), then the reciprocal of the condition number is
+ * computed as
+ * RCOND = 1 / ( norm(A) * norm(inv(A)) ).
+ *
+ * Arguments
+ * =========
+ *
+ * NORM (input) CHARACTER*1
+ * Specifies whether the 1-norm condition number or the
+ * infinity-norm condition number is required:
+ * = '1' or 'O': 1-norm;
+ * = 'I': Infinity-norm.
+ *
+ * UPLO (input) CHARACTER*1
+ * = 'U': A is upper triangular;
+ * = 'L': A is lower triangular.
+ *
+ * DIAG (input) CHARACTER*1
+ * = 'N': A is non-unit triangular;
+ * = 'U': A is unit triangular.
+ *
+ * N (input) INTEGER
+ * The order of the matrix A. N >= 0.
+ *
+ * A (input) COMPLEX*16 array, dimension (LDA,N)
+ * The triangular matrix A. If UPLO = 'U', the leading N-by-N
+ * upper triangular part of the array A contains the upper
+ * triangular matrix, and the strictly lower triangular part of
+ * A is not referenced. If UPLO = 'L', the leading N-by-N lower
+ * triangular part of the array A contains the lower triangular
+ * matrix, and the strictly upper triangular part of A is not
+ * referenced. If DIAG = 'U', the diagonal elements of A are
+ * also not referenced and are assumed to be 1.
+ *
+ * LDA (input) INTEGER
+ * The leading dimension of the array A. LDA >= max(1,N).
+ *
+ * RCOND (output) DOUBLE PRECISION
+ * The reciprocal of the condition number of the matrix A,
+ * computed as RCOND = 1/(norm(A) * norm(inv(A))).
+ *
+ * WORK (workspace) COMPLEX*16 array, dimension (2*N)
+ *
+ * RWORK (workspace) DOUBLE PRECISION array, dimension (N)
+ *
+ * INFO (output) INTEGER
+ * = 0: successful exit
+ * < 0: if INFO = -i, the i-th argument had an illegal value
+ *
+ * =====================================================================
+ *
+ * .. Parameters ..
+ DOUBLE PRECISION ONE, ZERO
+ PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
+ * ..
+ * .. Local Scalars ..
+ LOGICAL NOUNIT, ONENRM, UPPER
+ CHARACTER NORMIN
+ INTEGER IX, KASE, KASE1
+ DOUBLE PRECISION AINVNM, ANORM, SCALE, SMLNUM, XNORM
+ COMPLEX*16 ZDUM
+ * ..
+ * .. External Functions ..
+ LOGICAL LSAME
+ INTEGER IZAMAX
+ DOUBLE PRECISION DLAMCH, ZLANTR
+ EXTERNAL LSAME, IZAMAX, DLAMCH, ZLANTR
+ * ..
+ * .. External Subroutines ..
+ EXTERNAL XERBLA, ZDRSCL, ZLACON, ZLATRS
+ * ..
+ * .. Intrinsic Functions ..
+ INTRINSIC ABS, DBLE, DIMAG, MAX
+ * ..
+ * .. Statement Functions ..
+ DOUBLE PRECISION CABS1
+ * ..
+ * .. Statement Function definitions ..
+ CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) )
+ * ..
+ * .. Executable Statements ..
+ *
+ * Test the input parameters.
+ *
+ INFO = 0
+ UPPER = LSAME( UPLO, 'U' )
+ ONENRM = NORM.EQ.'1' .OR. LSAME( NORM, 'O' )
+ NOUNIT = LSAME( DIAG, 'N' )
+ *
+ IF( .NOT.ONENRM .AND. .NOT.LSAME( NORM, 'I' ) ) THEN
+ INFO = -1
+ ELSE IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
+ INFO = -2
+ ELSE IF( .NOT.NOUNIT .AND. .NOT.LSAME( DIAG, 'U' ) ) THEN
+ INFO = -3
+ ELSE IF( N.LT.0 ) THEN
+ INFO = -4
+ ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
+ INFO = -6
+ END IF
+ IF( INFO.NE.0 ) THEN
+ CALL XERBLA( 'ZTRCON', -INFO )
+ RETURN
+ END IF
+ *
+ * Quick return if possible
+ *
+ IF( N.EQ.0 ) THEN
+ RCOND = ONE
+ RETURN
+ END IF
+ *
+ RCOND = ZERO
+ SMLNUM = DLAMCH( 'Safe minimum' )*DBLE( MAX( 1, N ) )
+ *
+ * Compute the norm of the triangular matrix A.
+ *
+ ANORM = ZLANTR( NORM, UPLO, DIAG, N, N, A, LDA, RWORK )
+ *
+ * Continue only if ANORM > 0.
+ *
+ IF( ANORM.GT.ZERO ) THEN
+ *
+ * Estimate the norm of the inverse of A.
+ *
+ AINVNM = ZERO
+ NORMIN = 'N'
+ IF( ONENRM ) THEN
+ KASE1 = 1
+ ELSE
+ KASE1 = 2
+ END IF
+ KASE = 0
+ 10 CONTINUE
+ CALL ZLACON( N, WORK( N+1 ), WORK, AINVNM, KASE )
+ IF( KASE.NE.0 ) THEN
+ IF( KASE.EQ.KASE1 ) THEN
+ *
+ * Multiply by inv(A).
+ *
+ CALL ZLATRS( UPLO, 'No transpose', DIAG, NORMIN, N, A,
+ $ LDA, WORK, SCALE, RWORK, INFO )
+ ELSE
+ *
+ * Multiply by inv(A').
+ *
+ CALL ZLATRS( UPLO, 'Conjugate transpose', DIAG, NORMIN,
+ $ N, A, LDA, WORK, SCALE, RWORK, INFO )
+ END IF
+ NORMIN = 'Y'
+ *
+ * Multiply by 1/SCALE if doing so will not cause overflow.
+ *
+ IF( SCALE.NE.ONE ) THEN
+ IX = IZAMAX( N, WORK, 1 )
+ XNORM = CABS1( WORK( IX ) )
+ IF( SCALE.LT.XNORM*SMLNUM .OR. SCALE.EQ.ZERO )
+ $ GO TO 20
+ CALL ZDRSCL( N, SCALE, WORK, 1 )
+ END IF
+ GO TO 10
+ END IF
+ *
+ * Compute the estimate of the reciprocal condition number.
+ *
+ IF( AINVNM.NE.ZERO )
+ $ RCOND = ( ONE / ANORM ) / AINVNM
+ END IF
+ *
+ 20 CONTINUE
+ RETURN
+ *
+ * End of ZTRCON
+ *
+ END
*** ./liboctave/MatrixType.h.orig 2006-04-26 23:45:31.783156624 +0200
--- ./liboctave/MatrixType.h 2006-04-24 23:35:35.670059496 +0200
***************
*** 0 ****
--- 1,120 ----
+ /*
+
+ Copyright (C) 2006 David Bateman
+ Copyright (C) 2006 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.
+
+ */
+
+ #if !defined (octave_MatrixType_h)
+ #define octave_MatrixType_h
+
+ class Matrix;
+ class ComplexMatrix;
+
+ class
+ MatrixType
+ {
+ public:
+ enum matrix_type {
+ Unknown = 0,
+ Full,
+ Upper,
+ Lower,
+ Permuted_Upper,
+ Permuted_Lower,
+ Hermitian,
+ Rectangular
+ };
+
+ MatrixType (void) : typ (MatrixType::Unknown), nperm (0) { }
+
+ MatrixType (const MatrixType &a);
+
+ MatrixType (const Matrix &a);
+
+ MatrixType (const ComplexMatrix &a);
+
+ MatrixType (const matrix_type t);
+
+ MatrixType (const matrix_type t, const octave_idx_type np,
+ const octave_idx_type *p);
+
+ ~MatrixType (void);
+
+ MatrixType& operator = (const MatrixType& a);
+
+ int type (bool quiet = true);
+
+ int type (const Matrix &a);
+
+ int type (const ComplexMatrix &a);
+
+ bool is_upper_triangular (void) const
+ { return (typ == Upper || typ == Permuted_Upper); }
+
+ bool is_lower_triangular (void) const
+ { return (typ == Lower || typ == Permuted_Lower); }
+
+ bool is_hermitian (void) const
+ { return (typ == Hermitian); }
+
+ bool is_rectangular (void) const { return (typ == Rectangular); }
+
+ bool is_known (void) const { return (typ != Unknown); }
+
+ bool is_unknown (void) const { return (typ == Unknown); }
+
+ void info (void) const;
+
+ octave_idx_type * triangular_perm (void) const { return perm; }
+
+ void invalidate_type (void) { typ = Unknown; }
+
+ void mark_as_upper_triangular (void) { typ = Upper; }
+
+ void mark_as_lower_triangular (void) { typ = Lower; }
+
+ void mark_as_full (void) { typ = Full; }
+
+ void mark_as_rectangular (void) { typ = Rectangular; }
+
+ void mark_as_symmetric (void);
+
+ void mark_as_unsymmetric (void);
+
+ void mark_as_permuted (const octave_idx_type np, const octave_idx_type *p);
+
+ void mark_as_unpermuted (void);
+
+ MatrixType transpose (void) const;
+
+ private:
+ void type (int new_typ) { typ = static_cast<matrix_type>(new_typ); }
+
+ matrix_type typ;
+ octave_idx_type nperm;
+ octave_idx_type *perm;
+ };
+
+ #endif
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C++ ***
+ ;;; End: ***
+ */
*** ./liboctave/MatrixType.cc.orig 2006-04-26 23:45:34.319771000 +0200
--- ./liboctave/MatrixType.cc 2006-04-26 23:22:41.086534184 +0200
***************
*** 0 ****
--- 1,396 ----
+ /*
+
+ Copyright (C) 2006 David Bateman
+ Copyright (C) 2006 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
+
+ #include <vector>
+
+ #include "MatrixType.h"
+ #include "dMatrix.h"
+ #include "CMatrix.h"
+
+ // FIXME There is a large code duplication here
+
+ MatrixType::MatrixType (const MatrixType &a) : typ (a.typ),
+ nperm (a.nperm)
+ {
+ if (nperm != 0)
+ {
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = a.perm[i];
+ }
+ }
+
+ MatrixType::MatrixType (const Matrix &a)
+ {
+ octave_idx_type nrows = a.rows ();
+ octave_idx_type ncols = a.cols ();
+ nperm = 0;
+
+ if (ncols == nrows)
+ {
+ bool upper = true;
+ bool lower = true;
+ bool hermitian = true;
+
+ for (octave_idx_type j = 0; j < ncols; j++)
+ {
+ if (j < nrows)
+ {
+ if (a.elem (j,j) == 0.)
+ {
+ upper = false;
+ lower = false;
+ hermitian = false;
+ break;
+ }
+ if (a.elem (j,j) < 0.)
+ hermitian = false;
+ }
+ for (octave_idx_type i = 0; i < j; i++)
+ if (lower && a.elem (i,j) != 0.)
+ {
+ lower = false;
+ break;
+ }
+ for (octave_idx_type i = j+1; i < nrows; i++)
+ {
+ if (hermitian && a.elem (i, j) != a.elem (j, i))
+ hermitian = false;
+ if (upper && a.elem (i,j) != 0)
+ upper = false;
+ }
+ if (!upper && !lower && !hermitian)
+ break;
+ }
+
+ if (upper)
+ typ = MatrixType::Upper;
+ else if (lower)
+ typ = MatrixType::Lower;
+ else if (hermitian)
+ typ = MatrixType::Hermitian;
+ else if (ncols == nrows)
+ typ = MatrixType::Full;
+ }
+ else
+ typ = MatrixType::Rectangular;
+ }
+
+ MatrixType::MatrixType (const ComplexMatrix &a)
+ {
+ octave_idx_type nrows = a.rows ();
+ octave_idx_type ncols = a.cols ();
+ nperm = 0;
+
+ if (ncols == nrows)
+ {
+ bool upper = true;
+ bool lower = true;
+ bool hermitian = true;
+
+ for (octave_idx_type j = 0; j < ncols; j++)
+ {
+ if (j < ncols)
+ {
+ if (imag(a.elem (j,j)) == 0. &&
+ real(a.elem (j,j)) == 0.)
+ {
+ upper = false;
+ lower = false;
+ hermitian = false;
+ break;
+ }
+
+ if (imag(a.elem (j,j)) != 0. ||
+ real(a.elem (j,j)) < 0.)
+ hermitian = false;
+ }
+ for (octave_idx_type i = 0; i < j; i++)
+ if (lower && (real(a.elem (i,j)) != 0 || imag(a.elem (i,j)) != 0))
+ {
+ lower = false;
+ break;
+ }
+ for (octave_idx_type i = j+1; i < nrows; i++)
+ {
+ if (hermitian && a.elem (i, j) != conj(a.elem (j, i)))
+ hermitian = false;
+ if (upper && (real(a.elem (i,j)) != 0 ||
+ imag(a.elem (i,j)) != 0))
+ upper = false;
+ }
+ if (!upper && !lower && !hermitian)
+ break;
+ }
+
+ if (upper)
+ typ = MatrixType::Upper;
+ else if (lower)
+ typ = MatrixType::Lower;
+ else if (hermitian)
+ typ = MatrixType::Hermitian;
+ else if (ncols == nrows)
+ typ = MatrixType::Full;
+ }
+ else
+ typ = MatrixType::Rectangular;
+ }
+
+ MatrixType::MatrixType (const matrix_type t) : typ (MatrixType::Unknown),
+ nperm (0)
+ {
+ if (t == MatrixType::Full || t == MatrixType::Upper ||
+ t == MatrixType::Lower || t == MatrixType::Hermitian ||
+ t == MatrixType::Rectangular)
+ typ = t;
+ else
+ (*current_liboctave_warning_handler) ("Invalid matrix type");
+ }
+
+ MatrixType::MatrixType (const matrix_type t, const octave_idx_type np,
+ const octave_idx_type *p) : typ (MatrixType::Unknown),
+ nperm (0)
+ {
+ #if 0
+ if x(t == MatrixType::Permuted_Upper || t == MatrixType::Permuted_Lower)
+ {
+ typ = t;
+ nperm = np;
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = p[i];
+ }
+ else
+ (*current_liboctave_warning_handler) ("Invalid matrix type");
+ #else
+ (*current_liboctave_warning_handler)
+ ("Permuted triangular matrix not implemented");
+ #endif
+ }
+
+ MatrixType::~MatrixType (void)
+ {
+ if (nperm != 0)
+ {
+ delete [] perm;
+ }
+ }
+
+ MatrixType&
+ MatrixType::operator = (const MatrixType& a)
+ {
+ if (this != &a)
+ {
+ typ = a.typ;
+ nperm = a.nperm;
+
+ if (nperm != 0)
+ {
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = a.perm[i];
+ }
+ }
+
+ return *this;
+ }
+
+ int
+ MatrixType::type (bool quiet)
+ {
+ if (typ != MatrixType::Unknown)
+ {
+ if (!quiet)
+ (*current_liboctave_warning_handler)
+ ("Using Cached Matrix Matrix Type");
+
+ return typ;
+ }
+
+ typ = MatrixType::Unknown;
+
+ return typ;
+ }
+
+ int
+ MatrixType::type (const Matrix &a)
+ {
+ if (typ != MatrixType::Unknown)
+ {
+ // FIXME - Should we include something like the spumoni flag?
+ // (*current_liboctave_warning_handler)
+ // ("Using Cached Matrix Matrix Type");
+
+ return typ;
+ }
+
+ MatrixType tmp_typ (a);
+ typ = tmp_typ.typ;
+ nperm = tmp_typ.nperm;
+
+ if (nperm != 0)
+ {
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = tmp_typ.perm[i];
+ }
+
+ return typ;
+ }
+
+ int
+ MatrixType::type (const ComplexMatrix &a)
+ {
+ if (typ != MatrixType::Unknown)
+ {
+ // FIXME - Should we include something like the spumoni flag?
+ // (*current_liboctave_warning_handler)
+ // ("Using Cached Matrix Matrix Type");
+
+ return typ;
+ }
+
+ MatrixType tmp_typ (a);
+ typ = tmp_typ.typ;
+ nperm = tmp_typ.nperm;
+
+ if (nperm != 0)
+ {
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = tmp_typ.perm[i];
+ }
+
+ return typ;
+ }
+
+ void
+ MatrixType::info () const
+ {
+ if (typ == MatrixType::Unknown)
+ (*current_liboctave_warning_handler)
+ ("Unknown Matrix Type");
+ else if (typ == MatrixType::Upper)
+ (*current_liboctave_warning_handler)
+ ("Upper Triangular Matrix");
+ else if (typ == MatrixType::Lower)
+ (*current_liboctave_warning_handler)
+ ("Lower Triangular Matrix");
+ else if (typ == MatrixType::Permuted_Upper)
+ (*current_liboctave_warning_handler)
+ ("Permuted Upper Triangular Matrix");
+ else if (typ == MatrixType::Permuted_Lower)
+ (*current_liboctave_warning_handler)
+ ("Permuted Lower Triangular Matrix");
+ else if (typ == MatrixType::Hermitian)
+ (*current_liboctave_warning_handler)
+ ("Hermitian/Symmetric Matrix");
+ else if (typ == MatrixType::Rectangular)
+ (*current_liboctave_warning_handler)
+ ("Rectangular/Singular Matrix");
+ else if (typ == MatrixType::Full)
+ (*current_liboctave_warning_handler)
+ ("Full Matrix");
+ }
+
+ void
+ MatrixType::mark_as_symmetric (void)
+ {
+ if (typ == MatrixType::Full || typ == MatrixType::Hermitian ||
+ typ == MatrixType::Unknown)
+ typ = MatrixType::Hermitian;
+ else
+ (*current_liboctave_error_handler)
+ ("Can not mark current matrix type as symmetric");
+ }
+
+ void
+ MatrixType::mark_as_unsymmetric (void)
+ {
+ if (typ == MatrixType::Full || typ == MatrixType::Hermitian ||
+ typ == MatrixType::Unknown)
+ typ = MatrixType::Full;
+ }
+
+ void
+ MatrixType::mark_as_permuted (const octave_idx_type np, const octave_idx_type
*p)
+ {
+ #if 0
+ nperm = np;
+ perm = new octave_idx_type [nperm];
+ for (octave_idx_type i = 0; i < nperm; i++)
+ perm[i] = p[i];
+
+ if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+ typ = MatrixType::Permuted_Upper;
+ else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+ typ = MatrixType::Permuted_Lower;
+ else
+ (*current_liboctave_error_handler)
+ ("Can not mark current matrix type as symmetric");
+ #else
+ (*current_liboctave_warning_handler)
+ ("Permuted triangular matrix not implemented");
+ #endif
+ }
+
+ void
+ MatrixType::mark_as_unpermuted (void)
+ {
+ if (nperm)
+ {
+ nperm = 0;
+ delete [] perm;
+ }
+
+ if (typ == MatrixType::Upper || typ == MatrixType::Permuted_Upper)
+ typ = MatrixType::Upper;
+ else if (typ == MatrixType::Lower || typ == MatrixType::Permuted_Lower)
+ typ = MatrixType::Lower;
+ }
+
+ MatrixType
+ MatrixType::transpose (void) const
+ {
+ MatrixType retval (*this);
+ if (typ == MatrixType::Upper)
+ retval.typ = MatrixType::Lower;
+ else if (typ == MatrixType::Permuted_Upper)
+ retval.typ = MatrixType::Permuted_Lower;
+ else if (typ == MatrixType::Lower)
+ retval.typ = MatrixType::Upper;
+ else if (typ == MatrixType::Permuted_Lower)
+ retval.typ = MatrixType::Permuted_Upper;
+
+ return retval;
+ }
+
+ /*
+ ;;; Local Variables: ***
+ ;;; mode: C++ ***
+ ;;; End: ***
+ */
+
2006-04-26 David Bateman <address@hidden>
* CMatrix.cc (zpotrf, zpocon, zpotrs, ztrcon, ztrtrs):
External declaration of lapack triangular and Cholesky codes.
(ComplexMatrix::utsolve, ComplexMatrix::ltsolve,
ComplexMatrix::fsolve): New private solver codes for
upper, lower and LU/Cholesky solvers.
(ComplexMatrix::solve): New versions for cached matrix
type. Adapt old versions to call new versions
* CMatrix.h (utsolve, ltsolve, fsolve): Declaration of
new solvers.
(solve): New versions for cached matrix type.
* dMatrix.cc (dpotrf, dpocon, dpotrs, dtrcon, dtrtrs):
External declaration of lapack triangular and Cholesky codes.
(Matrix::utsolve, Matrix::ltsolve,
Matrix::fsolve): New private solver codes for
upper, lower and LU/Cholesky solvers.
(Matrix::solve): New versions for cached matrix
type. Adapt old versions to call new versions
* dMatrix.h (utsolve, ltsolve, fsolve): Declaration of
new solvers.
(solve): New versions for cached matrix type.
* MatrixType.cc: New file for class to cache matrix type.
* MAtrixType.h: ditto.
* Makefile.in (MATRIX_INC, MATRIX_SRC): Add MatrixType.h and
MatrixType.cc respectively.
2006-04-26 David Bateman <address@hidden>
* ov-base-mat.h: Add caching of matrix type, and code to supply
and copy matrix type.
* ov-bool-mat.h: Add caching to constructor.
* ov-re-mat.h: ditto.
* ov-cx-mat.h: ditto.
* ov.cc: Add to the BoolMatrix, Matrix and the ComplexMatrix
octave_value constructors, the ability to specify the matrix type.
* ov.h: Adapt declaration of above constructors.
* xdiv.cc (xdiv, xleftdiv): Pass the matrix type, simplfy since
the matrix solve function now calls lssolve if singular.
* xdiv.h (xdvi, xleftdiv): Update the declarations
* OPERATORS/op-cm-cm.cc, OPERATORS/op-cm-cs.cc,
OPERATORS/op-cm-m.cc, OPERATORS/op-cm-s.cc,
OPERATORS/op-cm-scm.cc, OPERATORS/op-cm-sm.cc,
OPERATORS/op-cs-cm.cc, OPERATORS/op-cs-m.cc,
OPERATORS/op-m-cm.cc, OPERATORS/op-m-cs.cc,
OPERATORS/op-m-m.cc, OPERATORS/op-m-s.cc,
OPERATORS/op-m-scm.cc, OPERATORS/op-m-sm.cc,
OPERATORS/op-s-cm.cc, OPERATORS/op-s-m.cc,
OPERATORS/op-scm-cm.cc, OPERATORS/op-scm-m.cc,
OPERATORS/op-sm-cm.cc, OPERATORS/op-sm-m.cc: Update use of
xdiv and xleftdiv functions to allow matrix type caching.
* DLD-FUNCTIONS/matrix_type.cc (Fmatrix_type): Update to allow typing
of Matrix, and ComplexMatrix types. Add new test code for this.
2006-04-26 David Bateman <address@hidden>
* lapack/dpocon.f, lapack/zpocon.f, lapack/dpotrs.f,
lapack/zpotrs.f, lapack/dtrcon.f, lapack/ztrcon.f,
lapack/dtrtrs.f, lapack/ztrtrs.f: New files.
- matrix type caching, triangualar and Cholesky solvers,
David Bateman <=