GRASS Programmer's Manual
6.4.2(2012)
|
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 }