New band matrix functions. Implemented new conjugate gradent band matrix

solver. Created tests for band matrix functions and solver. Simple band
matrix benchmarks are available. Documentation update.


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

include/gmath.h

@@ -139,38 +139,42 @@ 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);
+extern double **G_math_Asp_to_sband_matrix(G_math_spvector **, int, int);
+extern G_math_spvector **G_math_sband_matrix_to_Asp(double **, int, int, double);
+extern void G_math_print_spmatrix(G_math_spvector **, int);
+extern void G_math_Ax_sparse(G_math_spvector **, double *, double *, int );
 
 /*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);
+extern double **G_math_matrix_to_sband_matrix(double **, int, int);
+extern double **G_math_sband_matrix_to_matrix(double **, int, int);
+extern void G_math_Ax_sband(double ** A, double *x, double *y, 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 void G_math_solver_cholesky_sband(double **, double *, double *, int, int);
 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 );
 
+extern int G_math_solver_pcg(double **, double *, double *, int , int , double , int );
 extern int G_math_solver_cg(double **, double *, double *, int , int , double );
+extern int G_math_solver_cg_sband(double **, double *, double *, int, int, int, double);
 extern int G_math_solver_bicgstab(double **, double *, double *, int , int , double );
 extern int G_math_solver_sparse_jacobi(G_math_spvector **, double *, double *, int , int , double , double );
 extern int G_math_solver_sparse_gs(G_math_spvector **, double *, double *, int , int , double , double );
 extern int G_math_solver_sparse_pcg(G_math_spvector **, double *, double *, int , int , double , int );
 extern int G_math_solver_sparse_cg(G_math_spvector **, double *, double *, int , int , double );
 extern int G_math_solver_sparse_bicgstab(G_math_spvector **, double *, double *, int , int , double );
+
 /* solver algoithms and helper functions*/
 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_cholesky_sband_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);
+extern void G_math_cholesky_sband_substitution(double **, double *, double *, int, int);
 
 /*BLAS like level 1,2 and 3 functions*/
 
@@ -214,7 +218,6 @@ extern void G_math_daxpy(double *, double *, double , int );
 extern void G_math_saxpy(float *, float *, float , int );
 
 /*level 2 matrix - vector grass implementation with OpenMP thread support*/
-extern void G_math_Ax_sparse(G_math_spvector **, double *, double *, int );
 extern void G_math_d_Ax(double **, double *, double *, int , int );
 extern void G_math_f_Ax(float **, float *, float *, int , int );
 extern void G_math_d_x_dyad_y(double *, double *, double **, int, int );

lib/gmath/blas_level_2.c

@@ -27,38 +27,6 @@
 
 #define EPSILON 0.00000000000000001
 
-/*!
- * \brief Compute the matrix - vector product  
- * of sparse matrix **Asp and vector x.
- *
- * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
- *
- * y = A * x
- *
- *
- * \param Asp (G_math_spvector **) 
- * \param x (double) *)
- * \param y (double * )
- * \param rows (int)
- * \return (void)
- *
- * */
-void G_math_Ax_sparse(G_math_spvector ** Asp, double *x, double *y, int rows)
-{
-    int i, j;
-
-    double tmp;
-
-#pragma omp for schedule (static) private(i, j, tmp)
-    for (i = 0; i < rows; i++) {
-	tmp = 0;
-	for (j = 0; j < Asp[i]->cols; j++) {
-	    tmp += Asp[i]->values[j] * x[Asp[i]->index[j]];
-	}
-	y[i] = tmp;
-    }
-    return;
-}
 
 /*!
  * \brief Compute the matrix - vector product  

lib/gmath/gmathlib.dox

@@ -158,7 +158,6 @@ void G_math_saxpy(float *, float *, float , int )<br>
 
 \subsubsection blas_level_2 Level 2 matrix - vector GRASS implementation with OpenMP thread support
 
-void G_math_Ax_sparse(G_math_spvector **, double *, double *, int )<br>
 void G_math_d_Ax(double **, double *, double *, int , int )<br>
 void G_math_f_Ax(float **, float *, float *, int , int )<br>
 void G_math_d_x_dyad_y(double *, double *, double **, int, int )<br>
@@ -398,27 +397,28 @@ void G_math_free_spvector(G_math_spvector * )<br>
 int G_math_add_spvector(G_math_spvector **, G_math_spvector * , int )<br>
 G_math_spvector **G_math_A_to_Asp(double **, int, double)<br>
 double **G_math_Asp_to_A(G_math_spvector **, int)<br>
-double **G_math_Asp_to_band_matrix(G_math_spvector **, int, int)<br>
-G_math_spvector **G_math_band_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)<br>
+double **G_math_Asp_to_sband_matrix(G_math_spvector **, int, int)<br>
+G_math_spvector **G_math_sband_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)<br>
 void G_math_print_spmatrix(G_math_spvector ** Asp, int rows)<br>
+void G_math_Ax_sparse(G_math_spvector **, double *, double *, int )<br>
 
 <P>
 Conversion of band matrices
 <P>
-double **G_math_matrix_to_band_matrix(double **, int, int)<br>
-double **G_math_band_matrix_to_matrix(double **A, int rows, int bandwidth)<br>
-
+double **G_math_matrix_to_sband_matrix(double **, int, int)<br>
+double **G_math_sband_matrix_to_matrix(double **A, int rows, int bandwidth)<br>
+void G_math_Ax_sband(double ** A, double *x, double *y, int rows, int bandwidth)<br>
 <P>
 Direct linear equation solver
 <P>
 int G_math_solver_gauss(double **, double *, double *, int )<br>
 int G_math_solver_lu(double **, double *, double *, int )<br>
 int G_math_solver_cholesky(double **, double *, double *, int , int )<br>
-void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int bandwidth)<br>
-void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int bandwidth)<br>
+void G_math_solver_cholesky_sband(double **A, double *x, double *b, int rows, int bandwidth)<br>
+void G_math_solver_cholesky_sband(double **A, double *x, double *b, int rows, int bandwidth)<br>
 
 <P>
-Classic iterative linear equation solver for dense and sparse matrices
+Classic iterative linear equation solver for dense- and sparse-matrices
 <P>
 int G_math_solver_jacobi(double **, double *, double *, int , int , double , double )<br>
 int G_math_solver_gs(double **, double *, double *, int , int , double , double )<br>
@@ -426,10 +426,11 @@ int G_math_solver_sparse_jacobi(G_math_spvector **, double *, double *, int , in
 int G_math_solver_sparse_gs(G_math_spvector **, double *, double *, int , int , double , double )<br>
 
 <P>
-Krylov-Subspace iterative linear equation solver for dense and sparse matrices
+Krylov-Subspace iterative linear equation solver for dense-, band- and sparse-matrices
 <P>
 int G_math_solver_pcg(double **, double *, double *, int , int , double , int )<br>
 int G_math_solver_cg(double **, double *, double *, int , int , double )<br>
+int G_math_solver_cg_sband(double **, double *, double *, int, int, int, double)<br>
 int G_math_solver_bicgstab(double **, double *, double *, int , int , double )<br>
 int G_math_solver_sparse_pcg(G_math_spvector **, double *, double *, int , int , double , int )<br>
 int G_math_solver_sparse_cg(G_math_spvector **, double *, double *, int , int , double )<br>

lib/gmath/solvers_direct_cholesky_band.c

@@ -6,21 +6,21 @@
 #include <grass/glocale.h>
 
 /**
- * \brief Cholesky decomposition of a band matrix
+ * \brief Cholesky decomposition of a symmetric band matrix
  *
- * \param A (double**) the input band matrix
- * \param T (double**) the resulting lower tringular band matrix
+ * \param A (double**) the input symmetric band matrix
+ * \param T (double**) the resulting lower tringular sband matrix
  * \param rows (int) number of rows
- * \param bandwidth (int) the bandwidth of the band matrix
+ * \param bandwidth (int) the bandwidth of the symmetric band matrix
  *
  * */
 
