From maintainers-request at octave dot org Wed Dec 1 03:41:47 2004 Subject: Re: full/sparse/banded/triangular matrices... From: David Bateman To: "John W. Eaton" Cc: maintainers at octave dot org Date: Wed, 1 Dec 2004 10:38:48 +0100 According to John W. Eaton (on 12/01/04): > On 30-Nov-2004, David Bateman wrote: > > | 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. > > Yes, we will need to come up with a systematic way of handling all > these so they can be generated automatically (at least as much as is > possible). > > | 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 > > Wouldn't this last one only have to convert to Full if Scalar is > nonzero? In that case, we really do want the storage class to be a > dynamic attribute of the array class instead of a fixed type attribute > that would be required if we had separate "top-level" classes for each > storage type. Yes, I caught the Scalar zero case in the table I sent though. However, I don't see why this implicitly means that we are forced to having dynamic attributes... I see there as being two approaches, the better one that we both agree in the longer term is to have the storage classes as a dynamic property of the Array class. However, it seems to me that that is a really long term goal, but we can still get a lot of the benefits of such storage classes, at least in terms of computational speed with using seperate top level storage classes, of the form I mentioned. With using the differentiation at the top level class the dynamic conversion is handled at the top level. So in the short term where the triangular matrices are just a flag on the full (or sparse) matrix storage, the binop for add with a scalar might look like static octave_value oct_binop_add (const octave_value& a1, const octave_value& a2) { CAST_BINOP_ARGS (const octave_triangular_matrix &, const octave_scalar&); double d = v2.double_value (); if (d == 0) return octave_value (a1); else return octave_value (v1.matrix_value () + d, v1.attribute); } In fact in the long term the seperate top-level classes could also stay, and use different storage classes of the underlying Array class. In that case the dynamic change between classes would be handled in a similar manner. The advantage of this approach is that if the top level classes are defined like template octave_triangular_matrix : puble octave_matrix { ... }; and inherit from the underlying matrix class, then we only need to explicitly include the operations that preserve the attribute of the matrix class. This remains true even with additional storage classes in the Array type. > | So longer term I see the best situation as having Array evolve into > | something like > | > | template > | class > | Array > | { > | protected > | class FullRep > | { > | ... > | }; > | > | class SparseRep > | { > | ... > | }; > > Yes, I think this is a good direction. The rep classes don't have to > be nested, they are just handled that way now because it seemed good > to hide them. > > I'm not sure that we want to define math operations on the Array > object, so the operators would be separate. True, in fact if the operators are on derived classes of Array then it would be a single operator with special cases for each storage class. > To be able to take advantage of external libraries like lapack, we > would probably need to be able to expose some details about the > internals (we do that already, with fortran_vec). > > For efficiency, the operators would be defined on If we only have full, sparse and banded storage classes, then fortran_vec plus the size of the band for banded matrices is sufficient. The sparse class would need to expose the row and column indexing... > > | class BandRep > | { > | ... > | }; > | > | public: > | enum storage_type > | { > | DENSE, > | SPARSE, > | BANDED > | }; > | enum attribute_type > | { > | UNKNOWN, > | FULL, > | UPPER, > | LOWER, > | }; > > Do we really need these enums? Wouldn't you want to dispatch to the > rep class and have it do the right thing instead of using the enum > value? It seems redundant to have both and switches are bad form in > C++ when a dispatch would do the job (think about the trouble caused > by the enum when you add a new rep type that requires a new enum value). The reason I see for an enum is that imagine a case like a = zeros(n,n) ## Create n-by-n banded matrix with banded of +/-k for i = -k:k if (i <= 0) a ((0:(n+k-1))*(n+1) - k + 1) = foo (n+k); # Returns n+k < n values else a ((k:(n-1))*(n+1) + 1) = foo (n-k); # Returns n-k < n values endif endfor b = a \ randn(n,1); We don't know that "a" is banded. In fact we have no information about it, and its attribute is known. However, if the "\" operator attempts to determine the attribute of unknown matrices, then we can benefit from the acceleration for banded operations in any case. However this needs caching of the attribute. If the attribute wasn't cached then each time "\" was called on a full matrix its attribute would be checked, even if it had already been determined. I don't see how to treat this with some sort of enum.. 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