[GRASS-SVN] r42051 - in grass/trunk: include lib/gmath
lib/gmath/test
svn_grass at osgeo.org
svn_grass at osgeo.org
Wed Apr 28 06:41:22 EDT 2010
Author: huhabla
Date: 2010-04-28 06:41:14 -0400 (Wed, 28 Apr 2010)
New Revision: 42051
Added:
grass/trunk/lib/gmath/symmetric_band_matrix.c
Modified:
grass/trunk/include/gmath.h
grass/trunk/lib/gmath/blas_level_2.c
grass/trunk/lib/gmath/gmathlib.dox
grass/trunk/lib/gmath/solvers_direct_cholesky_band.c
grass/trunk/lib/gmath/solvers_krylov.c
grass/trunk/lib/gmath/sparse_matrix.c
grass/trunk/lib/gmath/test/bench_solver_direct.c
grass/trunk/lib/gmath/test/bench_solver_krylov.c
grass/trunk/lib/gmath/test/test_blas2.c
grass/trunk/lib/gmath/test/test_ccmath_wrapper.c
grass/trunk/lib/gmath/test/test_gmath_lib.h
grass/trunk/lib/gmath/test/test_main.c
grass/trunk/lib/gmath/test/test_matrix_conversion.c
grass/trunk/lib/gmath/test/test_solvers.c
grass/trunk/lib/gmath/test/test_tools.c
Log:
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.
Modified: grass/trunk/include/gmath.h
===================================================================
--- grass/trunk/include/gmath.h 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/include/gmath.h 2010-04-28 10:41:14 UTC (rev 42051)
@@ -139,38 +139,42 @@
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_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_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 );
Modified: grass/trunk/lib/gmath/blas_level_2.c
===================================================================
--- grass/trunk/lib/gmath/blas_level_2.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/blas_level_2.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -27,39 +27,7 @@
#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
* of matrix A and vector x.
Modified: grass/trunk/lib/gmath/gmathlib.dox
===================================================================
--- grass/trunk/lib/gmath/gmathlib.dox 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/gmathlib.dox 2010-04-28 10:41:14 UTC (rev 42051)
@@ -158,7 +158,6 @@
\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 @@
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_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>
Modified: grass/trunk/lib/gmath/solvers_direct_cholesky_band.c
===================================================================
--- grass/trunk/lib/gmath/solvers_direct_cholesky_band.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/solvers_direct_cholesky_band.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -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 @@
}
/**
- * \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 @@
}
/**
- * \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 @@
/*--------------------------------------------------------------------------------------*/
/* 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 @@
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 @@
/*--------------------------------------------------------------------------------------*/
/* 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 @@
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 @@
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;
-}
-
-
Modified: grass/trunk/lib/gmath/solvers_krylov.c
===================================================================
--- grass/trunk/lib/gmath/solvers_krylov.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/solvers_krylov.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -28,9 +28,9 @@
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 @@
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 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 @@
{
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 @@
{
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 @@
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 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_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 @@
{
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 @@
{
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 @@
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);
Modified: grass/trunk/lib/gmath/sparse_matrix.c
===================================================================
--- grass/trunk/lib/gmath/sparse_matrix.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/sparse_matrix.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -194,7 +194,7 @@
}
/*!
- * \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 @@
* \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 @@
/*!
- * \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 @@
}
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;
+}
Added: grass/trunk/lib/gmath/symmetric_band_matrix.c
===================================================================
--- grass/trunk/lib/gmath/symmetric_band_matrix.c (rev 0)
+++ grass/trunk/lib/gmath/symmetric_band_matrix.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -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;
+}
Modified: grass/trunk/lib/gmath/test/bench_solver_direct.c
===================================================================
--- grass/trunk/lib/gmath/test/bench_solver_direct.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/bench_solver_direct.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -94,6 +94,16 @@
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;
}
Modified: grass/trunk/lib/gmath/test/bench_solver_krylov.c
===================================================================
--- grass/trunk/lib/gmath/test/bench_solver_krylov.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/bench_solver_krylov.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -87,6 +87,17 @@
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);
Modified: grass/trunk/lib/gmath/test/test_blas2.c
===================================================================
--- grass/trunk/lib/gmath/test/test_blas2.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_blas2.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -65,12 +65,14 @@
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 @@
#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 @@
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 @@
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 @@
#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);
}
Modified: grass/trunk/lib/gmath/test/test_ccmath_wrapper.c
===================================================================
--- grass/trunk/lib/gmath/test/test_ccmath_wrapper.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_ccmath_wrapper.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -180,7 +180,7 @@
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 @@
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++;
}
Modified: grass/trunk/lib/gmath/test/test_gmath_lib.h
===================================================================
--- grass/trunk/lib/gmath/test/test_gmath_lib.h 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_gmath_lib.h 2010-04-28 10:41:14 UTC (rev 42051)
@@ -84,6 +84,7 @@
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);
Modified: grass/trunk/lib/gmath/test/test_main.c
===================================================================
--- grass/trunk/lib/gmath/test/test_main.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_main.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -102,7 +102,6 @@
G_gisinit(argv[0]);
module = G_define_module();
- module->keywords = _("test, gmath");
module->description
= _("Performs benchmarks, unit and integration tests for the gmath library");
Modified: grass/trunk/lib/gmath/test/test_matrix_conversion.c
===================================================================
--- grass/trunk/lib/gmath/test/test_matrix_conversion.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_matrix_conversion.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -109,7 +109,7 @@
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 @@
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 @@
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 @@
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;
Modified: grass/trunk/lib/gmath/test/test_solvers.c
===================================================================
--- grass/trunk/lib/gmath/test/test_solvers.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_solvers.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -287,6 +287,7 @@
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 @@
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 @@
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;
}
Modified: grass/trunk/lib/gmath/test/test_tools.c
===================================================================
--- grass/trunk/lib/gmath/test/test_tools.c 2010-04-27 17:33:41 UTC (rev 42050)
+++ grass/trunk/lib/gmath/test/test_tools.c 2010-04-28 10:41:14 UTC (rev 42051)
@@ -54,11 +54,9 @@
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 @@
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 */
/* ********************************************************************* */
More information about the grass-commit
mailing list