dolfin team mailing list archive

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

On Fri, 2006-06-30 at 16:41 +0200, Anders Logg wrote:
> On Fri, Jun 30, 2006 at 04:30:50PM +0200, Garth N. Wells wrote:
> > 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. 
> 
> I don't think we need to make templates (as usual...). I started
> looking at it and I've created an interface base class that specifies
> the interface for matrices that work with the uBlas Krylov solver: a
> size() function and a mult() function.
> 

:)

No trying to frustrate you, but I was planning to make uBlasMatrix a
templated class as the code in uBlasSparseMatrix and Densematrix is 95%
identical. I'd then add

  typedef uBlasMatrix<ublas_matrix> DenseMatrix; 
  typedef uBlasMatrix<ublas_sparse_matrix> uBlasSparseMatrix;
 

> > 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. 
> 
> Is it noticeable?
> 

It's equivalent to assigning a value to all vector elements.

> > 4. Add the class uBlasVirtualMatrix.
> 
> I called it uBlasKrylovMatrix and will rename VirtualMatrix to
> PETScKrylovMatrix. Then we can have a typedef for VirtualMatrix
> depending on the default matrix type.
> 

OK.

On another topic, I don't think that we should implement a uBlas sparse
direct solver (nor do I want to). It's way too complicated to make it
fast. There are however bindings for uBlas to UMFPACK which would be a
really simple way to have a good direct solver. UMFPACK is available as
a Debian package so it's easy to install and it's one the fastest LU
solvers around. What do you think?

BTW, can we fix the mailing list so I don't get two copies of
everything?

Garth

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