[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: solvetoep
|
From: |
David Bateman |
|
Subject: |
Re: solvetoep |
|
Date: |
Tue, 28 Nov 2006 23:00:23 +0100 |
|
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
Ok I took a further look at this patch and made a complete working
version based on your code. I addressed the three outstanding issues I
mentioned (reversal of output, singularity handler and toeplitz having
higher priority than positive definite), added tests and fixed a couple
of other issues with singularity handling in ltsolve and utsolve, and a
spurious warning message from matrix_type.
John, basically this is vastly superior solver for toeplitz matrices
than is currently used, and so I'd suggest including it as is. Gorazd
thinks he can possibly improve it further though..
D.
*** ./liboctave/CMatrix.cc.orig4 2006-11-15 20:52:45.000000000 +0100
--- ./liboctave/CMatrix.cc 2006-11-28 22:31:59.392060236 +0100
***************
*** 27,34 ****
#endif
#include <cfloat>
-
#include <iostream>
// FIXME
#ifdef HAVE_SYS_TYPES_H
--- 27,34 ----
#endif
#include <cfloat>
#include <iostream>
+ #include <vector>
// FIXME
#ifdef HAVE_SYS_TYPES_H
***************
*** 1539,1544 ****
--- 1539,1670 ----
}
ComplexMatrix
+ ComplexMatrix::toepsolve (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 ();
+ volatile int typ = mattype.type ();
+
+ info = 0;
+
+ if (nr == 0 || nc == 0 || nr != b.rows ())
+ (*current_liboctave_error_handler)
+ ("matrix dimension mismatch solution of linear equations");
+
+ else if (typ != MatrixType::Toeplitz)
+ (*current_liboctave_error_handler) ("incorrect matrix type");
+ else
+ {
+ octave_idx_type b_nc = b.cols();
+ retval.resize(nc, b_nc);
+
+ OCTAVE_LOCAL_BUFFER (Complex, c1, nr - 1);
+ OCTAVE_LOCAL_BUFFER (Complex, c2, nc - 1);
+
+ for (octave_idx_type j = 0; j < b_nc; j++)
+ {
+ Complex r1, r2, r3, r5, r6;
+
+ r1 = elem(0);
+
+ if (r1 == 0.)
+ {
+ info = -2;
+ goto singular;
+ }
+
+ retval.elem (0, j) = b.elem(0,j)/r1;
+
+ if ( nr == 1 )
+ continue;
+
+ for (octave_idx_type nsub = 2; nsub <= nr; nsub++)
+ {
+ r5 = elem(nsub-1);
+ r6 = elem(nr*(nsub-1));
+
+ if (nsub > 2)
+ {
+ c1[nsub-2] = r2;
+ for (octave_idx_type i = 1; i <= nsub - 2; i++)
+ {
+ r5 = r5 + elem(i) * c1[nsub-i-1];
+ r6 = r6 + elem(nr*i) * c2[i-1];
+ }
+ }
+
+ r2 = -r5/r1;
+ r3 = -r6/r1;
+ r1 = r1 + r5 * r3;
+
+ if (r1 == 0.)
+ {
+ info = -2;
+ goto singular;
+ }
+
+ if (nsub > 2)
+ {
+ r6 = c2[0];
+ c2[nsub-2] = 0;
+
+ for (octave_idx_type i = 2; i <= nsub - 1; i++)
+ {
+ r5 = c2[i-1];
+ c2[i-1] = c1[i-1] * r3 + r6;
+ c1[i-1] += r6 * r2;
+ r6 = r5;
+ }
+ }
+
+ c2[0] = r3;
+ r5 = 0;
+ for (octave_idx_type i = 1; i <= nsub - 1; i++)
+ r5 = r5 + elem(i) * retval.elem(nsub-i-1, j);
+
+ r6 = ( b(nsub-1) - r5 )/r1;
+
+ for (octave_idx_type i = 1; i <= nsub-1; i++)
+ retval.elem(i-1, j) += c2[i-1] * r6;
+
+ retval.elem(nsub-1, j) = r6;
+ }
+ }
+ }
+
+ singular:
+ if (calc_cond)
+ (*current_liboctave_error_handler)
+ ("calculation of condition number estimate not yet implemented");
+
+ if (info)
+ {
+ if (sing_handler)
+ sing_handler (rcond);
+ else
+ (*current_liboctave_error_handler)
+ ("matrix singular to machine precision");
+ mattype.mark_as_rectangular ();
+ }
+
+ return retval;
+ }
+
+ /*
+
+ %!test
+ %! a = 1i*toeplitz(1:5,2:6);
+ %! A = matrix_type(a,'full');
+ %! assert(a\ones(5,1),A\ones(5,1),1e-10)
+
+ */
+
+ ComplexMatrix
ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b,
octave_idx_type& info, double& rcond,
solve_singularity_handler sing_handler,
***************
*** 1636,1641 ****
--- 1762,1769 ----
(*current_liboctave_error_handler)
("unrecoverable error in dtrtrs");
}
+ else
+ mattype.mark_as_rectangular ();
}
}
else
***************
*** 1743,1748 ****
--- 1871,1878 ----
(*current_liboctave_error_handler)
("unrecoverable error in dtrtrs");
}
+ else
+ mattype.mark_as_rectangular ();
}
}
else
***************
*** 2033,2038 ****
--- 2163,2170 ----
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::Toeplitz)
+ retval = toepsolve (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)
*** ./liboctave/CMatrix.h.orig4 2006-10-27 03:45:54.000000000 +0200
--- ./liboctave/CMatrix.h 2006-11-28 16:26:45.827629583 +0100
***************
*** 153,158 ****
--- 153,164 ----
ComplexDET determinant (octave_idx_type& info, double& rcond, int calc_cond
= 1) const;
private:
+ // Toeplitz matrix solver
+ ComplexMatrix toepsolve (MatrixType &typ, const ComplexMatrix& b,
+ octave_idx_type& info, double& rcond,
+ solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
// Upper triangular matrix solvers
ComplexMatrix utsolve (MatrixType &typ, const ComplexMatrix& b,
octave_idx_type& info, double& rcond,
*** ./liboctave/dMatrix.h.orig4 2006-10-27 03:45:55.000000000 +0200
--- ./liboctave/dMatrix.h 2006-11-28 16:26:56.741099084 +0100
***************
*** 125,130 ****
--- 125,135 ----
DET determinant (octave_idx_type& info, double& rcond, int calc_cond = 1)
const;
private:
+ // Toeplitz matrix solver
+ Matrix toepsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
+ double& rcond, solve_singularity_handler sing_handler,
+ bool calc_cond = false) const;
+
// Upper triangular matrix solvers
Matrix utsolve (MatrixType &typ, const Matrix& b, octave_idx_type& info,
double& rcond, solve_singularity_handler sing_handler,
*** ./liboctave/dMatrix.cc.orig4 2006-11-15 20:52:45.000000000 +0100
--- ./liboctave/dMatrix.cc 2006-11-28 22:32:09.545515328 +0100
***************
*** 27,34 ****
#endif
#include <cfloat>
-
#include <iostream>
#include "Array-util.h"
#include "byte-swap.h"
--- 27,34 ----
#endif
#include <cfloat>
#include <iostream>
+ #include <vector>
#include "Array-util.h"
#include "byte-swap.h"
***************
*** 1201,1206 ****
--- 1201,1330 ----
}
Matrix
+ Matrix::toepsolve (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 ();
+ volatile int typ = mattype.type ();
+
+ info = 0;
+
+ if (nr == 0 || nc == 0 || nr != b.rows ())
+ (*current_liboctave_error_handler)
+ ("matrix dimension mismatch solution of linear equations");
+
+ else if (typ != MatrixType::Toeplitz)
+ (*current_liboctave_error_handler) ("incorrect matrix type");
+ else
+ {
+ octave_idx_type b_nc = b.cols();
+ retval.resize(nc, b_nc);
+
+ OCTAVE_LOCAL_BUFFER (double, c1, nr - 1);
+ OCTAVE_LOCAL_BUFFER (double, c2, nc - 1);
+
+ for (octave_idx_type j = 0; j < b_nc; j++)
+ {
+ double r1, r2, r3, r5, r6;
+
+ r1 = elem(0);
+
+ if (r1 == 0.)
+ {
+ info = -2;
+ goto singular;
+ }
+
+ retval.elem (0, j) = b.elem(0,j)/r1;
+
+ if ( nr == 1 )
+ continue;
+
+ for (octave_idx_type nsub = 2; nsub <= nr; nsub++)
+ {
+ r5 = elem(nsub-1);
+ r6 = elem(nr*(nsub-1));
+
+ if (nsub > 2)
+ {
+ c1[nsub-2] = r2;
+ for (octave_idx_type i = 1; i <= nsub - 2; i++)
+ {
+ r5 = r5 + elem(i) * c1[nsub-i-1];
+ r6 = r6 + elem(nr*i) * c2[i-1];
+ }
+ }
+
+ r2 = -r5/r1;
+ r3 = -r6/r1;
+ r1 = r1 + r5 * r3;
+
+ if (r1 == 0.)
+ {
+ info = -2;
+ goto singular;
+ }
+
+ if (nsub > 2)
+ {
+ r6 = c2[0];
+ c2[nsub-2] = 0;
+
+ for (octave_idx_type i = 2; i <= nsub - 1; i++)
+ {
+ r5 = c2[i-1];
+ c2[i-1] = c1[i-1] * r3 + r6;
+ c1[i-1] += r6 * r2;
+ r6 = r5;
+ }
+ }
+
+ c2[0] = r3;
+ r5 = 0;
+ for (octave_idx_type i = 1; i <= nsub - 1; i++)
+ r5 = r5 + elem(i) * retval.elem(nsub-i-1, j);
+
+ r6 = ( b(nsub-1) - r5 )/r1;
+
+ for (octave_idx_type i = 1; i <= nsub-1; i++)
+ retval.elem(i-1, j) += c2[i-1] * r6;
+
+ retval.elem(nsub-1, j) = r6;
+ }
+ }
+ }
+
+ singular:
+ if (calc_cond)
+ (*current_liboctave_error_handler)
+ ("calculation of condition number estimate not yet implemented");
+
+ if (info)
+ {
+ if (sing_handler)
+ sing_handler (rcond);
+ else
+ (*current_liboctave_error_handler)
+ ("matrix singular to machine precision");
+ mattype.mark_as_rectangular ();
+ }
+
+ return retval;
+ }
+
+ /*
+
+ %!test
+ %! a = toeplitz(1:5,2:6);
+ %! A = matrix_type(a,'full');
+ %! assert(a\ones(5,1),A\ones(5,1),1e-10)
+
+ */
+
+ Matrix
Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
double& rcond, solve_singularity_handler sing_handler,
bool calc_cond) const
***************
*** 1297,1302 ****
--- 1421,1428 ----
(*current_liboctave_error_handler)
("unrecoverable error in dtrtrs");
}
+ else
+ mattype.mark_as_rectangular ();
}
}
else
***************
*** 1403,1408 ****
--- 1529,1536 ----
(*current_liboctave_error_handler)
("unrecoverable error in dtrtrs");
}
+ else
+ mattype.mark_as_rectangular ();
}
}
else
***************
*** 1651,1656 ****
--- 1779,1786 ----
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::Toeplitz)
+ retval = toepsolve (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)
*** ./liboctave/MatrixType.h.orig4 2006-10-27 03:45:55.000000000 +0200
--- ./liboctave/MatrixType.h 2006-11-28 16:53:23.932498640 +0100
***************
*** 47,52 ****
--- 47,53 ----
Banded_Hermitian,
Tridiagonal,
Tridiagonal_Hermitian,
+ Toeplitz,
Rectangular
};
***************
*** 131,137 ****
void mark_as_lower_triangular (void) { typ = Lower; }
! void mark_as_tridiagonal (void) {typ = Tridiagonal; }
void mark_as_banded (const octave_idx_type ku, const octave_idx_type kl)
{ typ = Banded; upper_band = ku; lower_band = kl; }
--- 132,138 ----
void mark_as_lower_triangular (void) { typ = Lower; }
! void mark_as_tridiagonal (void) { typ = Tridiagonal; }
void mark_as_banded (const octave_idx_type ku, const octave_idx_type kl)
{ typ = Banded; upper_band = ku; lower_band = kl; }
***************
*** 152,157 ****
--- 153,160 ----
void mark_as_unpermuted (void);
+ void mark_as_toeplitz (void) { typ = Toeplitz; }
+
MatrixType transpose (void) const;
private:
*** ./liboctave/MatrixType.cc.orig4 2006-10-03 22:07:56.000000000 +0200
--- ./liboctave/MatrixType.cc 2006-11-28 22:36:49.942467213 +0100
***************
*** 103,112 ****
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;
--- 103,144 ----
typ = MatrixType::Upper;
else if (lower)
typ = MatrixType::Lower;
! else
! {
! if (hermitian)
! typ = MatrixType::Hermitian;
! else
! typ = MatrixType::Full;
!
! bool toeplitz = true;
!
! for (octave_idx_type i = 0; i < nrows; i++)
! {
! double val = a.elem(i,0);
!
! for (octave_idx_type j = 1; j < ncols - i; j++)
! if (val != a.elem(i+j, j))
! {
! toeplitz = false;
! goto NotToeplitz;
! }
! }
!
! for (octave_idx_type i = 1; i < ncols; i++)
! {
! double val = a.elem(0,i);
!
! for (octave_idx_type j = 1; j < nrows - i; j++)
! if (val != a.elem(j, i+j))
! {
! toeplitz = false;
! goto NotToeplitz;
! }
! }
! NotToeplitz:
! if (toeplitz)
! typ = MatrixType::Toeplitz;
! }
}
else
typ = MatrixType::Rectangular;
***************
*** 165,174 ****
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;
--- 197,238 ----
typ = MatrixType::Upper;
else if (lower)
typ = MatrixType::Lower;
! else
! {
! if (hermitian)
! typ = MatrixType::Hermitian;
! else
! typ = MatrixType::Full;
!
! bool toeplitz = true;
!
! for (octave_idx_type i = 0; i < nrows; i++)
! {
! Complex val = a.elem(i,0);
!
! for (octave_idx_type j = 1; j < ncols - i; j++)
! if (val != a.elem(i+j, j))
! {
! toeplitz = false;
! goto NotToeplitz;
! }
! }
!
! for (octave_idx_type i = 1; i < ncols; i++)
! {
! Complex val = a.elem(0,i);
!
! for (octave_idx_type j = 1; j < nrows - i; j++)
! if (val != a.elem(j, i+j))
! {
! toeplitz = false;
! goto NotToeplitz;
! }
! }
! NotToeplitz:
! if (toeplitz)
! typ = MatrixType::Toeplitz;
! }
}
else
typ = MatrixType::Rectangular;
***************
*** 834,840 ****
if (t == MatrixType::Full || t == MatrixType::Diagonal ||
t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper ||
t == MatrixType::Lower || t == MatrixType::Tridiagonal ||
! t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular)
typ = t;
else
(*current_liboctave_warning_handler) ("Invalid matrix type");
--- 898,905 ----
if (t == MatrixType::Full || t == MatrixType::Diagonal ||
t == MatrixType::Permuted_Diagonal || t == MatrixType::Upper ||
t == MatrixType::Lower || t == MatrixType::Tridiagonal ||
! t == MatrixType::Tridiagonal_Hermitian || t == MatrixType::Rectangular
||
! t == MatrixType::Toeplitz || t == MatrixType::Unknown)
typ = t;
else
(*current_liboctave_warning_handler) ("Invalid matrix type");
*** ./src/DLD-FUNCTIONS/matrix_type.cc.orig4 2006-05-19 07:32:18.000000000
+0200
--- ./src/DLD-FUNCTIONS/matrix_type.cc 2006-11-28 17:31:21.280963609 +0100
***************
*** 44,50 ****
@deftypefnx {Loadable Function} address@hidden =} matrix_type (@var{a},
'lower', @var{perm})\n\
@deftypefnx {Loadable Function} address@hidden =} matrix_type (@var{a},
'banded', @var{nl}, @var{nu})\n\
Identify the matrix type or mark a matrix as a particular type. This allows
rapid\n\
! for solutions of linear equations involving @var{a} to be performed. Called
with a\n\
single argument, @code{matrix_type} returns the type of the matrix and caches
it for\n\
future use. Called with more than one argument, @code{matrix_type} allows the
type\n\
of the matrix to be defined.\n\
--- 44,50 ----
@deftypefnx {Loadable Function} address@hidden =} matrix_type (@var{a},
'lower', @var{perm})\n\
@deftypefnx {Loadable Function} address@hidden =} matrix_type (@var{a},
'banded', @var{nl}, @var{nu})\n\
Identify the matrix type or mark a matrix as a particular type. This allows
rapid\n\
! solutions of linear equations involving @var{a} to be performed. Called with
a\n\
single argument, @code{matrix_type} returns the type of the matrix and caches
it for\n\
future use. Called with more than one argument, @code{matrix_type} allows the
type\n\
of the matrix to be defined.\n\
***************
*** 87,92 ****
--- 87,95 ----
with specialized code. In addition the matrix can be marked as positive
definite\n\
(Sparse matrices only)\n\
\n\
+ @item 'toeplitz'\n\
+ The matrix is assumed to be toeplitz.\n\
+ \n\
@item 'singular'\n\
The matrix is assumed to be singular and will be treated with a minimum norm
solution\n\
\n\
***************
*** 169,174 ****
--- 172,182 ----
retval = octave_value ("Tridiagonal Positive Definite");
else if (typ == MatrixType::Hermitian)
retval = octave_value ("Positive Definite");
+ // Don't allow sparse toeplitz for now
+ #if 0
+ else if (typ == MatrixType::Toeplitz)
+ retval = octave_value ("Toeplitz");
+ #endif
else if (typ == MatrixType::Rectangular)
{
if (args(0).rows() == args(0).columns())
***************
*** 241,246 ****
--- 249,259 ----
mattyp.mark_as_rectangular ();
else if (str_typ == "full")
mattyp.mark_as_full ();
+ // Don't allow sparse toeplitz for now
+ #if 0
+ else if (str_typ == "toeplitz")
+ mattyp.mark_as_toeplitz ();
+ #endif
else if (str_typ == "unknown")
mattyp.invalidate_type ();
else
***************
*** 342,347 ****
--- 355,362 ----
retval = octave_value ("Permuted Lower");
else if (typ == MatrixType::Hermitian)
retval = octave_value ("Positive Definite");
+ else if (typ == MatrixType::Toeplitz)
+ retval = octave_value ("Toeplitz");
else if (typ == MatrixType::Rectangular)
{
if (args(0).rows() == args(0).columns())
***************
*** 384,389 ****
--- 399,406 ----
mattyp.mark_as_rectangular ();
else if (str_typ == "full")
mattyp.mark_as_full ();
+ else if (str_typ == "toeplitz")
+ mattyp.mark_as_toeplitz ();
else if (str_typ == "unknown")
mattyp.invalidate_type ();
else
***************
*** 447,453 ****
## FIXME
## Disable tests for lower under-determined and upper over-determined
! ## matrices and this detection is disabled in MatrixType due to issues
## of non minimum norm solution being found.
%!assert(matrix_type(speye(10,10)),"Diagonal");
--- 464,470 ----
## FIXME
## Disable tests for lower under-determined and upper over-determined
! ## matrices as this detection is disabled in MatrixType due to issues
## of non minimum norm solution being found.
%!assert(matrix_type(speye(10,10)),"Diagonal");
***************
*** 487,492 ****
--- 504,510 ----
%!assert(matrix_type([speye(9,11);[1,sparse(1,8),1,0]]),"Lower");
%!assert(matrix_type([speye(9,11);[1,sparse(1,8),1,0]]([2,1,3:10],:)),"Permuted
Lower");
%!assert(matrix_type(spdiags(randn(10,4),[-2:1],10,9)),"Rectangular")
+ %!assert(matrix_type(sparse(toeplitz([1,0,0,0,2],[0,0,0,3,4]))),"Toeplitz")
%!assert(matrix_type(1i*speye(10,10)),"Diagonal");
%!assert(matrix_type(1i*speye(10,10)([2:10,1],:)),"Permuted Diagonal");
***************
*** 525,530 ****
--- 543,550 ----
%!assert(matrix_type([speye(9,11);[1i,sparse(1,8),1i,0]]),"Lower");
%!assert(matrix_type([speye(9,11);[1i,sparse(1,8),1i,0]]([2,1,3:10],:)),"Permuted
Lower");
%!assert(matrix_type(1i*spdiags(randn(10,4),[-2:1],10,9)),"Rectangular")
+ %!assert(matrix_type(1i*sparse(toeplitz([1,0,0,0,2],[0,0,0,3,4]))),"Toeplitz")
+
%!test
%! a = matrix_type(spdiags(randn(10,3),[-1,0,1],10,10),"Singular");
***************
*** 539,544 ****
--- 559,565 ----
%!test
%! a = matrix_type(ones(10,10),"Singular");
%! assert(matrix_type(a),"Singular");
+ %!assert(matrix_type(toeplitz([1,2],[3,4])),"Toeplitz")
%!assert(matrix_type(triu(1i*ones(10,10))),"Upper");
%!assert(matrix_type(triu(1i*ones(10,10),-1)),"Full");
***************
*** 549,554 ****
--- 570,576 ----
%!test
%! a = matrix_type(ones(10,10),"Singular");
%! assert(matrix_type(a),"Singular");
+ %!assert(matrix_type(1i*toeplitz([1,2],[3,4])),"Toeplitz")
*/
2006-11-28 David Bateman <address@hidden>
* MatrixType.h (matrix_type): Add Toeplitz to enum
(mark_as_toeplitz): New function.
* MatrixType.cc (MatrixType::MatrixType(const Matrix&)): Test for
toeplitz matrices. Don't give a warning for unknown matrix type.
(MatrixType::MatrixType(const ComplexMatrix&)): ditto.
* dMatrix.cc (Matrix::toepsolve(...)): New function to solve for
toeplitz matrices, based on code by Gorazd Brumen.
(Matrix::solve(...): Use it.
(Matrix::utsolve(...), Matrix::ltsolve(...)): Mark singular
matrices, so that the solver falls back to lssolve function.
* CMatrix.cc (ComplexMatrix::toepsolve(...)): ditto.
(ComplexMatrix::solve(...)): Use it.
(ComplexMatrix::utsolve(...), ComplexMatrix::ltsolve(...)): Mark
singular matrices, so that the solver falls back to lssolve function.
* dMatrix.h (toepsolve(...)): declaration of new function.
* CMatrix.h (toepsolve(...)): ditto.
2006-11-28 David Bateman <address@hidden>
* DLD-FUNCTIONS/matrix_type.cc: Allow full toeplitz matrices.
- solvetoep, Gorazd Brumen, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, Gorazd Brumen, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, Gorazd Brumen, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep,
David Bateman <=
- Re: solvetoep, Gorazd Brumen, 2006/11/28