GRASS Programmer's Manual  6.4.2(2012)
N_solvers_classic_iter.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:      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 static int sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
00027                                const char *type);
00028 static int jacobi(double **M, double *b, double *x, int rows, int maxit,
00029                   double sor, double error);
00030 static int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
00031                         double sor, double error);
00032 
00033 /* ******************************************************* *
00034  * ******** overrelaxed jacobian ************************* *
00035  * ******************************************************* */
00053 int N_solver_jacobi(N_les * L, int maxit, double sor, double error)
00054 {
00055 
00056     if (L->quad != 1) {
00057         G_warning(_("The linear equation system is not quadratic"));
00058         return -1;
00059     }
00060 
00061     if (L->type == N_NORMAL_LES)
00062         return jacobi(L->A, L->b, L->x, L->rows, maxit, sor, error);
00063     else
00064         return sparse_jacobi_gauss(L, maxit, sor, error,
00065                                    N_SOLVER_ITERATIVE_JACOBI);
00066 }
00067 
00068 
00069 /* ******************************************************* *
00070  * ********* overrelaxed gauss seidel ******************** *
00071  * ******************************************************* */
00090 int N_solver_SOR(N_les * L, int maxit, double sor, double error)
00091 {
00092 
00093     if (L->quad != 1) {
00094         G_warning(_("The linear equation system is not quadratic"));
00095         return -1;
00096     }
00097 
00098     if (L->type == N_NORMAL_LES)
00099         return gauss_seidel(L->A, L->b, L->x, L->rows, maxit, sor, error);
00100     else
00101         return sparse_jacobi_gauss(L, maxit, sor, error,
00102                                    N_SOLVER_ITERATIVE_SOR);
00103 }
00104 
00105 /* ******************************************************* *
00106  * ****** sparse jacobi and SOR algorithm **************** *
00107  * ******************************************************* */
00108 int
00109 sparse_jacobi_gauss(N_les * L, int maxit, double sor, double error,
00110                     const char *type)
00111 {
00112     int i, j, k, rows, finished = 0;
00113     double *Enew, *x, *b;
00114     double E, err = 0;
00115 
00116     x = L->x;
00117     b = L->b;
00118     rows = L->rows;
00119 
00120     Enew = vectmem(rows);
00121 
00122     for (k = 0; k < maxit; k++) {
00123         err = 0;
00124         {
00125             if (k == 0) {
00126                 for (j = 0; j < rows; j++) {
00127                     Enew[j] = x[j];
00128                 }
00129             }
00130             for (i = 0; i < rows; i++) {
00131                 E = 0;
00132                 if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0) {
00133                     for (j = 0; j < L->Asp[i]->cols; j++) {
00134                         E += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]];
00135                     }
00136                 }
00137                 else {
00138                     for (j = 0; j < L->Asp[i]->cols; j++) {
00139                         E += L->Asp[i]->values[j] * Enew[L->Asp[i]->index[j]];
00140                     }
00141                 }
00142                 Enew[i] = x[i] - sor * (E - b[i]) / L->Asp[i]->values[0];
00143             }
00144             for (j = 0; j < rows; j++) {
00145                 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
00146 
00147                 x[j] = Enew[j];
00148             }
00149         }
00150 
00151         if (strcmp(type, N_SOLVER_ITERATIVE_JACOBI) == 0)
00152             G_message(_("sparse Jacobi -- iteration %5i error %g\n"), k, err);
00153         else if (strcmp(type, N_SOLVER_ITERATIVE_SOR) == 0)
00154             G_message(_("sparse SOR -- iteration %5i error %g\n"), k, err);
00155 
00156         if (err < error) {
00157             finished = 1;
00158             break;
00159         }
00160     }
00161 
00162     G_free(Enew);
00163 
00164     return finished;
00165 }
00166 
00167 /* ******************************************************* *
00168  * ******* direct jacobi ********************************* *
00169  * ******************************************************* */
00170 int jacobi(double **M, double *b, double *x, int rows, int maxit, double sor,
00171            double error)
00172 {
00173     int i, j, k;
00174     double *Enew;
00175     double E, err = 0;
00176 
00177     Enew = vectmem(rows);
00178 
00179     for (j = 0; j < rows; j++) {
00180         Enew[j] = x[j];
00181     }
00182 
00183     for (k = 0; k < maxit; k++) {
00184         for (i = 0; i < rows; i++) {
00185             E = 0;
00186             for (j = 0; j < rows; j++) {
00187                 E += M[i][j] * x[j];
00188             }
00189             Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
00190         }
00191         err = 0;
00192         for (j = 0; j < rows; j++) {
00193             err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
00194 
00195             x[j] = Enew[j];
00196         }
00197         G_message(_("Jacobi -- iteration %5i error %g\n"), k, err);
00198         if (err < error)
00199             break;
00200     }
00201 
00202     return 1;
00203 }
00204 
00205 /* ******************************************************* *
00206  * ******* direct gauss seidel *************************** *
00207  * ******************************************************* */
00208 int gauss_seidel(double **M, double *b, double *x, int rows, int maxit,
00209                  double sor, double error)
00210 {
00211     int i, j, k;
00212     double *Enew;
00213     double E, err = 0;
00214 
00215     Enew = vectmem(rows);
00216 
00217     for (j = 0; j < rows; j++) {
00218         Enew[j] = x[j];
00219     }
00220 
00221     for (k = 0; k < maxit; k++) {
00222         for (i = 0; i < rows; i++) {
00223             E = 0;
00224             for (j = 0; j < rows; j++) {
00225                 E += M[i][j] * Enew[j];
00226             }
00227             Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
00228         }
00229         err = 0;
00230         for (j = 0; j < rows; j++) {
00231             err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
00232 
00233             x[j] = Enew[j];
00234         }
00235         G_message(_("SOR -- iteration %5i error %g\n"), k, err);
00236         if (err < error)
00237             break;
00238     }
00239 
00240     return 1;
00241 }
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines