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<T> 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<T> 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<T> evolve into

something like

template <class T>

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<T>::FullRep *fullrep;

typename Array<T>::SparseRep *sparserep;

typename Array<T>::BandedRep *bandedrep;

storage_type storage;

attribute_type attribute;

dim_vector dimensions;

idx_vector *idx;

int idx_count;

public:

...

Array (const Array<T>& 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

[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