According to Andy Adler <

[hidden email]> (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

[hidden email]
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