Dear developers,<br>after a loooong time of development i would like to commit some changes to the gmath and gdpe lib <br>and some dependent raster modules.<br><br>The following changes have been made:<br><br>* The gpde solver (gauss, jacobi, LU, Cholesky, CG, PCG and BiCGStab) are now located in the gmath library.<br>
Because of that several files in the gpde directory have been removed or updated.<br>* Tests and benchmarks for single and multi threaded computation are available.<br>* Blas level 1, 2 and 3 functions are implemented with tests, benchmarks and basic support for sparse matrices.<br>
The implemented BLAS functions are not optimized for high performance, but multi threaded with OpenMP.<br>* Some new solver uses the BLAS 1 and 2 functions of the grass implementation.<br>* 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.<br>
* r.gwflow and r3.gwflow are updated to use the gmath library.<br>* Updates to gmath.h and gisdefs.h<br><br>The patch file for the latest grass7 trunk is located here:<br><a href="http://www-pool.math.tu-berlin.de/~soeren/grass/files/soeren_gmath_gpde.diff.gz">http://www-pool.math.tu-berlin.de/~soeren/grass/files/soeren_gmath_gpde.diff.gz</a><br>
<br>A patched grass7 trunk is available here:<br><a href="http://www-pool.math.tu-berlin.de/~soeren/grass/files/grass_trunk_soeren_gpde_gmath.tar.bz2">http://www-pool.math.tu-berlin.de/~soeren/grass/files/grass_trunk_soeren_gpde_gmath.tar.bz2</a><br>
<br>To apply the patch unzip the diff file, change directory to the grass_trunk root directory and type:<br>cat soeren_gmath_gpde.diff | patch -p0<br><br>Why a grass BLAS implementation?<br>There was a discussion more than a year ago, to replace the dependency in <br>
the gmath library to the FORTRAN BLAS and LAPACK libraries with <br>the ATLAS BLAS/LAPACK implementation. <br>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,<br>
which uses the grass BLAS implementation as backup in case the ATLAS library is not available. Wrapper for ATLAS BLAS level <br>2 and 3 functions must be developed to gain any benefit from the ATLAS library.<br><br>Replacement of existing lu decomposition?<br>
Most of the new solvers for linear equation systems (Ax = b) are paralleled with OpenMP. This works<br>fine with the krylov-subspace methods (cg, pcg and bicgstab) which also work with sparse matrices,<br>but not with the gauss and lu decomposition solver. The speedup of booth solver is quite bad. <br>
I was unable to replace the existing lu decomposition, which is lend from numerical recipes,<br>with a high performance, parallel algorithm which works also with sparse matrices.<br>This is much to complicated for me. <br>
Because the new lu and gauss solver are simply paralleled with OpenMP they lack the support of pivot elements.<br>These new solver will not work with matrices with a bad condition (like the ones from the spline interpolation). <br>
So they are no replacement for the existing lu decomposition.<br><br><br><br>I have no write acces to the svn repository. If you think my updates are usefully, please commit them.<br><br>If you think these updates are not necessary i would like to suggest a radical different solution:<br>
<br>I have had a look on several solver packages in the last year. The most exiting packages to me are<br>the meschach implementation of real/complex vector/matrix computation including<br>most of the fastest linear equation/eigenvalue solver for dense and sparse matrices! <br>
and the Trilinos approach of the sandia national laboratories which is based on many free available serial and parallel solvers.<br><br>Meschach is available here:<br><a href="http://www.netlib.org/c/meschach/">http://www.netlib.org/c/meschach/</a><br>
<br>Trilinos is available here:<br><a href="http://trilinos.sandia.gov/">http://trilinos.sandia.gov/</a><br><br>Meschach is implemented in C and is able to replace the gmath library completely with much better implementations. <br>
I have integrated meschach into grass for testing purposes and, after some very basic testing, it works nicely. <br>But im unsure if the licence of meschach is compatible with the GPL? Can somebody please verify this? Im not an expert in licensing.<br>
<br>The Trilinos approach is quite different. It is a C++ framework with Python bindings which provides a simplified and<br>unified interface to most of the fastest serial and parallel solver available (BLAS, LAPACK, SuperLU, MUMPS, UMPFPACK, METIS....).<br>
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.<br>But because Trilinos is implemented in C++, all grass modules which should benefit from Trilinos must be implemented in C++ to.<br>
<br>So devs ... what do you think?<br><br>Best regards<br>Soeren<br>