Head section correction. New matrix conversion functions added.

Integrated the band matrix cholesky solver from lidar library
into gmath library. Added tests for band matrix solver and matrix
conversion functions. Updated the doxygen documentation.


git-svn-id: https://svn.osgeo.org/grass/grass/trunk@42049 15284696-431f-4ddb-bdfa-cd5b030d7da7

include/gmath.h

@@ -137,11 +137,21 @@ extern G_math_spvector **G_math_alloc_spmatrix(int );
 extern void G_math_free_spmatrix(G_math_spvector ** , int );
 extern void G_math_free_spvector(G_math_spvector * );
 extern int G_math_add_spvector(G_math_spvector **, G_math_spvector * , int );
+extern G_math_spvector **G_math_A_to_Asp(double **, int, double);
+extern double **G_math_Asp_to_A(G_math_spvector **, int);
+extern double **G_math_Asp_to_band_matrix(G_math_spvector **, int, int);
+extern G_math_spvector **G_math_band_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon);
+extern void G_math_print_spmatrix(G_math_spvector ** Asp, int rows);
+
+/*Symmetric band matrix handling */
+extern double **G_math_matrix_to_band_matrix(double **, int, int);
+extern double **G_math_band_matrix_to_matrix(double **A, int rows, int bandwidth);
 
 /*linear equation solver, most of them are multithreaded wih OpenMP*/
 extern int G_math_solver_gauss(double **, double *, double *, int );
 extern int G_math_solver_lu(double **, double *, double *, int );
 extern int G_math_solver_cholesky(double **, double *, double *, int , int );
+extern void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int bandwidth);
 extern int G_math_solver_jacobi(double **, double *, double *, int , int , double , double );
 extern int G_math_solver_gs(double **, double *, double *, int , int , double , double );
 extern int G_math_solver_pcg(double **, double *, double *, int , int , double , int );
@@ -157,8 +167,10 @@ extern int G_math_solver_sparse_bicgstab(G_math_spvector **, double *, double *,
 extern void G_math_gauss_elimination(double **, double *, int );
 extern void G_math_lu_decomposition(double **, double *, int );
 extern int G_math_cholesky_decomposition(double **, int , int );
+extern void G_math_cholesky_band_decomposition(double **, double **, int, int);
 extern void G_math_backward_substitution(double **, double *, double *, int );
 extern void G_math_forward_substitution(double **, double *, double *, int );
+extern void G_math_cholesky_band_substitution(double **, double *, double *, int, int);
 
 /*BLAS like level 1,2 and 3 functions*/
 

lib/gmath/ATLAS_wrapper_blas_level_1.c

@@ -1,14 +1,14 @@
 
 /*****************************************************************************
 *
-* MODULE:       Grass PDE Numerical Library
+* MODULE:       Grass numerical math interface
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
-* 		soerengebbert <at> gmx <dot> de
+* 		soerengebbert <at> googlemail <dot> com
 *               
-* PURPOSE:      blas level 1 like functions   
-* 		part of the gpde library
+* PURPOSE:      gras blas implementation
+* 		part of the gmath library
 *               
-* COPYRIGHT:    (C) 2000 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS

lib/gmath/blas_level_1.c

@@ -1,20 +1,20 @@
 
 /*****************************************************************************
- *
- * MODULE:       Grass Gmath Library
- * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2007
- * 		soerengebbert <at> gmx <dot> de
- *               
- * PURPOSE:      blas level 1 like functions   
- * 		part of the gmath library
- *               
- * COPYRIGHT:    (C) 2000 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *               License (>=v2). Read the file COPYING that comes with GRASS
- *               for details.
- *
- *****************************************************************************/
+*
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
+*               
+* PURPOSE:      gras blas implementation
+* 		part of the gmath library
+*               
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
+*
+*               This program is free software under the GNU General Public
+*               License (>=v2). Read the file COPYING that comes with GRASS
+*               for details.
+*
+*****************************************************************************/
 
 #include <math.h>
 #include <unistd.h>

lib/gmath/blas_level_2.c

@@ -1,20 +1,20 @@
 
 /*****************************************************************************
- *
- * MODULE:       Grass PDE Numerical Library
- * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2007
- * 		soerengebbert <at> gmx <dot> de
- *               
- * PURPOSE:      linear equation system solvers
- * 		part of the gpde library
- *               
- * COPYRIGHT:    (C) 2007 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *               License (>=v2). Read the file COPYING that comes with GRASS
- *               for details.
- *
- *****************************************************************************/
+*
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
+*               
+* PURPOSE:      gras blas implementation
+* 		part of the gmath library
+*               
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
+*
+*               This program is free software under the GNU General Public
+*               License (>=v2). Read the file COPYING that comes with GRASS
+*               for details.
+*
+*****************************************************************************/
 
 #include <math.h>
 #include <unistd.h>

lib/gmath/blas_level_3.c

@@ -1,14 +1,14 @@
 
 /*****************************************************************************
 *
-* MODULE:       Grass PDE Numerical Library
-* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2007
-* 		soerengebbert <at> gmx <dot> de
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
 *               
-* PURPOSE:      linear equation system solvers
-* 		part of the gpde library
+* PURPOSE:      gras blas implementation
+* 		part of the gmath library
 *               
-* COPYRIGHT:    (C) 2007 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS

lib/gmath/ccmath_grass_wrapper.c

@@ -1,3 +1,22 @@
+
+/*****************************************************************************
+*
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
+*               
+* PURPOSE:      ccmath library function wrapper
+* 		part of the gmath library
+*               
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
+*
+*               This program is free software under the GNU General Public
+*               License (>=v2). Read the file COPYING that comes with GRASS
+*               for details.
+*
+*****************************************************************************/
+
+
 #if defined(HAVE_CCMATH)
 #include <ccmath.h>
 #else

lib/gmath/solvers_classic_iter.c

@@ -1,14 +1,14 @@
 
 /*****************************************************************************
 *
-* MODULE:       Grass PDE Numerical Library
+* MODULE:       Grass numerical math interface
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
-* 		soerengebbert <at> gmx <dot> de
+* 		soerengebbert <at> googlemail <dot> com
 *               
 * PURPOSE:      linear equation system solvers
-* 		part of the gpde library
+* 		part of the gmath library
 *               
-* COPYRIGHT:    (C) 2007 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS

lib/gmath/solvers_direct.c

@@ -1,20 +1,20 @@
 
 /*****************************************************************************
- *
- * MODULE:       Grass PDE Numerical Library
- * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
- * 		soerengebbert <at> gmx <dot> de
- *               
- * PURPOSE:      direkt linear equation system solvers
- * 		part of the gpde library
- *               
- * COPYRIGHT:    (C) 2007 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *               License (>=v2). Read the file COPYING that comes with GRASS
- *               for details.
- *
- *****************************************************************************/
+*
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
+*               
+* PURPOSE:      linear equation system solvers
+* 		part of the gmath library
+*               
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
+*
+*               This program is free software under the GNU General Public
+*               License (>=v2). Read the file COPYING that comes with GRASS
+*               for details.
+*
+*****************************************************************************/
 
 #include <math.h>
 #include <unistd.h>

lib/gmath/solvers_direct_cholesky_band.c

@@ -0,0 +1,304 @@
+#include <stdlib.h>		
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+
+/**
+ * \brief Cholesky decomposition of a band matrix
+ *
+ * \param A (double**) the input band matrix
+ * \param T (double**) the resulting lower tringular band matrix
+ * \param rows (int) number of rows
+ * \param bandwidth (int) the bandwidth of the band matrix
+ *
+ * */
+
+void G_math_cholesky_band_decomposition(double **A, double **T, int rows, int bandwidth)
+{
+    int i, j, k, end;
+    double sum;
+
+    G_debug(2, "G_math_cholesky_band_decomposition(): n=%d  bandwidth=%d", rows, bandwidth);
+
+    for (i = 0; i < rows; i++) {
+	G_percent(i, rows, 2);
+	for (j = 0; j < bandwidth; j++) {
+	    sum = A[i][j];
+	    /* start = 1 */
+	    /* end = bandwidth - j or i + 1 */
+	    end = ((bandwidth - j) < (i + 1) ? (bandwidth - j) : (i + 1));
+	    for (k = 1; k < end; k++)
+		sum -= T[i - k][k] * T[i - k][j + k];
+	    if (j == 0) {
+		if (sum <= 0.0)
+		    G_fatal_error(_("Decomposition failed at row %i and col %i"), i, j);
+		T[i][0] = sqrt(sum);
+	    }
+	    else
+		T[i][j] = sum / T[i][0];
+	}
+    }
+
+    G_percent(i, rows, 2);
+    return;
+}
+
+/**
+ * \brief Cholesky band matrix solver for linear equation systems of type Ax = b 
+ *
+ * \param A (double**) the input band matrix
+ * \param x (double*) the resulting vector, result is written in here
+ * \param b (double*) the right hand side of Ax = b
+ * \param rows (int) number of rows
+ * \param bandwidth (int) the bandwidth of the band matrix
+ *
+ * */
+
+void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int bandwidth)
+{
+
+    double **T;
+
+    T = G_alloc_matrix(rows, bandwidth);
+
+    G_math_cholesky_band_decomposition(A, T, rows, bandwidth);	/* T computation                */
+    G_math_cholesky_band_substitution(T, x, b, rows, bandwidth);
+
+    G_free_matrix(T);
+
+    return;
+}
+
+/**
+ * \brief Forward and backward substitution of a lower tringular band matrix of A from system Ax = b
+ *
+ * \param T (double**) the lower band triangle band matrix
+ * \param x (double*) the resulting vector
+ * \param b (double*) the right hand side of Ax = b
+ * \param rows (int) number of rows
+ * \param bandwidth (int) the bandwidth of the band matrix
+ *
+ * */
+
+void G_math_cholesky_band_substitution(double **T, double *x, double *b, int rows, int bandwidth)
+{
+
+    int i, j, start, end;
+
+    /* Forward substitution */
+    x[0] = b[0] / T[0][0];
+    for (i = 1; i < rows; i++) {
+	x[i] = b[i];
+	/* start = 0 or i - bandwidth + 1 */
+	start = ((i - bandwidth + 1) < 0 ? 0 : (i - bandwidth + 1));
+	/* end = i */
+	for (j = start; j < i; j++)
+	    x[i] -= T[j][i - j] * x[j];
+	x[i] = x[i] / T[i][0];
+    }
+
+    /* Backward substitution */
+    x[rows - 1] = x[rows - 1] / T[rows - 1][0];
+    for (i = rows - 2; i >= 0; i--) {
+	/* start = i + 1 */
+	/* end = rows or bandwidth + i */
+	end = (rows < (bandwidth + i) ? rows : (bandwidth + i));
+	for (j = i + 1; j < end; j++)
+	    x[i] -= T[i][j - i] * x[j];
+	x[i] = x[i] / T[i][0];
+    }
+
+    return;
+}
+
+/*--------------------------------------------------------------------------------------*/
+/* Tcholetsky matrix invertion */
+
+void G_math_cholesky_band_invert(double **A, double *invAdiag, int rows, int bandwidth)
+{
+    double **T = NULL;
+    double *vect = NULL;
+    int i, j, k, start;
+    double sum;
+
+    T = G_alloc_matrix(rows, bandwidth);
+    vect = G_alloc_vector(rows);
+
+    /* T computation                */
+    G_math_cholesky_band_decomposition(A, T, rows, bandwidth);
+
+    /* T Diagonal invertion */
+    for (i = 0; i < rows; i++) {
+	T[i][0] = 1.0 / T[i][0];
+    }
+
+    /* A Diagonal invertion */
+    for (i = 0; i < rows; i++) {
+	vect[0] = T[i][0];
+	invAdiag[i] = vect[0] * vect[0];
+	for (j = i + 1; j < rows; j++) {
+	    sum = 0.0;
+	    /* start = i or j - bandwidth + 1 */
+	    start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
+	    /* end = j */
+	    for (k = start; k < j; k++) {
+		sum -= vect[k - i] * T[k][j - k];
+	    }
+	    vect[j - i] = sum * T[j][0];
+	    invAdiag[i] += vect[j - i] * vect[j - i];
+	}
+    }
+
+    G_free_matrix(T);
+    G_free_vector(vect);
+
+    return;
+}
+
+/*--------------------------------------------------------------------------------------*/
+/* Tcholetsky matrix solution and invertion */
+
+void G_math_solver_cholesky_band_invert(double **A, double *x, double *b, double *invAdiag,
+		   int rows, int bandwidth)
+{
+
+    double **T = NULL;
+    double *vect = NULL;
+    int i, j, k, start, end;
+    double sum;
+
+    T = G_alloc_matrix(rows, bandwidth);
+    vect = G_alloc_vector(rows);
+
+    /* T computation                */
+    G_math_cholesky_band_decomposition(A, T, rows, bandwidth);
+    G_math_cholesky_band_substitution(T, x, b, rows, bandwidth);
+
+    /* T Diagonal invertion */
+    for (i = 0; i < rows; i++) {
+	T[i][0] = 1.0 / T[i][0];
+    }
+
+    /* A Diagonal invertion */
+    for (i = 0; i < rows; i++) {
+	vect[0] = T[i][0];
+	invAdiag[i] = vect[0] * vect[0];
+	for (j = i + 1; j < rows; j++) {
+	    sum = 0.0;
+	    /* start = i or j - bandwidth + 1 */
+	    start = ((j - bandwidth + 1) < i ? i : (j - bandwidth + 1));
+	    /* end = j */
+	    for (k = start; k < j; k++) {
+		sum -= vect[k - i] * T[k][j - k];
+	    }
+	    vect[j - i] = sum * T[j][0];
+	    invAdiag[i] += vect[j - i] * vect[j - i];
+	}
+    }
+
+    G_free_matrix(T);
+    G_free_vector(vect);
+
+    return;
+}
+
+/**
+ * \brief Convert a symmetrix matrix into a band matrix
+ *
+ * \verbatim
+ Symmetric matrix with bandwidth of 3
+
+ 5 2 1 0
+ 2 5 2 1
+ 1 2 5 2
+ 0 1 2 5
+
+ will be converted into the band matrix
+ 
+ 5 2 1
+ 5 2 1
+ 5 2 0
+ 5 0 0
+
+  \endverbatim
+
+ * \param A (double**) the symmetric matrix
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \return B (double**) new created band matrix 
+ * */
+
+double **G_math_matrix_to_band_matrix(double **A, int rows, int bandwidth)
+{
+    int i, j;
+    double **B = G_alloc_matrix(rows, bandwidth);
+    double tmp;
+
+    for(i = 0; i < rows; i++) {
+       for(j = 0; j < bandwidth; j++) {
+          if(i + j < rows) {
+            tmp = A[i][i + j]; 
+            B[i][j] = tmp;
+          } else {
+            B[i][j] = 0.0;
+          }
+       }
+    }
+
+    return B;
+}
+
+
+/**
+ * \brief Convert a band matrix into a symmetric matrix
+ *
+ * \verbatim
+ Such a band matrix with banwidth 3
+ 
+ 5 2 1
+ 5 2 1
+ 5 2 0
+ 5 0 0
+
+ Will be converted into this symmetric matrix
+
+ 5 2 1 0
+ 2 5 2 1
+ 1 2 5 2
+ 0 1 2 5
+
+  \endverbatim
+ * \param A (double**) the band matrix
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \return B (double**) new created symmetric matrix 
+ * */
+
+double **G_math_band_matrix_to_matrix(double **A, int rows, int bandwidth)
+{
+    int i, j;
+    double **B = G_alloc_matrix(rows, rows);
+    double tmp;
+
+    for(i = 0; i < rows; i++) {
+       for(j = 0; j < bandwidth; j++) {
+          tmp = A[i][j];
+          if(i + j < rows) {
+            B[i][i + j] = tmp; 
+          } 
+       }
+    }
+
+    /*Symmetry*/
+    for(i = 0; i < rows; i++) {
+       for(j = i; j < rows; j++) {
+          B[j][i] = B[i][j]; 
+       }
+    }
+
+    return B;
+}
+
+

