From maintainers-request at octave dot org Wed Dec 8 07:01:47 2004 Subject: Re: permutations of matrix to triangular form? From: David Bateman To: Andy Adler Cc: David Bateman , maintainers@octave.org Date: Wed, 8 Dec 2004 13:59:31 +0100 According to Andy Adler (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 " at " 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 David dot Bateman at motorola dot com 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