-void G_math_cholesky_band_decomposition(double **A, double **T, int rows, int bandwidth)
+void G_math_cholesky_sband_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);
+    G_debug(2, "G_math_cholesky_sband_decomposition(): n=%d  bandwidth=%d", rows, bandwidth);
 
     for (i = 0; i < rows; i++) {
 	G_percent(i, rows, 2);
@@ -46,25 +46,25 @@ void G_math_cholesky_band_decomposition(double **A, double **T, int rows, int ba
 }
 
 /**
- * \brief Cholesky band matrix solver for linear equation systems of type Ax = b 
+ * \brief Cholesky symmetric band matrix solver for linear equation systems of type Ax = b 
  *
- * \param A (double**) the input band matrix
+ * \param A (double**) the input symmetric 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
+ * \param bandwidth (int) the bandwidth of the symmetric band matrix
  *
  * */
 
-void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int bandwidth)
+void G_math_solver_cholesky_sband(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_math_cholesky_sband_decomposition(A, T, rows, bandwidth);	/* T computation                */
+    G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
 
     G_free_matrix(T);
 
@@ -72,17 +72,17 @@ void G_math_solver_cholesky_band(double **A, double *x, double *b, int rows, int
 }
 
 /**
- * \brief Forward and backward substitution of a lower tringular band matrix of A from system Ax = b
+ * \brief Forward and backward substitution of a lower tringular symmetric band matrix of A from system Ax = b
  *
- * \param T (double**) the lower band triangle band matrix
+ * \param T (double**) the lower triangle symmetric 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
+ * \param bandwidth (int) the bandwidth of the symmetric band matrix
  *
  * */
 
-void G_math_cholesky_band_substitution(double **T, double *x, double *b, int rows, int bandwidth)
+void G_math_cholesky_sband_substitution(double **T, double *x, double *b, int rows, int bandwidth)
 {
 
     int i, j, start, end;
@@ -116,7 +116,7 @@ void G_math_cholesky_band_substitution(double **T, double *x, double *b, int row
 /*--------------------------------------------------------------------------------------*/
 /* Tcholetsky matrix invertion */
 
-void G_math_cholesky_band_invert(double **A, double *invAdiag, int rows, int bandwidth)
+void G_math_cholesky_sband_invert(double **A, double *invAdiag, int rows, int bandwidth)
 {
     double **T = NULL;
     double *vect = NULL;
@@ -127,7 +127,7 @@ void G_math_cholesky_band_invert(double **A, double *invAdiag, int rows, int ban
     vect = G_alloc_vector(rows);
 
     /* T computation                */
-    G_math_cholesky_band_decomposition(A, T, rows, bandwidth);
+    G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
 
     /* T Diagonal invertion */
     for (i = 0; i < rows; i++) {
@@ -160,7 +160,7 @@ void G_math_cholesky_band_invert(double **A, double *invAdiag, int rows, int ban
 /*--------------------------------------------------------------------------------------*/
 /* Tcholetsky matrix solution and invertion */
 
-void G_math_solver_cholesky_band_invert(double **A, double *x, double *b, double *invAdiag,
+void G_math_solver_cholesky_sband_invert(double **A, double *x, double *b, double *invAdiag,
 		   int rows, int bandwidth)
 {
 
@@ -173,8 +173,8 @@ void G_math_solver_cholesky_band_invert(double **A, double *x, double *b, double
     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);
+    G_math_cholesky_sband_decomposition(A, T, rows, bandwidth);
+    G_math_cholesky_sband_substitution(T, x, b, rows, bandwidth);
 
     /* T Diagonal invertion */
     for (i = 0; i < rows; i++) {
@@ -204,101 +204,3 @@ void G_math_solver_cholesky_band_invert(double **A, double *x, double *b, double
     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

@@ -28,9 +28,9 @@ static G_math_spvector **create_diag_precond_matrix(double **A,
 						    G_math_spvector ** Asp,
 						    int rows, int prec);
 static int solver_pcg(double **A, G_math_spvector ** Asp, double *x,
-		      double *b, int rows, int maxit, double err, int prec);
+		      double *b, int rows, int maxit, double err, int prec, int has_band, int bandwidth);
 static int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
-		     int rows, int maxit, double err);
+		     int rows, int maxit, double err, int has_band, int bandwidth);
 static int solver_bicgstab(double **A, G_math_spvector ** Asp, double *x,
 			   double *b, int rows, int maxit, double err);
 
@@ -61,10 +61,43 @@ int G_math_solver_pcg(double **A, double *x, double *b, int rows, int maxit,
 		      double err, int prec)
 {
 
-    return solver_pcg(A, NULL, x, b, rows, maxit, err, prec);
+    return solver_pcg(A, NULL, x, b, rows, maxit, err, prec, 0, 0);
 }
 
 /*!
+ * \brief The iterative preconditioned conjugate gradients solver for symmetric positive definite band matrices
+ *
+ * WARNING: The preconditioning of symmetric band matrices is not implemented yet
+ *
+ * This iterative solver works with symmetric positive definite band matrices.
+ *
+ * This solver solves the linear equation system:
+ *  A x = b
+ *
+ * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
+ * solver will abort the calculation and writes the current result into the vector x.
+ * The parameter <i>err</i> defines the error break criteria for the solver.
+ *
+ * \param A (double **) -- the positive definite band matrix
+ * \param x (double *) -- the value vector
+ * \param b (double *) -- the right hand side
+ * \param rows (int)
+ * \param bandwidth (int) -- bandwidth of matrix A
+ * \param maxit (int) -- the maximum number of iterations
+ * \param err (double) -- defines the error break criteria
+ * \param prec (int) -- the preconditioner which shoudl be used 1,2 or 3
+ * \return (int) -- 1 - success, 2 - not finisehd but success, 0 - matrix singular, -1 - could not solve the les
+ * 
+ * */
+int G_math_solver_pcg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit,
+		      double err, int prec)
+{
+    G_fatal_error("Preconditioning of band matrics is not implemented yet");
+    return solver_pcg(A, NULL, x, b, rows, maxit, err, prec, 1, bandwidth);
+}
+
+
+/*!
  * \brief The iterative preconditioned conjugate gradients solver for sparse symmetric positive definite matrices
  *
  * This iterative solver works with symmetric positive definite sparse matrices.
@@ -90,11 +123,11 @@ int G_math_solver_sparse_pcg(G_math_spvector ** Asp, double *x, double *b,
 			     int rows, int maxit, double err, int prec)
 {
 
-    return solver_pcg(NULL, Asp, x, b, rows, maxit, err, prec);
+    return solver_pcg(NULL, Asp, x, b, rows, maxit, err, prec, 0, 0);
 }
 
 int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
-	       int rows, int maxit, double err, int prec)
+	       int rows, int maxit, double err, int prec, int has_band, int bandwidth)
 {
     double *r, *z;
 
@@ -131,6 +164,8 @@ int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
     {
 	if (Asp)
 	    G_math_Ax_sparse(Asp, x, v, rows);
+	else if(has_band)
+	    G_math_Ax_sband(A, x, v, rows, bandwidth);
 	else
 	    G_math_d_Ax(A, x, v, rows, rows);
 
@@ -156,6 +191,8 @@ int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
 	{
 	    if (Asp)
 		G_math_Ax_sparse(Asp, p, v, rows);
+	    else if(has_band)
+		G_math_Ax_sband(A, p, v, rows, bandwidth);
 	    else
 		G_math_d_Ax(A, p, v, rows, rows);
 
@@ -180,6 +217,8 @@ int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
 	    if (m % 50 == 1) {
 		if (Asp)
 		    G_math_Ax_sparse(Asp, x, v, rows);
+		else if(has_band)
+		    G_math_Ax_sband(A, x, v, rows, bandwidth);
 		else
 		    G_math_d_Ax(A, x, v, rows, rows);
 
@@ -270,10 +309,38 @@ int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
 int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit,
 		     double err)
 {
-    return solver_cg(A, NULL, x, b, rows, maxit, err);
+    return solver_cg(A, NULL, x, b, rows, maxit, err, 0, 0);
 }
 
 /*!
+ * \brief The iterative conjugate gradients solver for symmetric positive definite band matrices
+ *
+ * This iterative solver works with symmetric positive definite band matrices.
+ *
+ * This solver solves the linear equation system:
+ *  A x = b
+ *
+ * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
+ * solver will abort the calculation and writes the current result into the vector x.
+ * The parameter <i>err</i> defines the error break criteria for the solver.
+ *
+ * \param A (double **) -- the symmetric positive definit band matrix
+ * \param x (double *) -- the value vector
+ * \param b (double *) -- the right hand side
+ * \param rows (int)
+ * \param bandwidth (int) -- the bandwidth of matrix A
+ * \param maxit (int) -- the maximum number of iterations
+ * \param err (double) -- defines the error break criteria
+ * \return (int) -- 1 - success, 2 - not finisehd but success, 0 - matrix singular, -1 - could not solve the les
+ * 
+ * */
+int G_math_solver_cg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit, double err)
+{
+    return solver_cg(A, NULL, x, b, rows, maxit, err, 1, bandwidth);
+}
+
+
+/*!
  * \brief The iterative conjugate gradients solver for sparse symmetric positive definite matrices
  *
  * This iterative solver works with symmetric positive definite sparse matrices.
@@ -297,12 +364,12 @@ int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit,
 int G_math_solver_sparse_cg(G_math_spvector ** Asp, double *x, double *b,
 			    int rows, int maxit, double err)
 {
-    return solver_cg(NULL, Asp, x, b, rows, maxit, err);
+    return solver_cg(NULL, Asp, x, b, rows, maxit, err, 0, 0);
 }
 
 
 int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
-	      int rows, int maxit, double err)
+	      int rows, int maxit, double err, int has_band, int bandwidth)
 {
     double *r;
 
@@ -332,6 +399,8 @@ int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
     {
 	if (Asp)
 	    G_math_Ax_sparse(Asp, x, v, rows);
+	else if(has_band)
+	    G_math_Ax_sband(A, x, v, rows, bandwidth);
 	else
 	    G_math_d_Ax(A, x, v, rows, rows);
 
@@ -356,6 +425,8 @@ int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
 	{
 	    if (Asp)
 		G_math_Ax_sparse(Asp, p, v, rows);
+	    else if(has_band)
+		G_math_Ax_sband(A, p, v, rows, bandwidth);
 	    else
 		G_math_d_Ax(A, p, v, rows, rows);
 
@@ -378,6 +449,8 @@ int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
 	    if (m % 50 == 1) {
 		if (Asp)
 		    G_math_Ax_sparse(Asp, x, v, rows);
+		else if(has_band)
+		    G_math_Ax_sband(A, x, v, rows, bandwidth);
 		else
 		    G_math_d_Ax(A, x, v, rows, rows);
 

lib/gmath/sparse_matrix.c

@@ -194,7 +194,7 @@ double **G_math_Asp_to_A(G_math_spvector ** Asp, int rows)
 }
 
 /*!
- * \brief Convert a symmetric sparse matrix into a band matrix
+ * \brief Convert a symmetric sparse matrix into a symmetric band matrix
  *
  \verbatim
  Symmetric matrix with bandwidth of 3
@@ -215,10 +215,10 @@ double **G_math_Asp_to_A(G_math_spvector ** Asp, int rows)
  * \param Asp (G_math_spvector **) 
  * \param rows (int)
  * \param bandwidth (int)
- * \return (double **) the resulting band matrix [rows][bandwidth]
+ * \return (double **) the resulting ymmetric band matrix [rows][bandwidth]
  *
  * */
-double **G_math_Asp_to_band_matrix(G_math_spvector ** Asp, int rows, int bandwidth)
+double **G_math_Asp_to_sband_matrix(G_math_spvector ** Asp, int rows, int bandwidth)
 {
     int i, j;
 
@@ -288,20 +288,20 @@ G_math_spvector **G_math_A_to_Asp(double **A, int rows, double epsilon)
 
 
 /*!
- * \brief Convert a band matrix into a sparse matrix
+ * \brief Convert a symmetric 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 A (double **) the symmetric 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)
+G_math_spvector **G_math_sband_matrix_to_Asp(double **A, int rows, int bandwidth, double epsilon)
 {
     int i, j;
 
@@ -342,3 +342,37 @@ G_math_spvector **G_math_band_matrix_to_Asp(double **A, int rows, int bandwidth,
     }
     return Asp;
 }
+
+
+/*!
+ * \brief Compute the matrix - vector product  
+ * of sparse matrix **Asp and vector x.
+ *
+ * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
+ *
+ * y = A * x
+ *
+ *
+ * \param Asp (G_math_spvector **) 
+ * \param x (double) *)
+ * \param y (double * )
+ * \param rows (int)
+ * \return (void)
+ *
+ * */
+void G_math_Ax_sparse(G_math_spvector ** Asp, double *x, double *y, int rows)
+{
+    int i, j;
+
+    double tmp;
+
+#pragma omp for schedule (static) private(i, j, tmp)
+    for (i = 0; i < rows; i++) {
+	tmp = 0;
+	for (j = 0; j < Asp[i]->cols; j++) {
+	    tmp += Asp[i]->values[j] * x[Asp[i]->index[j]];
+	}
+	y[i] = tmp;
+    }
+    return;
+}

lib/gmath/symmetric_band_matrix.c

@@ -0,0 +1,148 @@
+#include <stdlib.h>		
+#include <stdio.h>
+#include <math.h>
+#include <grass/gis.h>
+#include <grass/gmath.h>
+#include <grass/glocale.h>
+
+/**
+ * \brief Convert a symmetrix matrix into a symmetric 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 symmetric 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 symmetric band matrix 
+ * */
+
+double **G_math_matrix_to_sband_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 symmetric band matrix into a symmetric matrix
+ *
+ * \verbatim
+ Such a symmetric 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 symmetric band matrix
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \return B (double**) new created symmetric matrix 
+ * */
+
+double **G_math_sband_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;
+}
+
+
+/*!
+ * \brief Compute the matrix - vector product  
+ * of symmetric band matrix A and vector x.
+ *
+ * This function is multi-threaded with OpenMP and can be called within a parallel OpenMP region.
+ *
+ * y = A * x
+ *
+ *
+ * \param A (double **) 
+ * \param x (double) *)
+ * \param y (double * )
+ * \param rows (int)
+ * \param bandwidth (int)
+ * \return (void)
+ *
+ * */
+void G_math_Ax_sband(double ** A, double *x, double *y, int rows, int bandwidth)
+{
+    int i, j;
+
+    double tmp;
+
+#pragma omp for schedule (static) private(i, j, tmp)
+    for (i = 0; i < rows; i++) {
+	tmp = 0;
+	for (j = 0; j < bandwidth; j++) {
+	   if((i + j) < rows)
+	   	tmp += A[i][j]*x[i + j];
+	}
+	y[i] = tmp;
+    }
+#pragma omp for schedule (static) private(i, j, tmp)
+    for (i = 0; i < rows; i++) {
+	tmp = 0;
+	for (j = 1; j < bandwidth; j++) {
+	   	if(i + j < rows)
+			y[i + j] += A[i][j]*x[i];
+	}
+    }
+
+    return;
+}

lib/gmath/test/bench_solver_direct.c

@@ -94,6 +94,16 @@ int bench_solvers(int rows) {
     G_important_message("Computation time ccmath cholesky decomposition: %g\n", compute_time_difference(tstart, tend));
     G_math_free_les(les);
 
+    G_message("\t * benchmarking gmath cholesky band matrix decomposition solver with symmetric band matrix\n");
+
+    les = create_symmetric_band_les(rows);
+    gettimeofday(&tstart, NULL);
+    G_math_solver_cholesky_sband(les->A, les->x, les->b, les->rows, les->rows);
+    gettimeofday(&tend, NULL);
+    G_important_message("Computation time cholesky band matrix decomposition: %g\n", compute_time_difference(tstart, tend));
+    G_math_free_les(les);
+
+
     return 1;
 }
 

lib/gmath/test/bench_solver_krylov.c

@@ -87,6 +87,17 @@ int bench_solvers(int rows) {
     G_math_free_les(les);
     G_math_free_les(sples);
 
+    G_message("\t * benchmark cg solver with symmetric band matrix\n");
+
+    les = create_symmetric_band_les(rows);
+    
+    gettimeofday(&tstart, NULL);
+    G_math_solver_cg_sband(les->A, les->x, les->b, les->rows, les->rows, 250, 0.1e-9);
+    gettimeofday(&tend, NULL);
+    G_important_message("Computation time cg symmetric band matrix: %g\n", compute_time_difference(tstart, tend));
+    
+    G_math_free_les(les);
+
     G_message("\t * benchmark bicgstab solver with unsymmetric matrix\n");
 
     les = create_normal_unsymmetric_les(rows);

lib/gmath/test/test_blas2.c

@@ -65,12 +65,14 @@ int test_blas_level_2_double(void)
 
     int sum = 0;
     int rows = TEST_NUM_ROWS;
-    double **A, **B, **C, *x, *y, *z, a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0, f = 0.0, g = 0.0, h = 0.0, i = 0.0;
+    double **A, **B, **C, *x, *y, *z, a = 0.0, b = 0.0, c = 0.0, d = 0.0, e = 0.0, f = 0.0, g = 0.0, h = 0.0, i = 0.0, j =0.0;
 
     G_math_les *les;
     les = create_normal_unsymmetric_les(rows);
     G_math_les *sples;
-    sples = create_sparse_unsymmetric_les(rows);
+    sples = create_sparse_symmetric_les(rows);
+    G_math_les *bles;
+    bles = create_symmetric_band_les(rows);
 
     x = G_alloc_vector(rows);
     y = G_alloc_vector(rows);
@@ -86,28 +88,42 @@ int test_blas_level_2_double(void)
 
 #pragma omp parallel default(shared)
 {
-    G_math_Ax_sparse(sples->Asp, x, z, rows);
-    G_math_d_asum_norm(z, &a, rows);
+    G_message("Testing G_math_Ax_sparse");
+    G_math_Ax_sparse(sples->Asp, x, sples->b, rows);
+    G_math_print_les(sples);
+    G_math_d_asum_norm(sples->b, &a, rows);
 
+    G_message("Testing G_math_Ax_band");
+    G_math_Ax_sband(bles->A, x, bles->b, rows, rows);
+    G_math_print_les(bles);
+    G_math_d_asum_norm(bles->b, &j, rows);
+
+    G_message("Testing G_math_d_Ax");
     G_math_d_Ax(les->A, x, z, rows, rows);
     G_math_d_asum_norm(z, &b, rows);
 
+    G_message("Testing G_math_d_aAx_by");
     G_math_d_aAx_by(les->A, x, y, 1.0, 1.0, z, rows, rows);
     G_math_d_asum_norm(z, &c, rows);
 
+    G_message("Testing G_math_d_aAx_by");
     G_math_d_aAx_by(les->A, x, y, -1.0, 1.0, z, rows, rows);
     G_math_d_asum_norm(z, &d, rows);
 
+    G_message("Testing G_math_d_aAx_by");
     G_math_d_aAx_by(les->A, x, y, 1.0, 0.0, z, rows, rows);
     G_math_d_asum_norm(z, &e, rows);
 
+    G_message("Testing G_math_d_aAx_by");
     G_math_d_aAx_by(les->A, x, y, -1.0, -1.0, z, rows, rows);
     G_math_d_asum_norm(z, &f, rows);
 
+    G_message("Testing G_math_d_x_dyad_y");
     G_math_d_x_dyad_y(x, x, A, rows, rows);
     G_math_d_Ax(A, x, z, rows, rows);
     G_math_d_asum_norm(z, &g, rows);
 
+    G_message("Testing G_math_d_x_dyad_y");
     G_math_d_x_dyad_y(x, x, C, rows, rows);
     G_math_d_Ax(A, x, z, rows, rows);
     G_math_d_asum_norm(z, &h, rows);
@@ -120,6 +136,11 @@ int test_blas_level_2_double(void)
 	sum++;
     }
 
+    if(j - i > EPSILON) {
+    	G_message("Error in G_math_Ax_sband: %f != %f", i, j);
+	sum++;
+    }
+
     if(b - i > EPSILON) {
     	G_message("Error in G_math_d_Ax: %f != %f", i, b);
 	sum++;
@@ -163,6 +184,7 @@ int test_blas_level_2_double(void)
     G_free_vector(z);
 
     G_math_free_les(les);
+    G_math_free_les(bles);
     G_math_free_les(sples);
 
     if(A)
@@ -203,27 +225,34 @@ int test_blas_level_2_float(void)
 
 #pragma omp parallel default(shared)
 {
+    G_message("Testing G_math_f_Ax");
     G_math_f_Ax(les->A, x, z, rows, rows);
     G_math_f_asum_norm(z, &b, rows);
 
+    G_message("Testing G_math_f_aAx_by");
     G_math_f_aAx_by(les->A, x, y, 1.0, 1.0, z, rows, rows);
     G_math_f_asum_norm(z, &c, rows);
 
+    G_message("Testing G_math_f_aAx_by");
     G_math_f_aAx_by(les->A, x, y, -1.0, 1.0, z, rows, rows);
     G_math_f_asum_norm(z, &d, rows);
 
+    G_message("Testing G_math_f_aAx_by");
     G_math_f_aAx_by(les->A, x, y, 1.0, 0.0, z, rows, rows);
     G_math_f_asum_norm(z, &e, rows);
 
+    G_message("Testing G_math_f_aAx_by");
     G_math_f_aAx_by(les->A, x, y, -1.0, -1.0, z, rows, rows);
     G_math_f_asum_norm(z, &f, rows);
 
+    G_message("Testing G_math_f_x_dyad_y");
     G_math_f_x_dyad_y(x, x, A, rows, rows);
     G_math_f_Ax(A, x, z, rows, rows);
     G_math_f_asum_norm(z, &g, rows);
 
+    G_message("Testing G_math_f_x_dyad_y");
     G_math_f_x_dyad_y(x, x, C, rows, rows);
-     G_math_f_Ax(A, x, z, rows, rows);
+    G_math_f_Ax(A, x, z, rows, rows);
     G_math_f_asum_norm(z, &h, rows);
 }
 

lib/gmath/test/test_ccmath_wrapper.c

@@ -180,7 +180,7 @@ int test_ccmath_wrapper(void)
 	G_math_d_asum_norm(les->x, &val2, les->rows);
 	if ((val  - val2) > EPSILON_ITER)
 	{
-		G_warning("Error in G_math_eigv abs %2.20f != %i", val,
+		G_warning("Error in G_math_eigv abs %2.20f != %f", val,
 				val2);
 		sum++;
 	}
@@ -209,7 +209,7 @@ int test_ccmath_wrapper(void)
 	G_math_d_asum_norm(les->x, &val2, les->rows);
 	if ((val  - val2) > EPSILON_ITER)
 	{
-		G_warning("Error in G_math_eigen abs %2.20f != %i", val,
+		G_warning("Error in G_math_eigen abs %2.20f != %f", val,
 				val2);
 		sum++;
 	}

lib/gmath/test/test_gmath_lib.h

@@ -84,6 +84,7 @@ extern void fill_f_vector_scalar(float *x, float a, int rows);
 extern void fill_i_vector_scalar(int *x, int a, int rows);
 
 extern G_math_les *create_normal_symmetric_les(int rows);
+extern G_math_les *create_symmetric_band_les(int rows);
 extern G_math_les *create_normal_symmetric_pivot_les(int rows);
 extern G_math_les *create_normal_unsymmetric_les(int rows);
 extern G_math_les *create_sparse_symmetric_les(int rows);

lib/gmath/test/test_main.c

@@ -102,7 +102,6 @@ int main(int argc, char *argv[]) {
     G_gisinit(argv[0]);
 
     module = G_define_module();
-    module->keywords = _("test, gmath");
     module->description
             = _("Performs benchmarks, unit and integration tests for the gmath library");
 

lib/gmath/test/test_matrix_conversion.c

@@ -109,7 +109,7 @@ int test_matrix_conversion(void)
 
 	G_message("\t * Test matrix to band matrix conversion\n");
 
-        B = G_math_matrix_to_band_matrix(A, 5, 4);
+        B = G_math_matrix_to_sband_matrix(A, 5, 4);
 
 	print_matrix(B, 5, 4);
 
@@ -127,7 +127,7 @@ int test_matrix_conversion(void)
 
 	G_message("\t * Test sparse matrix to band matrix conversion\n");
 
-        D = G_math_Asp_to_band_matrix(Asp, 5, 4);
+        D = G_math_Asp_to_sband_matrix(Asp, 5, 4);
 
 	print_matrix(D, 5, 4);
 
@@ -143,7 +143,7 @@ int test_matrix_conversion(void)
 
 	G_message("\t * Test band matrix to matrix conversion\n");
 
-        E = G_math_band_matrix_to_matrix(D, 5, 4);
+        E = G_math_sband_matrix_to_matrix(D, 5, 4);
 
 	print_matrix(E, 5, 5);
 
@@ -159,7 +159,7 @@ int test_matrix_conversion(void)
 
 	G_message("\t * Test band matrix to sparse matrix conversion\n");
 
-        Asp2 = G_math_band_matrix_to_Asp(D, 5, 4, 0.0);
+        Asp2 = G_math_sband_matrix_to_Asp(D, 5, 4, 0.0);
 	G_math_print_spmatrix(Asp2, 5);
 
 	return sum;

lib/gmath/test/test_solvers.c

@@ -287,6 +287,7 @@ int test_solvers(void)
 	G_math_free_les(les);
 	G_math_free_les(sples);
 
+
 	G_message("\t * testing cg solver with symmetric bad conditioned matrix\n");
 
 	les = create_normal_symmetric_pivot_les(TEST_NUM_ROWS);
@@ -452,15 +453,15 @@ int test_solvers(void)
 	G_math_print_les(les);
 	G_math_free_les(les);
 
-	G_message("\t * testing cholesky band decomposition solver with symmetric matrix\n");
+	G_message("\t * testing cholesky band decomposition solver with symmetric band matrix 1\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);
+	les->A = G_math_matrix_to_sband_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);
+	/*cholesky*/G_math_solver_cholesky_sband(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)
 	{
@@ -471,6 +472,35 @@ int test_solvers(void)
 	G_math_print_les(les);
 	G_math_free_les(les);
 
+	G_message("\t * testing cholesky band decomposition solver with symmetric band matrix 2\n");
+        les = create_symmetric_band_les(TEST_NUM_ROWS);
+	G_math_print_les(les);
+
+	/*cholesky*/G_math_solver_cholesky_sband(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++;
+	}
+	G_math_print_les(les);
+	G_math_free_les(les);
+
+	G_message("\t * testing cg solver with symmetric band matrix\n");
+
+	les = create_symmetric_band_les(TEST_NUM_ROWS);
+
+	G_math_solver_cg_sband(les->A, les->x, les->b, les->rows, les->rows, 250, 0.1e-9);
+	G_math_d_asum_norm(les->x, &val, les->rows);
+	if ((val - (double)les->rows) > EPSILON_ITER)
+	{
+		G_warning("Error in G_math_solver_cg_sband abs %2.20f != %i", val, les->rows);
+		sum++;
+	}
+	G_math_print_les(les);
+	G_math_free_les(les);
+
 	return sum;
 }
 

lib/gmath/test/test_tools.c

@@ -54,11 +54,9 @@ G_math_les *create_normal_symmetric_les(int rows)
 		for (j = 0; j < size; j++)
 		{
 			if (j == i)
-				les->A[i][j] = (double)(1.0/(((double)i + 1.0) + ((double)j
-						+ 1.0)));
+				les->A[i][j] = (double)(1.0/(((double)i + 1.0) + ((double)j + 1.0)));
 			else
-				les->A[i][j] = (double)(1.0/((((double)i + 1.0) + ((double)j
-						+ 1.0) + 100.0)));
+				les->A[i][j] = (double)(1.0/((((double)i + 1.0) + ((double)j + 1.0) + 100.0)));
 
 			val += les->A[i][j];
 		}
@@ -69,6 +67,44 @@ G_math_les *create_normal_symmetric_les(int rows)
 	return les;
 }
 
+
+/* *************************************************************** */
+/* create a symmetric band matrix with values ** Hilbert matrix ** */
+/* *************************************************************** */
+G_math_les *create_symmetric_band_les(int rows)
+{
+	G_math_les *les;
+	int i, j;
+	int size =rows;
+	double val;
+
+        les = G_math_alloc_les(rows, G_MATH_NORMAL_LES);
+	for (i = 0; i < size; i++)
+	{
+		val = 0.0;
+		for (j = 0; j < size; j++)
+		{
+			if(i + j < size) {
+				les->A[i][j] = (double)(1.0/((((double)i + 1.0) + ((double)(i + j) + 1.0) + 100.0)));
+			} else if (j != i){
+			   	les->A[i][j] = 0.0;
+			}
+			if (j == i) {
+				les->A[i][0] = (double)(1.0/(((double)i + 1.0) + ((double)j + 1.0)));
+			} 
+
+			if (j == i) 
+				val += (double)(1.0/(((double)i + 1.0) + ((double)j + 1.0)));
+			else
+				val += (double)(1.0/((((double)i + 1.0) + ((double)j + 1.0) + 100.0)));
+		}
+		les->b[i] = val;
+		les->x[i] = 0.5;
+	}
+
+	return les;
+}
+
 /* ********************************************************************* */
 /* create a bad conditioned normal matrix with values ** Hilbert matrix  */
 /* ********************************************************************* */