23 #include "grass/N_pde.h"
27 static void scalar_product(
double *a,
double *
b,
double *scalar,
int rows);
28 static void sub_vectors(
double *source_a,
double *source_b,
double *result,
30 static void sub_vectors_scalar(
double *source_a,
double *source_b,
31 double *result,
double scalar_b,
int rows);
32 static void add_vectors(
double *source_a,
double *source_b,
double *result,
34 static void add_vectors_scalar(
double *source_a,
double *source_b,
35 double *result,
double scalar_b,
int rows);
36 static void add_vectors_scalar2(
double *source_a,
double *source_b,
37 double *result,
double scalar_a,
38 double scalar_b,
int rows);
39 static void scalar_vector_product(
double *a,
double *result,
double scalar,
41 static void sync_vectors(
double *source,
double *target,
int rows);
71 double a0 = 0, a1 = 0, mygamma, tmp = 0;
78 G_warning(_(
"The linear equation system is not quadratic"));
112 sub_vectors(b, v, r, rows);
117 #pragma omp for schedule (static) private(i) reduction(+:s)
118 for (i = 0; i < rows; i++) {
129 for (m = 0; m < maxit; m++) {
130 #pragma omp parallel default(shared)
140 #pragma omp for schedule (static) private(i) reduction(+:s)
141 for (i = 0; i < rows; i++) {
153 add_vectors_scalar(x, p, x, mygamma, rows);
162 sub_vectors(b, v, r, rows);
166 sub_vectors_scalar(r, v, r, mygamma, rows);
173 #pragma omp for schedule (static) private(i) reduction(+:s)
174 for (i = 0; i < rows; i++) {
186 if (a1 < 0 || a1 == 0 || a1 > 0) {
190 G_warning(_(
"Unable to solve the linear equation system"));
194 add_vectors_scalar(z, p, p, tmp, rows);
198 G_message(_(
"Sparse PCG -- iteration %i error %g\n"), m, a0);
200 G_message(_(
"PCG -- iteration %i error %g\n"), m, a0);
202 if (error_break == 1) {
249 double a0 = 0, a1 = 0, mygamma, tmp = 0;
255 G_warning(_(
"The linear equation system is not quadratic"));
261 G_warning(_(
"Matrix is not symmetric!"));
284 sub_vectors(b, v, r, rows);
285 sync_vectors(r, p, rows);
288 #pragma omp for schedule (static) private(i) reduction(+:s)
289 for (i = 0; i < rows; i++) {
300 for (m = 0; m < maxit; m++) {
301 #pragma omp parallel default(shared)
311 #pragma omp for schedule (static) private(i) reduction(+:s)
312 for (i = 0; i < rows; i++) {
324 add_vectors_scalar(x, p, x, mygamma, rows);
333 sub_vectors(b, v, r, rows);
336 sub_vectors_scalar(r, v, r, mygamma, rows);
340 #pragma omp for schedule (static) private(i) reduction(+:s)
341 for (i = 0; i < rows; i++) {
353 if (a1 < 0 || a1 == 0 || a1 > 0) {
357 G_warning(_(
"Unable to solve the linear equation system"));
361 add_vectors_scalar(r, p, p, tmp, rows);
365 G_message(_(
"Sparse CG -- iteration %i error %g\n"), m, a0);
367 G_message(_(
"CG -- iteration %i error %g\n"), m, a0);
369 if (error_break == 1) {
416 double s1 = 0.0, s2 = 0.0, s3 = 0.0;
417 double alpha = 0, beta = 0, omega, rr0 = 0, error;
423 G_warning(_(
"The linear equation system is not quadratic"));
446 sub_vectors(b, v, r, rows);
447 sync_vectors(r, r0, rows);
448 sync_vectors(r, p, rows);
456 for (m = 0; m < maxit; m++) {
458 #pragma omp parallel default(shared)
466 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3)
467 for (i = 0; i < rows; i++) {
477 if (error < 0 || error == 0 || error > 0) {
481 G_warning(_(
"Unable to solve the linear equation system"));
490 sub_vectors_scalar(r, v, s, alpha, rows);
497 #pragma omp for schedule (static) private(i) reduction(+:s1, s2)
498 for (i = 0; i < rows; i++) {
509 add_vectors_scalar2(p, s, r, alpha, omega, rows);
510 add_vectors(x, r, x, rows);
511 sub_vectors_scalar(s, t, r, omega, rows);
513 #pragma omp for schedule (static) private(i) reduction(+:s1)
514 for (i = 0; i < rows; i++) {
520 beta = alpha / omega * s1 / rr0;
524 sub_vectors_scalar(p, v, p, omega, rows);
525 add_vectors_scalar(r, p, p, beta, rows);
530 G_message(_(
"Sparse BiCGStab -- iteration %i error %g\n"), m,
533 G_message(_(
"BiCGStab -- iteration %i error %g\n"), m, error);
535 if (error_break == 1) {
568 void scalar_product(
double *a,
double *
b,
double *scalar,
int rows)
573 #pragma omp parallel for schedule (static) reduction(+:s)
574 for (i = 0; i < rows; i++) {
599 #pragma omp for schedule (static) private(i, j, tmp)
600 for (i = 0; i < L->
rows; i++) {
602 for (j = 0; j < L->
cols; j++) {
603 tmp += L->
A[i][j] * x[j];
627 #pragma omp for schedule (static) private(i, j, tmp)
628 for (i = 0; i < L->
rows; i++) {
630 for (j = 0; j < L->
Asp[i]->
cols; j++) {
653 add_vectors_scalar2(
double *a,
double *b,
double *result,
double scalar_a,
654 double scalar_b,
int rows)
658 #pragma omp for schedule (static)
659 for (i = 0; i < rows; i++) {
660 result[i] = scalar_a * a[i] + scalar_b * b[i];
680 add_vectors_scalar(
double *a,
double *b,
double *result,
double scalar_b,
685 #pragma omp for schedule (static)
686 for (i = 0; i < rows; i++) {
687 result[i] = a[i] + scalar_b * b[i];
707 sub_vectors_scalar(
double *a,
double *b,
double *result,
double scalar_b,
712 #pragma omp for schedule (static)
713 for (i = 0; i < rows; i++) {
714 result[i] = a[i] - scalar_b * b[i];
732 void add_vectors(
double *a,
double *b,
double *result,
int rows)
736 #pragma omp for schedule (static)
737 for (i = 0; i < rows; i++) {
738 result[i] = a[i] + b[i];
756 void sub_vectors(
double *a,
double *b,
double *result,
int rows)
760 #pragma omp for schedule (static)
761 for (i = 0; i < rows; i++) {
762 result[i] = a[i] - b[i];
780 void scalar_vector_product(
double *a,
double *result,
double scalar,
int rows)
784 #pragma omp for schedule (static)
785 for (i = 0; i < rows; i++) {
786 result[i] = scalar * a[i];
800 void sync_vectors(
double *source,
double *target,
int rows)
804 #pragma omp for schedule (static)
805 for (i = 0; i < rows; i++) {
806 target[i] = source[i];
832 G_warning(_(
"The linear equation system is not quadratic"));
836 G_debug(2,
"check_symmetry: Check if matrix is symmetric");
839 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
840 for (j = 0; j < L->
rows; j++) {
841 for (i = 1; i < L->
Asp[j]->
cols; i++) {
849 for (k = 1; k < L->
Asp[index]->
cols; k++) {
852 if ((value1 != value2)) {
853 if ((fabs((fabs(value1) - fabs(value2))) <
856 "check_symmetry: sparse matrix is unsymmetric, but within tolerance");
860 (
"Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf \ndifference = %12.18lf\nStop symmetry calculation.\n",
861 j, index, index, L->
Asp[index]->
index[k],
863 fabs(fabs(value1) - fabs(value2)));
873 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
874 for (i = 0; i < L->
rows; i++) {
875 for (j = i + 1; j < L->
rows; j++) {
876 if (L->
A[i][j] != L->
A[j][i]) {
877 if ((fabs(fabs(L->
A[i][j]) - fabs(L->
A[j][i])) <
880 "check_symmetry: matrix is unsymmetric, but within tolerance");
884 (
"Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf\ndifference = %12.18lf\nStop symmetry calculation.\n",
885 i, j, j, i, L->
A[i][j], L->
A[j][i],
886 fabs(fabs(L->
A[i][j]) - fabs(L->
A[j][i])));
920 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
921 for (i = 0; i < rows; i++) {
927 for (j = 0; j <
cols; j++)
928 sum += L->
A[i][j] * L->
A[i][j];
929 spvect->
values[0] = 1.0 / sqrt(sum);
933 for (j = 0; j <
cols; j++)
934 sum += fabs(L->
A[i][j]);
935 spvect->
values[0] = 1.0 / (sum);
938 spvect->
values[0] = 1.0 / L->
A[i][i];
941 spvect->
values[0] = 1.0 / L->
A[i][i];
945 spvect->
index[0] = i;
952 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
953 for (i = 0; i < rows; i++) {
959 for (j = 0; j < L->
Asp[i]->
cols; j++)
961 spvect->
values[0] = 1.0 / sqrt(sum);
965 for (j = 0; j < L->
Asp[i]->
cols; j++)
967 spvect->
values[0] = 1.0 / (sum);
976 spvect->
index[0] = i;
void G_free(void *buf)
Free allocated memory.
int N_solver_pcg(N_les *les, int maxit, double error, int prec)
The iterative preconditioned conjugate gradients solver for symmetric positive definite matrices...
The row vector of the sparse matrix.
void G_message(const char *msg,...)
Print a message to stderr.
N_spvector * N_alloc_spvector(int cols)
Allocate memory for a sparse vector.
#define N_DIAGONAL_PRECONDITION
double * vectmem(int rows)
Allocate vector memory.
void N_sparse_matrix_vector_product(N_les *les, double *source, double *result)
Calculates the matrix - vector product of sparse matrix L->Asp and vector x.
#define N_ROWSCALE_EUKLIDNORM_PRECONDITION
int N_solver_cg(N_les *les, int maxit, double error)
The iterative conjugate gradients solver for symmetric positive definite matrices.
N_les * N_alloc_les_A(int rows, int type)
Allocate memory for a quadratic linear equation system which includes the Matrix A.
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
int G_debug(int level, const char *msg,...)
Print debugging message.
int check_symmetry(N_les *L)
Check if the matrix in les is symmetric.
int N_solver_bicgstab(N_les *les, int maxit, double error)
The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices...
#define N_ROWSCALE_ABSSUMNORM_PRECONDITION
The linear equation system (les) structure.
N_les * N_create_diag_precond_matrix(N_les *les, int prec)
Compute a diagonal preconditioning matrix for krylov space solver.
void N_matrix_vector_product(N_les *les, double *source, double *result)
Calculates the matrix - vector product of matrix L->A and vector x.
int N_add_spvector_to_les(N_les *les, N_spvector *spvector, int row)
Adds a sparse vector to a sparse linear equation system at position row.