GRASS Programmer's Manual  6.4.2(2012)
N_solvers.c
Go to the documentation of this file.
00001 
00002 /*****************************************************************************
00003 *
00004 * MODULE:       Grass PDE Numerical Library
00005 * AUTHOR(S):    Soeren Gebbert, Berlin (GER) Dec 2006
00006 *               soerengebbert <at> gmx <dot> de
00007 *               
00008 * PURPOSE:      direkt linear equation system solvers
00009 *               part of the gpde library
00010 *               
00011 * COPYRIGHT:    (C) 2000 by the GRASS Development Team
00012 *
00013 *               This program is free software under the GNU General Public
00014 *               License (>=v2). Read the file COPYING that comes with GRASS
00015 *               for details.
00016 *
00017 *****************************************************************************/
00018 
00019 #include <math.h>
00020 #include <unistd.h>
00021 #include <stdio.h>
00022 #include <string.h>
00023 #include "grass/N_pde.h"
00024 #include "solvers_local_proto.h"
00025 
00026 /*prototypes */
00027 static void gauss_elimination(double **A, double *b, int rows);
00028 static void lu_decomposition(double **A, int rows);
00029 static int cholesky_decomposition(double **A, int rows);
00030 static void backward_solving(double **A, double *x, double *b, int rows);
00031 static void forward_solving(double **A, double *x, double *b, int rows);
00032 
00033 /***********************************************************
00034  * GAUSS elimination solver for Ax = b *********************
00035  * ********************************************************/
00046 int N_solver_gauss(N_les * les)
00047 {
00048 
00049     if (les->type != N_NORMAL_LES) {
00050         G_warning(_("The gauss elimination solver does not work with sparse matrices"));
00051         return 0;
00052     }
00053 
00054     if (les->quad != 1) {
00055         G_fatal_error(_("The linear equation system is not quadratic"));
00056         return 0;
00057     }
00058 
00059 
00060     G_message(_("Starting direct gauss elimination solver"));
00061 
00062     N_les_pivot_create(les);
00063     gauss_elimination(les->A, les->b, les->rows);
00064     backward_solving(les->A, les->x, les->b, les->rows);
00065 
00066     return 1;
00067 }
00068 
00069 /***********************************************************
00070  * LU solver for Ax = b ************************************
00071  * ********************************************************/
00082 int N_solver_lu(N_les * les)
00083 {
00084     int i;
00085     double *c, *tmpv;
00086 
00087     if (les->type != N_NORMAL_LES) {
00088         G_warning(_("The lu solver does not work with sparse matrices"));
00089         return 0;
00090     }
00091 
00092     if (les->quad != 1) {
00093         G_warning(_("The linear equation system is not quadratic"));
00094         return -1;
00095     }
00096 
00097 
00098     G_message(_("Starting direct lu decomposition solver"));
00099 
00100     tmpv = vectmem(les->rows);
00101     c = vectmem(les->rows);
00102 
00103     N_les_pivot_create(les);
00104     lu_decomposition(les->A, les->rows);
00105 
00106 #pragma omp parallel
00107     {
00108 
00109 #pragma omp for  schedule (static) private(i)
00110         for (i = 0; i < les->rows; i++) {
00111             tmpv[i] = les->A[i][i];
00112             les->A[i][i] = 1;
00113         }
00114 
00115 #pragma omp single
00116         {
00117             forward_solving(les->A, les->b, les->b, les->rows);
00118         }
00119 
00120 #pragma omp for  schedule (static) private(i)
00121         for (i = 0; i < les->rows; i++) {
00122             les->A[i][i] = tmpv[i];
00123         }
00124 
00125 #pragma omp single
00126         {
00127             backward_solving(les->A, les->x, les->b, les->rows);
00128         }
00129     }
00130 
00131     G_free(c);
00132     G_free(tmpv);
00133 
00134 
00135     return 1;
00136 }
00137 
00138 /***********************************************************
00139  * cholesky solver for Ax = b ******************************
00140  * ********************************************************/
00152 int N_solver_cholesky(N_les * les)
00153 {
00154     if (les->type != N_NORMAL_LES) {
00155         G_warning(_("The cholesky solver does not work with sparse matrices"));
00156         return 0;
00157     }
00158 
00159     if (les->quad != 1) {
00160         G_warning(_("The linear equation system is not quadratic"));
00161         return -1;
00162     }
00163 
00164     /* check for symmetry */
00165     if (check_symmetry(les) != 1) {
00166         G_warning(_("Matrix is not symmetric!"));
00167         return -3;
00168     }
00169 
00170     G_message(_("Starting cholesky decomposition solver"));
00171 
00172     if (cholesky_decomposition(les->A, les->rows) != 1) {
00173         G_warning(_("Unable to solve the linear equation system"));
00174         return -2;
00175     }
00176 
00177     forward_solving(les->A, les->b, les->b, les->rows);
00178     backward_solving(les->A, les->x, les->b, les->rows);
00179 
00180     return 1;
00181 }
00182 
00183 
00184 /***********************************************************
00185  * gauss elimination ***************************************
00186  * ********************************************************/
00198 void gauss_elimination(double **A, double *b, int rows)
00199 {
00200     int i, j, k;
00201     double tmpval = 0.0;
00202 
00203     for (k = 0; k < rows - 1; k++) {
00204 #pragma omp parallel for schedule (static) private(i, j, tmpval) shared(k, A, b, rows)
00205         for (i = k + 1; i < rows; i++) {
00206             tmpval = A[i][k] / A[k][k];
00207             b[i] = b[i] - tmpval * b[k];
00208             for (j = k + 1; j < rows; j++) {
00209                 A[i][j] = A[i][j] - tmpval * A[k][j];
00210             }
00211         }
00212     }
00213 
00214     return;
00215 }
00216 
00217 /***********************************************************
00218  * lu decomposition ****************************************
00219  * ********************************************************/
00230 void lu_decomposition(double **A, int rows)
00231 {
00232 
00233     int i, j, k;
00234 
00235     for (k = 0; k < rows - 1; k++) {
00236 #pragma omp parallel for schedule (static) private(i, j) shared(k, A, rows)
00237         for (i = k + 1; i < rows; i++) {
00238             A[i][k] = A[i][k] / A[k][k];
00239             for (j = k + 1; j < rows; j++) {
00240                 A[i][j] = A[i][j] - A[i][k] * A[k][j];
00241             }
00242         }
00243     }
00244 
00245     return;
00246 }
00247 
00248 /***********************************************************
00249  * cholesky decomposition **********************************
00250  * ********************************************************/
00262 int cholesky_decomposition(double **A, int rows)
00263 {
00264 
00265     int i, j, k;
00266     double sum_1 = 0.0;
00267     double sum_2 = 0.0;
00268     int error = 0;
00269 
00270 
00271     for (k = 0; k < rows; k++) {
00272 #pragma omp parallel private(i, j, sum_2) shared(A, rows, sum_1)
00273         {
00274 #pragma omp for schedule (static) private(j) reduction(+:sum_1)
00275             for (j = 0; j < k; j++) {
00276                 sum_1 += A[k][j] * A[k][j];
00277             }
00278 #pragma omp single
00279             {
00280                 if ((A[k][k] - sum_1) < 0) {
00281                     error++;    /*not allowed to exit the OpenMP region */
00282                 }
00283                 A[k][k] = sqrt(A[k][k] - sum_1);
00284                 sum_1 = 0.0;
00285             }
00286 #pragma omp for schedule (static) private(i, j, sum_2)
00287             for (i = k + 1; i < rows; i++) {
00288                 sum_2 = 0.0;
00289                 for (j = 0; j < k; j++) {
00290                     sum_2 += A[i][j] * A[k][j];
00291                 }
00292                 A[i][k] = (A[i][k] - sum_2) / A[k][k];
00293             }
00294         }
00295     }
00296 
00297     /*we need to copy the lower triangle matrix to the upper trianle */
00298 #pragma omp parallel for schedule (static) private(i, k) shared(A, rows)
00299     for (k = 0; k < rows; k++) {
00300         for (i = k + 1; i < rows; i++) {
00301             A[k][i] = A[i][k];
00302         }
00303     }
00304 
00305     if (error > 0) {
00306         G_warning("Matrix is not positive definite");
00307         return -1;
00308     }
00309 
00310     return 1;
00311 }
00312 
00313 
00314 /***********************************************************
00315  * backward solving ****************************************
00316  * ********************************************************/
00327 void backward_solving(double **A, double *x, double *b, int rows)
00328 {
00329     int i, j;
00330     double tmpval = 0.0;
00331 
00332     for (i = rows - 1; i >= 0; i--) {
00333         tmpval = 0;
00334         for (j = i + 1; j < rows; j++) {
00335             /*tmpval += A[i][j] * x[j]; */
00336             b[i] = b[i] - A[i][j] * x[j];
00337         }
00338         /*x[i] = (b[i] - tmpval) / A[i][i]; */
00339         x[i] = (b[i]) / A[i][i];
00340     }
00341 
00342     return;
00343 }
00344 
00345 
00346 /***********************************************************
00347  * forward solving *****************************************
00348  * ********************************************************/
00359 void forward_solving(double **A, double *x, double *b, int rows)
00360 {
00361     int i, j;
00362     double tmpval = 0.0;
00363 
00364     for (i = 0; i < rows; i++) {
00365         tmpval = 0;
00366         for (j = 0; j < i; j++) {
00367             tmpval += A[i][j] * x[j];
00368         }
00369         x[i] = (b[i] - tmpval) / A[i][i];
00370     }
00371 
00372     return;
00373 }
00374 
00375 
00376 /* ******************************************************* *
00377  * ***** solving a tridiagonal eq system ***************** *
00378  * ******************************************************* */
00379 void thomalg(double **M, double *V, int rows)
00380 {
00381     double *Vtmp;
00382     double *g;
00383     double b;
00384     int i;
00385 
00386     Vtmp = vectmem(rows);
00387     g = vectmem(rows);
00388 
00389     for (i = 0; i < rows; i++) {
00390         if (i == 0) {
00391             b = M[i][i];
00392             Vtmp[i] = V[i] / b;
00393         }
00394         else {
00395             b = M[i][i] - M[i][i - 1] * g[i - 1];
00396             Vtmp[i] = (V[i] - Vtmp[i - 1] * M[i][i - 1]) / b;
00397         }
00398         if (i < rows - 1) {
00399             g[i] = M[i][i + 1] / b;
00400         }
00401     }
00402 
00403     V[rows - 1] = Vtmp[rows - 1];
00404     for (i = rows - 2; i >= 0; i--) {
00405         V[i] = Vtmp[i] - g[i] * V[i + 1];
00406     }
00407 
00408     G_free(Vtmp);
00409     G_free(g);
00410 }
00411 
00412 
00413 /***********************************************************
00414  * vectmem *************************************************
00415  * ********************************************************/
00423 double *vectmem(int rows)
00424 {
00425     double *vector;
00426 
00427     vector = (double *)(G_calloc(rows, (sizeof(double))));
00428     return vector;
00429 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines