Re: solvetoep

[Gorazd Brumen]

Man, that was fast.

I will look into the papers that do the so-called "superfast" toeplitz
solution (that is for large systems only) and see if it can be
implemented. from the papers I see it is reasonably complicated +
there are some stability issues.

Keep you posted on my progress,

Gorazd


David Bateman wrote:
> 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.

-- 
Gorazd Brumen
Mail: address@hidden
WWW: http://www.gorazdbrumen.net/
PGP: Key at http://pgp.mit.edu, ID BCC93240