From maintainers-request at octave dot org Wed Dec  8 06:36:43 2004
Subject: Re: permutations of matrix to triangular form?
From: Andy Adler <adler at site dot uottawa dot ca>
To: David Bateman <David dot Bateman at motorola dot com>
cc: maintainers at octave dot org
Date: Wed, 8 Dec 2004 07:37:16 -0500 (EST)

On Wed, 8 Dec 2004, David Bateman wrote:

> Thinking yet again about implementing the triangular/banded solvers,
> I re-examined what matlab did as described at
>
> http://www.mathworks.com/access/helpdesk/help/techdoc/ref/mldivide.html
>
> Where it seems to test the type of the matrix (and presumably caches
> it) before solving with the "\" and "/" operators. Presumably
> something similar is done in a few other places like eig, etc.  Doing
> it matlabs way has the advantage that the banded and triangular types
> are not explicitly needed.
>
> Last night I was investigating how to do something similar with the
> sparse class, and basically just wrote a small piece of code to recognize
> the type of a matrix (I attach it here if you are interested, though it
> needs to be used with the sparse class I've been writing). The interest
> of this was to take Andy's back substitution code from his sparse inverse
> and use it directly in the "\" and "/" operations, and at the same time
> treat banded with LAPACK and diagonal matrices.

David,

This doesn't answer your question, but I think that this approach
(of having a small number of storage types (full,sparse) with 'hints'
about the solver) offers a nice set of possibilities.

Currently, your approach autodetects the matrix type, but these
properties could be made available to the user, such as

   a= sparse(diag(1:3));
   get(a, 'storage-type')
   a= set(a, 'storage-type', 'diagonal')

or some similar syntax. This would allow

1. Users who already know something about the matrix to write more
   efficient code

2. The possibilitiy of including more detailed information to the
   solver. (ie. iterative, precision thresholds )

   a= set(a, 'solver', 'iterative');

   Currently, doing this in Matlab with sparse code makes it look very
   unlike the simple solver cases. For example, in some finite element
   code I'm working on (http://cvs.sf.net/viewcvs.py/eidors3d/eidors3d/algorithms/np_2003/forward_solver.m)
   instead of
        V= E\I;
   one does
        %Permute the rows and columns to make the factors sparser
        E = E(pp,pp);
        In = I(pp,:);
        rr(pp)=1:max(size(pp));
        U = cholinc(E,tol/10);
        q_c =  U' \ In;
        Vn = U \ q_c;
        %De-permute the result for Cholesky
        V = Vn(rr,:);

I know you've suggested a similar idea before. I'm just pointing out that
it works with the framework you're working on now.

Another issue: There is the possibility of user's giving the wrong info
about the matrix.  My opinion would be that that is their problem.

--
Andy Adler <adler AT site.uOttawa.ca> 1(613)562-5800x6218