dolfin team mailing list archive

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

2007/10/8, Anders Logg <logg@xxxxxxxxx>:
> Martin Sandve Alnæs wrote:
> > 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.
>
> I'm not sure I understand. Don't you want to add something like
> GenericMatrix::type()?
>
> Then you can use it either to simply throw exceptions when the type is
> wrong or have a switch if you want to handle the other cases.
>
> /Anders

To check if "GenericVector * foo" is a MyVector *, the right way is to
use dynamic_cast<MyVector*>(foo) in C++, and isinstance(foo, MyVector)
in Python.

If we need a switch to handle subtypes, a typeswitch, the code doesn't
scale to new backends. If I add my own backend in my own code separate
from dolfin (like the pycc::CRS matrix will be), it shouldn't be
necessary to touch dolfin. All dolfin code that uses this typeswitch
will fail to handle my backend. This destroys much of the advantages
of inheritance, and should be avoided if possible.

Do you have any particular cases that has to be implemented with a typeswitch?

Quoting the Python manifest (try "import this" in python):
"Special cases aren't special enough to break the rules."

--
Martin