dolfin team mailing list archive

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

2007/10/6, 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

There's a huge difference between a _typeswitch_ in code that _uses_
the GenericTensor interface and a _typecheck_ in the _implementation_
of the interface. The typeswitch leads to dependencies all over the
code, the internal typecheck doesn't.

My suggested design has typechecks only within classes in a particular
backend, they don't need to know anything else than their own types.
That is, a GenericMatrix implementation only needs to know the type of
its matching GenericSparsityPattern implementation, and can use
dynamic_cast for the checks.

The same design could be applied to a mult function, the GenericMatrix
implementation only needs to know its matching GenericVector
implementation. To take this idea a bit further, GenericMatrix could
f.ex. have two functions "virtual GenericVector * create_row_vector()
const" and "create_col_vector" to construct vectors in its row and
column spaces. Only the implementation of MyMatrix and MyVector needs
to know the types.

This is pretty much like the design we came up with in ufc::form,
factory functions to construct "families" of objects that belong
together, does this design pattern have a name?

The main this doesn't solve, is how to create a vector when you don't
already have one. Some sort of singleton backend factory could solve
this though.

GenericVector * create_vector()
{
    return backend_factory()->create_vector();
}



select_backend("PetSC");
v = create_vector();
A = create_matrix();


or:

GenericVector * Vector()
{
    return backend_factory()->create_vector();
}

GenericVector * v = Vector();

--
Martin