[GRASS-dev] BLAS/LAPACK (Part II)

Mon Aug 20 17:37:39 EDT 2007

Hi Daniel,

2007/8/20, Daniel Bundala <bundala at gmail.com>:
> Guys,
>
> It is quite interesting, but I have had plans to replace v.generalize
> matrix code by "yours" library code. I have not studied G_matrix_*
> code carefully, but it seems to me that it is superior.

Unfortunately there are two libraries which handle
the solution of linear equation systems.
The gmath library with the G_math_* and G_vector_* functions
and the higher level gpde library with several multi threaded solvers
(N_solver_cg ...).
I have implemented the matrix vector functionality in the gpde library
again, because i wanted them multi threaded and i have had no idea
if the gmath functions are thread safe and easy to parallele.

>
> Firstly, Soeren wrote that the current code is multithreaded.

Yes, the code from the gpde library is multi threaded, but you can
link parallelled lapack and blas libraries to the gmath interface (scalapack).

> Secondly, someone mentioned, that it supports the sparse matrices.

The gpde library supports a simple sparse matrix implementation
and the matrix vector product functions.

> Support of sparse matrices would increase the efficiency of
> v.generalize since it uses only the sparse matrices.
> Thirdly, Soeren mentioned that the current code supports many methods
> my code doesnt support. To tell the truth, I have never heard about
> many of them (Well, I am still (young) student...)

Of what kind are your matrices? There are two very efficient
solver within the gpde lib:
1.) for sparse symmetric and positive definite matrices (conjugate gradients
with preconditioning cg/pcg)
http://en.wikipedia.org/wiki/Conjugate_gradient_method
and
2.)for sparse non symmetric, non definite matrices (stabilized bi
conjugate gradients)
http://en.wikipedia.org/wiki/Biconjugate_gradient_method

Those two methods are the most efficient linear equation solvers
available for large matrices.

>
> The only thing I am missing in the current code is the direct access
> to the elements of a matrix. But, this is quite dangerous and I really
> doubt whether this is a good API-desing.

The matrix implementation within the gpde library offers direct access
to the matrix entries and supports row shuffling by setting the
row pointer (important for for pivoting). So the programmer
have to assure thread safe access by himself.

>
> On the other hand, it is true that the current code is quite obscure,
> say. Also, it is tempting to replace fortran code by C code.
> Therefore, my suggestons are: clean library code and replace the
> current code by v.generalize code only if it is really faster. Some
> benchmarks are probably required, but I doubt that my code beats
> (optimized) library code.

The gpde library implementation is IMHO not faster than
the gmath and BLAS/LAPACK stuff. Eg: the gmath lu solver
is 30% faster than the gpde lu solver with pivoting.
But the gmath lu solver is code from the numerical recipes,
we have to rewrite this method. And the gpde lu solver runs on
multi processor machines.

I will present an implementation of several matrix-vector functions
in some days. And I'm open for any suggestions about API design. :)

Best regards
Soeren

>
> Daniel
>
> On 8/20/07, Markus Neteler <neteler at itc.it> wrote:
> > Hi Soeren,
> >
> > from my users point of view this sounds excellent. please go ahead...
> > You had already suggested it and there were apparently no
> > objections.
> >
> > thanks
> > Markus
> >
> > On Mon, Aug 20, 2007 at 05:20:18PM +0200, "Sören Gebbert" wrote:
> > > HI folks,
> > >
> > > -------- Original-Nachricht --------
> > > > Datum: Mon, 20 Aug 2007 12:29:14 +0300
> > > > Von: Wolf Bergenheim <wolf+grass at bergenheim.net>
> > > > An: GRASS Devel <grass5 at grass.itc.it>
> > > > CC: Daniel Bundala <daniel.bundala at oriel.ox.ac.uk>, Brad Douglas <rez at touchofmadness.com>
> > > > Betreff: Re: [GRASS-dev] BLAS/LAPACK (Part II)
> > >
> > > > On 17.08.2007 07:09, Brad Douglas wrote:
> > > > >
> > > > > What I propose is moving the matrix code from v.generalize
> > > >
> > > > +1
> > > >
> > > > > (in particular, matrix_inverse() ) to lib/gmath and simplifying the
> > > > existing
> > > > > MATRIX structure.
> > > > >
> > >
> > > I can easily integrate the matrix code from v.generailze into
> > > the gpde library, because the existing matrix structures are quite
> > > similar. Quadratic and sparse matrices are supported.
> > > The gpde library ships several vector-matrix and vector-vector
> > > functions with it, but currently as static functions within the krylov-space solvers. I can make them public (extern),
> > > so they can be accessed from out side of the krylov solvers.
> > >
> > > Many linear equation solvers are available
> > > within the gpde library:
> > > * direct solvers
> > > ** gauss elimination
> > > ** lu decomposition
> > > ** cholesky decomposition.
> > > * iterative solvers
> > > ** gauss seidel / SOR
> > > ** jacobi
> > > ** conjugate gradients (krylov space method)
> > > ** preconditioned conjugate gradients (krylov space method)
> > > ** biconjugate gradients stabilized (krylov space method)
> > >
> > > Everything is multithreaded with OpenMP (solver, matrix, vector operations and some array functions).
> > >
> > > And as you know, the lu code in gmath lib is a copy
> > > of the numerical recipes algorithm and not free.
> > >
> > > I would like to hear some suggestions.
> > >
> > > Best regards
> > > Soeren
> > >
> > >
> > > >
> > > > I think that would be a good idea, especially if you also want to use
> > > > that code. It is easier to maintain the code in one place.
> > > >
> > > > Brad do you know of any additional mathematics or similar things you'd
> > > > like to see in lib/gmath? Perhaps next year it could be a Summer of Code
> > > > project to add them ;)
> > > >
> > > > --Wolf
> > > >
> > > > --
> > > >
> > > > <:3 )---- Wolf Bergenheim ----( 8:>
> > > >
> > > > _______________________________________________
> > > > grass-dev mailing list
> > > > grass-dev at grass.itc.it
> > > > http://grass.itc.it/mailman/listinfo/grass-dev
> > >
> > > --
> > > Psssst! Schon vom neuen GMX MultiMessenger gehört?
> > > Der kanns mit allen: http://www.gmx.net/de/go/multimessenger
> > >
> > > _______________________________________________
> > > grass-dev mailing list
> > > grass-dev at grass.itc.it
> > > http://grass.itc.it/mailman/listinfo/grass-dev
> >
> > --
> > Markus Neteler  <neteler itc it>  http://mpa.itc.it/markus/
> > FBK-irst -  Centro per la Ricerca Scientifica e Tecnologica
> > MPBA - Predictive Models for Biol. & Environ. Data Analysis
> > Via Sommarive, 18        -       38050 Povo (Trento), Italy
> >
> > _______________________________________________
> > grass-dev mailing list
> > grass-dev at grass.itc.it
> > http://grass.itc.it/mailman/listinfo/grass-dev
> >
>
> _______________________________________________
> grass-dev mailing list
> grass-dev at grass.itc.it
> http://grass.itc.it/mailman/listinfo/grass-dev
>