lib/gmath/solvers_krylov.c

@@ -1,14 +1,14 @@
 
 /*****************************************************************************
 *
-* MODULE:       Grass PDE Numerical Library
+* MODULE:       Grass numerical math interface
 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
-* 		soerengebbert <at> gmx <dot> de
+* 		soerengebbert <at> googlemail <dot> com
 *               
 * PURPOSE:      linear equation system solvers
-* 		part of the gpde library
+* 		part of the gmath library
 *               
-* COPYRIGHT:    (C) 2000 by the GRASS Development Team
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
 *
 *               This program is free software under the GNU General Public
 *               License (>=v2). Read the file COPYING that comes with GRASS

lib/gmath/sparse_matrix.c

@@ -1,20 +1,20 @@
 
 /*****************************************************************************
- *
- * MODULE:       Grass Gmath Library
- * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
- * 		soerengebbert <at> gmx <dot> de
- *               
- * PURPOSE:      functions to manage linear equation systems
- * 		part of the gmath library
- *               
- * COPYRIGHT:    (C) 2000 by the GRASS Development Team
- *
- *               This program is free software under the GNU General Public
- *               License (>=v2). Read the file COPYING that comes with GRASS
- *               for details.
- *
- *****************************************************************************/
+*
+* MODULE:       Grass numerical math interface
+* AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+* 		soerengebbert <at> googlemail <dot> com
+*               
+* PURPOSE:      linear equation system solvers
+* 		part of the gmath library
+*               
+* COPYRIGHT:    (C) 2010 by the GRASS Development Team
+*
+*               This program is free software under the GNU General Public
+*               License (>=v2). Read the file COPYING that comes with GRASS
+*               for details.
+*
+*****************************************************************************/
 
 #include <stdlib.h>
 #include <math.h>
