Re: matrix_type problems

[Vic Norton]

On 12/3/06, at 3:29 PM +0100, David Bateman wrote:
> 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.

In fact, David, I do want to force the matrix to be of type 'Upper' if that 
will speed up performance.

Here is my situation. I have an iterative routine that updates QR 
factorizations of submatrices of an m-by-n matrix A. At each iteration I have 
an index vector J of length nJ <= min( m, n ), and I need to compute
   e = A( : , J ) \ u;
This can be done with 2 * nJ * nJ * ( m - nJ/3 ) flops. (flop = floating point 
operation.) The vector u is fixed, but J and nJ change from one iteration to 
the next, usually by the addition or deletion of a single index. The matrix A( 
: , nJ ) always has rank nJ.

In my algorithm I maintain a QR factorization of A( : , J ) in two arrays, QQ 
and RR, of sufficient size. Thus A( : , J ) = QQ( : , 1:nJ ) * RR( 1:nJ , 1:nJ 
), so that e can be computed as
   uQ = QQ( : , 1:nJ )' * u;   e = RR( 1:nJ , 1:nJ ) \ uQ;
This computation of e can be done with nJ * ( 2 * m + nJ ) flops, if one takes 
advantage of the triangular structure of RR( 1:nJ , 1:nJ ).

I do want to take advantage of the triangular stucture of RR. Since a typical 
update of the QR factorization takes at most 6 * nJ * ( m + nJ ) flops, the QR 
method should be much faster than the direct, e = A( : , J ) \ u, computation, 
especially for large nJ.

My question is this. Should there be a speed advantage if I do
   RR( 1:nJ , 1:nJ ) = matrix_type( RR( 1:nJ , 1:nJ ) , "upper" );
   uQ = QQ( : , 1:nJ )' * u;   e = RR( 1:nJ , 1:nJ ) \ uQ;
on each iteration rather than doing just the last line? Or could I accomplish 
the same thing by declaring
   RR = matrix_type( RR , "upper" );
before the iterations start? Would this line force all subsequent  RR( 1:nJ , 
1:nJ )'s to be treated as upper triangular when e is computed?

Thanks for your advice.

Regards,

Vic