Re: permutations of matrix to triangular form?

[David Bateman]

According to Andy Adler <address@hidden> (on 12/08/04):
> 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.

Not my idea yet, as all I've done is implement detecting of the matrix
type as per the matlab mldivide spec and in the same precedence. Including
that in the sparse solvers themselves is the next step

> 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')

The matlab way is that diagonal and banded matrices are only available
in the sparse type, as why would you store them any other way. Finding
a diagonal, banded or unpermuted triangular matrix is a releatively
low cost operation, and so allowing the user to specify the type just
introduces the possiblity that they'll get it wrong without a huge
speed gain. I'd hesitate to do that.

> 2. The possibilitiy of including more detailed information to the
>    solver. (ie. iterative, precision thresholds )
> 
>    a= set(a, 'solver', 'iterative');

This is what John was against in a previous thread as the setting the 
solver type can be disassociated from the operation itself and increases
the potential for confusion.

>    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 agree this is ugly, but it is safer. Perhaps the way around this is
to introduce mldivide/mrdivide functions (we'll have to do it anyway for
matlab style "@" classes), but with an additional argument that can force
a particular solver to be used. In this manner, we remain compatiable,
don't disassociate the command for the solver type from the operator
itself and make the above simplier. Which is simplier


    E = set(E, 'solver', 'cholesky');
    V= E\I;

or

    V = mldivide (E, I, 'cholesky')

This of course means that "mldivide" would have to be explicitly used
in such cases.

> 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.

Its too easy in your scheme to have the setting of the solver far away
from its use, with potential to introduce horrible bugs that are difficult
to find. One debugging aid in that case is the spparms('spmoni') flag
that can be used to detect the solver type used. Then the problem becomes
one of educating users. I what a simple life, so I feel that its probably
still bad to have such a scheme.

Cheers
David

-- 
David Bateman                                address@hidden
Motorola CRM                                 +33 1 69 35 48 04 (Ph) 
Parc Les Algorithmes, Commune de St Aubin    +33 1 69 35 77 01 (Fax) 
91193 Gif-Sur-Yvette FRANCE

The information contained in this communication has been classified as: 

[x] General Business Information 
[ ] Motorola Internal Use Only 
[ ] Motorola Confidential Proprietary