Glas List: Re: [glas] multi-dimensional arrays?

From: Peter Gottschling (pgottsch_at_[hidden])
Date: 2006-12-18 09:28:48

Hi Laurent,

Concerning nested matrix and vector types, we need to ask ourselves 
first if we need special algorithms to handle it or if a least in some 
cases the operators on the values already use the right algorithm.

Second, is the use of operators efficient enough performance-wise?  In 
a matrix multiplication c(i, k)+= a(i, j) * b(j, k) is correct for any 
type with consistent operators but not very efficient if a, b, and c 
are matrices. Unless we have really sophisticated expression templates.

Next question, do we need special concepts?  Say e.g. 
RecursiveTransposable for a recursive function transpose

How would it look like?

concept RecursiveTransposable<typename Matrix>
{
    where RecursiveTransposable<Matrix::value_type>;
    // ...
};

would recursively descent and break if it reaches a scalar value that 
has no value_type.

Conversely, something like:

auto // ???
concept RecursiveTransposable<typename Matrix>
{
     Matrix transpose(Matrix );
}

Now we say this is defined for arithmetic scalars

template<typename T>
     where Arithmetic<T>
concept_map RecursiveTransposable<T> {}

And now enable it from bottom up:

template<typename T>
     where RecursiveTransposable<T>
concept_map RecursiveTransposable<matrix<T> > {}

Extending this to block-matrices,
(A B)
(C D)
where A, B, C, and D are matrices of different type. Say we have a type 
for heterogeneous blocks

template <typename MA, typename MB, typename MC, typename MD>
hb_matrix { ... };

We can define:
template <typename MA, typename MB, typename MC, typename MD>
     where RecursiveTransposable<MA>
         &&  RecursiveTransposable<MB>
         &&  RecursiveTransposable<MC>
         &&  RecursiveTransposable<MD>
concept_map RecursiveTransposable<hb_matrix<MA, MB, MC, MD> > {}

However, this are only some preliminary thoughts.  We need to look, 
which algorithms we want to implement recursively.  And it looks like  
that there can be a great difference between a syntactically correct 
implementation and a high-performance implementation.

Cheers,
Peter

On 15.12.2006, at 05:49, plagne.laurent_at_[hidden] wrote:

> Hi Karl,
>
> IMHO, the matrix< matrix<..> > can be more or less difficult to 
> implement
> depending on the constraints you accept. As an example, if you decide 
> that the
> library must remain open to third party codes and handle different 
> kind of
> matrix data storage : this is pretty easy to do with simple matrix type
> (matrix<double>) because the return type of A(i,j) will be a simple 
> ref to
> double (or const ref). In the 2D case, matrix< matrix<double> > A, the 
> accessor
> operator A(i,j) will be naturally a ref to matrix<double>. And you 
> loose the
> compatibility with third party data elements...
> One can avoid this problem but it is not a trivial issue.
>
> Another point in nested matrix issue is that it can help (again IMHO) 
> to
> validate the very impressive work the GLAS team has done with vector 
> space
> concepts (HI Peter). For example, the 0 and 1 matrix counter-part are 
> to be
> efficiently treated.
>
> Last point is maybe not a so general issue : In my case, I have to 
> handle sparse
> matrix with a complicated but highly structured sparsity pattern. This 
> pattern
> can be express with a reduced set of simple sparsity patterns : e.g A
> =Dense<Sparse<Diagonal..> > > + Sparse<Tridiagonal<Dense<..> > >
>
>
> Best
>
> Laurent
>
>
>> I think this is an important point. Toon and I have thought about this
>> for a long time. We will certainly talk about this in detail. 
>> Providing
>> the container matrix< matrix<T> > is not hard. I recall that 
>> algorithms
>> are application/interpretation dependent, and we had a lively 
>> discussion
>> on the list on this topic. So, I agree that 'some' support is needed,
>> although it is not yet clear to me what 'some' here means.
>>
>> Best,
>
>
>
> _______________________________________________
> glas mailing list
> glas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/glas
>
------------
Peter Gottschling, Ph.D.
Research Associate
Open Systems Laboratory
Indiana University
135 Lindley Hall
Bloomington, IN 47405
Tel.: [phone]   Fax: [phone]
http://www.osl.iu.edu/~pgottsch