octave-maintainers
[Top][All Lists]
Advanced

[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 18:10:57 +0100
User-agent: Thunderbird 1.5.0.7 (X11/20060921)

Gorazd Brumen wrote:
> Hi David, all,
>
> OK. I am a bit inexperienced and have only recently realized the
> following points:
>
> a) The code is a C++ rewrittement of the netlib solver. Perhaps it
> would be better to use netlib fortran code and then do a
> C++ wrapper around it, since then we dont have to "worry" about the
> netlib development more. If some better algorithm comes into the
> netlib - voila, our code changes automagically.
>   
This might be better, though some benchmarks should be done. If your is
faster I'd stick with your code, as it removes an external dependency..

> b) The code is O(n^2) which is better than the general O(n^3),
> it is from netlib (supposedly good, have no idea how good the netlib
> code is) and uses a lot less space than the general inverse 2*n, instead
> of the full n^2 through
>
> toeplitz (c,r)\b
>
> The matter of fact is that there exist even better routines, which
> do the job in O(n log ^2 (n)) time, but with a big coefficient and
> are better than the proposed solver for matrices of sizes > 256.
> I have not yet seen that code, but have read the papers about it.
>   

Maybe make the code polymorphic in that case. Treat one way if n < 256
and another elsewise.. Though benchmarks are again needed to find the
right change over point.


> c) I dont quite know what you mean by sparse.
>   
Matrices with lots of zeros that aren't stored.. See the "sparse" function..

> Please, suggest what is to be done.
>   
Based on the code you sent, I made a quick first attempt at a patch for
full toeplitz matrices. The toepsolve functions need to test for
singular matrices and flag them in the info variable and it might be
good if there was code to calculate the condition number of the toeplitz
as well (thus the warning in the build of the toepsolve functions).
These also need to be expanded to handle the sparse matrices as well.

In any case there appears to be a reversal of the return vector
somewhere in the code, and a spurious warning. However, this almost
works.... See for example.

octave:1> a = toeplitz(1:4,5:8)
warning: toeplitz: column wins diagonal conflict
a =

   1   6   7   8
   2   1   6   7
   3   2   1   6
   4   3   2   1

octave:2> matrix_type(a)
ans = Toeplitz
octave:3> a \ ones(4,1)
ans =

   0.052023
   0.034682
   0.023121
   0.202312

octave:4> A = matrix_type(a,'full')
warning: Invalid matrix type
A =

   1   6   7   8
   2   1   6   7
   3   2   1   6
   4   3   2   1

octave:5> matrix_type(A)
ans = Full
octave:6> A \ ones(4,1)
ans =

   0.202312
   0.023121
   0.034682
   0.052023

<snip>

octave:14> n = 1024; a = toeplitz(1:n,2:(n+1));
warning: toeplitz: column wins diagonal conflict
octave:15> matrix_type(a)
ans = Toeplitz
octave:16> tic; a \ ones(n,1); toc
Elapsed time is 0.143250 seconds.
octave:17> A = matrix_type(a,'full');
warning: Invalid matrix type
octave:18> tic; A \ ones(n,1); toc
Elapsed time is 3.116371 seconds.

It also seems that this is a big speed improvement for these cases.. In
any case would you like to take a look at reversing the output so that
it is the same for an LU solution of this, fixing the singular flag and
condition number calculation, adding sparse versions, benchmarking, etc,
then I believe this patch is a good start...

Regards
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

*** ./liboctave/CMatrix.cc.orig4        2006-11-15 20:52:45.000000000 +0100
--- ./liboctave/CMatrix.cc      2006-11-28 17:57:28.927993977 +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,1632 ----
  }
  
  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 ();
+ 
+   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);
+         retval.elem (0, j) = b.elem(0,j)/r1;
+   
+         if ( nr == 1 ) 
+           continue;
+ 
+         for (octave_idx_type nsub = 2; nsub <= nr; nsub++)
+           {
+             r5 = elem(nr*(nsub-1));
+             r6 = elem(nsub-1);
+ 
+             if (nsub > 2) 
+               {
+                 c1[nsub-2] = r2;
+                 for (octave_idx_type i = 1; i <= nsub - 2; i++) 
+                   {
+                     r5 = r5 + elem(nr*i) * c1[nsub-i-1];
+                     r6 = r6 + elem(i) * c2[i-1];
+                   }
+               }
+ 
+             r2 = -r5/r1;
+             r3 = -r6/r1;
+             r1 = r1 + r5 * r3;
+     
+             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(nr*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;
+           }
+       }
+     }
+ 
+   return retval;
+ }
+ 
+ ComplexMatrix
  ComplexMatrix::utsolve (MatrixType &mattype, const ComplexMatrix& b, 
                        octave_idx_type& info, double& rcond, 
                        solve_singularity_handler sing_handler,
***************
*** 2033,2038 ****
--- 2121,2128 ----
      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 17:57:22.735186103 +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,1294 ----
  }
  
  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 ();
+ 
+   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);
+         retval.elem (0, j) = b.elem(0,j)/r1;
+   
+         if ( nr == 1 ) 
+           continue;
+ 
+         for (octave_idx_type nsub = 2; nsub <= nr; nsub++)
+           {
+             r5 = elem(nr*(nsub-1));
+             r6 = elem(nsub-1);
+ 
+             if (nsub > 2) 
+               {
+                 c1[nsub-2] = r2;
+                 for (octave_idx_type i = 1; i <= nsub - 2; i++) 
+                   {
+                     r5 = r5 + elem(nr*i) * c1[nsub-i-1];
+                     r6 = r6 + elem(i) * c2[i-1];
+                   }
+               }
+ 
+             r2 = -r5/r1;
+             r3 = -r6/r1;
+             r1 = r1 + r5 * r3;
+     
+             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(nr*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;
+           }
+       }
+     }
+ 
+   return retval;
+ }
+ 
+ Matrix
  Matrix::utsolve (MatrixType &mattype, const Matrix& b, octave_idx_type& info,
                double& rcond, solve_singularity_handler sing_handler,
                bool calc_cond) const
***************
*** 1651,1656 ****
--- 1739,1746 ----
      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 18:00:39.516981029 +0100
***************
*** 105,112 ****
        typ = MatrixType::Lower;
        else if (hermitian)
        typ = MatrixType::Hermitian;
!       else if (ncols == nrows)
!       typ = MatrixType::Full;
      }
    else
      typ = MatrixType::Rectangular;
--- 105,143 ----
        typ = MatrixType::Lower;
        else if (hermitian)
        typ = MatrixType::Hermitian;
!       else 
!       {
!         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::Full;
!       }
      }
    else
      typ = MatrixType::Rectangular;
***************
*** 168,174 ****
        else if (hermitian)
        typ = MatrixType::Hermitian;
        else if (ncols == nrows)
!       typ = MatrixType::Full;
      }
    else
      typ = MatrixType::Rectangular;
--- 199,236 ----
        else if (hermitian)
        typ = MatrixType::Hermitian;
        else if (ncols == nrows)
!       {
!         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::Full;
!       }
      }
    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");
--- 896,903 ----
    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)
      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")
  
  */
  

reply via email to

[Prev in Thread] Current Thread [Next in Thread]