[GRASS-dev] gmath and gpde update

Soeren Gebbert soerengebbert at googlemail.com
Sat Jan 10 10:36:14 EST 2009


Dear developers,
after a loooong time of development i would like to commit some changes to
the gmath and gdpe lib
and some dependent raster modules.

The following changes have been made:

* The gpde solver (gauss, jacobi, LU, Cholesky, CG, PCG and BiCGStab) are
now located in the gmath library.
  Because of that several files in the gpde directory have been removed or
updated.
* Tests and benchmarks for single and multi threaded computation are
available.
* Blas level 1, 2 and 3 functions are implemented with tests, benchmarks and
basic support for sparse matrices.
  The implemented BLAS functions are not optimized for high performance, but
multi threaded with OpenMP.
* Some new solver uses the BLAS 1 and 2 functions of the grass
implementation.
* An ATLAS wrapper for BLAS level 1 functions is available, BUT level 2 and
3 function are still work in progress and not finished yet.
* r.gwflow and r3.gwflow are updated to use the gmath library.
* Updates to gmath.h and gisdefs.h

The patch file for the latest grass7 trunk is located here:
http://www-pool.math.tu-berlin.de/~soeren/grass/files/soeren_gmath_gpde.diff.gz

A patched grass7 trunk is available here:
http://www-pool.math.tu-berlin.de/~soeren/grass/files/grass_trunk_soeren_gpde_gmath.tar.bz2

To apply the patch unzip the diff file, change directory to the grass_trunk
root directory and type:
cat soeren_gmath_gpde.diff | patch -p0

Why a grass BLAS implementation?
There was a discussion more than a year ago, to replace the dependency in
the gmath library  to the FORTRAN BLAS and LAPACK libraries with
the ATLAS BLAS/LAPACK implementation.
Because of that i have implemented a simple grass BLAS level 1, 2 and 3
library and a wrapper for ATLAS BLAS level 1 functions,
which uses the grass BLAS implementation as backup in case the ATLAS library
is not available. Wrapper for ATLAS BLAS level
2 and 3 functions must be developed to gain any benefit from the ATLAS
library.

Replacement of existing lu decomposition?
Most of the new solvers for linear equation systems (Ax = b) are paralleled
with OpenMP. This works
fine with the krylov-subspace methods (cg, pcg and bicgstab) which also work
with sparse matrices,
but not with the gauss and lu decomposition solver. The speedup of booth
solver is quite bad.
I was unable to replace the existing lu decomposition, which is lend from
numerical recipes,
with a high performance, parallel algorithm which works also with sparse
matrices.
This is much to complicated for me.
Because the new lu and gauss solver are simply paralleled with OpenMP they
lack the support of pivot elements.
These new solver will not work with matrices with a bad condition (like the
ones from the spline interpolation).
So they are no replacement for the existing lu decomposition.



I have no write acces to the svn repository. If you think my updates are
usefully, please commit them.

If you think these updates are not necessary i would like to suggest a
radical different solution:

I have had a look on several solver packages in the last year. The most
exiting packages to me are
the meschach implementation of real/complex vector/matrix computation
including
most of the fastest linear equation/eigenvalue solver for dense and sparse
matrices!
and the Trilinos approach of the sandia national laboratories which is based
on many free available serial and parallel solvers.

Meschach is available here:
http://www.netlib.org/c/meschach/

Trilinos is available here:
http://trilinos.sandia.gov/

Meschach is implemented in C and is able to replace the gmath library
completely with much better implementations.
I have integrated meschach into grass for testing purposes and, after some
very basic testing, it works nicely.
But im unsure if the licence of meschach is compatible with the GPL? Can
somebody please verify this? Im not an expert in licensing.

The Trilinos approach is quite different. It is a C++ framework with Python
bindings which provides a simplified and
unified interface to most of the fastest serial and parallel solver
available (BLAS, LAPACK, SuperLU, MUMPS, UMPFPACK, METIS....).
It was never easier for me to start a parallel ilu/gmres computation of a
huge sparse matrix on a cluster with some simple Python instructions.
But because Trilinos is implemented in C++, all grass modules which should
benefit from Trilinos must be implemented in C++ to.

So devs ... what do you think?

Best regards
Soeren
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-dev/attachments/20090110/9f382726/attachment.html


More information about the grass-dev mailing list