dolfin team mailing list archive

[Anders Logg <logg@xxxxxxxxx>]

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.

> There are some member function of PETScSparseMatrix which are used in
> ode which I intentionally didn't implement for uBlasSparseMatrix. An
> example is PETScSparseMatrix::mult(const real x[], uint row), which
> looks nasty to me. On top of that, the same operation is really compact
> using the uBlas data types. Mutliplying row i of (sparse) matrix A with
> the vector x is just
> 
>   real a = inner_prod(row(A,i), x); 
> 
> Garth

The mult() function is somewhat special and should probably be
removed. It's just there for convenience. The multi-adaptive solvers
need to compute a specific component f_i of the right-hand side and if
the right-hand side is something like f(u) = Au then f_i = A.mult(u, i).

/Anders