|
|
Glas :Re: [glas] multi-dimensional arrays? |
From: plagne.laurent_at_[hidden]
Date: 2006-12-22 12:15:53
Hi Peter, Hi Karl, sorry for the delay.
Since my own lib (with a new name LEGOLAS++) handles nested matrices and vectors
I can reply when I had to address the issues you have mentioned. Up o now, I did
not address the difficult and important matrix<->matrix (BLAS3) issue in the
general case so I do not pretend to hold all the answers...
>>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.
All the matrices hold a specific sparsity pattern (Dense, Sparse, Banded,
Tridiag..) which defines a corresponding default matrix-vector algorithm. As an
example, A Dense matrix will rely (by default) on a basic matrix-vector algo
which is modified to take in account the assign mode (Y+=A*X, or Y=A*X,..) :
struct DenseMatrixVectorProduct{
template<class ASSIGN_MODE>
struct Engine{
template <class MATRIX,class VECTOR, class VECTOR_INOUT>
static inline void apply(const MATRIX & A ,
const VECTOR & X ,
VECTOR_INOUT & Y)
{
ASSIGN_MODE::initialize(Y);
const int NR=matrix.nrows();
const int NC=matrix.ncols();
for (int i=0; i< NR ; i++){
for (int j=0; j< NC ; j++){
ASSIGN_MODE::accumulate(Y[i],A(i,j)*X[j]);
}
}
}
return ;
}
};
In the case of Y=A*X :
ASSIGN_MODE::initialize(Y)<=> Y=0.0
ASSIGN_MODE::accumulate(Y[i],A(i,j)*X[j])<=>Y[i]+=A(i,j)*X[j]
In the case of Y+=A*X :
ASSIGN_MODE::initialize(Y)<=> DoNothing
ASSIGN_MODE::accumulate(Y[i],A(i,j)*X[j])<=>Y[i]+=A(i,j)*X[j]
In the case of Y-=A*X :
ASSIGN_MODE::initialize(Y)<=> DoNothing
ASSIGN_MODE::accumulate(Y[i],A(i,j)*X[j])<=>Y[i]-=A(i,j)*X[j]
If A(i,j) is a submatrix, it will invoke automatically a sub-algorithm
linked with its sparsity pattern.
In addition, the user can provide its own algorithm at every level. For example,
if the matrix is Dense AND contains double, float or complex elements AND is
stored AND is stored via compatible 1D (RowMajor or ColMajor) vector : then a
xgemm algo can be used. Up to now, I did not use Traits to link the better
possible algorithm available corresponding to a given matrix type. The
important issue was for me that the user deal only (in a first approach) with
the sparsity patterns and still be able to tune the code (via typedef) at a
later time. My main concern was that this choice remains always possible and
that is why the user can also provides its own container usually prescribed by
a third part lib (e.g. ATLAS).
Since the vector case is much simpler in my case (only dense vector) a Traits
class
is used to automatically invoke the best lib available for a given vectorial
operation
(X+=aY, dot, X=0,..) in the 1D case.
For the transpose issue, I have tried different strategies but I doubt that I
found a perfect one :
every MVP algorithms come with a pair Transpose class method which is called
recursively
via a transpose function. For example :
Y=transpose(A)*X lead to use
struct DenseMatrixVectorProduct{
....
struct Transpose{
template<class ASSIGN_MODE>
struct Engine{
template <class MATRIX,class VECTOR, class VECTOR_INOUT>
static inline void apply(const MATRIX & A ,
const VECTOR & X ,
VECTOR_INOUT & Y)
{
ASSIGN_MODE::initialize(Y);
const int NR=matrix.nrows();
const int NC=matrix.ncols();
for (int i=0; i< NR ; i++){
for (int j=0; j< NC ; j++){
ASSIGN_MODE::accumulate(Y[j],transpose(A(i,j))*X[i]);
}
}
}
return ;
}
};
I choose this pretty heavy method because it lead to access the matrix elements
in the same order than the direct product. On can think that the traversal
order is closely linked with the matrix storage which is not modified.
>Next question, do we need special concepts? Say e.g.
>RecursiveTransposable for a recursive function transpose
To me, every matrix can be transposed... So no additional concept
>would recursively descent and break if it reaches a scalar value that
>has no value_type.
In LEGOLAS++, every matrix elements type (especially scalars) comes with a
companion ElementDataDriver class which is used to ensure that every elements
will model the ElementMatrix Concept. This concept includes the definition of
transpose function, dot product, zero,...
Again Peter, I think that the matrix-matrix issue addressed by GLAS is much more
difficult.
>> Extending this to block-matrices,
Strictly speaking I do not deal with block-matrices but only with Matrix with a
ElementMatrix concept which is modeled by matrices.
The same thing for vectors.
Hi Karl :
I do not deal with blocked algorithm used in optimized xgemm implementation for
Dense matrices. While this issue can be seen as an implementation policy, the
sparsity pattern of a matrix (chosen be the user) is NOT an implementation
detail which should be abstracted away. On the contrary I try to design my lib
to allow a very detailed description of this sparsity
pattern.
To me, a matrix and a 2D vectors are very different types. Matrix are mainly
linear operator acting on Vectors.
The container of a dense matrix can eventually be a 2D array (vector) but the
semantic level of a matrix and the corresponding container should be strictly
separated since a matrix does not imply a container. Permutation matrices,
finite different stencils etc.. are matrix and usually are not stored at all...
Best,
Laurent
Quoting Peter Gottschling <pgottsch_at_[hidden]>:
> 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.: +1-812-855-3608 Fax: +1-812-856-0853
> http://www.osl.iu.edu/~pgottsch
>
> _______________________________________________
> glas mailing list
> glas_at_[hidden]
> http://lists.boost.org/mailman/listinfo.cgi/glas
>