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: 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 }