23 #include "grass/N_pde.h"
27 static void gauss_elimination(
double **A,
double *
b,
int rows);
28 static void lu_decomposition(
double **A,
int rows);
29 static int cholesky_decomposition(
double **A,
int rows);
30 static void backward_solving(
double **A,
double *x,
double *
b,
int rows);
31 static void forward_solving(
double **A,
double *x,
double *
b,
int rows);
50 G_warning(_(
"The gauss elimination solver does not work with sparse matrices"));
55 G_fatal_error(_(
"The linear equation system is not quadratic"));
60 G_message(_(
"Starting direct gauss elimination solver"));
63 gauss_elimination(les->
A, les->
b, les->
rows);
64 backward_solving(les->
A, les->
x, les->
b, les->
rows);
88 G_warning(_(
"The lu solver does not work with sparse matrices"));
93 G_warning(_(
"The linear equation system is not quadratic"));
98 G_message(_(
"Starting direct lu decomposition solver"));
104 lu_decomposition(les->
A, les->
rows);
109 #pragma omp for schedule (static) private(i)
110 for (i = 0; i < les->
rows; i++) {
111 tmpv[i] = les->
A[i][i];
117 forward_solving(les->
A, les->
b, les->
b, les->
rows);
120 #pragma omp for schedule (static) private(i)
121 for (i = 0; i < les->
rows; i++) {
122 les->
A[i][i] = tmpv[i];
127 backward_solving(les->
A, les->
x, les->
b, les->
rows);
155 G_warning(_(
"The cholesky solver does not work with sparse matrices"));
159 if (les->
quad != 1) {
160 G_warning(_(
"The linear equation system is not quadratic"));
166 G_warning(_(
"Matrix is not symmetric!"));
170 G_message(_(
"Starting cholesky decomposition solver"));
172 if (cholesky_decomposition(les->
A, les->
rows) != 1) {
173 G_warning(_(
"Unable to solve the linear equation system"));
177 forward_solving(les->
A, les->
b, les->
b, les->
rows);
178 backward_solving(les->
A, les->
x, les->
b, les->
rows);
198 void gauss_elimination(
double **A,
double *
b,
int rows)
203 for (k = 0; k < rows - 1; k++) {
204 #pragma omp parallel for schedule (static) private(i, j, tmpval) shared(k, A, b, rows)
205 for (i = k + 1; i < rows; i++) {
206 tmpval = A[i][k] / A[k][k];
207 b[i] = b[i] - tmpval * b[k];
208 for (j = k + 1; j < rows; j++) {
209 A[i][j] = A[i][j] - tmpval * A[k][j];
230 void lu_decomposition(
double **A,
int rows)
235 for (k = 0; k < rows - 1; k++) {
236 #pragma omp parallel for schedule (static) private(i, j) shared(k, A, rows)
237 for (i = k + 1; i < rows; i++) {
238 A[i][k] = A[i][k] / A[k][k];
239 for (j = k + 1; j < rows; j++) {
240 A[i][j] = A[i][j] - A[i][k] * A[k][j];
262 int cholesky_decomposition(
double **A,
int rows)
271 for (k = 0; k < rows; k++) {
272 #pragma omp parallel private(i, j, sum_2) shared(A, rows, sum_1)
274 #pragma omp for schedule (static) private(j) reduction(+:sum_1)
275 for (j = 0; j < k; j++) {
276 sum_1 += A[k][j] * A[k][j];
280 if ((A[k][k] - sum_1) < 0) {
283 A[k][k] = sqrt(A[k][k] - sum_1);
286 #pragma omp for schedule (static) private(i, j, sum_2)
287 for (i = k + 1; i < rows; i++) {
289 for (j = 0; j < k; j++) {
290 sum_2 += A[i][j] * A[k][j];
292 A[i][k] = (A[i][k] - sum_2) / A[k][k];
298 #pragma omp parallel for schedule (static) private(i, k) shared(A, rows)
299 for (k = 0; k < rows; k++) {
300 for (i = k + 1; i < rows; i++) {
306 G_warning(
"Matrix is not positive definite");
327 void backward_solving(
double **A,
double *x,
double *b,
int rows)
332 for (i = rows - 1; i >= 0; i--) {
334 for (j = i + 1; j < rows; j++) {
336 b[i] = b[i] - A[i][j] * x[j];
339 x[i] = (b[i]) / A[i][i];
359 void forward_solving(
double **A,
double *x,
double *b,
int rows)
364 for (i = 0; i < rows; i++) {
366 for (j = 0; j < i; j++) {
367 tmpval += A[i][j] * x[j];
369 x[i] = (b[i] - tmpval) / A[i][i];
389 for (i = 0; i < rows; i++) {
395 b = M[i][i] - M[i][i - 1] * g[i - 1];
396 Vtmp[i] = (V[i] - Vtmp[i - 1] * M[i][i - 1]) / b;
399 g[i] = M[i][i + 1] /
b;
403 V[rows - 1] = Vtmp[rows - 1];
404 for (i = rows - 2; i >= 0; i--) {
405 V[i] = Vtmp[i] - g[i] * V[i + 1];
427 vector = (
double *)(G_calloc(rows, (
sizeof(
double))));
void G_free(void *buf)
Free allocated memory.
def error(msg)
Display an error message using g.message -e
int N_les_pivot_create(N_les *les)
Optimize the structure of the linear equation system with a common pivoting strategy.
int G_warning(const char *msg,...)
Print a warning message to stderr.
void G_message(const char *msg,...)
Print a message to stderr.
double * vectmem(int rows)
Allocate vector memory.
int N_solver_gauss(N_les *les)
The gauss elimination solver for quardatic matrices.
int N_solver_cholesky(N_les *les)
The choleksy decomposition solver for quardatic, symmetric positiv definite matrices.
int check_symmetry(N_les *L)
Check if the matrix in les is symmetric.
void thomalg(double **M, double *V, int rows)
int G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
The linear equation system (les) structure.
int N_solver_lu(N_les *les)
The LU solver for quardatic matrices.