dolfin team mailing list archive

["Martin Sandve Alnæs" <martinal@xxxxxxxxx>]

2007/10/5, Anders Logg <logg@xxxxxxxxx>:
> Martin Sandve Alnæs wrote:
>  > Ola and I did some thinking on the SparsityPattern vs Epetra_CRSGraph
>  > issue, ie the case where you want the Assembler to use an external
>  > class for building the sparsity pattern. I think we came up with a
>  > nice solution, in fact it looks much like the GenericTensor::create()
>  > function we added earlier (as well as the ufc::form::create_*
>  > functions).
>  >
>  > Basically, a GenericTensor knows the type of the
>  > GenericSparsityPattern it expects in init(pattern). Thus we can do
>  > this:
>  >
>  > in application:
>  >   A = Matrix()
>  >   assemble(A, ...)
>  >
>  > in assemble:
>  >   pattern = A.createPattern()
>  >   initPattern(pattern, ...)
>
> Is initPattern() a member function of Assembler?
>
> And that function computes the sparsity pattern from the mesh, dofmaps
> etc? Then it looks to me that it may only do so using the public
> interface of GenericSparsityPattern. Is that enough to build an
> Epetra_CRSGraph?


initPattern() is meant to be a member of Assembler, which runs over
the mesh and calls dofmap->tabulate_dofs(cell...) like during assembly
and passes rows and cols to insertBlock below, like the add function
in GenericTensor but without the values. However, I didn't think about
rank 1 and 0 tensor here, only matrices need this mesh iteration. And
I'm not on top of the parallell stuff...


>  >   A.init(pattern)
>  >   ...
>  >
>  > with
>  >
>  > class GenericSparsityPattern
>  > {
>  >   virtual void insertBlock(rows, cols, ...)  = 0;
>  >   ...
>  > }
>
> What does insertBlock do? Is this all we need?


It adds the "element tensor" with given global row and col indices to
the sparsity
pattern. We probably need some more minor functions like setting dimensions.


>  > class GenericTensor
>  > {
>  >   GenericSparsityPattern * createPattern() const = 0;
>  >
>  >   // similar to what we added earlier:
>  >   GenericTensor * create() const = 0;
>  > }
>  >
>  > Then, for Epetra, we will have the classes:
>  >
>  > class EpetraMatrix: public GenericMatrix
>  > class EpetraMatrixPattern: public GenericSparsityPattern
>  > class EpetraVector: public GenericVector
>  > class EpetraVectorPattern: public GenericSparsityPattern
>  >
>  > Other implementations can reuse the current SparsityPattern if they
> want to.
>  >
>  > In my opinion, the only ugly part is that
>  > GenericTensor::init(GenericSparsityPattern & pattern) won't really
>  > accept any subclass, but should have a type check and throw an
>  > exception. But we can't get around that with a virtual init
>  > function.
>
> The same problem shows up if we want to add mult() to the
> GenericMatrix interface (or any other function that involves
> interaction between the two different GenericFoo classes).
>
> So should we add type() as a method to all GenericFoo classes?
>
> If the GenericFoo classes should be extendable by a user (who is not
> allowed to edit GenericFoo), then type() could return a unique uint
> (not enum) for each backend (uBlas = 1, PETSc = 2, Epetra = 3). Then a
> user can add new backends 4, 5, etc.


What should we use that for? It smells of type switches, which means
code that uses this will have to be extended if a new backend is
added. It's better to have the functionality we need in the interface,
so the backend is purely specified by the GenericTensor object.

In PyCC, Ola is currently rewriting CRS to be a direct implementation
of GenericMatrix. That's a typical task that users of other external
libraries might want to do, to use dolfin to assemble into their own
matrix formats without touching dolfin. If everything is contained in
the generic interface, that's possible.


>  > And at last there's still the parallellity troll waiting, probably
>  > requiring some more changes. But I think the distribution pattern can
>  > be merged nicely into the sparsity pattern similar to the Petra
>  > design.
>
> ok.
>
> /Anders

But don't trust me on this...

--
Martin