@@ -194,6 +194,52 @@ double **G_math_Asp_to_A(G_math_spvector ** Asp, int rows)
 }
 
 /*!
+ * \brief Convert a symmetric sparse matrix into a band matrix
+ *
+ \verbatim
+ Symmetric matrix with bandwidth of 3
+
+ 5 2 1 0
+ 2 5 2 1
+ 1 2 5 2
+ 0 1 2 5
+
+ will be converted into the band matrix
+ 
+ 5 2 1
+ 5 2 1
+ 5 2 0
+ 5 0 0
+
+ \endverbatim
+ * \param Asp (G_math_spvector **) 
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \return (double **) the resulting band matrix [rows][bandwidth]
+ *
+ * */
+double **G_math_Asp_to_band_matrix(G_math_spvector ** Asp, int rows, int bandwidth)
+{
+    int i, j;
+
+    double **A = NULL;
+
+    A = G_alloc_matrix(rows, bandwidth);
+
+    for (i = 0; i < rows; i++) {
+	for (j = 0; j < Asp[i]->cols; j++) {
+	   if(Asp[i]->index[j] == i) {
+	      A[i][0] = Asp[i]->values[j];
+	   } else if (Asp[i]->index[j] > i) {
+	      A[i][Asp[i]->index[j] - i] = Asp[i]->values[j];
+	   }
+	}
+    }
+    return A;
+}
+
+
+/*!
  * \brief Convert a quadratic matrix into a sparse matrix
  *
  * This function is multi-threaded with OpenMP. It creates its own parallel OpenMP region.
@@ -238,3 +284,61 @@ G_math_spvector **G_math_A_to_Asp(double **A, int rows, double epsilon)
     }
     return Asp;
 }
+
+
+
+/*!
+ * \brief Convert a band matrix into a sparse matrix
+ *
+ * WARNING:
+ * This function is experimental, do not use.
+ * Only the upper triangle matrix of the band strcuture is copied.
+ *
+ * \param A (double **) the band matrix
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \param epsilon (double) -- non-zero values are greater then epsilon
+ * \return (G_math_spvector **)
+ *
+ * */
+G_math_spvector **G_math_band_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)
+{
+    int i, j;
+
+    int nonull, count = 0;
+
+    G_math_spvector **Asp = NULL;
+
+    Asp = G_math_alloc_spmatrix(rows);
+
+    for (i = 0; i < rows; i++) {
+	nonull = 0;
+	/*Count the number of non zero entries */
+	for (j = 0; j < bandwidth; j++) {
+	    if (A[i][j] > epsilon)
+		nonull++;
+	}
+
+	/*Allocate the sparse vector and insert values */
+
+	G_math_spvector *v = G_math_alloc_spvector(nonull);
+
+	count = 0;
+	if (A[i][0] > epsilon) {
+	    v->index[count] = i;
+	    v->values[count] = A[i][0];
+	    count++;
+	}
+
+	for (j = 1; j < bandwidth; j++) {
+	    if (A[i][j] > epsilon && i + j < rows) {
+		v->index[count] = i + j;
+		v->values[count] = A[i][j];
+		count++;
+	    }
+	}
+	/*Add vector to sparse matrix */
+	G_math_add_spvector(Asp, v, i);
+    }
+    return Asp;
+}

