23#include <grass/gmath.h>
24#include <grass/glocale.h>
44 G_message(_(
"Starting direct gauss elimination solver"));
72 G_message(_(
"Starting direct lu decomposition solver"));
82#pragma omp for schedule(static) private(i)
83 for (i = 0; i < rows; i++) {
93#pragma omp for schedule(static) private(i)
94 for (i = 0; i < rows; i++) {
128 G_message(_(
"Starting cholesky decomposition solver"));
131 G_warning(_(
"Unable to solve the linear equation system"));
159 for (k = 0; k < rows - 1; k++) {
160#pragma omp parallel for schedule(static) private(i, j, tmpval) \
161 shared(k, A, b, rows)
162 for (i = k + 1; i < rows; i++) {
163 tmpval = A[i][k] / A[k][k];
164 b[i] =
b[i] - tmpval *
b[k];
165 for (j = k + 1; j < rows; j++) {
166 A[i][j] = A[i][j] - tmpval * A[k][j];
191 for (k = 0; k < rows - 1; k++) {
192#pragma omp parallel for schedule(static) private(i, j) shared(k, A, rows)
193 for (i = k + 1; i < rows; i++) {
194 A[i][k] = A[i][k] / A[k][k];
195 for (j = k + 1; j < rows; j++) {
196 A[i][j] = A[i][j] - A[i][k] * A[k][j];
220 int i = 0, j = 0, k = 0;
233 for (k = 0; k < rows; k++) {
234#pragma omp parallel for schedule (static) private(i, j, sum_2) shared(A, k) reduction(+:sum_1)
235 for (j = 0; j < k; j++) {
236 sum_1 += A[k][j] * A[k][j];
239 if (0 > (A[k][k] - sum_1)) {
240 G_warning(
"Matrix is not positive definite. break.");
243 A[k][k] = sqrt(A[k][k] - sum_1);
246 if ((k + bandwidth) > rows) {
250 colsize = k + bandwidth;
253#pragma omp parallel for schedule(static) private(i, j, sum_2) \
254 shared(A, k, sum_1, colsize)
256 for (i = k + 1; i < colsize; i++) {
258 for (j = 0; j < k; j++) {
259 sum_2 += A[i][j] * A[k][j];
261 A[i][k] = (A[i][k] - sum_2) / A[k][k];
265#pragma omp parallel for schedule(static) private(i, k) shared(A, rows)
266 for (k = 0; k < rows; k++) {
267 for (i = k + 1; i < rows; i++) {
289 for (i = rows - 1; i >= 0; i--) {
290 for (j = i + 1; j < rows; j++) {
291 b[i] =
b[i] - A[i][j] *
x[j];
293 x[i] = (
b[i]) / A[i][i];
315 for (i = 0; i < rows; i++) {
317 for (j = 0; j < i; j++) {
318 tmpval += A[i][j] *
x[j];
320 x[i] = (
b[i] - tmpval) / A[i][i];
void G_free(void *buf)
Free allocated memory.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
void G_message(const char *msg,...)
Print a message to stderr.
void G_warning(const char *msg,...)
Print a warning message to stderr.
int G_math_solver_lu(double **A, double *x, double *b, int rows)
The LU solver for quardatic matrices.
void G_math_lu_decomposition(double **A, double *b UNUSED, int rows)
lu decomposition
int G_math_solver_gauss(double **A, double *x, double *b, int rows)
The gauss elimination solver for quardatic matrices.
void G_math_forward_substitution(double **A, double *x, double *b, int rows)
forward substitution
int G_math_solver_cholesky(double **A, double *x, double *b, int bandwidth, int rows)
The choleksy decomposition solver for quardatic, symmetric positive definite matrices.
void G_math_backward_substitution(double **A, double *x, double *b, int rows)
backward substitution
void G_math_gauss_elimination(double **A, double *b, int rows)
Gauss elimination.
int G_math_cholesky_decomposition(double **A, int rows, int bandwidth)
cholesky decomposition for symmetric, positive definite matrices with bandwidth optimization