[GRASS-dev] Licence problems in gmath LU decomposition

[Sören Gebbert]

I have created a new patch for the rst library to test the gpde solvers
with the rst matrices. I have added a new function to  rst/interp_float/segmen2d.c (solve_matrix).
If you want to test the different gpde solvers (only cg and bicgstab will work),
just remove or add a comment within the new function.
The created linear eqution system can be printed to stdout,
just remove the comment before N_print_les(les).

This is a quick and dirty hack and only for testing purpose. :)

You have to patch the rst library and v.surf.rst to get it work.

This is the v.surf.rst patch:

Index: Makefile
===================================================================
RCS file: /home/grass/grassrepository/grass6/vector/v.surf.rst/Makefile,v
retrieving revision 1.7
diff -u -r1.7 Makefile
--- Makefile    20 Feb 2004 17:43:33 -0000      1.7
+++ Makefile    2 May 2007 01:16:20 -0000
@@ -5,7 +5,7 @@
  EXTRA_CLEAN_DIRS=doxygenhtml

  LIBES     = $(INTERPDATALIB) $(QTREELIB) $(INTERPFLLIB) $(BITMAPLIB) $(LINKMLIB) \
-            $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(SITESLIB) $(GISLIB) $(DATETIMELIB) $(GMATHLIB)
+            $(VECTLIB) $(VECTLIB_REAL) $(DBMILIB) $(SITESLIB) $(GISLIB) $(DATETIMELIB) $(GMATHLIB) $(GPDELIB)
  DEPENDENCIES = $(INTERPDATADEP) $(QTREEDEP) $(INTERPFLDEP) $(BITMAPDEP) $(LINKMDEP) \
              $(VECTDEP) $(DBMIDEP) $(GISDEP) $(DATETIMEDEP) $(GMATHDEP)
  EXTRA_INC = $(VECT_INC)





The rst library patch is attached.

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 -r2.3 matrix.c
--- interp_float/matrix.c	29 Mar 2006 20:05:08 -0000	2.3
+++ interp_float/matrix.c	2 May 2007 01:14:15 -0000
@@ -154,12 +154,6 @@
         matrix[i][j] = A[m];
       }
     }
-
-    if (G_ludcmp(matrix,n_points+1,indx,&d)<=0) { /* find the inverse of the mat
-rix */
-      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 -r2.5 segmen2d.c
--- interp_float/segmen2d.c	9 Feb 2006 03:08:58 -0000	2.5
+++ interp_float/segmen2d.c	2 May 2007 01:14:17 -0000
@@ -8,10 +8,12 @@
 #include <math.h>
 #include <grass/gis.h>
 #include <grass/glocale.h>
+#include <grass/N_pde.h>
 
 #include <grass/interpf.h>
 
 static double  smallest_segment( struct multtree *, int);
+void solve_matrix(struct quaddata *data, double *b, double **matrix, int *indx, int ofs);
 
 int IL_interp_segments_2d (
     struct interp_params *params,
@@ -249,31 +251,21 @@
      }
       }/* segment area test */
 
-        if (!params->cv){
-        if (params->matrix_create(params,data->points,data->n_points,
+        if (!params->cv) {
+        	if (params->matrix_create(params,data->points,data->n_points,
                                        matrix,indx) < 0)  return -1;
-        }
-        else if (segtest == 1)
-	{
-                if (params->matrix_create(params,data->points,data->n_points-1,
-                                        matrix,indx) < 0)  return -1;
-        }
 
-        if (!params->cv) {
-        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);
-        params->check_points(params,data,b,ertot,zmin,dnorm,skip_point);
+		  /*solve the equation system*/
+		  solve_matrix(data, b, matrix, indx, 1);
+        	  params->check_points(params,data,b,ertot,zmin,dnorm,skip_point);
         } else if (segtest == 1)
         {
-                for (i = 0; i < data->n_points-1; i++)
-                  b[i+1] = data->points[i].z;
-                  b[0] = 0.;
-                  G_lubksb(matrix,data->n_points,indx,b);
+                if (params->matrix_create(params,data->points,data->n_points-1,
+                                        matrix,indx) < 0)  return -1;
+		  /*solve the equation system*/
+		  solve_matrix(data, b, matrix, indx, 0);
                   params->check_points(params,data,b,ertot,zmin,dnorm,skip_point);
         }
-
         } /*cv loop */
         
 	if(!params->cv)
@@ -340,4 +332,44 @@
 }
 
 
+/*Solve linear equation system with different solvers*/
+void solve_matrix(struct quaddata *data, double *b, double **matrix, int *indx, int ofs)
+{
+  int i, j;
+  double d;
 
+  N_les * les  = N_alloc_les(data->n_points + ofs, N_NORMAL_LES);
+
+  for (i = 0; i < data->n_points - 1 - ofs; i++)
+    b[i+1] = data->points[i].z;
+  b[0] = 0.;
+
+  /*create the gpde lineare equation system and transfer the data*/
+  for(i = 0; i < data->n_points + ofs; i++) {
+    les->b[i] = b[i];
+    for(j = 0; j < data->n_points + ofs; j++) {
+      les->A[i][j] = matrix[i][j];
+    }
+  }
+  /*solve it with the bicgstab algorithm*/
+  //N_print_les(les);
+  //N_solver_gauss(les);
+  //N_solver_lu(les);
+  //N_solver_jacobi(les, 100, 1, 0.00001);
+  //N_solver_SOR(les, 100, 1, 0.00001);
+  //N_solver_bicgstab(les, 100000, 0.000000001);
+  N_solver_cg(les, 100000, 0.000000001);
+
+  /*transfer the result*/
+  for(i = 0; i < data->n_points + ofs; i++) 
+  b[i] = les->x[i];
+
+  /*release memory*/
+  N_free_les(les);
+
+  /* for testing, uncomment this if you want to test the lu decomposition */
+  //G_ludcmp(matrix,data->n_points + ofs,indx,&d);
+  //G_lubksb(matrix,data->n_points + ofs,indx,b);
+  
+return;
+}