dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

Martin Sandve Alnæs wrote:
> 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.

Yes, I agree. But it looks to me like your suggestion (for sparsity
pattern) would introduce the same type check.

> 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.

Yes, good. But to avoid the type switches (if we want to avoid them,
see below) then maybe we should keep a single class SparsityPattern
and force every backend to use it. This means we move the
Epetra-specific stuff to EpetraMatrix::init(). The init function
receives a SparsityPattern which should contain everything needed to
build an Epetra-specific sparsity pattern.

This means we need to copy some data, but I don't see how this can be
avoided if we absolutely don't want any type switches.

Avoiding type switches completely also means that we can't add
something like mult() to GenericMatrix. I suspect we will be tempted
to (re)introduce type switches (and a full envelope-letter design) at
some point to have a more complete interface.

/Anders


>
>
>>  > 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...
>