[GRASS-dev] Licence problems in gmath LU decomposition
Helena Mitasova
hmitaso at unity.ncsu.edu
Tue May 1 17:36:42 EDT 2007
Soeren -
by default rst does not have 0 on the diagonal - it has the value of
the smoothing parameter
which is 0.1 or 0.01 (I would have to check, I just found it is
missing in the man page).
This keeps the solution stable.
It runs with zero smoothing too if you have sufficiently high tension,
but it gets unstable for small tension and zero smoothing. That is
why the defaults and
overall internal adjustments of the parameters are set to avoid it
and the program
gives warnings about the overshoots. There is more here (see e.g.
fig. 2)
http://skagit.meas.ncsu.edu/~helena/gmslab/papers/IEEEGRSL2005.pdf
But the G_ludcmp() needs to be replaced - it was discussed long time
ago and
I somehow thought that it was done and the function does not include
the num. recipes
program.
Thanks for the patch, I planned to look at the solver because of the
performance issues,
so this helps, although the strange performance behavior since
somewhere GRASS5*
may be more related to recently reported stunning degradation of
performance for
rasters between 4.* and 5.* than in lineq. solver. I always though it
is slower due to the
floating point introduction and did not realize that it affects
integer rasters too,
Helena
-------------- next part --------------
A non-text attachment was scrubbed...
Name: LINEQS.c
Type: application/octet-stream
Size: 3030 bytes
Desc: not available
Url : http://lists.osgeo.org/pipermail/grass-dev/attachments/20070501/ad796268/LINEQS.obj
-------------- next part --------------
Helena Mitasova
Dept. of Marine, Earth and Atm. Sciences
1125 Jordan Hall, NCSU Box 8208,
Raleigh NC 27695
http://skagit.meas.ncsu.edu/~helena/
On May 1, 2007, at 3:44 PM, S?ren Gebbert wrote:
> Dear developers,
> i just noticed that the code of is an exact copy of
> the numerical recipes algorithm. :(
> We need to replace it.
>
> I have implemented a new simple parallel LU decomposition in the
> gpde library (since several months in CVS).
> But this algorithm lacks the ability of pivoting (it would be very
> hard to implement this as a parallel version).
> So it is not usable within the rst library in grass and maybe
> others? This is a big problem.
>
> The rst library creates linear equation systems with a zero entry
> in the diagonal matrix. All the
> algorithms implemented in gpde library to solve linear equation
> systems are not able to
> handle this (LU, Gauss, jacobi, SOR) except the bicgstab algorithm.
> The bicgstab is a state of the art iterative linear equation solver
> algorithm for non-symmetric matrices.
> But it will fail on very bad conditioned matrices, so it is not the
> best choice if a direct solve can be used.
>
> I have created a very simple patch for the rst library (regular
> spline with tension) to use the bicgstab
> algorithm form the gpde lib instead of G_ludcmp(). Just for testing.
> It seems to work and it is not slower than G_ludcmp(). The patch is
> attached.
> It would be better if the rst lib would create matrices with a
> better condition and without
> zero entries at the diagonal matrix. Then we could use other gpde
> algorithms than bicgstab.
>
> The main question is, what to do with G_ludcmp()? There is a free
> implementation of the LU
> decomposition algorithm with pivoting in the GSL library, which we
> can use. Or we use the bicgstab
> algorithm in cases where the simple LU algorithm from the gpde
> library will fail?
>
> Any suggestions are welcome
> Best regards
> Soeren
> Index: interp_float/matrix.c
> ===================================================================
> RCS file: /home/grass/grassrepository/grass6/lib/rst/interp_float/
> matrix.c,v
> retrieving revision 2.3
> diff -u -u -r2.3 matrix.c
> --- interp_float/matrix.c 29 Mar 2006 20:05:08 -0000 2.3
> +++ interp_float/matrix.c 1 May 2007 19:41:01 -0000
> @@ -154,12 +154,13 @@
> matrix[i][j] = A[m];
> }
> }
> -
> - if (G_ludcmp(matrix,n_points+1,indx,&d)<=0) { /* find the
> inverse of the mat
> -rix */
> +/*
> + if (G_ludcmp(matrix,n_points+1,indx,&d)<=0) {
> +
> fprintf(stderr,"G_ludcmp() failed! n=%d\n",n_points);
> return -1;
> }
> + */
> /*
> G_free_vector(A);
> */
> Index: interp_float/segmen2d.c
> ===================================================================
> RCS file: /home/grass/grassrepository/grass6/lib/rst/interp_float/
> segmen2d.c,v
> retrieving revision 2.5
> diff -u -u -r2.5 segmen2d.c
> --- interp_float/segmen2d.c 9 Feb 2006 03:08:58 -0000 2.5
> +++ interp_float/segmen2d.c 1 May 2007 19:41:02 -0000
> @@ -8,6 +8,7 @@
> #include <math.h>
> #include <grass/gis.h>
> #include <grass/glocale.h>
> +#include <grass/N_pde.h>
>
> #include <grass/interpf.h>
>
> @@ -263,7 +264,30 @@
> for (i = 0; i < data->n_points; i++)
> b[i+1] = data->points[i].z;
> b[0] = 0.;
> - G_lubksb(matrix,data->n_points+1,indx,b);
> +
> + /*create the gpde lineare equation system and transfer the data*/
> + N_les * les = N_alloc_les(data->n_points + 1, N_NORMAL_LES);
> + for(i = 0; i < data->n_points + 1; i++) {
> + les->b[i] = b[i];
> + for(j = 0; j < data->n_points + 1; j++) {
> + les->A[i][j] = matrix[i][j];
> + }
> + }
> + /*solve it with the bicgstab algorithm*/
> + N_solver_bicgstab(les, 100000, 0.000000001);
> + double d;
> +
> + /* for testing */
> + //G_ludcmp(matrix,data->n_points + 1,indx,&d);
> + //G_lubksb(matrix,data->n_points + 1,indx,b);
> +
> + /*transfer the result*/
> + for(i = 0; i < data->n_points + 1; i++)
> + b[i] = les->x[i];
> +
> + /*release memory*/
> + N_free_les(les);
> +
> params->check_points
> (params,data,b,ertot,zmin,dnorm,skip_point);
> } else if (segtest == 1)
> {
> @@ -273,7 +297,6 @@
> G_lubksb(matrix,data->n_points,indx,b);
> params->check_points
> (params,data,b,ertot,zmin,dnorm,skip_point);
> }
> -
> } /*cv loop */
>
> if(!params->cv)
> _______________________________________________
> grass-dev mailing list
> grass-dev at grass.itc.it
> http://grass.itc.it/mailman/listinfo/grass-dev
More information about the grass-dev
mailing list