[GRASS-SVN] r38889 - in grass/branches/develbranch_6/lib/gpde: .
test
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 27 16:52:10 EDT 2009
Author: huhabla
Date: 2009-08-27 16:52:10 -0400 (Thu, 27 Aug 2009)
New Revision: 38889
Added:
grass/branches/develbranch_6/lib/gpde/test/test.gpde.lib.html
Removed:
grass/branches/develbranch_6/lib/gpde/N_les_pivot.c
grass/branches/develbranch_6/lib/gpde/N_solvers.c
grass/branches/develbranch_6/lib/gpde/N_solvers_classic_iter.c
grass/branches/develbranch_6/lib/gpde/N_solvers_krylov.c
grass/branches/develbranch_6/lib/gpde/solvers_local_proto.h
grass/branches/develbranch_6/lib/gpde/test/test_solvers.c
Modified:
grass/branches/develbranch_6/lib/gpde/N_geom.c
grass/branches/develbranch_6/lib/gpde/N_gradient.c
grass/branches/develbranch_6/lib/gpde/N_gradient_calc.c
grass/branches/develbranch_6/lib/gpde/N_gwflow.c
grass/branches/develbranch_6/lib/gpde/N_gwflow.h
grass/branches/develbranch_6/lib/gpde/N_heatflow.h
grass/branches/develbranch_6/lib/gpde/N_les.c
grass/branches/develbranch_6/lib/gpde/N_les_assemble.c
grass/branches/develbranch_6/lib/gpde/N_pde.h
grass/branches/develbranch_6/lib/gpde/N_solute_transport.c
grass/branches/develbranch_6/lib/gpde/N_solute_transport.h
grass/branches/develbranch_6/lib/gpde/N_tools.c
grass/branches/develbranch_6/lib/gpde/gpdelib.dox
grass/branches/develbranch_6/lib/gpde/test/Makefile
grass/branches/develbranch_6/lib/gpde/test/test_arrays.c
grass/branches/develbranch_6/lib/gpde/test/test_assemble.c
grass/branches/develbranch_6/lib/gpde/test/test_geom.c
grass/branches/develbranch_6/lib/gpde/test/test_gpde_lib.h
grass/branches/develbranch_6/lib/gpde/test/test_gradient.c
grass/branches/develbranch_6/lib/gpde/test/test_gwflow.c
grass/branches/develbranch_6/lib/gpde/test/test_les.c
grass/branches/develbranch_6/lib/gpde/test/test_main.c
grass/branches/develbranch_6/lib/gpde/test/test_solute_transport.c
grass/branches/develbranch_6/lib/gpde/test/test_tools.c
Log:
Moving the linear equation solver and tests into the gmath library
Modified: grass/branches/develbranch_6/lib/gpde/N_geom.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_geom.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_geom.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -73,6 +73,7 @@
N_geom_data *N_init_geom_data_3d(G3D_Region * region3d, N_geom_data * geodata)
{
N_geom_data *geom = geodata;
+
struct Cell_head region2d;
#pragma omp critical
@@ -115,9 +116,13 @@
N_geom_data * geodata)
{
N_geom_data *geom = geodata;
+
struct Cell_head backup;
+
double meters;
+
short ll = 0;
+
int i;
@@ -158,10 +163,12 @@
"N_init_geom_data_2d: calculating the areas for non parametric projection");
geom->planimetric = 0;
- if (geom->area != NULL)
+ if (geom->area != NULL) {
G_free(geom->area);
- else
geom->area = G_calloc(geom->rows, sizeof(double));
+ } else {
+ geom->area = G_calloc(geom->rows, sizeof(double));
+ }
/*fill the area vector */
for (i = 0; i < geom->rows; i++) {
Modified: grass/branches/develbranch_6/lib/gpde/N_gradient.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_gradient.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_gradient.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -116,6 +116,7 @@
N_gradient_2d * gradient, int col, int row)
{
double NC = 0, SC = 0, WC = 0, EC = 0;
+
N_gradient_2d *grad = gradient;
@@ -250,6 +251,7 @@
int depth)
{
double NC, SC, WC, EC, TC, BC;
+
N_gradient_3d *grad = gradient;
NC = N_get_array_3d_d_value(field->y_array, col, row, depth);
@@ -631,6 +633,7 @@
N_gradient_neighbours_y * y)
{
N_gradient_neighbours_2d *grad;
+
int fail = 0;
G_debug(5,
@@ -706,9 +709,13 @@
int row)
{
double NWN, NEN, WC, EC, SWS, SES;
+
double NWW, NEE, NC, SC, SWW, SEE;
+
N_gradient_neighbours_2d *grad = NULL;
+
N_gradient_neighbours_x *grad_x = NULL;
+
N_gradient_neighbours_y *grad_y = NULL;
@@ -832,6 +839,7 @@
N_gradient_neighbours_z * zb)
{
N_gradient_neighbours_3d *grad;
+
int fail = 0;
G_debug(5,
Modified: grass/branches/develbranch_6/lib/gpde/N_gradient_calc.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_gradient_calc.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_gradient_calc.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -29,8 +29,11 @@
void N_calc_gradient_field_2d_stats(N_gradient_field_2d * field)
{
double minx, miny;
+
double maxx, maxy;
+
double sumx, sumy;
+
int nonullx, nonully;
G_debug(3,
@@ -110,8 +113,11 @@
gradfield)
{
int i, j;
+
int rows, cols;
+
double dx, dy, p1, p2, r1, r2, mean, grad, res;
+
N_gradient_field_2d *field = gradfield;
@@ -245,10 +251,15 @@
N_array_2d * y_comp)
{
int i, j;
+
int rows, cols;
+
double vx, vy;
+
N_array_2d *x = x_comp;
+
N_array_2d *y = y_comp;
+
N_gradient_2d grad;
@@ -300,8 +311,11 @@
void N_calc_gradient_field_3d_stats(N_gradient_field_3d * field)
{
double minx, miny, minz;
+
double maxx, maxy, maxz;
+
double sumx, sumy, sumz;
+
int nonullx, nonully, nonullz;
G_debug(3,
@@ -392,8 +406,11 @@
gradfield)
{
int i, j, k;
+
int cols, rows, depths;
+
double dx, dy, dz, p1, p2, r1, r2, mean, grad, res;
+
N_gradient_field_3d *field = gradfield;
@@ -582,11 +599,17 @@
N_array_3d * z_comp)
{
int i, j, k;
+
int rows, cols, depths;
+
double vx, vy, vz;
+
N_array_3d *x = x_comp;
+
N_array_3d *y = y_comp;
+
N_array_3d *z = z_comp;
+
N_gradient_3d grad;
Modified: grass/branches/develbranch_6/lib/gpde/N_gwflow.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_gwflow.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_gwflow.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -22,7 +22,7 @@
/* ***************** N_gwflow_data3d ***************************** */
/* *************************************************************** */
/*!
- * \brief Alllocate memory for the groundwater calculation data structure in 3 dimensions
+ * \brief Allocate memory for the groundwater calculation data structure in 3 dimensions
*
* The groundwater calculation data structure will be allocated including
* all appendant 3d and 2d arrays. The offset for the 3d arrays is one
@@ -264,14 +264,23 @@
int row, int depth)
{
double hc_e = 0, hc_w = 0, hc_n = 0, hc_s = 0, hc_t = 0, hc_b = 0;
+
double dx, dy, dz, Ax, Ay, Az;
+
double hc_x, hc_y, hc_z;
+
double hc_xw, hc_yn, hc_zt;
+
double hc_xe, hc_ys, hc_zb;
+
double hc_start;
+
double Ss, r, nf, q;
+
double C, W, E, N, S, T, B, V;
+
N_data_star *mat_pos;
+
N_gwflow_data3d *data;
/*cast the void pointer to the right data structure */
@@ -371,22 +380,39 @@
int row)
{
double T_e = 0, T_w = 0, T_n = 0, T_s = 0;
+
double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
+
double dx, dy, Az;
+
double hc_x, hc_y;
+
double z, top;
+
double hc_xw, hc_yn;
+
double z_xw, z_yn;
+
double hc_xe, hc_ys;
+
double z_xe, z_ys;
+
double hc, hc_start;
+
double Ss, r, nf, q;
+
double C, W, E, N, S, V;
+
N_gwflow_data2d *data;
+
N_data_star *mat_pos;
+
double river_vect = 0; /*entry in vector */
+
double river_mat = 0; /*entry in matrix */
+
double drain_vect = 0; /*entry in vector */
+
double drain_mat = 0; /*entry in matrix */
/*cast the void pointer to the right data structure */
Modified: grass/branches/develbranch_6/lib/gpde/N_gwflow.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_gwflow.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_gwflow.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -104,5 +104,6 @@
extern N_gwflow_data2d *N_alloc_gwflow_data2d(int cols, int rows, int river,
int drain);
extern void N_free_gwflow_data3d(N_gwflow_data3d * data);
+
extern void N_free_gwflow_data2d(N_gwflow_data2d * data);
#endif
Modified: grass/branches/develbranch_6/lib/gpde/N_heatflow.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_heatflow.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_heatflow.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -62,6 +62,8 @@
extern N_heatflow_data3d *N_alloc_heatflow_data3d(int depths, int rows,
int cols);
extern N_heatflow_data2d *N_alloc_heatflow_data2d(int rows, int cols);
+
extern void N_free_heatflow_data3d(N_heatflow_data3d * data);
+
extern void N_free_heatflow_data2d(N_heatflow_data2d * data);
#endif
Modified: grass/branches/develbranch_6/lib/gpde/N_les.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_les.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_les.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -17,30 +17,10 @@
*****************************************************************************/
#include "grass/N_pde.h"
+#include "grass/gmath.h"
#include <stdlib.h>
-/*!
- * \brief Allocate memory for a sparse vector
- *
- * \param cols int
- * \return N_spvector *
- *
- * */
-N_spvector *N_alloc_spvector(int cols)
-{
- N_spvector *spvector;
- G_debug(4, "Allocate memory for a sparse vector with %i cols\n", cols);
-
- spvector = (N_spvector *) G_calloc(1, sizeof(N_spvector));
-
- spvector->cols = cols;
- spvector->index = (int *)G_calloc(cols, sizeof(int));
- spvector->values = (double *)G_calloc(cols, sizeof(double));
-
- return spvector;
-}
-
/*!
* \brief Allocate memory for a (not) quadratic linear equation system which includes the Matrix A, vector x and vector b
*
@@ -198,6 +178,7 @@
N_les *N_alloc_les_param(int cols, int rows, int type, int parts)
{
N_les *les;
+
int i;
if (type == N_SPARSE_LES)
@@ -234,60 +215,19 @@
les->quad = 0;
if (type == N_SPARSE_LES) {
- les->Asp = (N_spvector **) G_calloc(rows, sizeof(N_spvector *));
+ les->Asp = G_math_alloc_spmatrix(rows);
les->type = N_SPARSE_LES;
}
else {
- les->A = (double **)G_calloc(rows, sizeof(double *));
- for (i = 0; i < rows; i++) {
- les->A[i] = (double *)G_calloc(cols, sizeof(double));
- }
+ les->A = G_alloc_matrix(rows, cols);
les->type = N_NORMAL_LES;
}
return les;
}
-
/*!
- * \brief Adds a sparse vector to a sparse linear equation system at position row
*
- * Return 1 for success and -1 for failure
- *
- * \param les N_les *
- * \param spvector N_spvector *
- * \param row int
- * \return int
- *
- * */
-int N_add_spvector_to_les(N_les * les, N_spvector * spvector, int row)
-{
-
-
- if (les != NULL) {
- if (les->type != N_SPARSE_LES)
- return -1;
-
- if (les->rows > row) {
- G_debug(5,
- "Add sparse vector %p to the sparse linear equation system at row %i\n",
- spvector, row);
- les->Asp[row] = spvector;
- }
- else
- return -1;
-
- }
- else {
- return -1;
- }
-
-
- return 1;
-}
-
-/*!
- *
* \brief prints the linear equation system to stdout
*
* <p>
@@ -354,29 +294,6 @@
}
/*!
- * \brief Release the memory of the sparse vector
- *
- * \param spvector N_spvector *
- * \return void
- *
- * */
-void N_free_spvector(N_spvector * spvector)
-{
- if (spvector) {
- if (spvector->values)
- G_free(spvector->values);
- if (spvector->index)
- G_free(spvector->index);
- G_free(spvector);
-
- spvector = NULL;
- }
-
- return;
-}
-
-
-/*!
* \brief Release the memory of the linear equation system
*
* \param les N_les *
@@ -403,21 +320,13 @@
if (les->type == N_SPARSE_LES) {
if (les->Asp) {
- for (i = 0; i < les->rows; i++)
- if (les->Asp[i])
- N_free_spvector(les->Asp[i]);
-
- G_free(les->Asp);
+ G_math_free_spmatrix(les->Asp, les->rows);
}
}
else {
if (les->A) {
- for (i = 0; i < les->rows; i++)
- if (les->A[i])
- G_free(les->A[i]);
-
- G_free(les->A);
+ G_free_matrix(les->A);
}
}
Modified: grass/branches/develbranch_6/lib/gpde/N_les_assemble.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_les_assemble.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_les_assemble.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -23,15 +23,17 @@
/* local protos */
static int make_les_entry_2d(int i, int j, int offset_i, int offset_j,
int count, int pos, N_les * les,
- N_spvector * spvect, N_array_2d * cell_count,
- N_array_2d * status, N_array_2d * start_val,
- double entry, int cell_type);
+ G_math_spvector * spvect,
+ N_array_2d * cell_count, N_array_2d * status,
+ N_array_2d * start_val, double entry,
+ int cell_type);
static int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
int offset_k, int count, int pos, N_les * les,
- N_spvector * spvect, N_array_3d * cell_count,
- N_array_3d * status, N_array_3d * start_val,
- double entry, int cell_type);
+ G_math_spvector * spvect,
+ N_array_3d * cell_count, N_array_3d * status,
+ N_array_3d * start_val, double entry,
+ int cell_type);
/* *************************************************************** *
* ********************** N_alloc_5star ************************** *
@@ -561,9 +563,13 @@
int cell_type)
{
int i, j, count = 0, pos = 0;
+
int cell_type_count = 0;
+
int **index_ij;
+
N_array_2d *cell_count;
+
N_les *les = NULL;
G_debug(2,
@@ -658,11 +664,11 @@
N_data_star *items = call->callback(data, geom, i, j);
/* we need a sparse vector pointer anytime */
- N_spvector *spvect = NULL;
+ G_math_spvector *spvect = NULL;
/*allocate a sprase vector */
if (les_type == N_SPARSE_LES) {
- spvect = N_alloc_spvector(items->count);
+ spvect = G_math_alloc_spvector(items->count);
}
/* initial conditions */
les->x[count] = N_get_array_2d_d_value(start_val, i, j);
@@ -737,7 +743,7 @@
/*How many entries in the les */
if (les->type == N_SPARSE_LES) {
spvect->cols = pos + 1;
- N_add_spvector_to_les(les, spvect, count);
+ G_math_add_spvector(les->Asp, spvect, count);
}
if (items)
@@ -785,9 +791,13 @@
N_array_2d * status, N_array_2d * start_val)
{
int rows, cols;
+
int count = 0;
+
int i, j, x, y, stat;
+
double *dvect1;
+
double *dvect2;
G_debug(2,
@@ -818,11 +828,11 @@
#pragma omp parallel default(shared)
{
- /*performe the matrix vector product */
+ /*performe the matrix vector product and */
if (les->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(les, dvect1, dvect2);
+ G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
else
- N_matrix_vector_product(les, dvect1, dvect2);
+ G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
#pragma omp for schedule (static) private(i)
for (i = 0; i < les->cols; i++)
les->b[i] = les->b[i] - dvect2[i];
@@ -876,12 +886,14 @@
/* **** make an entry in the les (2d) ***************************** */
/* **************************************************************** */
int make_les_entry_2d(int i, int j, int offset_i, int offset_j, int count,
- int pos, N_les * les, N_spvector * spvect,
+ int pos, N_les * les, G_math_spvector * spvect,
N_array_2d * cell_count, N_array_2d * status,
N_array_2d * start_val, double entry, int cell_type)
{
int K;
+
int di = offset_i;
+
int dj = offset_j;
K = N_get_array_2d_c_value(cell_count, i + di, j + dj) -
@@ -1013,9 +1025,13 @@
int cell_type)
{
int i, j, k, count = 0, pos = 0;
+
int cell_type_count = 0;
+
N_array_3d *cell_count;
+
N_les *les = NULL;
+
int **index_ij;
G_debug(2,
@@ -1121,11 +1137,11 @@
/*create the entries for the */
N_data_star *items = call->callback(data, geom, i, j, k);
- N_spvector *spvect = NULL;
+ G_math_spvector *spvect = NULL;
/*allocate a sprase vector */
if (les_type == N_SPARSE_LES)
- spvect = N_alloc_spvector(items->count);
+ spvect = G_math_alloc_spvector(items->count);
/* initial conditions */
les->x[count] = N_get_array_3d_d_value(start_val, i, j, k);
@@ -1191,7 +1207,7 @@
/*How many entries in the les */
if (les->type == N_SPARSE_LES) {
spvect->cols = pos + 1;
- N_add_spvector_to_les(les, spvect, count);
+ G_math_add_spvector(les->Asp, spvect, count);
}
if (items)
@@ -1238,9 +1254,13 @@
N_array_3d * status, N_array_3d * start_val)
{
int rows, cols, depths;
+
int count = 0;
+
int i, j, x, y, z, stat;
+
double *dvect1;
+
double *dvect2;
G_debug(2,
@@ -1277,9 +1297,9 @@
{
/*performe the matrix vector product and */
if (les->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(les, dvect1, dvect2);
+ G_math_Ax_sparse(les->Asp, dvect1, dvect2, les->rows);
else
- N_matrix_vector_product(les, dvect1, dvect2);
+ G_math_d_Ax(les->A, dvect1, dvect2, les->rows, les->cols);
#pragma omp for schedule (static) private(i)
for (i = 0; i < les->cols; i++)
les->b[i] = les->b[i] - dvect2[i];
@@ -1335,17 +1355,20 @@
/* **************************************************************** */
int make_les_entry_3d(int i, int j, int k, int offset_i, int offset_j,
int offset_k, int count, int pos, N_les * les,
- N_spvector * spvect, N_array_3d * cell_count,
+ G_math_spvector * spvect, N_array_3d * cell_count,
N_array_3d * status, N_array_3d * start_val,
double entry, int cell_type)
{
int K;
+
int di = offset_i;
+
int dj = offset_j;
+
int dk = offset_k;
- K = N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) -
- N_get_array_3d_d_value(cell_count, i, j, k);
+ K = (int)N_get_array_3d_d_value(cell_count, i + di, j + dj, k + dk) -
+ (int)N_get_array_3d_d_value(cell_count, i, j, k);
if (cell_type == N_CELL_ACTIVE) {
if ((int)N_get_array_3d_d_value(status, i + di, j + dj, k + dk) >
Deleted: grass/branches/develbranch_6/lib/gpde/N_les_pivot.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_les_pivot.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_les_pivot.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,92 +0,0 @@
-
-/*****************************************************************************
-*
-* MODULE: Grass PDE Numerical Library
-* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
-* soerengebbert <at> gmx <dot> de
-*
-* PURPOSE: linear equation system pivoting strategy
-* part of the gpde 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.
-*
-*****************************************************************************/
-
-#include <math.h>
-#include <unistd.h>
-#include <stdio.h>
-#include <string.h>
-#include "grass/N_pde.h"
-#include "solvers_local_proto.h"
-
-
-#define TINY 1.0e-20
-
-
-
-/*!
- * \brief Optimize the structure of the linear equation system with a common pivoting strategy
- *
- * Create a optimized linear equation system for
- * direct solvers: gauss and lu decomposition.
- *
- * The rows are permuted based on the pivot elements.
- *
- * This algorithm will modify the provided linear equation system
- * and should only be used with the gauss elimination and lu decomposition solver.
- *
- * \param les * N_les -- the linear equation system
- * \return int - the number of swapped rows
- *
- *
- * */
-int N_les_pivot_create(N_les * les)
-{
- int num = 0; /*number of changed rows */
- int i, j, k;
- double max;
- int number = 0;
- double tmpval = 0.0, s = 0.0;
- double *link = NULL;
-
- G_debug(2, "N_les_pivot_create: swap rows if needed");
- for (i = 0; i < les->rows; i++) {
- max = fabs(les->A[i][i]);
- number = i;
- for (j = i; j < les->rows; j++) {
- s = 0.0;
- for (k = i; k < les->rows; k++) {
- s += fabs(les->A[j][i]);
- }
- /*search for the pivot element */
- if (max < fabs(les->A[j][i]) / s) {
- max = fabs(les->A[j][i]);
- number = j;
- }
- }
- if (max == 0) {
- max = TINY;
- G_warning("Matrix is singular");
- }
- /*if an pivot element was found, swap the les entries */
- if (number != i) {
-
- G_debug(4, "swap row %i with row %i", i, number);
-
- tmpval = les->b[number];
- les->b[number] = les->b[i];
- les->b[i] = tmpval;
-
- link = les->A[number];
- les->A[number] = les->A[i];
- les->A[i] = link;
- num++;
- }
- }
-
- return num;
-}
Modified: grass/branches/develbranch_6/lib/gpde/N_pde.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_pde.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_pde.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -18,6 +18,7 @@
#include <grass/gis.h>
#include <grass/G3d.h>
#include <grass/glocale.h>
+#include <grass/gmath.h>
#ifndef _N_PDE_H_
#define _N_PDE_H_
@@ -74,17 +75,6 @@
/* *************************************************************** */
/*!
- * \brief The row vector of the sparse matrix
- * */
-typedef struct
-{
- int cols; /*Number of entries */
- double *values; /*The non null values of the row */
- int *index; /*the index number */
-} N_spvector;
-
-
-/*!
* \brief The linear equation system (les) structure
*
* This structure manages the Ax = b system.
@@ -98,26 +88,33 @@
double *x; /*the value vector */
double *b; /*the right side of Ax = b */
double **A; /*the normal quadratic matrix */
- N_spvector **Asp; /*the sparse matrix */
+ G_math_spvector **Asp; /*the sparse matrix */
int rows; /*number of rows */
int cols; /*number of cols */
int quad; /*is the matrix quadratic (1-quadratic, 0 not) */
int type; /*the type of the les, normal == 0, sparse == 1 */
} N_les;
-extern N_spvector *N_alloc_spvector(int cols);
extern N_les *N_alloc_les_param(int cols, int rows, int type, int param);
+
extern N_les *N_alloc_les(int rows, int type);
+
extern N_les *N_alloc_les_A(int rows, int type);
+
extern N_les *N_alloc_les_Ax(int rows, int type);
+
extern N_les *N_alloc_les_Ax_b(int rows, int type);
+
extern N_les *N_alloc_nquad_les(int cols, int rows, int type);
+
extern N_les *N_alloc_nquad_les_A(int cols, int rows, int type);
+
extern N_les *N_alloc_nquad_les_Ax(int cols, int rows, int type);
+
extern N_les *N_alloc_nquad_les_Ax_b(int cols, int rows, int type);
+
extern void N_print_les(N_les * les);
-extern int N_add_spvector_to_les(N_les * les, N_spvector * vector, int row);
-extern void N_free_spvector(N_spvector * vector);
+
extern void N_free_les(N_les * les);
/* *************************************************************** */
@@ -146,32 +143,16 @@
} N_geom_data;
extern N_geom_data *N_alloc_geom_data(void);
+
extern void N_free_geom_data(N_geom_data * geodata);
+
extern N_geom_data *N_init_geom_data_3d(G3D_Region * region3d,
N_geom_data * geodata);
extern N_geom_data *N_init_geom_data_2d(struct Cell_head *region,
N_geom_data * geodata);
extern double N_get_geom_data_area_of_cell(N_geom_data * geom, int row);
-
/* *************************************************************** */
-/* *************** LINEARE EQUATION SOLVER PART ****************** */
-/* *************************************************************** */
-extern int N_solver_gauss(N_les * les);
-extern int N_solver_lu(N_les * les);
-extern int N_solver_cholesky(N_les * les);
-extern int N_solver_jacobi(N_les * L, int maxit, double sor, double error);
-extern int N_solver_SOR(N_les * L, int maxit, double sor, double error);
-extern int N_solver_cg(N_les * les, int maxit, double error);
-extern int N_solver_pcg(N_les * les, int maxit, double error, int prec);
-extern int N_solver_bicgstab(N_les * les, int maxit, double error);
-extern void N_matrix_vector_product(N_les * les, double *source,
- double *result);
-extern void N_sparse_matrix_vector_product(N_les * les, double *source,
- double *result);
-extern N_les *N_create_diag_precond_matrix(N_les * les, int prec);
-
-/* *************************************************************** */
/* *************** READING RASTER AND VOLUME DATA **************** */
/* *************************************************************** */
@@ -187,13 +168,19 @@
} N_array_2d;
extern N_array_2d *N_alloc_array_2d(int cols, int rows, int offset, int type);
+
extern void N_free_array_2d(N_array_2d * data_array);
+
extern int N_get_array_2d_type(N_array_2d * array2d);
+
extern void N_get_array_2d_value(N_array_2d * array2d, int col, int row,
void *value);
extern CELL N_get_array_2d_c_value(N_array_2d * array2d, int col, int row);
+
extern FCELL N_get_array_2d_f_value(N_array_2d * array2d, int col, int row);
+
extern DCELL N_get_array_2d_d_value(N_array_2d * array2d, int col, int row);
+
extern void N_put_array_2d_value(N_array_2d * array2d, int col, int row,
char *value);
extern void N_put_array_2d_c_value(N_array_2d * array2d, int col, int row,
@@ -203,17 +190,25 @@
extern void N_put_array_2d_d_value(N_array_2d * array2d, int col, int row,
DCELL value);
extern int N_is_array_2d_value_null(N_array_2d * array2d, int col, int row);
+
extern void N_put_array_2d_value_null(N_array_2d * array2d, int col, int row);
+
extern void N_print_array_2d(N_array_2d * data);
+
extern void N_print_array_2d_info(N_array_2d * data);
+
extern void N_copy_array_2d(N_array_2d * source, N_array_2d * target);
+
extern double N_norm_array_2d(N_array_2d * array1, N_array_2d * array2,
int type);
extern N_array_2d *N_math_array_2d(N_array_2d * array1, N_array_2d * array2,
N_array_2d * result, int type);
extern int N_convert_array_2d_null_to_zero(N_array_2d * a);
+
extern N_array_2d *N_read_rast_to_array_2d(char *name, N_array_2d * array);
+
extern void N_write_array_2d_to_rast(N_array_2d * array, char *name);
+
extern void N_calc_array_2d_stats(N_array_2d * a, double *min, double *max,
double *sum, int *nonzero, int withoffset);
@@ -230,7 +225,9 @@
extern N_array_3d *N_alloc_array_3d(int cols, int rows, int depths,
int offset, int type);
extern void N_free_array_3d(N_array_3d * data_array);
+
extern int N_get_array_3d_type(N_array_3d * array3d);
+
extern void N_get_array_3d_value(N_array_3d * array3d, int col, int row,
int depth, void *value);
extern float N_get_array_3d_f_value(N_array_3d * array3d, int col, int row,
@@ -248,13 +245,17 @@
extern void N_put_array_3d_value_null(N_array_3d * array3d, int col, int row,
int depth);
extern void N_print_array_3d(N_array_3d * data);
+
extern void N_print_array_3d_info(N_array_3d * data);
+
extern void N_copy_array_3d(N_array_3d * source, N_array_3d * target);
+
extern double N_norm_array_3d(N_array_3d * array1, N_array_3d * array2,
int type);
extern N_array_3d *N_math_array_3d(N_array_3d * array1, N_array_3d * array2,
N_array_3d * result, int type);
extern int N_convert_array_3d_null_to_zero(N_array_3d * a);
+
extern N_array_3d *N_read_rast3d_to_array_3d(char *name, N_array_3d * array,
int mask);
extern void N_write_array_3d_to_rast3d(N_array_3d * array, char *name,
@@ -371,11 +372,17 @@
extern void N_set_les_callback_2d_func(N_les_callback_2d * data,
N_data_star * (*callback_func_2d) ());
extern N_les_callback_3d *N_alloc_les_callback_3d(void);
+
extern N_les_callback_2d *N_alloc_les_callback_2d(void);
+
extern N_data_star *N_alloc_5star(void);
+
extern N_data_star *N_alloc_7star(void);
+
extern N_data_star *N_alloc_9star(void);
+
extern N_data_star *N_alloc_27star(void);
+
extern N_data_star *N_create_5star(double C, double W, double E, double N,
double S, double V);
extern N_data_star *N_create_7star(double C, double W, double E, double N,
@@ -431,6 +438,7 @@
int cell_Type);
extern int N_les_pivot_create(N_les * les);
+
int N_les_integrate_dirichlet_2d(N_les * les, N_geom_data * geom,
N_array_2d * status, N_array_2d * start_val);
int N_les_integrate_dirichlet_3d(N_les * les, N_geom_data * geom,
@@ -459,12 +467,19 @@
/* *************************************************************** */
extern double N_calc_arith_mean(double a, double b);
+
extern double N_calc_arith_mean_n(double *a, int size);
+
extern double N_calc_geom_mean(double a, double b);
+
extern double N_calc_geom_mean_n(double *a, int size);
+
extern double N_calc_harmonic_mean(double a, double b);
+
extern double N_calc_harmonic_mean_n(double *a, int size);
+
extern double N_calc_quad_mean(double a, double b);
+
extern double N_calc_quad_mean_n(double *a, int size);
/* *************************************************************** */
@@ -472,6 +487,7 @@
/* *************************************************************** */
extern double N_full_upwinding(double sprod, double distance, double D);
+
extern double N_exp_upwinding(double sprod, double distance, double D);
@@ -644,25 +660,33 @@
extern N_gradient_2d *N_alloc_gradient_2d(void);
+
extern void N_free_gradient_2d(N_gradient_2d * grad);
+
extern N_gradient_2d *N_create_gradient_2d(double NC, double SC, double WC,
double EC);
extern int N_copy_gradient_2d(N_gradient_2d * source, N_gradient_2d * target);
+
extern N_gradient_2d *N_get_gradient_2d(N_gradient_field_2d * field,
N_gradient_2d * gradient, int col,
int row);
extern N_gradient_3d *N_alloc_gradient_3d(void);
+
extern void N_free_gradient_3d(N_gradient_3d * grad);
+
extern N_gradient_3d *N_create_gradient_3d(double NC, double SC, double WC,
double EC, double TC, double BC);
extern int N_copy_gradient_3d(N_gradient_3d * source, N_gradient_3d * target);
+
extern N_gradient_3d *N_get_gradient_3d(N_gradient_field_3d * field,
N_gradient_3d * gradient, int col,
int row, int depth);
extern N_gradient_neighbours_x *N_alloc_gradient_neighbours_x(void);
+
extern void N_free_gradient_neighbours_x(N_gradient_neighbours_x * grad);
+
extern N_gradient_neighbours_x *N_create_gradient_neighbours_x(double NWN,
double NEN,
double WC,
@@ -673,7 +697,9 @@
N_gradient_neighbours_x * target);
extern N_gradient_neighbours_y *N_alloc_gradient_neighbours_y(void);
+
extern void N_free_gradient_neighbours_y(N_gradient_neighbours_y * grad);
+
extern N_gradient_neighbours_y *N_create_gradient_neighbours_y(double NWW,
double NEE,
double NC,
@@ -684,7 +710,9 @@
N_gradient_neighbours_y * target);
extern N_gradient_neighbours_z *N_alloc_gradient_neighbours_z(void);
+
extern void N_free_gradient_neighbours_z(N_gradient_neighbours_z * grad);
+
extern N_gradient_neighbours_z *N_create_gradient_neighbours_z(double NWZ,
double NZ,
double NEZ,
@@ -698,7 +726,9 @@
N_gradient_neighbours_z * target);
extern N_gradient_neighbours_2d *N_alloc_gradient_neighbours_2d(void);
+
extern void N_free_gradient_neighbours_2d(N_gradient_neighbours_2d * grad);
+
extern N_gradient_neighbours_2d
* N_create_gradient_neighbours_2d(N_gradient_neighbours_x * x,
N_gradient_neighbours_y * y);
@@ -711,7 +741,9 @@
extern N_gradient_neighbours_3d *N_alloc_gradient_neighbours_3d(void);
+
extern void N_free_gradient_neighbours_3d(N_gradient_neighbours_3d * grad);
+
extern N_gradient_neighbours_3d
* N_create_gradient_neighbours_3d(N_gradient_neighbours_x * xt,
N_gradient_neighbours_x * xc,
@@ -725,11 +757,14 @@
N_gradient_neighbours_3d * target);
extern void N_print_gradient_field_2d_info(N_gradient_field_2d * field);
+
extern void N_calc_gradient_field_2d_stats(N_gradient_field_2d * field);
extern N_gradient_field_2d *N_alloc_gradient_field_2d(int cols, int rows);
+
extern void N_free_gradient_field_2d(N_gradient_field_2d * field);
+
extern int N_copy_gradient_field_2d(N_gradient_field_2d * source,
N_gradient_field_2d * target);
extern N_gradient_field_2d *N_compute_gradient_field_2d(N_array_2d * pot,
@@ -743,11 +778,13 @@
N_array_2d * y_comp);
extern void N_print_gradient_field_3d_info(N_gradient_field_3d * field);
+
extern void N_calc_gradient_field_3d_stats(N_gradient_field_3d * field);
extern N_gradient_field_3d *N_alloc_gradient_field_3d(int cols, int rows,
int depths);
extern void N_free_gradient_field_3d(N_gradient_field_3d * field);
+
extern int N_copy_gradient_field_3d(N_gradient_field_3d * source,
N_gradient_field_3d * target);
extern N_gradient_field_3d *N_compute_gradient_field_3d(N_array_3d * pot,
Modified: grass/branches/develbranch_6/lib/gpde/N_solute_transport.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_solute_transport.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_solute_transport.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -30,21 +30,35 @@
int row, int depth)
{
double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0, Df_t = 0, Df_b = 0;
+
double dx, dy, dz, Az;
+
double diff_x, diff_y, diff_z;
+
double diff_xw, diff_yn;
+
double diff_xe, diff_ys;
+
double diff_zt, diff_zb;
+
double cin = 0, cg, cg_start;
+
double R, nf, cs, q;
+
double C, W, E, N, S, T, B, V;
+
double vw = 0, ve = 0, vn = 0, vs = 0, vt = 0, vb = 0;
+
double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0, Ds_t = 0, Ds_b = 0;
+
double Dw = 0, De = 0, Dn = 0, Ds = 0, Dt = 0, Db = 0;
+
double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5, rt = 0.5, rb = 0.5;
N_solute_transport_data3d *data = NULL;
+
N_data_star *mat_pos;
+
N_gradient_3d grad;
/*cast the void pointer to the right data structure */
@@ -186,27 +200,47 @@
int row)
{
double Df_e = 0, Df_w = 0, Df_n = 0, Df_s = 0;
+
double z_e = 0, z_w = 0, z_n = 0, z_s = 0;
+
double dx, dy, Az;
+
double diff_x, diff_y;
+
double disp_x, disp_y;
+
double z;
+
double diff_xw, diff_yn;
+
double disp_xw, disp_yn;
+
double z_xw, z_yn;
+
double diff_xe, diff_ys;
+
double disp_xe, disp_ys;
+
double z_xe, z_ys;
+
double cin = 0, cg, cg_start;
+
double R, nf, cs, q;
+
double C, W, E, N, S, V, NE, NW, SW, SE;
+
double vw = 0, ve = 0, vn = 0, vs = 0;
+
double Ds_w = 0, Ds_e = 0, Ds_n = 0, Ds_s = 0;
+
double Dw = 0, De = 0, Dn = 0, Ds = 0;
+
double rw = 0.5, re = 0.5, rn = 0.5, rs = 0.5;
N_solute_transport_data2d *data = NULL;
+
N_data_star *mat_pos;
+
N_gradient_2d grad;
/*cast the void pointer to the right data structure */
@@ -578,8 +612,11 @@
void N_calc_solute_transport_transmission_2d(N_solute_transport_data2d * data)
{
int i, j, count = 1;
+
int cols, rows;
+
double c;
+
N_gradient_2d grad;
cols = data->grad->cols;
@@ -651,9 +688,13 @@
void N_calc_solute_transport_disptensor_2d(N_solute_transport_data2d * data)
{
int i, j;
+
int cols, rows;
+
double vx, vy, vv;
+
double disp_xx, disp_yy, disp_xy;
+
N_gradient_2d grad;
cols = data->grad->cols;
@@ -710,9 +751,13 @@
void N_calc_solute_transport_disptensor_3d(N_solute_transport_data3d * data)
{
int i, j, k;
+
int cols, rows, depths;
+
double vx, vy, vz, vv;
+
double disp_xx, disp_yy, disp_zz, disp_xy, disp_xz, disp_yz;
+
N_gradient_3d grad;
cols = data->grad->cols;
Modified: grass/branches/develbranch_6/lib/gpde/N_solute_transport.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_solute_transport.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_solute_transport.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -95,6 +95,7 @@
extern N_solute_transport_data2d *N_alloc_solute_transport_data2d(int cols,
int rows);
extern void N_free_solute_transport_data3d(N_solute_transport_data3d * data);
+
extern void N_free_solute_transport_data2d(N_solute_transport_data2d * data);
/*compute the dispersivity tensor */
Deleted: grass/branches/develbranch_6/lib/gpde/N_solvers.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_solvers.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_solvers.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,429 +0,0 @@
-
-/*****************************************************************************
-*
-* 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) 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 <math.h>
-#include <unistd.h>
-#include <stdio.h>
-#include <string.h>
-#include "grass/N_pde.h"
-#include "solvers_local_proto.h"
-
-/*prototypes */
-static void gauss_elimination(double **A, double *b, int rows);
-static void lu_decomposition(double **A, int rows);
-static int cholesky_decomposition(double **A, int rows);
-static void backward_solving(double **A, double *x, double *b, int rows);
-static void forward_solving(double **A, double *x, double *b, int rows);
-
-/***********************************************************
- * GAUSS elimination solver for Ax = b *********************
- * ********************************************************/
-/*!
- * \brief The gauss elimination solver for quardatic matrices
- *
- * This solver does not support sparse matrices
- * The matrix A will be overwritten.
- * The result is written to the vector x in the N_les structure
- *
- * \param les N_les *
- * \return int -- 1 success, 0 solver does not work with sprase matrices, -1 matrix is not quadratic, -2 unable to solve the les
- * */
-int N_solver_gauss(N_les * les)
-{
-
- if (les->type != N_NORMAL_LES) {
- G_warning(_("The gauss elimination solver does not work with sparse matrices"));
- return 0;
- }
-
- if (les->quad != 1) {
- G_fatal_error(_("The linear equation system is not quadratic"));
- return 0;
- }
-
-
- G_message(_("Starting direct gauss elimination solver"));
-
- N_les_pivot_create(les);
- gauss_elimination(les->A, les->b, les->rows);
- backward_solving(les->A, les->x, les->b, les->rows);
-
- return 1;
-}
-
-/***********************************************************
- * LU solver for Ax = b ************************************
- * ********************************************************/
-/*!
- * \brief The LU solver for quardatic matrices
- *
- * This solver does not support sparse matrices
- * The matrix A will be overwritten.
- * The result is written to the vector x in the N_les structure
- *
- * \param les N_les *
- * \return int -- 1 success, 0 solver does not work with sprase matrices, -1 matrix is not quadratic, -2 unable to solve the les
- * */
-int N_solver_lu(N_les * les)
-{
- int i;
- double *c, *tmpv;
-
- if (les->type != N_NORMAL_LES) {
- G_warning(_("The lu solver does not work with sparse matrices"));
- return 0;
- }
-
- if (les->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
-
- G_message(_("Starting direct lu decomposition solver"));
-
- tmpv = vectmem(les->rows);
- c = vectmem(les->rows);
-
- N_les_pivot_create(les);
- lu_decomposition(les->A, les->rows);
-
-#pragma omp parallel
- {
-
-#pragma omp for schedule (static) private(i)
- for (i = 0; i < les->rows; i++) {
- tmpv[i] = les->A[i][i];
- les->A[i][i] = 1;
- }
-
-#pragma omp single
- {
- forward_solving(les->A, les->b, les->b, les->rows);
- }
-
-#pragma omp for schedule (static) private(i)
- for (i = 0; i < les->rows; i++) {
- les->A[i][i] = tmpv[i];
- }
-
-#pragma omp single
- {
- backward_solving(les->A, les->x, les->b, les->rows);
- }
- }
-
- G_free(c);
- G_free(tmpv);
-
-
- return 1;
-}
-
-/***********************************************************
- * cholesky solver for Ax = b ******************************
- * ********************************************************/
-/*!
- * \brief The choleksy decomposition solver for quardatic, symmetric
- * positiv definite matrices
- *
- * This solver does not support sparse matrices
- * The matrix A will be overwritten.
- * The result is written to the vector x in the N_les structure
- *
- * \param les N_les *
- * \return int -- 1 success, 0 solver does not work with sprase matrices, -1 matrix is not quadratic, -2 unable to solve the les, -3 matrix is not symmetric
- * */
-int N_solver_cholesky(N_les * les)
-{
- if (les->type != N_NORMAL_LES) {
- G_warning(_("The cholesky solver does not work with sparse matrices"));
- return 0;
- }
-
- if (les->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
- /* check for symmetry */
- if (check_symmetry(les) != 1) {
- G_warning(_("Matrix is not symmetric!"));
- return -3;
- }
-
- G_message(_("Starting cholesky decomposition solver"));
-
- if (cholesky_decomposition(les->A, les->rows) != 1) {
- G_warning(_("Unable to solve the linear equation system"));
- return -2;
- }
-
- forward_solving(les->A, les->b, les->b, les->rows);
- backward_solving(les->A, les->x, les->b, les->rows);
-
- return 1;
-}
-
-
-/***********************************************************
- * gauss elimination ***************************************
- * ********************************************************/
-/*!
- * \brief Gauss elimination
- *
- * The matrix will be overwritten with the decomposite form
- *
- * \param A double **
- * \param b double *
- * \param rows int
- * \return void
- *
- * */
-void gauss_elimination(double **A, double *b, int rows)
-{
- int i, j, k;
- double tmpval = 0.0;
-
- for (k = 0; k < rows - 1; k++) {
-#pragma omp parallel for schedule (static) private(i, j, tmpval) shared(k, A, b, rows)
- for (i = k + 1; i < rows; i++) {
- tmpval = A[i][k] / A[k][k];
- b[i] = b[i] - tmpval * b[k];
- for (j = k + 1; j < rows; j++) {
- A[i][j] = A[i][j] - tmpval * A[k][j];
- }
- }
- }
-
- return;
-}
-
-/***********************************************************
- * lu decomposition ****************************************
- * ********************************************************/
-/*!
- * \brief lu decomposition
- *
- * The matrix will be overwritten with the decomposite form
- *
- * \param A double **
- * \param rows int
- * \return void
- *
- * */
-void lu_decomposition(double **A, int rows)
-{
-
- int i, j, k;
-
- for (k = 0; k < rows - 1; k++) {
-#pragma omp parallel for schedule (static) private(i, j) shared(k, A, rows)
- for (i = k + 1; i < rows; i++) {
- A[i][k] = A[i][k] / A[k][k];
- for (j = k + 1; j < rows; j++) {
- A[i][j] = A[i][j] - A[i][k] * A[k][j];
- }
- }
- }
-
- return;
-}
-
-/***********************************************************
- * cholesky decomposition **********************************
- * ********************************************************/
-/*!
- * \brief cholesky decomposition for symmetric, positiv definite matrices
- *
- * The provided matrix will be overwritten with the lower and
- * upper triangle matrix A = LL^T
- *
- * \param A double **
- * \param rows int
- * \return void
- *
- * */
-int cholesky_decomposition(double **A, int rows)
-{
-
- int i, j, k;
- double sum_1 = 0.0;
- double sum_2 = 0.0;
- int error = 0;
-
-
- for (k = 0; k < rows; k++) {
-#pragma omp parallel private(i, j, sum_2) shared(A, rows, sum_1)
- {
-#pragma omp for schedule (static) private(j) reduction(+:sum_1)
- for (j = 0; j < k; j++) {
- sum_1 += A[k][j] * A[k][j];
- }
-#pragma omp single
- {
- if ((A[k][k] - sum_1) < 0) {
- error++; /*not allowed to exit the OpenMP region */
- }
- A[k][k] = sqrt(A[k][k] - sum_1);
- sum_1 = 0.0;
- }
-#pragma omp for schedule (static) private(i, j, sum_2)
- for (i = k + 1; i < rows; i++) {
- sum_2 = 0.0;
- for (j = 0; j < k; j++) {
- sum_2 += A[i][j] * A[k][j];
- }
- A[i][k] = (A[i][k] - sum_2) / A[k][k];
- }
- }
- }
-
- /*we need to copy the lower triangle matrix to the upper trianle */
-#pragma omp parallel for schedule (static) private(i, k) shared(A, rows)
- for (k = 0; k < rows; k++) {
- for (i = k + 1; i < rows; i++) {
- A[k][i] = A[i][k];
- }
- }
-
- if (error > 0) {
- G_warning("Matrix is not positive definite");
- return -1;
- }
-
- return 1;
-}
-
-
-/***********************************************************
- * backward solving ****************************************
- * ********************************************************/
-/*!
- * \brief backward solving
- *
- * \param A double **
- * \param x double *
- * \param b double *
- * \param rows int
- * \return void
- *
- * */
-void backward_solving(double **A, double *x, double *b, int rows)
-{
- int i, j;
- double tmpval = 0.0;
-
- for (i = rows - 1; i >= 0; i--) {
- tmpval = 0;
- for (j = i + 1; j < rows; j++) {
- /*tmpval += A[i][j] * x[j]; */
- b[i] = b[i] - A[i][j] * x[j];
- }
- /*x[i] = (b[i] - tmpval) / A[i][i]; */
- x[i] = (b[i]) / A[i][i];
- }
-
- return;
-}
-
-
-/***********************************************************
- * forward solving *****************************************
- * ********************************************************/
-/*!
- * \brief forward solving
- *
- * \param A double **
- * \param x double *
- * \param b double *
- * \param rows int
- * \return void
- *
- * */
-void forward_solving(double **A, double *x, double *b, int rows)
-{
- int i, j;
- double tmpval = 0.0;
-
- for (i = 0; i < rows; i++) {
- tmpval = 0;
- for (j = 0; j < i; j++) {
- tmpval += A[i][j] * x[j];
- }
- x[i] = (b[i] - tmpval) / A[i][i];
- }
-
- return;
-}
-
-
-/* ******************************************************* *
- * ***** solving a tridiagonal eq system ***************** *
- * ******************************************************* */
-void thomalg(double **M, double *V, int rows)
-{
- double *Vtmp;
- double *g;
- double b;
- int i;
-
- Vtmp = vectmem(rows);
- g = vectmem(rows);
-
- for (i = 0; i < rows; i++) {
- if (i == 0) {
- b = M[i][i];
- Vtmp[i] = V[i] / b;
- }
- else {
- b = M[i][i] - M[i][i - 1] * g[i - 1];
- Vtmp[i] = (V[i] - Vtmp[i - 1] * M[i][i - 1]) / b;
- }
- if (i < rows - 1) {
- g[i] = M[i][i + 1] / b;
- }
- }
-
- V[rows - 1] = Vtmp[rows - 1];
- for (i = rows - 2; i >= 0; i--) {
- V[i] = Vtmp[i] - g[i] * V[i + 1];
- }
-
- G_free(Vtmp);
- G_free(g);
-}
-
-
-/***********************************************************
- * vectmem *************************************************
- * ********************************************************/
-/*!
- * \brief Allocate vector memory
- *
- * \param rows int
- * \return double *
- *
- * */
-double *vectmem(int rows)
-{
- double *vector;
-
- vector = (double *)(G_calloc(rows, (sizeof(double))));
- return vector;
-}
Deleted: grass/branches/develbranch_6/lib/gpde/N_solvers_classic_iter.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_solvers_classic_iter.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_solvers_classic_iter.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,241 +0,0 @@
-
-/*****************************************************************************
-*
-* MODULE: Grass PDE Numerical Library
-* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
-* soerengebbert <at> gmx <dot> de
-*
-* PURPOSE: linear equation system solvers
-* part of the gpde 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.
-*
-*****************************************************************************/
-
-#include <math.h>
-#include <unistd.h>
-#include <stdio.h>
-#include <string.h>
-#include "grass/N_pde.h"
-#include "solvers_local_proto.h"
-
-static int sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
- const char *type);
-static int jacobi(double **M, double *b, double *x, int rows, int maxit,
- double sor, double error);
-static int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
- double sor, double error);
-
-/* ******************************************************* *
- * ******** overrelaxed jacobian ************************* *
- * ******************************************************* */
-/*!
- * \brief The iterative jacobian solver for regular matrices
- *
- * The result is written to the vector L->x of the les.
- * This iterative solver works with sparse matrices and regular quadratic matrices.
- *
- * 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 L->x.
- * The parameter <i>err</i> defines the error break criteria for the solver.
- *
- * \param L N_les * -- the linear equatuin system
- * \param maxit int -- the maximum number of iterations
- * \param sor double -- defines the successive overrelaxion parameter [0:1]
- * \param error double -- defines the error break criteria
- * \return int -- 1=success, -1=could not solve the les
- *
- * */
-int N_solver_jacobi(N_les * L, int maxit, double sor, double error)
-{
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
- if (L->type == N_NORMAL_LES)
- return jacobi(L->A, L->b, L->x, L->rows, maxit, sor, error);
- else
- return sparse_jacobi_gauss(L, maxit, sor, error,
- N_SOLVER_ITERATIVE_JACOBI);
-}
-
-
-/* ******************************************************* *
- * ********* overrelaxed gauss seidel ******************** *
- * ******************************************************* */
-/*!
- * \brief The iterative overrelaxed gauss seidel solver for regular matrices
- *
- * The result is written to the vector L->x of the les.
- * This iterative solver works with sparse matrices and regular quadratic matrices.
- *
- * 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 L->x.
- * The parameter <i>err</i> defines the error break criteria for the solver.
- *
- * \param L N_les * -- the linear equatuin system
- * \param maxit int -- the maximum number of iterations
- * \param sor double -- defines the successive overrelaxion parameter [0:1]
- * \param error double -- defines the error break criteria
- * \return int -- 1=success, -1=could not solve the les
- *
- * */
-
-int N_solver_SOR(N_les * L, int maxit, double sor, double error)
-{
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
- if (L->type == N_NORMAL_LES)
- return gauss_seidel(L->A, L->b, L->x, L->rows, maxit, sor, error);
- else
- return sparse_jacobi_gauss(L, maxit, sor, error,
- N_SOLVER_ITERATIVE_SOR);
-}
-
-/* ******************************************************* *
- * ****** sparse jacobi and SOR algorithm **************** *
- * ******************************************************* */
-int
-sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
- const char *type)
-{
- int i, j, k, rows, finished = 0;
- double *Enew, *x, *b;
- double E, err = 0;
-
- x = L->x;
- b = L->b;
- rows = L->rows;
-
- Enew = vectmem(rows);
-
- for (k = 0; k < maxit; k++) {
- err = 0;
- {
- if (k == 0) {
- for (j = 0; j < rows; j++) {
- Enew[j] = x[j];
- }
- }
- for (i = 0; i < rows; i++) {
- E = 0;
- if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0) {
- for (j = 0; j < L->Asp[i]->cols; j++) {
- E += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]];
- }
- }
- else {
- for (j = 0; j < L->Asp[i]->cols; j++) {
- E += L->Asp[i]->values[j] * Enew[L->Asp[i]->index[j]];
- }
- }
- Enew[i] = x[i] - sor * (E - b[i]) / L->Asp[i]->values[0];
- }
- for (j = 0; j < rows; j++) {
- err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
-
- x[j] = Enew[j];
- }
- }
-
- if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0)
- G_message(_("sparse Jacobi -- iteration %5i error %g\n"), k, err);
- else if (strcmp(type, N_SOLVER_ITERATIVE_SOR) == 0)
- G_message(_("sparse SOR -- iteration %5i error %g\n"), k, err);
-
- if (err < error) {
- finished = 1;
- break;
- }
- }
-
- G_free(Enew);
-
- return finished;
-}
-
-/* ******************************************************* *
- * ******* direct jacobi ********************************* *
- * ******************************************************* */
-int jacobi(double **M, double *b, double *x, int rows, int maxit, double sor,
- double error)
-{
- int i, j, k;
- double *Enew;
- double E, err = 0;
-
- Enew = vectmem(rows);
-
- for (j = 0; j < rows; j++) {
- Enew[j] = x[j];
- }
-
- for (k = 0; k < maxit; k++) {
- for (i = 0; i < rows; i++) {
- E = 0;
- for (j = 0; j < rows; j++) {
- E += M[i][j] * x[j];
- }
- Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
- }
- err = 0;
- for (j = 0; j < rows; j++) {
- err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
-
- x[j] = Enew[j];
- }
- G_message(_("Jacobi -- iteration %5i error %g\n"), k, err);
- if (err < error)
- break;
- }
-
- return 1;
-}
-
-/* ******************************************************* *
- * ******* direct gauss seidel *************************** *
- * ******************************************************* */
-int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
- double sor, double error)
-{
- int i, j, k;
- double *Enew;
- double E, err = 0;
-
- Enew = vectmem(rows);
-
- for (j = 0; j < rows; j++) {
- Enew[j] = x[j];
- }
-
- for (k = 0; k < maxit; k++) {
- for (i = 0; i < rows; i++) {
- E = 0;
- for (j = 0; j < rows; j++) {
- E += M[i][j] * Enew[j];
- }
- Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
- }
- err = 0;
- for (j = 0; j < rows; j++) {
- err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
-
- x[j] = Enew[j];
- }
- G_message(_("SOR -- iteration %5i error %g\n"), k, err);
- if (err < error)
- break;
- }
-
- return 1;
-}
Deleted: grass/branches/develbranch_6/lib/gpde/N_solvers_krylov.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_solvers_krylov.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_solvers_krylov.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,982 +0,0 @@
-
-/*****************************************************************************
-*
-* MODULE: Grass PDE Numerical Library
-* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
-* soerengebbert <at> gmx <dot> de
-*
-* PURPOSE: linear equation system solvers
-* part of the gpde 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.
-*
-*****************************************************************************/
-
-#include <math.h>
-#include <unistd.h>
-#include <stdio.h>
-#include <string.h>
-#include "grass/N_pde.h"
-#include "solvers_local_proto.h"
-
-/*local protos */
-static void scalar_product(double *a, double *b, double *scalar, int rows);
-static void sub_vectors(double *source_a, double *source_b, double *result,
- int row);
-static void sub_vectors_scalar(double *source_a, double *source_b,
- double *result, double scalar_b, int rows);
-static void add_vectors(double *source_a, double *source_b, double *result,
- int rows);
-static void add_vectors_scalar(double *source_a, double *source_b,
- double *result, double scalar_b, int rows);
-static void add_vectors_scalar2(double *source_a, double *source_b,
- double *result, double scalar_a,
- double scalar_b, int rows);
-static void scalar_vector_product(double *a, double *result, double scalar,
- int rows);
-static void sync_vectors(double *source, double *target, int rows);
-
-
-/* ******************************************************* *
- * *** preconditioned conjugate gradients **************** *
- * ******************************************************* */
-/*!
- * \brief The iterative preconditioned conjugate gradients solver for symmetric positive definite matrices
- *
- * This iterative solver works with symmetric positive definite sparse matrices and
- * symmetric positive definite regular quadratic matrices.
- *
- * 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 L->x.
- * The parameter <i>err</i> defines the error break criteria for the solver.
- *
- * \param L N_les * -- the linear equatuin system
- * \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 N_DIAGONAL_PRECONDITION, N_ROWSUM_PRECONDITION
- * \return int -- 1 - success, 2 - not finisehd but success, 0 - matrix singular, -1 - could not solve the les
- *
- * */
-int N_solver_pcg(N_les * L, int maxit, double err, int prec)
-{
- double *r, *z;
- double *p;
- double *v;
- double *x, *b;
- double s = 0.0;
- double a0 = 0, a1 = 0, mygamma, tmp = 0;
- int m, rows, i;
- int finished = 2;
- int error_break;
- N_les *M;
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
- /* check for symmetry */
- if (check_symmetry(L) != 1) {
- G_warning(_("Matrix is not symmetric!"));
- }
-
- x = L->x;
- b = L->b;
- rows = L->rows;
-
- r = vectmem(rows);
- p = vectmem(rows);
- v = vectmem(rows);
- z = vectmem(rows);
-
- error_break = 0;
-
- /*compute the preconditioning matrix */
- M = N_create_diag_precond_matrix(L, prec);
-
- /*
- * residual calculation
- */
-#pragma omp parallel
- {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, x, v);
- else
- N_matrix_vector_product(L, x, v);
-
- sub_vectors(b, v, r, rows);
- /*performe the preconditioning */
- N_sparse_matrix_vector_product(M, r, p);
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += p[i] * r[i];
- }
- }
-
- a0 = s;
- s = 0.0;
-
- /* ******************* */
- /* start the iteration */
- /* ******************* */
- for (m = 0; m < maxit; m++) {
-#pragma omp parallel default(shared)
- {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, p, v);
- else
- N_matrix_vector_product(L, p, v);
-
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += v[i] * p[i];
- }
-
- /* barrier */
-#pragma omp single
- {
- tmp = s;
- mygamma = a0 / tmp;
- s = 0.0;
- }
-
- add_vectors_scalar(x, p, x, mygamma, rows);
-
- if (m % 50 == 1) {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, x, v);
- else
- N_matrix_vector_product(L, x, v);
-
- sub_vectors(b, v, r, rows);
-
- }
- else {
- sub_vectors_scalar(r, v, r, mygamma, rows);
- }
-
- /*performe the preconditioning */
- N_sparse_matrix_vector_product(M, r, z);
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += z[i] * r[i];
- }
-
- /* barrier */
-#pragma omp single
- {
- a1 = s;
- tmp = a1 / a0;
- a0 = a1;
- s = 0.0;
-
- if (a1 < 0 || a1 == 0 || a1 > 0) {
- ;
- }
- else {
- G_warning(_("Unable to solve the linear equation system"));
- error_break = 1;
- }
- }
- add_vectors_scalar(z, p, p, tmp, rows);
- }
-
- if (L->type == N_SPARSE_LES)
- G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0);
- else
- G_message(_("PCG -- iteration %i error %g\n"), m, a0);
-
- if (error_break == 1) {
- finished = -1;
- break;
- }
-
-
- if (a0 < err) {
- finished = 1;
- break;
- }
- }
-
- G_free(r);
- G_free(p);
- G_free(v);
- G_free(z);
-
- return finished;
-}
-
-
-/* ******************************************************* *
- * ****************** conjugate gradients **************** *
- * ******************************************************* */
-/*!
- * \brief The iterative conjugate gradients solver for symmetric positive definite matrices
- *
- * This iterative solver works with symmetric positive definite sparse matrices and
- * symmetric positive definite regular quadratic matrices.
- *
- * 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 L->x.
- * The parameter <i>err</i> defines the error break criteria for the solver.
- *
- * \param L N_les * -- the linear equatuin system
- * \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 N_solver_cg(N_les * L, int maxit, double err)
-{
- double *r;
- double *p;
- double *v;
- double *x, *b;
- double s = 0.0;
- double a0 = 0, a1 = 0, mygamma, tmp = 0;
- int m, rows, i;
- int finished = 2;
- int error_break;
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
- /* check for symmetry */
- if (check_symmetry(L) != 1) {
- G_warning(_("Matrix is not symmetric!"));
- }
-
- x = L->x;
- b = L->b;
- rows = L->rows;
-
- r = vectmem(rows);
- p = vectmem(rows);
- v = vectmem(rows);
-
- error_break = 0;
- /*
- * residual calculation
- */
-#pragma omp parallel
- {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, x, v);
- else
- N_matrix_vector_product(L, x, v);
-
- sub_vectors(b, v, r, rows);
- sync_vectors(r, p, rows);
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += r[i] * r[i];
- }
- }
-
- a0 = s;
- s = 0.0;
-
- /* ******************* */
- /* start the iteration */
- /* ******************* */
- for (m = 0; m < maxit; m++) {
-#pragma omp parallel default(shared)
- {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, p, v);
- else
- N_matrix_vector_product(L, p, v);
-
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += v[i] * p[i];
- }
-
- /* barrier */
-#pragma omp single
- {
- tmp = s;
- mygamma = a0 / tmp;
- s = 0.0;
- }
-
- add_vectors_scalar(x, p, x, mygamma, rows);
-
- if (m % 50 == 1) {
- /* matrix vector multiplication */
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, x, v);
- else
- N_matrix_vector_product(L, x, v);
-
- sub_vectors(b, v, r, rows);
- }
- else {
- sub_vectors_scalar(r, v, r, mygamma, rows);
- }
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += r[i] * r[i];
- }
-
- /* barrier */
-#pragma omp single
- {
- a1 = s;
- tmp = a1 / a0;
- a0 = a1;
- s = 0.0;
-
- if (a1 < 0 || a1 == 0 || a1 > 0) {
- ;
- }
- else {
- G_warning(_("Unable to solve the linear equation system"));
- error_break = 1;
- }
- }
- add_vectors_scalar(r, p, p, tmp, rows);
- }
-
- if (L->type == N_SPARSE_LES)
- G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0);
- else
- G_message(_("CG -- iteration %i error %g\n"), m, a0);
-
- if (error_break == 1) {
- finished = -1;
- break;
- }
-
- if (a0 < err) {
- finished = 1;
- break;
- }
- }
-
- G_free(r);
- G_free(p);
- G_free(v);
-
- return finished;
-}
-
-/* ******************************************************* *
- * ************ biconjugate gradients ******************** *
- * ******************************************************* */
-/*!
- * \brief The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices
- *
- * The result is written to the vector L->x of the les.
- * This iterative solver works with sparse matrices and regular quadratic matrices.
- *
- * 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 L->x.
- * The parameter <i>err</i> defines the error break criteria for the solver.
- *
- * \param L N_les * -- the linear equatuin system
- * \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 N_solver_bicgstab(N_les * L, int maxit, double err)
-{
- double *r;
- double *r0;
- double *p;
- double *v;
- double *s;
- double *t;
- double *x, *b;
- double s1 = 0.0, s2 = 0.0, s3 = 0.0;
- double alpha = 0, beta = 0, omega, rr0 = 0, error;
- int m, rows, i;
- int finished = 2;
- int error_break;
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return -1;
- }
-
-
- x = L->x;
- b = L->b;
- rows = L->rows;
- r = vectmem(rows);
- r0 = vectmem(rows);
- p = vectmem(rows);
- v = vectmem(rows);
- s = vectmem(rows);
- t = vectmem(rows);
-
- error_break = 0;
-
-#pragma omp parallel
- {
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, x, v);
- else
- N_matrix_vector_product(L, x, v);
- sub_vectors(b, v, r, rows);
- sync_vectors(r, r0, rows);
- sync_vectors(r, p, rows);
- }
-
- s1 = s2 = s3 = 0.0;
-
- /* ******************* */
- /* start the iteration */
- /* ******************* */
- for (m = 0; m < maxit; m++) {
-
-#pragma omp parallel default(shared)
- {
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, p, v);
- else
- N_matrix_vector_product(L, p, v);
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3)
- for (i = 0; i < rows; i++) {
- s1 += r[i] * r[i];
- s2 += r[i] * r0[i];
- s3 += v[i] * r0[i];
- }
-
-#pragma omp single
- {
- error = s1;
-
- if (error < 0 || error == 0 || error > 0) {
- ;
- }
- else {
- G_warning(_("Unable to solve the linear equation system"));
- error_break = 1;
- }
-
- rr0 = s2;
- alpha = rr0 / s3;
- s1 = s2 = s3 = 0.0;
- }
-
- sub_vectors_scalar(r, v, s, alpha, rows);
- if (L->type == N_SPARSE_LES)
- N_sparse_matrix_vector_product(L, s, t);
- else
- N_matrix_vector_product(L, s, t);
-
- /* scalar product */
-#pragma omp for schedule (static) private(i) reduction(+:s1, s2)
- for (i = 0; i < rows; i++) {
- s1 += t[i] * s[i];
- s2 += t[i] * t[i];
- }
-
-#pragma omp single
- {
- omega = s1 / s2;
- s1 = s2 = 0.0;
- }
-
- add_vectors_scalar2(p, s, r, alpha, omega, rows);
- add_vectors(x, r, x, rows);
- sub_vectors_scalar(s, t, r, omega, rows);
-
-#pragma omp for schedule (static) private(i) reduction(+:s1)
- for (i = 0; i < rows; i++) {
- s1 += r[i] * r0[i];
- }
-
-#pragma omp single
- {
- beta = alpha / omega * s1 / rr0;
- s1 = s2 = s3 = 0.0;
- }
-
- sub_vectors_scalar(p, v, p, omega, rows);
- add_vectors_scalar(r, p, p, beta, rows);
- }
-
-
- if (L->type == N_SPARSE_LES)
- G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m,
- error);
- else
- G_message(_("BiCGStab -- iteration %i error %g\n"), m, error);
-
- if (error_break == 1) {
- finished = -1;
- break;
- }
-
- if (error < err) {
- finished = 1;
- break;
- }
- }
-
- G_free(r);
- G_free(r0);
- G_free(p);
- G_free(v);
- G_free(s);
- G_free(t);
-
- return finished;
-}
-
-/*!
- * \brief Calculates the scalar product of vector a and b
- *
- * The result is written to variable scalar
- *
- * \param a double *
- * \param b double *
- * \param scalar double *
- * \param rows int
- * \return void
- *
- * */
-void scalar_product(double *a, double *b, double *scalar, int rows)
-{
- int i;
- double s = 0.0;
-
-#pragma omp parallel for schedule (static) reduction(+:s)
- for (i = 0; i < rows; i++) {
- s += a[i] * b[i];
- }
-
- *scalar = s;
- return;
-}
-
-/*!
- * \brief Calculates the matrix - vector product of matrix L->A and vector x
- *
- * The result is written to vector named result. This function only works with
- * regular quadratic matrices.
- *
- * \param L N_les *
- * \param x double *
- * \param result double *
- * \return void
- *
- * */
-void N_matrix_vector_product(N_les * L, double *x, double *result)
-{
- int i, j;
- double tmp;
-
-#pragma omp for schedule (static) private(i, j, tmp)
- for (i = 0; i < L->rows; i++) {
- tmp = 0;
- for (j = 0; j < L->cols; j++) {
- tmp += L->A[i][j] * x[j];
- }
- result[i] = tmp;
- }
- return;
-}
-
-/*!
- * \brief Calculates the matrix - vector product of sparse matrix L->Asp and vector x
- *
- * The result is written to vector named result. This function only works with
- * sparse matrices matrices.
- *
- * \param L N_les *
- * \param x double *
- * \param result double *
- * \return void
- *
- * */
-void N_sparse_matrix_vector_product(N_les * L, double *x, double *result)
-{
- int i, j;
- double tmp;
-
-#pragma omp for schedule (static) private(i, j, tmp)
- for (i = 0; i < L->rows; i++) {
- tmp = 0;
- for (j = 0; j < L->Asp[i]->cols; j++) {
- tmp += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]];
- }
- result[i] = tmp;
- }
- return;
-}
-
-/*!
- * \brief Multipiles the vector a and b with the scalars scalar_a and scalar_b and adds them
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param b double *
- * \param result double *
- * \param scalar_a double
- * \param scalar_b double
- * \param rows int
- *
- * */
-void
-add_vectors_scalar2(double *a, double *b, double *result, double scalar_a,
- double scalar_b, int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = scalar_a * a[i] + scalar_b * b[i];
- }
-
- return;
-}
-
-/*!
- * \brief Multipiles the vector b with the scalar scalar_b and adds a to b
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param b double *
- * \param result double *
- * \param scalar_b double
- * \param rows int
- *
- * */
-void
-add_vectors_scalar(double *a, double *b, double *result, double scalar_b,
- int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = a[i] + scalar_b * b[i];
- }
-
- return;
-}
-
-/*!
- * \brief Multipiles the vector b with the scalar scalar_b and substracts b from a
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param b double *
- * \param result double *
- * \param scalar_b double
- * \param rows int
- *
- * */
-void
-sub_vectors_scalar(double *a, double *b, double *result, double scalar_b,
- int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = a[i] - scalar_b * b[i];
- }
-
- return;
-}
-
-/*!
- * \brief Adds a to b
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param b double *
- * \param result double *
- * \param rows int
- *
- * */
-void add_vectors(double *a, double *b, double *result, int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = a[i] + b[i];
- }
-
- return;
-}
-
-/*!
- * \brief Substracts b from a
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param b double *
- * \param result double *
- * \param rows int
- *
- * */
-void sub_vectors(double *a, double *b, double *result, int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = a[i] - b[i];
- }
-
- return;
-}
-
-/*!
- * \brief Multipiles the vector a with the scalar named scalar
- *
- *
- * The result is written to the vector named result.
- *
- * \param a double *
- * \param result double *
- * \param scalar double *
- * \param rows int
- *
- * */
-void scalar_vector_product(double *a, double *result, double scalar, int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- result[i] = scalar * a[i];
- }
-
- return;
-}
-
-/*!
- * \brief Copies the source vector to the target vector
- *
- * \param source double *
- * \param target double *
- * \param rows int
- *
- * */
-void sync_vectors(double *source, double *target, int rows)
-{
- int i;
-
-#pragma omp for schedule (static)
- for (i = 0; i < rows; i++) {
- target[i] = source[i];
- }
-
- return;
-}
-
-/* ******************************************************* *
- * ***** Check if matrix is symmetric ******************** *
- * ******************************************************* */
-/*!
- * \brief Check if the matrix in les is symmetric
- *
- * \param L N_les* -- the linear equation system
- * \return int -- 1 = symmetric, 0 = unsymmetric
- *
- * */
-
-int check_symmetry(N_les * L)
-{
- int i, j, k;
- double value1 = 0;
- double value2 = 0;
- int index;
- int symm = 0;
-
- if (L->quad != 1) {
- G_warning(_("The linear equation system is not quadratic"));
- return 0;
- }
-
- G_debug(2, "check_symmetry: Check if matrix is symmetric");
-
- if (L->type == N_SPARSE_LES) {
-#pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
- for (j = 0; j < L->rows; j++) {
- for (i = 1; i < L->Asp[j]->cols; i++) {
- value1 = 0;
- value2 = 0;
- index = 0;
-
- index = L->Asp[j]->index[i];
- value1 = L->Asp[j]->values[i];
-
- for (k = 1; k < L->Asp[index]->cols; k++) {
- if (L->Asp[index]->index[k] == j) {
- value2 = L->Asp[index]->values[k];
- if ((value1 != value2)) {
- if ((fabs((fabs(value1) - fabs(value2))) <
- SYMM_TOLERANCE)) {
- G_debug(5,
- "check_symmetry: sparse matrix is unsymmetric, but within tolerance");
- }
- else {
- G_warning
- ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf \ndifference = %12.18lf\nStop symmetry calculation.\n",
- j, index, index, L->Asp[index]->index[k],
- value1, value2,
- fabs(fabs(value1) - fabs(value2)));
- symm++;
- }
- }
- }
- }
- }
- }
- }
- else {
-#pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
- for (i = 0; i < L->rows; i++) {
- for (j = i + 1; j < L->rows; j++) {
- if (L->A[i][j] != L->A[j][i]) {
- if ((fabs(fabs(L->A[i][j]) - fabs(L->A[j][i])) <
- SYMM_TOLERANCE)) {
- G_debug(5,
- "check_symmetry: matrix is unsymmetric, but within tolerance");
- }
- else {
- G_warning
- ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf\ndifference = %12.18lf\nStop symmetry calculation.\n",
- i, j, j, i, L->A[i][j], L->A[j][i],
- fabs(fabs(L->A[i][j]) - fabs(L->A[j][i])));
- symm++;
- }
- }
- }
- }
- }
-
- if (symm > 0)
- return 0;
-
- return 1;
-}
-
-
-/*!
- * \brief Compute a diagonal preconditioning matrix for krylov space solver
- *
- * \param L N_les*
- * \pram prec int -- the preconditioner which should be choosen N_DIAGONAL_PRECONDITION, N_ROWSUM_PRECONDITION
- * \return M N_les* -- the preconditioning matrix
- *
- * */
-N_les *N_create_diag_precond_matrix(N_les * L, int prec)
-{
- N_les *L_new;
- int rows = L->rows;
- int cols = L->cols;
- int i, j;
- double sum;
-
- L_new = N_alloc_les_A(rows, N_SPARSE_LES);
-
- if (L->type == N_NORMAL_LES) {
-#pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
- for (i = 0; i < rows; i++) {
- N_spvector *spvect = N_alloc_spvector(1);
-
- switch (prec) {
- case N_ROWSCALE_EUKLIDNORM_PRECONDITION:
- sum = 0;
- for (j = 0; j < cols; j++)
- sum += L->A[i][j] * L->A[i][j];
- spvect->values[0] = 1.0 / sqrt(sum);
- break;
- case N_ROWSCALE_ABSSUMNORM_PRECONDITION:
- sum = 0;
- for (j = 0; j < cols; j++)
- sum += fabs(L->A[i][j]);
- spvect->values[0] = 1.0 / (sum);
- break;
- case N_DIAGONAL_PRECONDITION:
- spvect->values[0] = 1.0 / L->A[i][i];
- break;
- default:
- spvect->values[0] = 1.0 / L->A[i][i];
- }
-
-
- spvect->index[0] = i;
- spvect->cols = 1;;
- N_add_spvector_to_les(L_new, spvect, i);
-
- }
- }
- else {
-#pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
- for (i = 0; i < rows; i++) {
- N_spvector *spvect = N_alloc_spvector(1);
-
- switch (prec) {
- case N_ROWSCALE_EUKLIDNORM_PRECONDITION:
- sum = 0;
- for (j = 0; j < L->Asp[i]->cols; j++)
- sum += L->Asp[i]->values[j] * L->Asp[i]->values[j];
- spvect->values[0] = 1.0 / sqrt(sum);
- break;
- case N_ROWSCALE_ABSSUMNORM_PRECONDITION:
- sum = 0;
- for (j = 0; j < L->Asp[i]->cols; j++)
- sum += fabs(L->Asp[i]->values[j]);
- spvect->values[0] = 1.0 / (sum);
- break;
- case N_DIAGONAL_PRECONDITION:
- spvect->values[0] = 1.0 / L->Asp[i]->values[0];
- break;
- default:
- spvect->values[0] = 1.0 / L->Asp[i]->values[0];
- }
-
- spvect->index[0] = i;
- spvect->cols = 1;;
- N_add_spvector_to_les(L_new, spvect, i);
- }
- }
- return L_new;
-}
Modified: grass/branches/develbranch_6/lib/gpde/N_tools.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/N_tools.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/N_tools.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -53,6 +53,7 @@
double N_calc_arith_mean_n(double *a, int size)
{
double val = 0.0;
+
int i;
for (i = 0; i < size; i++)
@@ -96,6 +97,7 @@
double N_calc_geom_mean_n(double *a, int size)
{
double val = 1;
+
int i;
for (i = 0; i < size; i++)
@@ -140,6 +142,7 @@
double N_calc_harmonic_mean_n(double *a, int size)
{
double val = 0;
+
int i;
for (i = 0; i < size; i++)
@@ -189,6 +192,7 @@
double N_calc_quad_mean_n(double *a, int size)
{
double val = 0;
+
int i;
for (i = 0; i < size; i++)
Modified: grass/branches/develbranch_6/lib/gpde/gpdelib.dox
===================================================================
--- grass/branches/develbranch_6/lib/gpde/gpdelib.dox 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/gpdelib.dox 2009-08-27 20:52:10 UTC (rev 38889)
@@ -11,10 +11,9 @@
solution of PDE's within GRASS.
Therefore it provides functions to create and solve linear equation systems.
The library design is thread safe and supports threaded parallelism with OpenMP.
-Most of the available solvers (expect the gauss seidel and jacobi solver)
-and the assembling of the linear equation systems are parallelized with OpenMP.
+The assembling of the linear equation systems is parallelized with OpenMP.
The creation of a linear equation system can be done by using quadratic or sparse matrices.
-Sparse and quadratic matrices are supported by the iterative equation solvers,
+Sparse and quadratic matrices are supported by the iterative equation solvers of the gmath library,
the direct equation solvers only support regular quadratic matrices.
<p>
To enable parallel computing, you will need a compiler with OpenMP capabilities
@@ -31,9 +30,6 @@
The basis for the discretisation is the geometrical structure of raster and volume maps.
<br>
<br>
-There are plans to implement heat flow in two and three dimensions.
-<br>
-<br>
Currently only planimetric projections are supported for numerical calculations.
<br>
<br>
@@ -52,99 +48,14 @@
to manage the data
2.) Assemble the linear equation system
- use the sparse matrix structure to save lots of memory and cpu time
-3.) Solve the linear equation system with the gauss, lu, jacobi, sor, cg or bicgstab method
- - always prefer the iterative krylov-space methods for solving
- - if the linear equation systems are in a bad condition try the jacobian solver
+3.) Solve the linear equation system with the gauss, lu, cg, pcg or bicgstab method
+ - always prefer the iterative krylov-space (cg, pcg and bicgstab) methods for solving
- the available direct solvers don't support sparse matrices
4.) Write the result back as raster/volume map
\endverbatim
-The following code shows how to implement this principle with the GRASS PDE lib:
+The modules r.gwflow and r3.gwflow are examples how to use the gpde library.
-\verbatim
-/* *************************************************************** */
-/* ****** calculate 2d dimensional groundwater flow ************** */
-/* *************************************************************** */
-int main()
-{
- int i, j;
- N_gwflow_data2d *data = NULL; /* data structure with the groundwater data */
- N_geom_data *geom = NULL; /* geometry of the calculation area (region) */
- N_les *les = NULL; /* the linear equation system structure */
- N_les_callback_2d *call = NULL; /* the callback for the groundwater flow calulation */
-
- /* allocate the callback structure */
- call = N_alloc_les_callback_2d();
-
- /* set the callback to groundwater flow calculation */
- N_set_les_callback_2d_func(call, (*N_callback_gwflow_2d));
-
- /* allocate the groundwater data structure with one million cells */
- data = N_alloc_gwflow_data2d(1000, 1000);
-
- /* get the current region */
- G_get_set_window(®ion);
-
- /* allocate and initiate the geometry structure for geometry and area calculation
- needed by the groundwater calculation callback */
- geom = N_init_geom_data_2d(®ion, geom);
-
-
- /*fill the data arrays with raster maps*/
- N_read_rast_to_array_2d("phead_map_name", data->phead);
- N_read_rast_to_array_2d("phead_map_name", data->phead_start);
- N_read_rast_to_array_2d("status_map_name", data->status);
- N_read_rast_to_array_2d("hc_x_map_name", data->hc_x);
- N_read_rast_to_array_2d("hc_y_map_name", data->hc_y);
- N_read_rast_to_array_2d("q_map_name", data->q);
- N_read_rast_to_array_2d("s_map_name", data->s);
- N_read_rast_to_array_2d("top_map_name", data->top);
- N_read_rast_to_array_2d("bottom_map_name", data->bottom);
- N_read_rast_to_array_2d("r_map_name", data->r);
-
- /* Set the inactive cells to zero, to assure a no flow boundary */
- /* notice: the data arrays are of type DCELL, so the put_*_d_value functions are used */
- for (y = 0; y < geom->rows; y++) {
- for (x = 0; x < geom->cols; x++) {
- stat = (int)N_get_array_2d_d_value(data->status, x, y);
- if (stat == N_CELL_INACTIVE) { /*only inactive cells */
- N_put_array_2d_d_value(data->hc_x, x, y, 0);
- N_put_array_2d_d_value(data->hc_y, x, y, 0);
- N_put_array_2d_d_value(data->s, x, y, 0);
- N_put_array_2d_d_value(data->q, x, y, 0);
- }
- }
- }
-
- /*Assemble the sparse matrix */
- les = N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start, (void *)data, call);
-
- /*Solve the linear equation system with the cg method*/
- N_solver_cg(les, 100000, 0.1e-8);
-
- /* now copy the result from the x vector of the linear equation system
- into a data array, this function is not provided by the gpde lib,
- you have to write your own: take a look at r.gwflow
- */
- copy_x_vect_to_data_array(les->x, data->phead, geom);
-
- /*write the x vector of the les back into a raster map*/
- N_write_array_2d_to_rast(data->phead, "output_map_name");
-
- /*free the memory*/
- N_free_les(les);
- N_free_gwflow_data2d(data);
- N_free_geom_data(geom);
- G_free(call);
-
- return 0;
-}
-
-\endverbatim
-
-The assembling of lineare equation systems is based on a 5 point star for raster maps and a
-7 point star for volume maps, implemented as finite volume/differences discretization.
-
<p>
All raster and volume maps which are used to create a linear equation system
must be loaded completely into the memory.<br>
@@ -397,45 +308,6 @@
<p>
N_data_star *#N_callback_template_2d(void *data, N_geom_data * geom, int row, int col);
-
-\section solvers Available solvers
-
-\subsection direct_solvers Direct solvers
-
-Two direct solvers are implemented, the gauss elimination solver and the lu decomposition solver.
-The direct solvers only work with regular quadratic matrices, not with sparse matrices.
-
-<p>
-int #N_solver_gauss (N_les * les);
-
-<p>
-int #N_solver_lu (N_les * les);
-
-
-\subsection iterative_solvers Iterative solvers
-
-The iterative solvers work with regular quadartic and sparse matrices.
-<p>
-To solve symmetric and positive definite linear equation systems the iterative
-conjugated gradient method with additional diagonal preconditioning are implemented.
-
-int #N_solver_cg(N_les * les, int maxit, double error);
-<p>
-int #N_solver_pcg(N_les * les, int maxit, double error);
-
-<p>
-To solve unsymmetric non definite linear equation system the iterative BiCGSatb method is implemented
-
-<p>
-int #N_solver_bicgstab(N_les * les, int maxit, double error);
-
-<p>
-Additionally jacobi and Gauss-Seidel iterative solvers with relaxation are implemented
-<p>
-int #N_solver_jacobi (N_les * L, int maxit, double sor, double error);
-<p>
-int #N_solver_SOR (N_les * L, int maxit, double sor, double error);
-
\section available_pdes Implemented PDE's
Groundwater flow in 2 and 3 dimensions are implemented.
Deleted: grass/branches/develbranch_6/lib/gpde/solvers_local_proto.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/solvers_local_proto.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/solvers_local_proto.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,22 +0,0 @@
-
-/*****************************************************************************
-*
-* MODULE: Grass PDE Numerical Library
-* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
-* soerengebbert <at> gmx <dot> de
-*
-* PURPOSE: local prototypes for linear equation system solvers
-* part of the gpde 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.
-*
-*****************************************************************************/
-
-#define SYMM_TOLERANCE 1.0e-18
-
-double *vectmem(int size);
-int check_symmetry(N_les * les);
Modified: grass/branches/develbranch_6/lib/gpde/test/Makefile
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/Makefile 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/Makefile 2009-08-27 20:52:10 UTC (rev 38889)
@@ -2,8 +2,8 @@
PGM=test.gpde.lib
-LIBES = $(GISLIB) $(G3DLIB) $(GPDELIB)
-DEPENDENCIES = $(GISDEP) $(G3DDEP) $(GPDEDEP)
+LIBES = $(GISLIB) $(G3DLIB) $(GPDELIB) $(GMATHLIB)
+DEPENDENCIES = $(GISDEP) $(G3DDEP) $(GPDEDEP) $(GMATDEP)
include $(MODULE_TOPDIR)/include/Make/Module.make
Added: grass/branches/develbranch_6/lib/gpde/test/test.gpde.lib.html
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test.gpde.lib.html (rev 0)
+++ grass/branches/develbranch_6/lib/gpde/test/test.gpde.lib.html 2009-08-27 20:52:10 UTC (rev 38889)
@@ -0,0 +1,3 @@
+
+
+Take a look at the module command line help for more information.
Modified: grass/branches/develbranch_6/lib/gpde/test/test_arrays.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_arrays.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_arrays.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -332,12 +332,10 @@
data11 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, CELL_TYPE);
data2 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, FCELL_TYPE);
N_print_array_2d_info(data2);
- data22 =
- N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, FCELL_TYPE);
+ data22 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, FCELL_TYPE);
data3 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE);
N_print_array_2d_info(data3);
- data33 =
- N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE);
+ data33 = N_alloc_array_2d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, 1, DCELL_TYPE);
/*Fill the first arrays with data */
@@ -529,8 +527,7 @@
N_math_array_2d(data2, data22, tmp, N_ARRAY_SUM);
res = N_convert_array_2d_null_to_zero(tmp);
if (res == 0) {
- G_warning
- ("test_array_2d: error in N_convert_array_2d_null_to_zero ");
+ G_warning("test_array_2d: error in N_convert_array_2d_null_to_zero ");
sum++;
}
N_free_array_2d(tmp);
@@ -627,19 +624,19 @@
/*Alloacte memory for all arrays */
data1 =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 2, FCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
+ FCELL_TYPE);
N_print_array_3d_info(data1);
data11 =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 2, FCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
+ FCELL_TYPE);
data2 =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 2, DCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
+ DCELL_TYPE);
N_print_array_3d_info(data2);
data22 =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 2, DCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 2,
+ DCELL_TYPE);
/*Fill the first arrays with data */
Modified: grass/branches/develbranch_6/lib/gpde/test/test_assemble.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_assemble.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_assemble.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -73,7 +73,7 @@
N_put_array_2d_c_value(data, i, j, 2);
}
else {
- N_put_array_2d_c_value(data, i, j, 1);
+ N_put_array_2d_c_value(data, i, j, 1);
}
}
}
@@ -98,7 +98,7 @@
N_put_array_2d_d_value(data, i, j, 50);
}
else {
- N_put_array_2d_d_value(data, i, j, 1);
+ N_put_array_2d_d_value(data, i, j, 1);
}
}
}
@@ -114,8 +114,8 @@
int i, j, k;
data =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 1, FCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 1,
+ FCELL_TYPE);
#pragma omp parallel for private (i, j, k) shared (data)
for (k = 0; k < TEST_N_NUM_DEPTHS; k++)
@@ -127,7 +127,7 @@
}
else {
- N_put_array_3d_f_value(data, i, j, k, 1.0);
+ N_put_array_3d_f_value(data, i, j, k, 1.0);
}
}
}
@@ -145,8 +145,8 @@
int i, j, k;
data =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 1, DCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 1,
+ DCELL_TYPE);
#pragma omp parallel for private (i, j, k) shared (data)
for (k = 0; k < TEST_N_NUM_DEPTHS; k++)
@@ -159,7 +159,7 @@
}
else {
- N_put_array_3d_f_value(data, i, j, k, 1);
+ N_put_array_3d_f_value(data, i, j, k, 1);
}
}
}
@@ -199,29 +199,19 @@
geom->cols = TEST_N_NUM_COLS;
/*Assemble the matrix */
- les =
- N_assemble_les_3d(N_SPARSE_LES, geom, status, start_val, NULL, call);
+ les = N_assemble_les_3d(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_3d_active(N_SPARSE_LES, geom, status, start_val, NULL,
- call);
+ les = N_assemble_les_3d_active(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, status, start_val,
- NULL, call);
+ les = N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_les_integrate_dirichlet_3d(les, geom, status, start_val);
N_free_les(les);
- les =
- N_assemble_les_3d(N_NORMAL_LES, geom, status, start_val, NULL, call);
+ les = N_assemble_les_3d(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_3d_active(N_NORMAL_LES, geom, status, start_val, NULL,
- call);
+ les = N_assemble_les_3d_active(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, status, start_val,
- NULL, call);
+ les = N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_les_integrate_dirichlet_3d(les, geom, status, start_val);
N_free_les(les);
@@ -261,29 +251,19 @@
geom->cols = TEST_N_NUM_COLS;
/*Assemble the matrix */
- les =
- N_assemble_les_2d(N_SPARSE_LES, geom, status, start_val, NULL, call);
+ les = N_assemble_les_2d(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_2d_active(N_SPARSE_LES, geom, status, start_val, NULL,
- call);
+ les = N_assemble_les_2d_active(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, status, start_val,
- NULL, call);
+ les = N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, status, start_val, NULL, call);
N_les_integrate_dirichlet_2d(les, geom, status, start_val);
N_free_les(les);
- les =
- N_assemble_les_2d(N_NORMAL_LES, geom, status, start_val, NULL, call);
+ les = N_assemble_les_2d(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_2d_active(N_NORMAL_LES, geom, status, start_val, NULL,
- call);
+ les = N_assemble_les_2d_active(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_free_les(les);
- les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, status, start_val,
- NULL, call);
+ les = N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, status, start_val, NULL, call);
N_les_integrate_dirichlet_2d(les, geom, status, start_val);
N_free_les(les);
Modified: grass/branches/develbranch_6/lib/gpde/test/test_geom.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_geom.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_geom.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -111,8 +111,7 @@
area += N_get_geom_data_area_of_cell(geom, i);
if (area == 0) {
- G_warning
- ("Wrong area calculation in N_get_geom_data_area_of_cell");
+ G_warning("Wrong area calculation in N_get_geom_data_area_of_cell");
sum++;
}
}
@@ -149,8 +148,7 @@
area += N_get_geom_data_area_of_cell(geom, i);
if (area == 0) {
- G_warning
- ("Wrong area calculation in N_get_geom_data_area_of_cell");
+ G_warning("Wrong area calculation in N_get_geom_data_area_of_cell");
sum++;
}
}
@@ -161,8 +159,7 @@
area += N_get_geom_data_area_of_cell(geom, i);
if (area == 0) {
- G_warning
- ("Wrong area calculation in N_get_geom_data_area_of_cell");
+ G_warning("Wrong area calculation in N_get_geom_data_area_of_cell");
sum++;
}
}
Modified: grass/branches/develbranch_6/lib/gpde/test/test_gpde_lib.h
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_gpde_lib.h 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_gpde_lib.h 2009-08-27 20:52:10 UTC (rev 38889)
@@ -24,7 +24,7 @@
#define TEST_N_NUM_DEPTHS 10
/* Array test functions */
-extern int unit_test_arrays(void);
+extern int unit_test_arrays (void);
/* matrix assembling */
extern int unit_test_assemble(void);
@@ -32,9 +32,6 @@
/* gradient creation and handling tests */
extern int unit_test_gradient(void);
-/* direct and iterative solvers */
-extern int unit_test_solvers(void);
-
/* test the meth tools of gpde */
extern int unit_test_tools(void);
@@ -42,10 +39,10 @@
extern int unit_test_geom_data(void);
/* les creation */
-extern int unit_test_les_creation(void);
+extern int unit_test_les_creation (void);
-/*gwflow */
-extern int integration_test_gwflow(void);
+/*gwflow*/
+extern int integration_test_gwflow (void);
/* solute transport */
extern int integration_test_solute_transport(void);
Modified: grass/branches/develbranch_6/lib/gpde/test/test_gradient.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_gradient.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_gradient.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -68,7 +68,7 @@
#pragma omp parallel for private (i, j) shared (data)
for (j = 0; j < TEST_N_NUM_ROWS; j++) {
for (i = 0; i < TEST_N_NUM_COLS; i++) {
- N_put_array_2d_c_value(data, i, j, 1);
+ N_put_array_2d_c_value(data, i, j, 1);
}
}
return data;
@@ -87,7 +87,7 @@
#pragma omp parallel for private (i, j) shared (data)
for (j = 0; j < TEST_N_NUM_ROWS; j++) {
for (i = 0; i < TEST_N_NUM_COLS; i++) {
- N_put_array_2d_d_value(data, i, j, (double)i * j);
+ N_put_array_2d_d_value(data, i, j, (double)i*j);
}
}
return data;
@@ -102,17 +102,17 @@
int i, j, k;
data =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 1, FCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 1,
+ FCELL_TYPE);
#pragma omp parallel for private (i, j, k) shared (data)
for (k = 0; k < TEST_N_NUM_DEPTHS; k++) {
for (j = 0; j < TEST_N_NUM_ROWS; j++) {
for (i = 0; i < TEST_N_NUM_COLS; i++) {
- N_put_array_3d_f_value(data, i, j, k, 1.0);
+ N_put_array_3d_f_value(data, i, j, k, 1.0);
+ }
}
}
- }
return data;
@@ -127,16 +127,16 @@
int i, j, k;
data =
- N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS,
- 1, DCELL_TYPE);
+ N_alloc_array_3d(TEST_N_NUM_COLS, TEST_N_NUM_ROWS, TEST_N_NUM_DEPTHS, 1,
+ DCELL_TYPE);
#pragma omp parallel for private (i, j, k) shared (data)
for (k = 0; k < TEST_N_NUM_DEPTHS; k++)
for (j = 0; j < TEST_N_NUM_ROWS; j++) {
for (i = 0; i < TEST_N_NUM_COLS; i++) {
- N_put_array_3d_f_value(data, i, j, k, (double)i * j * k);
+ N_put_array_3d_f_value(data, i, j, k, (float)i*j*k);
+ }
}
- }
return data;
@@ -174,10 +174,9 @@
pot = create_potential_array_3d();
field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL);
- field =
- N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
+ field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
- /*compute stats */
+ /*compute stats*/
N_calc_gradient_field_3d_stats(field);
N_print_gradient_field_3d_info(field);
@@ -187,7 +186,7 @@
N_free_array_3d(pot);
relax = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
- pot = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
+ pot = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
xcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
ycomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
zcomp = N_alloc_array_3d(3, 3, 3, 0, DCELL_TYPE);
@@ -267,8 +266,7 @@
geom->cols = 3;
field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, NULL);
- field =
- N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
+ field = N_compute_gradient_field_3d(pot, relax, relax, relax, geom, field);
N_print_gradient_field_3d_info(field);
grad = N_get_gradient_3d(field, NULL, 0, 0, 0);
@@ -302,10 +300,10 @@
N_compute_gradient_field_components_3d(field, xcomp, ycomp, zcomp);
/*
- N_print_array_3d(xcomp);
- N_print_array_3d(ycomp);
- N_print_array_3d(zcomp);
- */
+ N_print_array_3d(xcomp);
+ N_print_array_3d(ycomp);
+ N_print_array_3d(zcomp);
+ */
N_free_gradient_field_3d(field);
G_free(geom);
@@ -354,7 +352,7 @@
field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
- /*compute stats */
+ /*compute stats*/
N_calc_gradient_field_2d_stats(field);
N_print_gradient_field_2d_info(field);
@@ -364,7 +362,7 @@
N_free_array_2d(pot);
relax = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
- pot = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
+ pot = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
xcomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
ycomp = N_alloc_array_2d(3, 3, 0, DCELL_TYPE);
@@ -401,55 +399,46 @@
field = N_compute_gradient_field_2d(pot, relax, relax, geom, field);
N_print_gradient_field_2d_info(field);
- /*common gradient calculation */
+ /*common gradient calculation*/
grad = N_get_gradient_2d(field, NULL, 0, 0);
- G_message
- ("Gradient 2d: pos 0,0 NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1\n",
- grad->NC, grad->SC, grad->WC, grad->EC);
+ G_message("Gradient 2d: pos 0,0 NC %g == 0 ; SC %g == 2 ; WC %g == 0 ; EC %g == -1\n",
+ grad->NC, grad->SC, grad->WC, grad->EC);
grad = N_get_gradient_2d(field, grad, 1, 0);
- G_message
- ("Gradient 2d: pos 1,0 NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4\n",
- grad->NC, grad->SC, grad->WC, grad->EC);
+ G_message("Gradient 2d: pos 1,0 NC %g == 0 ; SC %g == 5 ; WC %g == -1 ; EC %g == -4\n",
+ grad->NC, grad->SC, grad->WC, grad->EC);
N_free_gradient_2d(grad);
grad = N_get_gradient_2d(field, NULL, 1, 1);
- G_message
- ("Gradient 2d: pos 1,1 NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3\n",
- grad->NC, grad->SC, grad->WC, grad->EC);
+ G_message("Gradient 2d: pos 1,1 NC %g == 5 ; SC %g == 8 ; WC %g == -4 ; EC %g == -3\n",
+ grad->NC, grad->SC, grad->WC, grad->EC);
grad = N_get_gradient_2d(field, grad, 1, 2);
- G_message
- ("Gradient 2d: pos 1,2 NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10\n",
- grad->NC, grad->SC, grad->WC, grad->EC);
+ G_message("Gradient 2d: pos 1,2 NC %g == 8 ; SC %g == 0 ; WC %g == -7 ; EC %g == -10\n",
+ grad->NC, grad->SC, grad->WC, grad->EC);
N_free_gradient_2d(grad);
grad = N_get_gradient_2d(field, NULL, 2, 2);
- G_message
- ("Gradient 2d: pos 2,2 NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0\n",
- grad->NC, grad->SC, grad->WC, grad->EC);
+ G_message("Gradient 2d: pos 2,2 NC %g ==15 ; SC %g == 0 ; WC %g == -10 ; EC %g == 0\n",
+ grad->NC, grad->SC, grad->WC, grad->EC);
N_free_gradient_2d(grad);
N_compute_gradient_field_components_2d(field, xcomp, ycomp);
- /*gradient neighbour calculation */
+ /*gradient neighbour calculation*/
grad_2d = N_get_gradient_neighbours_2d(field, NULL, 1, 1);
grad_2d = N_get_gradient_neighbours_2d(field, grad_2d, 1, 1);
- G_message
- ("N_gradient_neighbours_x; pos 1,1 NWN %g NEN %g WC %g EC %g SWS %g SES %g\n",
- grad_2d->x->NWN, grad_2d->x->NEN, grad_2d->x->WC, grad_2d->x->EC,
- grad_2d->x->SWS, grad_2d->x->SES);
- G_message
- ("N_gradient_neighbours_y: pos 1,1 NWW %g NEE %g NC %g SC %g SWW %g SEE %g\n",
- grad_2d->y->NWW, grad_2d->y->NEE, grad_2d->y->NC, grad_2d->y->SC,
- grad_2d->y->SWW, grad_2d->y->SEE);
-
+ G_message("N_gradient_neighbours_x; pos 1,1 NWN %g NEN %g WC %g EC %g SWS %g SES %g\n",
+ grad_2d->x->NWN, grad_2d->x->NEN, grad_2d->x->WC, grad_2d->x->EC, grad_2d->x->SWS, grad_2d->x->SES);
+ G_message("N_gradient_neighbours_y: pos 1,1 NWW %g NEE %g NC %g SC %g SWW %g SEE %g\n",
+ grad_2d->y->NWW, grad_2d->y->NEE, grad_2d->y->NC, grad_2d->y->SC, grad_2d->y->SWW, grad_2d->y->SEE);
+
N_free_gradient_neighbours_2d(grad_2d);
/*
- N_print_array_2d(xcomp);
- N_print_array_2d(ycomp);
- */
+ N_print_array_2d(xcomp);
+ N_print_array_2d(ycomp);
+ */
N_free_gradient_field_2d(field);
G_free(geom);
Modified: grass/branches/develbranch_6/lib/gpde/test/test_gwflow.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_gwflow.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_gwflow.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -15,6 +15,9 @@
*
*****************************************************************************/
+#include <grass/gmath.h>
+
+
#include <grass/gis.h>
#include <grass/glocale.h>
#include <grass/N_pde.h>
@@ -108,9 +111,7 @@
int i, j;
N_gwflow_data2d *data;
- data =
- N_alloc_gwflow_data2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1,
- 1);
+ data = N_alloc_gwflow_data2d(TEST_N_NUM_COLS_LOCAL, TEST_N_NUM_ROWS_LOCAL, 1, 1);
#pragma omp parallel for private (i, j) shared (data)
for (j = 0; j < TEST_N_NUM_ROWS_LOCAL; j++) {
@@ -179,37 +180,37 @@
/*CG*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_cg(les, 100, 0.1e-8);
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
- /*PCG N_DIAGONAL_PRECONDITION */ les =
+ /*PCG N_DIAGONAL_PRECONDITION*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les =
+ /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les =
+ /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
/*CG*/ les =
- N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_cg(les, 100, 0.1e-8);
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
@@ -217,41 +218,41 @@
/*CG*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
-
- N_solver_cg(les, 100, 0.1e-8);
+
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
- /*PCG N_DIAGONAL_PRECONDITION */ les =
+ /*PCG N_DIAGONAL_PRECONDITION*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
-
- N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION */ les =
+ /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
-
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION */ les =
+ /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
-
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
/*CG*/ les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_cg(les, 100, 0.1e-8);
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
@@ -259,30 +260,30 @@
/*BICG*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
- N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
@@ -290,45 +291,45 @@
/*GUASS*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_gauss(les);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*LU*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_lu(les);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*GUASS*/ les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_gauss(les);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*LU*/ les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_lu(les);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
- /*Cholesky */ les =
+ /*Cholesky*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_cholesky(les);
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
N_print_les(les);
N_free_les(les);
- /*Cholesky */ les =
- N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ /*Cholesky*/ les =
+ N_assemble_les_3d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_3d(les, geom, data->status, data->phead_start);
- N_solver_cholesky(les);
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
N_print_les(les);
N_free_les(les);
@@ -365,53 +366,85 @@
geom->cols = TEST_N_NUM_COLS_LOCAL;
- /*Assemble the matrix */
+ /*Assemble the matrix */
/*
*/
/*CG*/ les =
N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_cg(les, 100, 0.1e-8);
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
+ /*PCG N_DIAGONAL_PRECONDITION*/ les =
+ N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+ N_print_les(les);
+ N_free_les(les);
+
+ /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION*/ les =
+ N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
+ N_print_les(les);
+ N_free_les(les);
+
+ /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION*/ les =
+ N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+ G_math_solver_sparse_pcg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
+ N_print_les(les);
+ N_free_les(les);
+
+
/*CG*/ les =
- N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_cg(les, 100, 0.1e-8);
+ G_math_solver_sparse_cg(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
+
/*CG*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_cg(les, 100, 0.1e-8);
+
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
+ /*PCG N_DIAGONAL_PRECONDITION*/ les =
+ N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+ N_print_les(les);
+ N_free_les(les);
- /*CG*/ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
- N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_cg(les, 100, 0.1e-8);
+ /*PCG N_ROWSCALE_EUKLIDNORM_PRECONDITION*/ les =
+ N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG*/ les =
+ /*PCG N_ROWSCALE_ABSSUMNORM_PRECONDITION*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+
+ G_math_solver_pcg(les->A, les->x, les->b, les->rows, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
N_print_les(les);
N_free_les(les);
- /*PCG*/ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ /*CG*/ les =
+ N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
+ G_math_solver_cg(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
@@ -419,83 +452,79 @@
/*BICG*/ les =
N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
- N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status,
- data->phead_start, (void *)data, call);
- N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
-
/*BICG*/ les =
- N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ N_assemble_les_2d_dirichlet(N_SPARSE_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
- /*GAUSS*/ les =
+
+ /*GUASS*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_gauss(les);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
- /*GAUSS*/ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
- N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_gauss(les);
+ /*LU*/ les =
+ N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
-
- /*LU*/ les =
- N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
+ /*GUASS*/ les =
+ N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_lu(les);
+ N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*LU*/ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_lu(les);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
- /*Cholesky */ les =
+ /*Cholesky*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->phead_start,
(void *)data, call);
- N_solver_cholesky(les);
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
N_print_les(les);
N_free_les(les);
- /*Cholesky */ les =
- N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status,
- data->phead_start, (void *)data, call);
+ /*Cholesky*/ les =
+ N_assemble_les_2d_dirichlet(N_NORMAL_LES, geom, data->status, data->phead_start,
+ (void *)data, call);
N_les_integrate_dirichlet_2d(les, geom, data->status, data->phead_start);
- N_solver_cholesky(les);
+ G_math_solver_cholesky(les->A, les->x, les->b, les->rows, les->rows);
N_print_les(les);
N_free_les(les);
-
-
-
-
+
N_free_gwflow_data2d(data);
G_free(geom);
G_free(call);
Modified: grass/branches/develbranch_6/lib/gpde/test/test_les.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_les.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_les.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -17,6 +17,7 @@
#include <grass/gis.h>
#include <grass/glocale.h>
+#include <grass/gmath.h>
#include <grass/N_pde.h>
#include "test_gpde_lib.h"
@@ -49,7 +50,7 @@
/* *************************************************************** */
int test_les(void)
{
- N_spvector *spvector = NULL;
+ G_math_spvector *spvector = NULL;
N_les *les = NULL;
N_les *sples = NULL;
int i, j;
@@ -106,7 +107,7 @@
#pragma omp parallel for private(i, j) shared(sples, spvector)
for (i = 0; i < TEST_N_NUM_ROWS; i++) {
- spvector = N_alloc_spvector(TEST_N_NUM_ROWS);
+ spvector = G_math_alloc_spvector(TEST_N_NUM_ROWS);
for (j = 0; j < TEST_N_NUM_ROWS; j++)
if (i != j)
@@ -115,7 +116,7 @@
spvector->index[0] = i;
spvector->values[0] = -1e2 - i;
- N_add_spvector_to_les(sples, spvector, i);
+ G_math_add_spvector(sples->Asp, spvector, i);
sples->x[i] = 273.15 + i;
sples->b[i] = 1e2 - i;
}
@@ -139,7 +140,7 @@
}
for (i = 0; i < TEST_N_NUM_ROWS; i++) {
- spvector = N_alloc_spvector(TEST_N_NUM_ROWS);
+ spvector = G_math_alloc_spvector(TEST_N_NUM_ROWS);
for (j = 0; j < TEST_N_NUM_ROWS; j++)
if (i != j)
@@ -148,7 +149,7 @@
spvector->index[0] = i;
spvector->values[0] = -1e2 - i;
- N_add_spvector_to_les(sples, spvector, i);
+ G_math_add_spvector(sples->Asp, spvector, i);
sples->x[i] = 273.15 + i;
sples->b[i] = 1e2 - i;
}
Modified: grass/branches/develbranch_6/lib/gpde/test/test_main.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_main.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_main.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -34,7 +34,7 @@
paramType param; /*Parameters */
/*- prototypes --------------------------------------------------------------*/
-static void set_params(void); /*Fill the paramType structure */
+static void set_params(void); /*Fill the paramType structure */
/* ************************************************************************* */
/* Set up the arguments we are expecting ********************************** */
@@ -45,7 +45,7 @@
param.unit->key = "unit";
param.unit->type = TYPE_STRING;
param.unit->required = NO;
- param.unit->options = "array,assemble,geom,gradient,les,solver,tools";
+ param.unit->options = "array,assemble,geom,gradient,les,tools";
param.unit->description = _("Choose the unit tests to run");
param.integration = G_define_option();
@@ -100,7 +100,6 @@
returnstat += unit_test_gradient();
returnstat += unit_test_geom_data();
returnstat += unit_test_les_creation();
- returnstat += unit_test_solvers();
returnstat += unit_test_tools();
}
@@ -133,9 +132,6 @@
if (strcmp(param.unit->answers[i], "les") == 0)
returnstat += unit_test_les_creation();
- if (strcmp(param.unit->answers[i], "solver") == 0)
- returnstat += unit_test_solvers();
-
if (strcmp(param.unit->answers[i], "tools") == 0)
returnstat += unit_test_tools();
@@ -150,11 +146,10 @@
if (strcmp(param.integration->answers[i], "gwflow") == 0)
returnstat += integration_test_gwflow();
- if (strcmp(param.integration->answers[i], "heatflow") == 0) ; /*nothing to do for now */
+ if (strcmp(param.integration->answers[i], "heatflow") == 0); /*nothing to do for now */
- if (strcmp(param.integration->answers[i], "transport") ==
- 0)
- returnstat += integration_test_solute_transport();
+ if (strcmp(param.integration->answers[i], "transport") == 0)
+ returnstat += integration_test_solute_transport();
i++;
}
@@ -162,10 +157,10 @@
}
}
- if (returnstat != 0)
- G_warning("Errors detected while testing the gpde lib");
+ if(returnstat != 0)
+ G_warning("Errors detected while testing the gpde lib");
else
- G_message("\n-- gpde lib tests finished successfully --");
+ G_message("\n-- gpde lib tests finished successfully --");
return (returnstat);
}
Modified: grass/branches/develbranch_6/lib/gpde/test/test_solute_transport.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_solute_transport.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_solute_transport.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -51,7 +51,8 @@
if (sum > 0)
G_warning(_("\n-- solute_transport integration tests failure --"));
else
- G_message(_("\n-- solute_transport integration tests finished successfully --"));
+ G_message(_
+ ("\n-- solute_transport integration tests finished successfully --"));
return sum;
}
@@ -143,9 +144,9 @@
N_put_array_2d_d_value(data->cs, i, j, 1.0);
}
}
- /*dispersivity length */
- data->al = 0.2;
- data->at = 0.02;
+ /*dispersivity length*/
+ data->al = 0.2;
+ data->at = 0.02;
@@ -168,7 +169,7 @@
N_set_les_callback_3d_func(call, (*N_callback_solute_transport_3d)); /*solute_transport 3d */
data = create_solute_transport_data_3d();
-
+
N_calc_solute_transport_disptensor_3d(data);
data->dt = 86400;
@@ -187,59 +188,32 @@
/*Assemble the matrix */
/*
*/
- /*Jacobi */ les =
- N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_jacobi(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
- /*jacobi */ les =
- N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_jacobi(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
- /*SOR*/ les =
- N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_SOR(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
- /*SOR*/ les =
- N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_SOR(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
/*BICG*/ les =
N_assemble_les_3d(N_SPARSE_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*BICG*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
/*GUASS*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_gauss(les);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*LU*/ les =
N_assemble_les_3d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_lu(les);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
@@ -302,65 +276,37 @@
N_free_gradient_field_2d(field);
N_compute_gradient_field_2d(pot, relax, relax, geom, data->grad);
- /*The dispersivity tensor */
+ /*The dispersivity tensor*/
N_calc_solute_transport_disptensor_2d(data);
/*Assemble the matrix */
/*
*/
- /*Jacobi */ les =
+ /*BICG*/ les =
N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_jacobi(les, 100, 1, 0.1e-8);
+ G_math_solver_sparse_bicgstab(les->Asp, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
- /*jacobi */ les =
- N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_jacobi(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
- /*SOR*/ les =
- N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_SOR(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
- /*SOR*/ les =
- N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_SOR(les, 100, 1, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
/*BICG*/ les =
- N_assemble_les_2d(N_SPARSE_LES, geom, data->status, data->c_start,
- (void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
- N_print_les(les);
- N_free_les(les);
-
- /*BICG*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_bicgstab(les, 100, 0.1e-8);
+ G_math_solver_bicgstab(les->A, les->x, les->b, les->rows, 100, 0.1e-8);
N_print_les(les);
N_free_les(les);
- /*GAUSS*/ les =
+ /*GUASS*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_gauss(les);
+ G_math_solver_gauss(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
/*LU*/ les =
N_assemble_les_2d(N_NORMAL_LES, geom, data->status, data->c_start,
(void *)data, call);
- N_solver_lu(les);
+ G_math_solver_lu(les->A, les->x, les->b, les->rows);
N_print_les(les);
N_free_les(les);
Deleted: grass/branches/develbranch_6/lib/gpde/test/test_solvers.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_solvers.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_solvers.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -1,235 +0,0 @@
-
-/*****************************************************************************
-*
-* 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/N_pde.h>
-#include "test_gpde_lib.h"
-
-/* prototypes */
-static int test_solvers(void);
-static N_les *create_normal_les(int rows);
-static N_les *create_sparse_les(int rows);
-
-/* ************************************************************************* */
-/* Performe the solver unit tests ****************************************** */
-/* ************************************************************************* */
-int unit_test_solvers(void)
-{
- int sum = 0;
-
- G_message(_("\n++ Running solver unit tests ++"));
-
- sum += test_solvers();
-
- if (sum > 0)
- G_warning(_("\n-- Solver unit tests failure --"));
- else
- G_message(_("\n-- Solver unit tests finished successfully --"));
-
- return sum;
-}
-
-/* *************************************************************** */
-/* create a normal matrix with values ** Hilbert matrix ********** */
-/* *************************************************************** */
-N_les *create_normal_les(int rows)
-{
- N_les *les;
- int i, j;
- int size = rows;
- double val;
-
- les = N_alloc_les(rows, N_NORMAL_LES);
- for (i = 0; i < size; i++) {
- val = 0.0;
- for (j = 0; j < size; j++) {
- les->A[i][j] = (double)(1.0 / (((double)i + 1.0) +
- ((double)j + 1.0) - 1.0));
- val += les->A[i][j];
- }
- les->b[i] = val;
- }
-
- return les;
-}
-
-/* *************************************************************** */
-/* create a sparse matrix with values ** Hilbert matrix ********** */
-/* *************************************************************** */
-N_les *create_sparse_les(int rows)
-{
- N_les *les;
- N_spvector *spvector;
- int i, j;
- double val;
-
- les = N_alloc_les(rows, N_SPARSE_LES);
-
- for (i = 0; i < rows; i++) {
- spvector = N_alloc_spvector(rows);
- val = 0;
-
- for (j = 0; j < rows; j++) {
- spvector->values[j] =
- (double)(1.0 / (((double)i + 1.0) + ((double)j + 1.0) - 1.0));
- spvector->index[j] = j;
- val += spvector->values[j];
- }
-
- N_add_spvector_to_les(les, spvector, i);
- les->b[i] = val;
- }
-
-
- return les;
-}
-
-
-/* *************************************************************** */
-/* Test all implemented solvers for sparse and normal matrices *** */
-/* *************************************************************** */
-int test_solvers(void)
-{
- N_les *les;
- N_les *sples;
-
- G_message("\t * testing jacobi solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_jacobi(les, 100, 1, 0.1e-4);
- /*N_print_les(les); */
- N_solver_jacobi(sples, 100, 1, 0.1e-4);
- /*N_print_les(sples); */
-
- N_free_les(les);
- N_free_les(sples);
-
-
- G_message("\t * testing SOR solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_SOR(les, 100, 1, 0.1e-4);
- /*N_print_les(les); */
- N_solver_SOR(sples, 100, 1, 0.1e-4);
- /*N_print_les(sples); */
-
- N_free_les(les);
- N_free_les(sples);
-
- G_message("\t * testing cg solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_cg(les, 100, 0.1e-8);
- /*N_print_les(les); */
- N_solver_cg(sples, 100, 0.1e-8);
- /*N_print_les(sples); */
-
- N_free_les(les);
- N_free_les(sples);
-
- G_message("\t * testing pcg solver with N_DIAGONAL_PRECONDITION\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_pcg(les, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
- N_print_les(les);
- N_solver_pcg(sples, 100, 0.1e-8, N_DIAGONAL_PRECONDITION);
- N_print_les(sples);
-
- N_free_les(les);
- N_free_les(sples);
-
- G_message
- ("\t * testing pcg solver with N_ROWSCALE_EUKLIDNORM_PRECONDITION\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
- N_print_les(les);
- N_solver_pcg(sples, 100, 0.1e-8, N_ROWSCALE_EUKLIDNORM_PRECONDITION);
- N_print_les(sples);
-
- N_free_les(les);
- N_free_les(sples);
-
- G_message
- ("\t * testing pcg solver with N_ROWSCALE_ABSSUMNORM_PRECONDITION\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_pcg(les, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
- N_print_les(les);
- N_solver_pcg(sples, 100, 0.1e-8, N_ROWSCALE_ABSSUMNORM_PRECONDITION);
- N_print_les(sples);
-
- N_free_les(les);
- N_free_les(sples);
-
-
- G_message("\t * testing bicgstab solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
- sples = create_sparse_les(TEST_N_NUM_ROWS);
-
- N_solver_bicgstab(les, 100, 0.1e-8);
- /*N_print_les(les); */
- N_solver_bicgstab(sples, 100, 0.1e-8);
- /*N_print_les(sples); */
-
- N_free_les(les);
- N_free_les(sples);
-
- G_message("\t * testing gauss elimination solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
-
- /*GAUSS*/ N_solver_gauss(les);
- N_print_les(les);
-
- N_free_les(les);
-
- G_message("\t * testing lu decomposition solver\n");
-
- les = create_normal_les(TEST_N_NUM_ROWS);
-
- /*LU*/ N_solver_lu(les);
- N_print_les(les);
-
- N_free_les(les);
-
- les = create_normal_les(TEST_N_NUM_ROWS);
-
- /*cholesky */ N_solver_cholesky(les);
- N_print_les(les);
-
- N_free_les(les);
-
-
- return 0;
-}
Modified: grass/branches/develbranch_6/lib/gpde/test/test_tools.c
===================================================================
--- grass/branches/develbranch_6/lib/gpde/test/test_tools.c 2009-08-27 20:50:27 UTC (rev 38888)
+++ grass/branches/develbranch_6/lib/gpde/test/test_tools.c 2009-08-27 20:52:10 UTC (rev 38889)
@@ -48,149 +48,139 @@
/* *************************************************************** */
int test_mean_calc(void)
{
- double a, b, mean_n, mean, vector, distance, D, weight;
- double v[2];
- double array[10] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
- int i;
- int sum = 0;
- char buff1[10];
+ double a, b, mean_n, mean, vector, distance, D, weight;
+ double v[2];
+ double array[10] = {0,0,0,0,0,0,0,0,0,0};
+ int i;
+ int sum = 0;
+ char buff1[10];
- for (i = 0; i < 10; i++)
- array[i] += i;
+ for(i = 0; i < 10; i++)
+ array[i] += i;
- a = 1.0 / 3.0;
- b = 3.0;
- v[0] = a;
- v[1] = b;
+ a = 1.0/3.0;
+ b = 3.0;
+ v[0] = a; v[1] = b;
- /*arith mean */
- mean = N_calc_arith_mean(a, b);
- G_message("N_calc_arith_mean: calc a %g and b %g = %12.18lf", a, b, mean);
- mean_n = N_calc_arith_mean_n(v, 2);
- G_message("N_calc_arith_mean_n: calc a %g and b %g = %12.18lf", v[0],
- v[1], mean_n);
- if (mean != mean_n)
- sum++;
+ /*arith mean*/
+ mean = N_calc_arith_mean(a, b);
+ G_message("N_calc_arith_mean: calc a %g and b %g = %12.18lf", a, b, mean);
+ mean_n = N_calc_arith_mean_n(v, 2);
+ G_message("N_calc_arith_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1], mean_n);
+ if(mean != mean_n)
+ sum++;
- /*geom mean */
- mean = N_calc_geom_mean(a, b);
- G_message("N_calc_geom_mean: calc a %g and b %g = %12.18lf", a, b, mean);
- mean_n = N_calc_geom_mean_n(v, 2);
- G_message("N_calc_geom_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1],
- mean_n);
- if (mean != mean_n)
- sum++;
+ /*geom mean*/
+ mean = N_calc_geom_mean(a, b);
+ G_message("N_calc_geom_mean: calc a %g and b %g = %12.18lf", a, b, mean);
+ mean_n = N_calc_geom_mean_n(v, 2);
+ G_message("N_calc_geom_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1], mean_n);
+ if(mean != mean_n)
+ sum++;
- /*harmonic mean */
- mean = N_calc_harmonic_mean(a, b);
- G_message("N_calc_harmonic_mean: calc a %g and b %g = %12.18lf", a, b,
- mean);
- mean_n = N_calc_harmonic_mean_n(v, 2);
- G_message("N_calc_harmonic_mean_n: calc a %g and b %g = %12.18lf", v[0],
- v[1], mean_n);
- if (mean != mean_n)
- sum++;
- /*null test */
- a = 2;
- b = 0;
- v[0] = a;
- v[1] = b;
- mean = N_calc_harmonic_mean(a, b);
- G_message("N_calc_harmonic_mean: calc a %g and b %g = %12.18lf", a, b,
- mean);
- mean_n = N_calc_harmonic_mean_n(v, 2);
- G_message("N_calc_harmonic_mean_n: calc a %g and b %g = %12.18lf", v[0],
- v[1], mean_n);
- if (mean != mean_n)
- sum++;
+ /*harmonic mean*/
+ mean = N_calc_harmonic_mean(a, b);
+ G_message("N_calc_harmonic_mean: calc a %g and b %g = %12.18lf", a, b, mean);
+ mean_n = N_calc_harmonic_mean_n(v, 2);
+ G_message("N_calc_harmonic_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1], mean_n);
+ if(mean != mean_n)
+ sum++;
+ /*null test*/
+ a = 2;
+ b = 0;
+ v[0] = a; v[1] = b;
+ mean = N_calc_harmonic_mean(a, b);
+ G_message("N_calc_harmonic_mean: calc a %g and b %g = %12.18lf", a, b, mean);
+ mean_n = N_calc_harmonic_mean_n(v, 2);
+ G_message("N_calc_harmonic_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1], mean_n);
+ if(mean != mean_n)
+ sum++;
+
+ /*quadratic mean*/
+ a = 1.0/3.0;
+ b = 3.0;
+ v[0] = a; v[1] = b;
- /*quadratic mean */
- a = 1.0 / 3.0;
- b = 3.0;
- v[0] = a;
- v[1] = b;
+ mean = N_calc_quad_mean(a, b);
+ G_message("N_calc_quad_mean: calc a %g and b %g = %12.18lf", a, b, mean);
+ mean_n = N_calc_quad_mean_n(v, 2);
+ G_message("N_calc_quad_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1], mean_n);
+ if(mean != mean_n)
+ sum++;
- mean = N_calc_quad_mean(a, b);
- G_message("N_calc_quad_mean: calc a %g and b %g = %12.18lf", a, b, mean);
- mean_n = N_calc_quad_mean_n(v, 2);
- G_message("N_calc_quad_mean_n: calc a %g and b %g = %12.18lf", v[0], v[1],
- mean_n);
- if (mean != mean_n)
- sum++;
+ /*Test the full upwind stabailization*/
+ vector= -0.000001;
+ distance = 20;
+ D = 0.000001;
- /*Test the full upwind stabailization */
- vector = -0.000001;
- distance = 20;
- D = 0.000001;
+ weight = N_full_upwinding(vector, distance, D);
+ G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
- weight = N_full_upwinding(vector, distance, D);
- G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
-
- if (weight != 0) {
+ if(weight != 0)
+ {
G_warning("Error detected in N_full_upwinding");
sum++;
- }
+ }
- vector = 0.000001;
+ vector= 0.000001;
- weight = N_full_upwinding(vector, distance, D);
- G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
- if (weight != 1) {
+ weight = N_full_upwinding(vector, distance, D);
+ G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
+ if(weight != 1)
+ {
G_warning("Error detected in N_full_upwinding");
sum++;
- }
+ }
- D = 0.0;
+ D = 0.0;
- weight = N_full_upwinding(vector, distance, D);
- G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
- if (weight != 0.5) {
+ weight = N_full_upwinding(vector, distance, D);
+ G_message("N_full_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
+ if(weight != 0.5)
+ {
G_warning("Error detected in N_full_upwinding");
sum++;
- }
+ }
- /*Test the exponential upwind stabailization */
- vector = -0.000001;
- distance = 20;
- D = 0.000001;
+ /*Test the exponential upwind stabailization*/
+ vector= -0.000001;
+ distance = 20;
+ D = 0.000001;
- weight = N_exp_upwinding(vector, distance, D);
- G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
- sprintf(buff1, "%1.2lf", weight);
- sscanf(buff1, "%lf", &weight);
+ weight = N_exp_upwinding(vector, distance, D);
+ G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
+ sprintf(buff1, "%1.2lf", weight);
+ sscanf(buff1, "%lf", &weight);
- if (weight != 0.05) {
+ if(weight != 0.05)
+ {
G_warning("Error detected in N_exp_upwinding");
sum++;
- }
+ }
- vector = 0.000001;
+ vector= 0.000001;
- weight = N_exp_upwinding(vector, distance, D);
- G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
- sprintf(buff1, "%1.2lf", weight);
- sscanf(buff1, "%lf", &weight);
+ weight = N_exp_upwinding(vector, distance, D);
+ G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
+ sprintf(buff1, "%1.2lf", weight);
+ sscanf(buff1, "%lf", &weight);
- if (weight != 0.95) {
+ if(weight != 0.95)
+ {
G_warning("Error detected in N_exp_upwinding");
sum++;
- }
+ }
- D = 0.0;
+ D = 0.0;
- weight = N_exp_upwinding(vector, distance, D);
- G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n",
- vector, distance, D, weight);
- if (weight != 0.5) {
+ weight = N_exp_upwinding(vector, distance, D);
+ G_message("N_exp_upwinding: vector %g distance %g D %g weight %g\n", vector, distance, D, weight);
+ if(weight != 0.5)
+ {
G_warning("Error detected in N_exp_upwinding");
sum++;
- }
+ }
- return sum;
+ return sum;
}
More information about the grass-commit
mailing list