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 free.fr
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
>>glaslists.boost.org
>>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.: +1-812-855-3608 Fax: +1-812-856-0853
>http://www.osl.iu.edu
/~pgottsch
>
>_______________________________________________
>glas mailing list
>glaslists.boost.org
>http
://lists.boost.org/mailman/listinfo.cgi/glas
>
>
Disclaimer: http
://www.kuleuven.be/cwis/email_disclaimer.htm
_______________________________________________
glas mailing list
glaslists.boost.org
http
://lists.boost.org/mailman/listinfo.cgi/glas |