dolfin team mailing list archive

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

On Mon, 2006-07-03 at 12:20 +0200, Anders Logg wrote:
> On Fri, Jun 30, 2006 at 05:51:36PM +0200, Garth N. Wells wrote:
> > 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.
> 
> Sounds good.
> 
> But with DenseMatrix replaced with uBlasDenseMatrix?
> 

Yes. This will be useful because some libraries (like UMFPACK) require
compressed column storage, so we can then have something like 

  class uBlasSparseColMatrix : 
                public uBlasMatrix<ublas_sparse_col_major_matrix> {}



> > > 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.
> 
> Yes, sounds good. Should there also be PETScSparseMatrix and
> PETScDenseMatrix so the structure is the same?
> 
> The difference would be in the implementation of the two different
> matrix types. One is handled by a base class template, the other one
> with a suitable argument to the constructor.
> 

This can be done, but do we want to encourage the use of PETSc for dense
matrices?

Garth

> > > 
> > > > > > 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.
> 
> Very strange. I'll double-check when I send this email that it is only
> addressed to the list.
> 
> I'll see if I can setup the virtual matrix stuff for uBlas matrices
> this afternoon and then go ahead and make the ODE solvers PETSc
> independent.
> 
> /Anders
> 
> 
> 
> > 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
> > 
> _______________________________________________
> DOLFIN-dev mailing list
> DOLFIN-dev@xxxxxxxxxx
> http://www.fenics.org/mailman/listinfo/dolfin-dev