dolfin team mailing list archive

["Garth N. Wells" <g.n.wells@xxxxxxxxxx>]

On Sun, 2006-05-07 at 17:18 +0200, Anders Logg wrote:
> On Sun, May 07, 2006 at 03:59:09PM +0200, Garth N. Wells wrote:
> > On Sun, 2006-05-07 at 14:10 +0200, Anders Logg wrote:
> > > When adding the option to compile without PETSc, I took the
> > > opportunity change the precomputation of quadrature weights for the
> > > ODE solvers to use DenseMatrix instead of the standard PETSc
> > > Matrix. (The matrices are small and dense anyway.)
> > > 
> > > It seems that DenseMatrix::invert() sometimes gives unexpected results
> > > and the reason seems to be a memory access error somewhere in ublas. I
> > > don't know if this is a problem with ublas or the way we use
> > > it. Garth, could you look again at DenseMatrix::invert() and see if
> > > you can spot something obvious?
> > > 
> > > I have attached the error log from running valgrind and here is the code that
> > > reproduces the error:
> > > 
> > >    DenseMatrix A(3, 3);
> > >    A(0, 0) = A(1, 1) = A(2, 2) = 1.0;
> > >    A.invert();
> > > 
> > 
> > When you create a matrix, it is not initialised to zero by Boost. You
> > can do
> > 
> >     A.clear(); // the Boost function
> > 
> > or
> 
> Thanks| Works perfectly now. I added clear() to DenseMatrix and
> DenseVector because I don't see a point in not initializing a dense
> matrix. It's probably a very small additional overhead and so much
> more user friendly.
> 

Have you done any benchmarks? If you configure with
--enable-optimization, you'll see about a factor 10 speed-up with
DenseMatrix and DenseVector.

> >     A = 0; // DOLFIN wrapper for A.clear() (provides compatibility with
> > Matrix, but I prefer A.clear() because A = x only works for x = zero.)
> 
> Assignment to a scalar should work for any value of the scalar for a
> dense matrix. It's only for the sparse where we need the special case.
> (Then we don't want to initialize all elements as nonzero.)
> 

I agree. I would like to remove the function "operator= real 0" from
Matrix and use something else to imply initialisation. clear() is good,
but I recall that it's used for something else in Matrix. I want to keep
the meaning of functions that are common to DenseMatrix and DenseVector
identical. 

> > I left it up to the user to initialise the matrix because often all
> > values will be set by the user, so initialisation is just extra overhead
> > from some problems.
> 
> I would argue that we put in the clear() automatically. There's a
> trade-off of course, but if you don't make any decisions at all for
> the user you end up with PETSc and a lot of overhead in writing the
> code.

Is that a veiled comment on PETSc? :-). On this point (protecting the
user), what if we remove or make protected the functions that provide
access to individual elements of Matrix and Vector. The problem when
setting values is that 

  MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
  MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

are called every time and element is set which is really slow. PETSc is
not designed for element-by-element access.

Another point - now that we have a Boost dependency, how about using
more of its features? It has some nice smart pointer classes so we could
get rid of a lot of "new" and "delete", would be useful for
sub-functions that share a vector. Also, uBlas has classes "c_vector"
and "c_matrix" in which the data is stored in a regular C++ array and
provides vector functionality (the maximum size must be known at compile
time). It has very low overhead compared to a regular array. This could
be useful in BilinearForm to have a nice-looking matrix-like
representation and, as Johan J suggested before, could be used for the
affine maps. 


Garth

> 
> /Anders
> 
> 
> > Garth
> > 
> > > How about adding solve(DenseVector& x, const DenseVector& b) to DenseMatrix?
> > > 
> > > /Anders
> > > _______________________________________________
> > > DOLFIN-dev mailing list
> > > DOLFIN-dev@xxxxxxxxxx
> > > http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
> 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/cgi-bin/mailman/listinfo/dolfin-dev
-- 
Dr. Garth N. Wells
Faculty of Civil Engineering and Geosciences
Delft University of Technology
Stevinweg 1
2628 CN Delft
Netherlands

tel.     +31 15 278 7922
fax.     +31 15 278 6383
e-mail   g.n.wells@xxxxxxxxxx
url      http://www.mechanics.citg.tudelft.nl/~garth