From maintainers-request at octave dot org Tue Nov 30 16:51:25 2004 Subject: full/sparse/banded/triangular matrices... From: David Bateman To: maintainers at octave dot org Date: Tue, 30 Nov 2004 23:51:42 +0100 In thinking a little bit about what to do with Andy's spinv function, I found myself asking lots of questions about how to address triangular matrices, since the heart of his function is the back subsitution on a sparse upper triangular matrix. It seemed to me that this routine has applications beyond just this function, for example in a cholesky based solver. To achieve something like this it would make sense to have sparse triangular matrices, perhaps in a manner like that in ov-re-tri.cc in octave-forge. The risk is an explosion in the number of mixed operators. This is even more true if banded matrices are addressed at the same time. Why do we want sparse/triangular/banded matrices? What is the ideal implementation of full, sparse, triangular and banded matrices? What should the classes in liboctave look like? Is there some way to get to that vision in easier steps? Why seems to be clear. Reduced memory usage and faster computation. However, the two issues aren't necessarily related. There is nothing to stop us storing the zero values, as a fast way of at least getting some of the benefits of banded or triangular matrices in terms of more rapid computation. The LAPACK routines for banded routines can easily be used in this case just by selecting the LDA parameter appropriately and there are special version of the LAPACK triangular back substitute functions for a full matrix with only the upper/lower part defined. In the longer run there should be ideally three different storage classes, Full, Sparse and Banded. There is only a factor of two saving in memory use in a triangular storage class and I'm not sure its worth having a special storage class for this, but rather flag for the Full and Sparse classes as being triangular (an attribute in matlab terminology). This means that all of the operators on sparse/full matrices also directly apply to triangular matrices, with only the things like "/", "\", etc needing special casing for acceleration. There would also need to be a means of preserving, modifying or destroying the attribute during other operations. For example UpperTriangular .* Full -> UpperTriangular UpperTriangular .' -> LowerTriangular UpperTriangular + Scalar -> Full At at the outset something similar might be done for a banded matrix. This might be either done by defining classes for triangular/ banded matrices in the same way as in ov-re-tri.cc in octave-forge, and expanding the number of specialized operators on these class to define the preservation of the attribute, or the attribute might be explicitly included in the Array class and the operators in MArray, dNDArray, etc modified to preserve the attributes. In the long run it would be better to do it in the Array class as then even applications that don't like to liboctinterp.so can have access to the benefits of these classes. But in the short term it is probably easier and certainly less distruptive to do it in the same manner as ov-re-tri.cc. So longer term I see the best situation as having Array evolve into something like template class Array { protected class FullRep { ... }; class SparseRep { ... }; class BandRep { ... }; public: enum storage_type { DENSE, SPARSE, BANDED }; enum attribute_type { UNKNOWN, FULL, UPPER, LOWER, }; protected: typename Array::FullRep *fullrep; typename Array::SparseRep *sparserep; typename Array::BandedRep *bandedrep; storage_type storage; attribute_type attribute; dim_vector dimensions; idx_vector *idx; int idx_count; public: ... Array (const Array& a) : attribute (a.attribute), dimensions (a.dimensions), idx (0), idx_count (0) { switch (a.storage) { case SPARSE: sparserep = a.sparserep; sparserep->count++; break; case BANDED: bandedrep = a.bandedrep; bandedrep->count++; break; case DENSE: default: fullrep = a.fullrep; fullrep->count++; break; } } ... }; with the consequent changes to all of the operators in the derived classes. This would allow such crazy things as banded cell arrays as being possible. However such changes are likely to be very disruptive as it imposes changes at the very lowest level of how octave currently does things, and so is probably not something that should be done before 3.0. Should the other approach of having the classes class octave_triangular_matrix : public octave_matrix { ... }; class octave_triangular_complex_matrix : public octave_complex_matrix { ...}; class octave_banded_matrix : public octave_matrix { ... }; class octave_banded_complex_matrix : public octave_complex_matrix { ... }; class octave_triangular_sparse_matrix : public octave_sparse_matrix { ... }; class octave_triangular_sparse_complex_matrix : public octave_sparse_complex_matrix { ... }; be implemented? I believe yes, as the implementation is relatively straight forward. I didn't include banded sparse matrices as there seems to be no point. If we missed a few of the attribute preserving operations in the op-*.cc files then it just falls back to the current case of converting to a full/sparse version.. Also even the majority of the attribute preserving functions use the existing code, with just the returned octave_value being adapted appropriately, so there is a low risk of introduces major errors. This also means that only cases where the attribute is preserved need to be treated, and a conversion to a full matrix can fall back to the full matrix code. Then it just comes done to making tables like Unary Ops OP INPUT UP LO BND --- -- -- --- ! UP LO BND + UP LO BND - UP LO BND transpose LO UP BND hermitian LO UP BND Binary OPS UP/UP UP/LO UP/BND UP/FULL UP/SCAL ---------- ----- ----- ------ ------- ------- + UP FULL FULL FULL FULL - UP FULL FULL FULL FULL * UP FULL FULL FULL UP / ** ** ** FULL UP ^ NA NA NA NA UP \ ** ** ** ** ** < UP FULL FULL FULL FULL > UP FULL FULL FULL FULL == FULL FULL FULL FULL FULL >= FULL FULL FULL FULL FULL <= FULL FULL FULL FULL FULL .* UP UP UP/BND UP UP .^ FULL FULL FULL FULL UP (!=0) ** Accelerated special case to identify what to do for the op-*.cc files, and the rest is straight forward. If we want to get a little bit more complicated we could include the acceleration for things like BND+BND or UP+UP, that don't calculate the zero terms for a bit of extra speed. Though its not these operations that are going to dominant any calculation, but rather the "/" and "\" operators, so I'm not sure I see any advantage. As a final simplification the constructors of the octave_values might be modified like enum attribute_type { UNKNOWN, FULL, BANDED UPPER, LOWER, }; octave_value::octave_value (const NDArray& a, attribute_type typ = UNKNOWN, int band_size = -1) { switch (typ) { case BANDED: rep = new octave_banded_matrix (a, band_size); break; case UPPER: rep = new octave_triangular_matrix (a, octave_triangular_matrix::UPPER); break; case LOWER: rep = new octave_triangular_matrix (a, octave_triangular_matrix::LOWER); break; default rep = new octave_matrix (a); break; } rep->count = 1; maybe_mutate (); } and similarly for all of the other constructors, so that the process of writing the op-*.cc files might then be simplified to a cut-and-paste with a new macro to call the constructor with the correct matrix type. The default args allow the exist code to work as is. Finally the functions chol, lu, triu, tril would need to be modified to return the correct matrix type, and if a matrix is flagged as unknown the "/" and "\" operators should probably attempt to determine the matrix type, so the attribute should be cached for the octave_matrix and octave_complex_matrix type. Anyway, these are just my late night ruminations on the topic. It all seems relatively straight forward, and seems a natural addition at the same time as sparse matrices in the octave core. But its up to John to see what he whats to do. In any case I hope this rather long mail might kick off a small discussion. Regards 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