[GRASS-dev] Licence problems in gmath LU decomposition

[Sören Gebbert]

Dear developers,
i just noticed that the code of G_ludcmp() 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
-------------- next part --------------
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)