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