lib/gmath/test/test_gmath_lib.h

@@ -103,6 +103,9 @@ extern G_math_f_les *create_normal_unsymmetric_f_nquad_les_A(int rows, int cols)
 /* direct and iterative solvers */
 extern int unit_test_solvers(void);
 
+/* Test the matrix conversion dense -> band ->sparse and vis versa */
+extern int unit_test_matrix_conversion(void);
+
 /* ccmath wrapper tests*/
 int unit_test_ccmath_wrapper(void);
 

lib/gmath/test/test_main.c

@@ -43,7 +43,7 @@ void set_params(void) {
     param.unit->key = "unit";
     param.unit->type = TYPE_STRING;
     param.unit->required = NO;
-    param.unit->options = "blas1,blas2,blas3,solver,ccmath";
+    param.unit->options = "blas1,blas2,blas3,solver,ccmath,matconv";
     param.unit->description = _("Choose the unit tests to run");
 
     param.integration = G_define_option();
@@ -122,6 +122,7 @@ int main(int argc, char *argv[]) {
         returnstat += unit_test_blas_level_2();
         returnstat += unit_test_blas_level_3();
         returnstat += unit_test_solvers();
+        returnstat += unit_test_matrix_conversion();
         returnstat += unit_test_ccmath_wrapper();
 
     }
@@ -153,6 +154,9 @@ int main(int argc, char *argv[]) {
                     if (strcmp(param.unit->answers[i], "ccmath") == 0)
                         returnstat += unit_test_ccmath_wrapper();
 
+                    if (strcmp(param.unit->answers[i], "matconv") == 0)
+                        returnstat += unit_test_matrix_conversion();
+
                     i++;
                 }
         }

