dolfin team mailing list archive

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

On Fri, 2006-06-30 at 16:00 +0200, Anders Logg wrote:
> On Fri, Jun 30, 2006 at 03:08:45PM +0200, Garth N. Wells wrote:
> > On Fri, 2006-06-30 at 12:29 +0200, Anders Logg wrote:
> > > On Fri, Jun 30, 2006 at 08:48:59AM +0200, Anders Logg wrote:
> > > > On Fri, Jun 30, 2006 at 12:32:25AM +0200, Garth N. Wells wrote:
> > > > > On Thu, 2006-06-29 at 20:20 +0200, Anders Logg wrote: 
> > > > > > I've started trying to remove the PETSc ifdefs from the ODE solvers,
> > > > > > but ran into some problems. The ODE solver uses the calls
> > > > > > 
> > > > > >     Vector x;
> > > > > >     real* xx = x.array();
> > > > > >     // operate on xx
> > > > > >     x.restore(xx);
> > > > > > 
> > > > > > This is supported only by PETSc vectors.
> > > > > > 
> > > > > > There are a couple of different solutions:
> > > > > > 
> > > > > > 1. Don't use Vector::array().
> > > > > > 
> > > > > > 2. Add Vector::array() to uBlas vectors.
> > > > > > 
> > > > > > 3. Do 2 and add it to GenericVector as well.
> > > > > > 
> > > > > > Any thoughts? I have a feeling 1 is the best solution, but it's also
> > > > > > the one that requires the most work (for me to rewrite the ODE
> > > > > > solvers).
> > > > > > 
> > > > > 
> > > > > I think that option 1 is best. I don't like options 2 & 3 as they take
> > > > > away the advantages of using uBlas. What exactly do you do with the
> > > > > Vector::array() in ode?
> > > > 
> > > > The solution is (for some of the solvers) stored as a Vector u, but when
> > > > calling the user-supplied right-hand side f(u, t), the argument is
> > > > given as a plain C array since the user will often want to do a lot of
> > > > indexing of the vector in the definition of the right-hand side and I
> > > > suspect indexing is slow for PETSc vectors.
> > > > 
> > > > > How fast are the PETSc functions VecGetValues
> > > > > and VecSetValue, and how much time is involved in calling
> > > > > VecAssemblyBegin() and VecAssemblyEnd() (which are currently called
> > > > > after every set)? I use Vector::array in the Cahn-Hilliard solver to do
> > > > > things like take the log of each term. Not sure how to do this
> > > > > efficiently without using array which will work for both PETSc and uBlas
> > > > > types. 
> > > > 
> > > > I'll make a small test to see how slow it really is.
> > > 
> > > Here are some results. The timings are for setting and then accessing
> > > all elements in a vector of size 10000:
> > > 
> > > Setting values in PETSc vector:   0.0708
> > > Accessing values in PETSc vector: 0.00031
> > > Setting values in uBlas vector:   1.42e-05
> > > Accessing values in uBlas vector: 4.6e-06
> > > Setting values in STL vector:     1.41e-05
> > > Accessing values in STL vector:   4.7e-06
> > > Setting values in C array:        1.4e-05
> > > Accessing values in C array:      1.41e-05
> > > 
> > > Some conclusions:
> > > 
> > >  - Setting and accessing values element-wise in a PETSc array
> > >    is terribly slow
> > >
> > >  - uBlas, STL vector and C array are equally fast, don't know
> > >    why access for C array is slower
> > > 
> > 
> > When I was first testing uBlas I saw that C arrays were often a bit
> > slower. I didn't understand why either.
> > 
> > 
> > > It seems the natural is to modify the interface for the ODE solver from
> > > 
> > >     void f(const real u[], real t, real y[]);
> > > 
> > > to
> > > 
> > >     void f(const DenseVector& u, real t, DenseVector& y);
> > > 
> > > Since the ODE solvers are currently not parallel anyway, I might as
> > > well use DenseVector in combination with the uBlas solvers. This
> > > should speed up the ODE solvers for small systems as well.
> > > 
> > 
> > For the future, it might be possible to create a uBlas vector that
> > addresses the array returned by PETScVector::array(). It's important
> > that we find a solution for manipulating array elements efficiently and
> > cleanly.
> 
> DenseVector seems like a good option. It's the fastest and it works
> with some of the solvers.
> 
> Maybe the array() should be special for the PETSc vector since that's
> the only place where it's needed?
> 
> > > Garth, we need to extend the new uBlas solvers you added so that they
> > > work with a VirtualMatrix, that is, something which provides the
> > > function
> > > 
> > >     void mult(const DenseVector& x, DenseVector& y) const;
> > > 
> > > Would that be much work to add?
> > > 
> > 
> > Should be OK.
> 
> ok, I might try to fix it.
> 

I've had a quick look, and the only things that needs to be done are

1. Template the uBlasKrylovSolver solve functions,

  template <class M >  
  uint solve(const M& A, DenseVector& x, const DenseVector& b);

I wanted to do this anyway as the solvers can also work with dense
matrices. 

2. Add a function uBlasSparseMatrix::mult_plus(const x, y) to compute y
+= Ax

3. Use A.mult(x, y) in place of axpy_prod(A, x, y, true) and
A.mult_plus(x,y) in place of axpy_prod(A, x, y, false) in
uBlasKrylovSolver.

4. Make the uBlasKrylovSolver::mult() functions templates so that they
can accept ublas_vector. The alternative is to use Densevector in place
of ublas_vector in uBlasKrylovSolver. I didn't do this originally as the
DenseVector constructor calls clear() which is not needed and does
involve some overhead. 
 
4. Add the class uBlasVirtualMatrix.

5. Generalise the preconditioner. 

Garth

> /Anders