dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

On Sat, Mar 10, 2007 at 04:23:55PM +0100, Martin Sandve Alnæs wrote:
> 2007/3/9, Garth N. Wells <g.n.wells@xxxxxxxxxx>:
> > Anders Logg wrote:
> > > On Fri, Mar 09, 2007 at 07:10:31PM +0100, Garth N. Wells wrote:
> > >> Anders Logg wrote:
> > >>> On Fri, Mar 09, 2007 at 03:07:43PM +0100, Garth N. Wells wrote:
> > >>>> Why do we need a class AssemblyMatrix? GenericMatrix/GenericVector
> > >>>> should provide the necessary  virtual functions for assembly.
> > >>>>
> > >>>> Garth
> > >>> Yes, definitely. The assembler assumes it gets a GenericTensor that
> > >>> has a function
> > >>>
> > >>>     void add(real* block, uint* size, uint** entries)
> > >>>
> > >>> (entries is probably a bad choice of name here, let's change it)
> > >>> This makes it possible to implement just one assembly function and
> > >>> keep the Assembler class small.
> > >>>
> > >> Then GenericVector and GenericMatrix should then derive from GenericTensor?
> > >
> > > Yes. So far it is just two functions that are needed, init() and add().
> > >
> > >>> Anyway, the class AssemblyMatrix is just a temporary implementation of
> > >>> something that implements the GenericTensor interface. My suggestion
> > >>> is that we overload the add() function above in both GenericMatrix and
> > >>> GenericVector. The overloaded functions will map to the correct add()
> > >>> function in the respective classes. For GenericMatrix, we'll have
> > >>>
> > >>>     void add(real* block, uint* size, uint** entries)
> > >>>     {
> > >>>         add(block, size[0], entries[0], size[1], entries[1]);
> > >>>     }
> > >>>
> > >>> So, even if AssemblyMatrix could be of some use (it seems to be a fast
> > >>> option for assembly, the STL vector of maps beat PETSc in at least one
> > >>> benchmark I made)
> > >> Don't draw conclusions too quickly. STL vector of maps is fast to
> > >> assemble, but is close to useless for linear algebra. The current
> > >> implementation for the assembly uBLAS matrices use something like STL
> > >> vector of maps  (fast to assemble but slow for linear algebra) and
> > >> converts it to compressed row storage (fast for linear algebra). The
> > >> conversion can take some time. A STL vector of integer maps should be
> > >> created and then used to initialise a particular sparse storage format
> > >> (before first use only). Then we wouldn't have to convert matrix formats.
> > >
> > > Do you mean this should be done in DofMap?
> > >
> >
> > Yes. The STL vector of maps is cheap to construct, and can be used to
> > initialise the non-zero layout of a sparse matrix.
> >
> 
> You should use vector<set<int>>, not vector<map<int>>, if the goal is
> only to build the sparsity pattern, as this gives lower memory
> overhead since the set is a map with no value.

Yes, we use vector<set<int>> for the sparsity pattern, but in the
class AssemblyMatrix (which is something I use temporarily until we
make GenericMatrix and GenericTensor subclasses of GenericTensor) I
have used vector<map<int>> to represent the sparse matrix.

> Anders: I have an implementation of this in
> pycc/MatSparse/CRSGraph.cpp and pycc/FEMAssembler/CRSAssembler.cpp,
> feel free to use any suitable parts of it.

We (Garth) has an implementation of this in the class OldDofMap in
DOLFIN that should be merged with the "new" class DofMap for UFC-based
assembly. It looks very similar to what you have in
CRSAssembler/CRSGraph:

for (CellIterator cell(*mesh); !cell.end(); ++cell)
{
  element[0]->nodemap(row_dof, *cell, *mesh);
  for (uint i = 0; i < element[0]->spacedim(); ++i)
  {
    element[1]->nodemap(column_dof, *cell, *mesh);
    for (uint j = 0; j < element[1]->spacedim(); ++j)
      matrix_sparsity_pattern[ row_dof[i] ].insert( column_dof[j] );
  }
}

Here, matrix_sparsity_pattern is a vector<set<int>> and the nodemap
function is the old (non-UFC) name of the tabulate_dofs() function.

Is this the same as your implementation? (Other than splitting the
insertion into a separate call insert_block()?)

> In my brief tests, using vector<set<int>> and using it for allocating
> the pycc::CRS matrix is about as fast as Epetra_CRSGraph/Epetra_CRS.

ok. It will be interesting to compare the two codes in a set of
benchmarks.

/Anders


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