dolfin team mailing list archive

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

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

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




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