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/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 you do:
  Foo* foo = dynamic_cast<Foo*>(bar);
foo will be NULL if bar wasn't a Foo* after all, or the correctly
typed pointer if it was.


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


Ok.


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


I see, but do we really need or want to be able to do that? I'd think
that people usually will only want to use one backend at a time, but I
could be wrong.


> But in summary, I think your suggestion looks good.

Great.


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


I'll ask Ola to do it since he's working on implementing pycc::CRS as
a subclass of dolfin::GenericMatrix right now.

--
Martin