lib/gmath/test/test_matrix_conversion.c

@@ -0,0 +1,167 @@
+/*****************************************************************************
+ *
+ * MODULE:       Grass PDE Numerical Library
+ * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
+ * 		soerengebbert <at> gmx <dot> de
+ *               
+ * PURPOSE:      Unit tests for les solving
+ *
+ * COPYRIGHT:    (C) 2000 by the GRASS Development Team
+ *
+ *               This program is free software under the GNU General Public
+ *               License (>=v2). Read the file COPYING that comes with GRASS
+ *               for details.
+ *
+ *****************************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <grass/glocale.h>
+#include <grass/gmath.h>
+#include "test_gmath_lib.h"
+
+#define EPSILON_DIRECT 1.0E-10
+#define EPSILON_ITER 1.0E-4
+
+/* prototypes */
+static int test_matrix_conversion(void);
+
+/* ************************************************************************* */
+/* Performe the solver unit tests ****************************************** */
+/* ************************************************************************* */
+int unit_test_matrix_conversion(void)
+{
+	int sum = 0;
+
+	G_message(_("\n++ Running matrix conversion unit tests ++"));
+
+	sum += test_matrix_conversion();
+
+	if (sum > 0)
+		G_warning(_("\n-- Matrix conversion unit tests failure --"));
+	else
+		G_message(_("\n-- Matrix conversion unit tests finished successfully --"));
+
+	return sum;
+}
+
+/* ************************************************************************* */
+/* ************************************************************************* */
+/* ************************************************************************* */
+
+static void print_matrix(double **A, int rows, int cols)
+{
+   int i, j;
+   for(i = 0; i < rows; i++) {
+      for(j = 0; j < cols; j++) {
+	printf("%g ", A[i][j]);
+      }
+      printf("\n");
+   }
+}
+
+/* *************************************************************** */
+/* Test all implemented solvers for sparse and normal matrix *** */
+/* *************************************************************** */
+int test_matrix_conversion(void)
+{
+	int sum = 0;
+	int i, j;
+	double val = 0.0;
+	double **A;
+	double **B;
+	double **C;
+	double **D;
+	double **E;
+	double **F;
+	G_math_spvector **Asp;
+	G_math_spvector **Asp2;
+
+	A = G_alloc_matrix(5,5);
+	F = G_alloc_matrix(5,5);
+	G_message("\t * Creating symmetric matrix\n");
+
+	A[0][0] = 8;
+	A[1][1] = 7;
+	A[2][2] = 6;
+	A[3][3] = 5;
+	A[4][4] = 4;
+	
+	A[0][1] = 3;
+	A[0][2] = 0;
+	A[0][3] = 1;
+	A[0][4] = 0;
+	
+	A[1][2] = 3;
+	A[1][3] = 0;
+	A[1][4] = 1;
+	
+	A[2][3] = 3;
+	
+	A[3][4] = 3;
+
+	for(i = 0; i < 5; i++)
+	   for(j = 0; j < 5; j++)
+	      A[j][i] = A[i][j];
+
+	print_matrix(A, 5, 5);
+
+	G_message("\t * Test matrix to band matrix conversion\n");
+
+        B = G_math_matrix_to_band_matrix(A, 5, 4);
+
+	print_matrix(B, 5, 4);
+
+	G_message("\t * Test matrix to sparse matrix conversion\n");
+
+        Asp = G_math_A_to_Asp(A, 5, 0.0);
+
+	G_math_print_spmatrix(Asp, 5);
+
+	G_message("\t * Test sparse matrix to matrix conversion\n");
+
+        C = G_math_Asp_to_A(Asp, 5);
+
+	print_matrix(C, 5, 5);
+
+	G_message("\t * Test sparse matrix to band matrix conversion\n");
+
+        D = G_math_Asp_to_band_matrix(Asp, 5, 4);
+
+	print_matrix(D, 5, 4);
+
+	/*Check the band matrix results*/
+	G_math_d_aA_B(B, D, -1, F, 5, 4);
+	G_math_d_asum_norm(F[0], &val, 20);
+	   
+	if (val != 0)
+	{
+		G_warning("Error in band matrix conversion");
+		sum++;
+	}
+
+	G_message("\t * Test band matrix to matrix conversion\n");
+
+        E = G_math_band_matrix_to_matrix(D, 5, 4);
+
+	print_matrix(E, 5, 5);
+
+	/*Check the matrix results*/
+	G_math_d_aA_B(A, E, -1, F, 5, 5);
+	G_math_d_asum_norm(F[0], &val, 25);
+
+	if (val != 0)
+	{
+		G_warning("Error in matrix conversion");
+		sum++;
+	}
+
+	G_message("\t * Test band matrix to sparse matrix conversion\n");
+
+        Asp2 = G_math_band_matrix_to_Asp(D, 5, 4, 0.0);
+	G_math_print_spmatrix(Asp2, 5);
+
+	return sum;
+}
+

