dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sun, Oct 24, 2004 at 12:37:54PM +0300, dkouroun@xxxxxxxxx wrote:

> Dear developers, Currently I am writing expression-templates code
> which I am going to publish , that allow for evaluation of
> variational forms. You do not need any any external language for
> that. Expression templates result to code that runs as fast as
> hand-optimizing code. In fact I know 2 other ways to do the same
> thing.  One is only applicable to tetrahedra(3D)-triangles(2D) and
> the other is more general. I expect them though to be slower than
> expression templates.  How much slower? I will find out in the next
> few days! The good thing is
> that either way, one is finally able to write a variational formulation like:
> 
> Variable<Simplex, P2, _3D> u,v;
> f = u*v + k*grad(u)*grad(v); //where k is constant.

Looks interesting. Yes, templates is one possible approach. Generating
code is another one. We will see which approach turns out to be the
fastest and it will be good to have some competition.

> In order for it to be more general though, one should be able to use,
> elements of arbitrary order. Then I have a problem with integration in
> triangle-tetrahedra, because the integration points-weights that are available,
> as far as I know, are 7th order for triangles and 5th order for
> tetrahedra. Do you know any generic routine which generates points-weights
> for quadrature in triangles-tetrahedra for arbitrary order? 

Yes, FIAT does this. It will hopefully be available on the FEniCS web
site soon. Meanwhile, take a look at

    http://people.cs.uchicago.edu/~kirby/index.cgi?page=fiat

> As it concerns the Matrix-Vector Library it depends on the goal of the whole
> project. If you aim at parallel stuff, then you also have to decide which kind 
> of parallel stuff you want to support. Then one has the following 2 options:
> 1:) For SMPs the best way, at my opinion,to go through is to create your own
> expression-template library with support for OpenMP directives. 
> 2:) For distributed parallel stuff, since dolfin's is a FEM oriented software
> and FEM use sparse  kernels, the is no point developing or wrapping existing 
> (dense)Matrix-Vector parallel(MPI) libraries. On the other hand distributed
> parallel FEM codes should implement Domain-Decomposition rather than parallel 
> sparse linear solvers.

We're trying to come to a decision concerning the linear algebra. Any
input is appreciated. There will be a separate thread about this.
 
> Anders let me ask you something! When solving non-linear problems,
> or when using an iterative solver(GMRES, BiCGSTAB ...) we are
> computing a lot of times the residual (A*x-b) where A is a sparse
> matrix. The matrix vector multiplication routine is a very slow one
> and its performance is also bounded by bandwidth. Do you have any
> idea of an alternative way? Some people say, that first a
> renumbering should be performed to minimize the sparse-matrix
> bandwidth and then use optimized(Goto, Atlas BLAS) kernels to
> perform the multiplication blocks by blocks! This can be applied of
> course if your matrix has dense sub-blocks. Others convert the CSR
> format to BSR(block compressed sparse row) format. What do you think
> about that? Do you have something in mind?

It will all depend on the decision we make on the linear algebra, but
I don't think the current implementation of matrix-vector-multiply is
slow. I think we compared it to a number of library (blas) routines
and the comparison was favorable.

/Anders