dolfin team mailing list archive

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

On Fri, 2006-06-30 at 17:28 +0200, Anders Logg wrote:
> On Fri, Jun 30, 2006 at 05:01:57PM +0200, Garth N. Wells wrote:
> > 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;
> 
> My problem with templates is that it easily propagates. If you make
> uBlasMatrix a template, then it might force us to template also
> functions that work on a uBlasMatrix, like uBlasKrylovSolver. It
> doesn't propagate to FEM, since uBlasMatrix will be a subclass of
> GenericMatrix (I hope).
> 
> I don't mind templates for data structures like uBlasMatrix if it
> helps reduce the code and we can make it so that it doesn't propagate.
>
> But I think we can find a good solution that does both things: using
> templates to reduce the code and have one uBlas implementation and
> at the same time avoid propagating the templates.
> 
> I suggest creating a new template uBlasMatrix as you suggest, but
> making it a subclass of GenericMatrix (so it works with assembly and
> possibly other algorithms) and also a subclass of uBlasKrylovMatrix.
> 
> uBlasKrylovSolver then accepts as argument a uBlasKrylovMatrix, which
> could then be one of the subclasses DenseMatrix and uBlasSparseMatrix.
> 
> Would that work?
> 

I thought about this propagation problem, and the increased compile time
with templates. Rather than using typedefs, we can have

   class uBlasSparseMatrix : public uBlasMatrix<ublas_sparse_matrix> {}

   class DenseMatrix : public uBlasMatrix<ublas_matrix> {}

This way, uBlasSparseMatrix and DenseMatrix are not templates.

> I also suggest in the same spirit that we rename VirtualMatrix to
> PETScKrylovMatrix.
> 
> Furthermore, to make it really structured, we could have the following
> types:
> 
>     uBlasSparseMatrix    sparse ublas, typedef from uBlasMatrix template
>     uBlasDenseMatrix     dense ublas, typedef from uBlasMatrix template
>     uBlasVector          there is only one type of uBlas vector
> 
>     PETScSparseMatrix    sparse petsc
>    (PETScDenseMatrix     not needed?)
>     PETScVector          there is only one type of PETSc vector
> 
>     typedef PETScSparseMatrix SparseMatrix;
>     typedef PETScVector SparseVector;
> 
>     typedef uBlasDenseMatrix DenseMatrix;
>     typedef uBlasVector DenseVector;
> 
>     typedef SparseMatrix Matrix;
>     typedef SparseVector Vector;
> 
> In addition, we have a suitable ifdef if PETSc is not available to
> replace PETSc with sparse uBlas.

I was thinking along the same lines. It will provide a nice structure
for any linear algebra packages which may come along in the future. What
about "PETScSparseMatrix" becoming "PETScMatrix"? For PETSc matrices,
the only difference (in terms of the interface) between dense and sparse
matrices is the argument passed to the PETSc function MatCreate. This
way, the "real" action will be in

   PETScMatrix
   uBlasMatrix

and other derived matrix classes/typedefs will be specialisations of
these two.

> 
> > > > 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.
> 
> Yes, but I suspect it's not noticeable compared to the other stuff
> that is done in the uBlasKrylovSolver. It would be easy to test.
>
> > > > 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?
> 
> You've received two copies on the last couple of mails because I've
> just pressed reply and the mails were sent to you. Then I noticed I
> didn't post it to the list so I had to do it all over again. Let's
> hope I remember this time and you only receive one reply... :-)
>
> It might be that the mails I've received from you were addressed to me
> and then cc to the list and not just to the list? Or there used to be
> a reply-to field that is now missing.
> 

I'm definitely receiving double of everything.

Garth

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