lib/gmath/test/test_solvers.c

@@ -431,7 +431,7 @@ int test_solvers(void)
 	G_math_d_asum_norm(les->x, &val, les->rows);
 	if ((val - (double)les->rows) > EPSILON_DIRECT)
 	{
-		G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
+		G_warning("Error in G_math_solver_lu abs %2.20f != %i", val,
 				les->rows);
 		sum++;
 	}
@@ -445,7 +445,26 @@ int test_solvers(void)
 	G_math_d_asum_norm(les->x, &val, les->rows);
 	if ((val - (double)les->rows) > EPSILON_DIRECT)
 	{
-		G_warning("Error in G_math_solver_gauss abs %2.20f != %i", val,
+		G_warning("Error in G_math_solver_solver_cholesky abs %2.20f != %i", val,
+				les->rows);
+		sum++;
+	}
+	G_math_print_les(les);
+	G_math_free_les(les);
+
+	G_message("\t * testing cholesky band decomposition solver with symmetric matrix\n");
+	les = create_normal_symmetric_les(TEST_NUM_ROWS);
+	G_math_print_les(les);
+	/* Create a band matrix*/
+	G_message("\t * Creating symmetric band matrix\n");
+	les->A = G_math_matrix_to_band_matrix(les->A, les->rows, les->rows);
+	G_math_print_les(les);
+
+	/*cholesky*/G_math_solver_cholesky_band(les->A, les->x, les->b, les->rows,les->rows);
+	G_math_d_asum_norm(les->x, &val, les->rows);
+	if ((val - (double)les->rows) > EPSILON_DIRECT)
+	{
+		G_warning("Error in G_math_solver_solver_cholesky_band abs %2.20f != %i", val,
 				les->rows);
 		sum++;
 	}