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

From: Karl Meerbergen (Karl.Meerbergen_at_[hidden])
Date: 2006-12-18 10:19:23

Hi,

If I understand well, the reason for having nested matrices is to 
organise blocking in matrices and vectors for exploiting memory hierarchy.

1) Should the user see this structure? Does the user access the data via 
the nesting or can he address directly? For example, in BLAS and LAPACK 
implementations, this is fully transparant.

2) Former discussions on the mailing list have shown different 
interpretations of nested vectors and matrices.

For example, for some, vector< vector<T> > was interpreted as a vector 
organized in blocks, for others it was a matrix. Therefore, I would 
prefer dedicated names such as blocked_vector<T> and matrix_as_vector<T> 
or something like that.

The internals in the matvec prod and matmat prod should use some form of 
blocking, but should this be visible inside the vector and matrix classes?

Karl

Peter Gottschling wrote:

>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
>
>_______________________________________________
>glas mailing list
>glas_at_[hidden]
>http://lists.boost.org/mailman/listinfo.cgi/glas
>  
>

Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm