← Back to team overview

dolfin team mailing list archive

Re: About SparsityPattern

 

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

ok, good. This is what I was looking for. I didn't realize that we have something that corresponds to isinstance() in C++.

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.

I agree. But I thought your suggestion involved adding a type() function. I think you can have type() and type switches and still scale (by having each new backend "register" the new backend) but the longer we can avoid it the better.

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

There are quite a few (all involving operating on objects with different backends), like

  x += y;

where x is a PETSc vector and y is a uBlas vector.

But in summary, I think your suggestion looks good.

Just set up something like ~/martinal/hg/dolfin on fenics.org and we can pull your changes (preferrably in small portions).

/Anders


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

--
Martin
_______________________________________________
DOLFIN-dev mailing list
DOLFIN-dev@xxxxxxxxxx
http://www.fenics.org/mailman/listinfo/dolfin-dev


Follow ups

References