[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: solvetoep
|
From: |
Gorazd Brumen |
|
Subject: |
Re: solvetoep |
|
Date: |
Tue, 28 Nov 2006 20:32:56 +0100 |
|
User-agent: |
Thunderbird 1.5.0.8 (X11/20061116) |
Hi David,
I am willing to have a look at it and work on it and I think I
can provide sth. within the next days.
Here a couple of questions:
a) What will be the means of communication? Should I post on the
mailing list a a patch against the cvs tree?
b) I am not as wizardy as you are about this code, so it will take
me a couple of days at least to get it running (or as in Blow: "I can do
it, but it will take me at least a year").
c) As I can see, you havent made any changes to the cvs tree yet
(I dont know exactly how this thing goes), so, should I start from
your patch against the cvs and work from there or is your
patch just a reference to start from?
Gorazd
David Bateman wrote:
> 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
>
>
>
> ------------------------------------------------------------------------
>
> *** ./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")
>
> */
>
--
Gorazd Brumen
Mail: address@hidden
WWW: http://www.gorazdbrumen.net/
PGP: Key at http://pgp.mit.edu, ID BCC93240
- 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 <=
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, David Bateman, 2006/11/28
- Re: solvetoep, Gorazd Brumen, 2006/11/28