[GRASS-SVN] r42025 - in grass/trunk: include lib/gmath
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Apr 24 19:05:11 EDT 2010
Author: huhabla
Date: 2010-04-24 19:05:10 -0400 (Sat, 24 Apr 2010)
New Revision: 42025
Modified:
grass/trunk/include/gmath.h
grass/trunk/lib/gmath/solvers_direct.c
Log:
Removed pivot creation. Better function naming.
Modified: grass/trunk/include/gmath.h
===================================================================
--- grass/trunk/include/gmath.h 2010-04-24 22:57:31 UTC (rev 42024)
+++ grass/trunk/include/gmath.h 2010-04-24 23:05:10 UTC (rev 42025)
@@ -145,6 +145,7 @@
extern int G_math_solver_jacobi(double **, double *, double *, int , int , double , double );
extern int G_math_solver_gs(double **, double *, double *, int , int , double , double );
extern int G_math_solver_pcg(double **, double *, double *, int , int , double , int );
+
extern int G_math_solver_cg(double **, double *, double *, int , int , double );
extern int G_math_solver_bicgstab(double **, double *, double *, int , int , double );
extern int G_math_solver_sparse_jacobi(G_math_spvector **, double *, double *, int , int , double , double );
@@ -156,11 +157,9 @@
extern void G_math_gauss_elimination(double **, double *, int );
extern void G_math_lu_decomposition(double **, double *, int );
extern int G_math_cholesky_decomposition(double **, int , int );
-extern void G_math_backward_solving(double **, double *, double *, int );
-extern void G_math_forward_solving(double **, double *, double *, int );
-extern int G_math_pivot_create(double **, double *, int , int );
+extern void G_math_backward_substitution(double **, double *, double *, int );
+extern void G_math_forward_substitution(double **, double *, double *, int );
-
/*BLAS like level 1,2 and 3 functions*/
/*level 1 vector - vector grass implementation with OpenMP thread support*/
Modified: grass/trunk/lib/gmath/solvers_direct.c
===================================================================
--- grass/trunk/lib/gmath/solvers_direct.c 2010-04-24 22:57:31 UTC (rev 42024)
+++ grass/trunk/lib/gmath/solvers_direct.c 2010-04-24 23:05:10 UTC (rev 42025)
@@ -45,7 +45,7 @@
G_message(_("Starting direct gauss elimination solver"));
G_math_gauss_elimination(A, b, rows);
- G_math_backward_solving(A, x, b, rows);
+ G_math_backward_substitution(A, x, b, rows);
return 1;
}
@@ -88,7 +88,7 @@
#pragma omp single
{
- G_math_forward_solving(A, b, b, rows);
+ G_math_forward_substitution(A, b, b, rows);
}
#pragma omp for schedule (static) private(i)
@@ -98,7 +98,7 @@
#pragma omp single
{
- G_math_backward_solving(A, x, b, rows);
+ G_math_backward_substitution(A, x, b, rows);
}
}
@@ -134,8 +134,8 @@
return -2;
}
- G_math_forward_solving(A, b, b, rows);
- G_math_backward_solving(A, x, b, rows);
+ G_math_forward_substitution(A, b, b, rows);
+ G_math_backward_substitution(A, x, b, rows);
return 1;
}
@@ -158,10 +158,6 @@
double tmpval = 0.0;
- /*compute the pivot -- commented out, because its meaningless
- to compute it only nth times. */
- /*G_math_pivot_create(A, b, rows, 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++) {
@@ -194,10 +190,6 @@
int i, j, k;
- /*compute the pivot -- commented out, because its meaningless
- to compute it only nth times. */
- /*G_math_pivot_create(A, b, rows, 0); */
-
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++) {
@@ -284,7 +276,7 @@
}
/*!
- * \brief backward solving
+ * \brief backward substitution
*
* \param A double **
* \param x double *
@@ -293,7 +285,7 @@
* \return void
*
* */
-void G_math_backward_solving(double **A, double *x, double *b, int rows)
+void G_math_backward_substitution(double **A, double *x, double *b, int rows)
{
int i, j;
@@ -308,7 +300,7 @@
}
/*!
- * \brief forward solving
+ * \brief forward substitution
*
* \param A double **
* \param x double *
@@ -317,7 +309,7 @@
* \return void
*
* */
-void G_math_forward_solving(double **A, double *x, double *b, int rows)
+void G_math_forward_substitution(double **A, double *x, double *b, int rows)
{
int i, j;
@@ -333,84 +325,3 @@
return;
}
-
-
-/*!
- * \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 A double ** - a quadratic matrix
- * \param b double * - the right hand vector, if not available set it to NULL
- * \param rows int
- * \param start int -- the row
- * \return int - the number of swapped rows
- *
- *
- * */
-int G_math_pivot_create(double **A, double *b, int rows, int start)
-{
- 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;
-
- link = G_alloc_vector(rows);
-
- G_debug(2, "G_math_pivot_create: swap rows if needed");
- for (i = start; i < rows; i++) {
- s = 0.0;
- for (k = i + 1; k < rows; k++) {
- s += fabs(A[i][k]);
- }
- max = fabs(A[i][i]) / s;
- number = i;
- for (j = i + 1; j < rows; j++) {
- s = 0.0;
- for (k = j; k < rows; k++) {
- s += fabs(A[j][k]);
- }
- /*search for the pivot element */
- if (max < fabs(A[j][i]) / s) {
- max = fabs(A[j][i] / s);
- 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);
-
- if (b != NULL) {
- tmpval = b[number];
- b[number] = b[i];
- b[i] = tmpval;
- }
- G_math_d_copy(A[number], link, rows);
- G_math_d_copy(A[i], A[number], rows);
- G_math_d_copy(link, A[i], rows);
- num++;
- }
- }
-
- G_free_vector(link);
-
- return num;
-}
More information about the grass-commit
mailing list