dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

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

/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