Re: matrix_type problems

[David Bateman]

Vic Norton wrote:
> The idea of a matrix_type function makes a lot of sense to me, but what we 
> have now doesn't seem to work in Octave 2.9.9.
> 
> For example:
> 
>    octave> A = eye(5)
>    A =
>    
>       1   0   0   0   0
>       0   1   0   0   0
>       0   0   1   0   0
>       0   0   0   1   0
>       0   0   0   0   1
> 
>    octave> matrix_type(A)
>    ans = Upper
>    octave> A(2,1) = 1
>    A =
>    
>       1   0   0   0   0
>       1   1   0   0   0
>       0   0   1   0   0
>       0   0   0   1   0
>       0   0   0   0   1
>    
>    octave> matrix_type(A)
>    ans = Upper
> 
> 
> It would be great if octave could take advantage of upper triangular 
> structure, but the above exercise doesn't give me much confidence.
> 
> I would also like to see a simple test of whether a matrix is upper 
> triangular. Obviously
> 
>    octave> strcmp(matrix_type(A), "Upper")
>    ans =  1
> 
> doesn't work, because matrix_type (with one argument) doesn't work.
> 
> Regards,
> 
> Vic

Ouch, seems I forgot to invalidate the matrix type on an assign for full
matrices.. Try the following instead

A = speye(5)
matrix_type(A)
A(2,1) = 1
matrix_type(A)

and you'll see things work fine. Octave does take advantage of the
triangular structure, that remains true. Rather what you've done is
spotted a bug, that should only occur in the circumstance where the
matrix type is resolved to something prior to an assignment that changes
it type. That might happen if you explicitly call matrix_type as you've
done, but

A = eye(5); x = A \ ones(5,1); A(2,1) = 1;

will also have the same effect as the \ operator will resolve the matrix
type. Anyway attached is the trivial fix.. After applying this patch I get

octave:1> a = eye(5)
a =

   1   0   0   0   0
   0   1   0   0   0
   0   0   1   0   0
   0   0   0   1   0
   0   0   0   0   1

octave:2> matrix_type(a)
ans = Upper
octave:3> a(2,1) = 1
a =

   1   0   0   0   0
   1   1   0   0   0
   0   0   1   0   0
   0   0   0   1   0
   0   0   0   0   1

octave:4> matrix_type(a)
ans = Lower

As for the test that a matrix is upper triangular, A =
matrix_type(A,'Upper') will force the matrix A to be of type 'Upper',
without having to have the type probed. This gives a slight gain in
performance and is valid for the output of such things as a Cholesky
factorization. However, this isn't want you wanted. To probe whether a
matrix is upper triangular something like
strcmp(matrix_stype(A),'Upper') would be more appropriate.

Cheers
David

2006-12-03  David Bateman  <address@hidden>

            * ov-base-mat.cc (void octave_base_matrix<MT>::assign
            (const octave_value_list&, const MT&)): Invalidate matrix
            type on assignment

Index: src/ov-base-mat.cc
===================================================================
RCS file: /cvs/octave/src/ov-base-mat.cc,v
retrieving revision 1.37
diff -c -r1.37 ov-base-mat.cc
*** src/ov-base-mat.cc  14 Jul 2006 20:29:35 -0000      1.37
--- src/ov-base-mat.cc  3 Dec 2006 19:59:44 -0000
***************
*** 199,204 ****
--- 199,208 ----
      matrix.set_index (idx(i).index_vector ());
  
    ::assign (matrix, rhs, MT::resize_fill_value ());
+ 
+ 
+   // Invalidate the matrix type
+   typ.invalidate_type ();
  }
  
  template <class MT>