Re: full/sparse/banded/triangular matrices...

[David Bateman]

According to John W. Eaton <address@hidden> (on 12/01/04):
> On 30-Nov-2004, David Bateman <address@hidden> 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<T> 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<T> 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 <class T> 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<T> type.

> | So longer term I see the best situation as having Array<T> evolve into
> | something like
> | 
> | template <class T>
> | 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<T> 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                                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