24#include <grass/gmath.h>
25#include <grass/glocale.h>
27static G_math_spvector **create_diag_precond_matrix(
double **A,
28 G_math_spvector **Asp,
30static int solver_pcg(
double **A, G_math_spvector **Asp,
double *
x,
double *
b,
31 int rows,
int maxit,
double err,
int prec,
int has_band,
33static int solver_cg(
double **A, G_math_spvector **Asp,
double *
x,
double *
b,
34 int rows,
int maxit,
double err,
int has_band,
36static int solver_bicgstab(
double **A, G_math_spvector **Asp,
double *
x,
37 double *
b,
int rows,
int maxit,
double err);
69 return solver_pcg(A,
NULL,
x,
b, rows, maxit,
err, prec, 0, 0);
102 int bandwidth,
int maxit,
double err,
int prec)
104 G_fatal_error(
"Preconditioning of band matrics is not implemented yet");
105 return solver_pcg(A,
NULL,
x,
b, rows, maxit,
err, prec, 1, bandwidth);
134 int rows,
int maxit,
double err,
int prec)
137 return solver_pcg(
NULL, Asp,
x,
b, rows, maxit,
err, prec, 0, 0);
140int solver_pcg(
double **A, G_math_spvector **Asp,
double *
x,
double *
b,
141 int rows,
int maxit,
double err,
int prec,
int has_band,
152 double a0 = 0, a1 = 0, mygamma, tmp = 0;
170 M = create_diag_precond_matrix(A, Asp, rows, prec);
189#pragma omp for schedule(static) private(i) reduction(+ : s)
190 for (i = 0; i < rows; i++) {
201 for (m = 0; m < maxit; m++) {
202#pragma omp parallel default(shared)
212#pragma omp for schedule(static) private(i) reduction(+ : s)
213 for (i = 0; i < rows; i++) {
245#pragma omp for schedule(static) private(i) reduction(+ : s)
246 for (i = 0; i < rows; i++) {
258 if (a1 < 0 || a1 == 0 || a1 > 0) {
262 G_warning(_(
"Unable to solve the linear equation system"));
270 G_message(_(
"Sparse PCG -- iteration %i error %g\n"), m, a0);
272 G_message(_(
"PCG -- iteration %i error %g\n"), m, a0);
274 if (error_break == 1) {
322 return solver_cg(A,
NULL,
x,
b, rows, maxit,
err, 0, 0);
351 int bandwidth,
int maxit,
double err)
353 return solver_cg(A,
NULL,
x,
b, rows, maxit,
err, 1, bandwidth);
381 int rows,
int maxit,
double err)
383 return solver_cg(
NULL, Asp,
x,
b, rows, maxit,
err, 0, 0);
386int solver_cg(
double **A, G_math_spvector **Asp,
double *
x,
double *
b,
int rows,
387 int maxit,
double err,
int has_band,
int bandwidth)
397 double a0 = 0, a1 = 0, mygamma, tmp = 0;
426#pragma omp for schedule(static) private(i) reduction(+ : s)
427 for (i = 0; i < rows; i++) {
438 for (m = 0; m < maxit; m++) {
439#pragma omp parallel default(shared)
449#pragma omp for schedule(static) private(i) reduction(+ : s)
450 for (i = 0; i < rows; i++) {
479#pragma omp for schedule(static) private(i) reduction(+ : s)
480 for (i = 0; i < rows; i++) {
492 if (a1 < 0 || a1 == 0 || a1 > 0) {
496 G_warning(_(
"Unable to solve the linear equation system"));
504 G_message(_(
"Sparse CG -- iteration %i error %g\n"), m, a0);
506 G_message(_(
"CG -- iteration %i error %g\n"), m, a0);
508 if (error_break == 1) {
551 int maxit,
double err)
553 return solver_bicgstab(A,
NULL,
x,
b, rows, maxit,
err);
581 int rows,
int maxit,
double err)
583 return solver_bicgstab(
NULL, Asp,
x,
b, rows, maxit,
err);
586int solver_bicgstab(
double **A, G_math_spvector **Asp,
double *
x,
double *
b,
587 int rows,
int maxit,
double err)
601 double s1 = 0.0, s2 = 0.0, s3 = 0.0;
603 double alpha = 0, beta = 0, omega, rr0 = 0, error;
637 for (m = 0; m < maxit; m++) {
639#pragma omp parallel default(shared)
647#pragma omp for schedule(static) private(i) reduction(+ : s1, s2, s3)
648 for (i = 0; i < rows; i++) {
658 if (error < 0 || error == 0 || error > 0) {
662 G_warning(_(
"Unable to solve the linear equation system"));
678#pragma omp for schedule(static) private(i) reduction(+ : s1, s2)
679 for (i = 0; i < rows; i++) {
694#pragma omp for schedule(static) private(i) reduction(+ : s1)
695 for (i = 0; i < rows; i++) {
701 beta = alpha / omega * s1 / rr0;
710 G_message(_(
"Sparse BiCGStab -- iteration %i error %g\n"), m,
713 G_message(_(
"BiCGStab -- iteration %i error %g\n"), m, error);
715 if (error_break == 1) {
746G_math_spvector **create_diag_precond_matrix(
double **A, G_math_spvector **Asp,
749 G_math_spvector **Msp;
751 unsigned int i, j, cols = (
unsigned int)rows;
760#pragma omp parallel for schedule(static) private(i, j, sum) \
761 shared(A, Msp, rows, cols, prec)
762 for (i = 0; i < (
unsigned int)rows; i++) {
766 case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
768 for (j = 0; j < cols; j++)
769 sum += A[i][j] * A[i][j];
770 spvect->values[0] = 1.0 / sqrt(sum);
772 case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
774 for (j = 0; j < cols; j++)
775 sum += fabs(A[i][j]);
776 spvect->values[0] = 1.0 / (sum);
778 case G_MATH_DIAGONAL_PRECONDITION:
780 spvect->values[0] = 1.0 / A[i][i];
784 spvect->index[0] = i;
791#pragma omp parallel for schedule(static) private(i, j, sum) \
792 shared(Asp, Msp, rows, cols, prec)
793 for (i = 0; i < (
unsigned int)rows; i++) {
797 case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
799 for (j = 0; j < Asp[i]->cols; j++)
800 sum += Asp[i]->values[j] * Asp[i]->values[j];
801 spvect->values[0] = 1.0 / sqrt(sum);
803 case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
805 for (j = 0; j < Asp[i]->cols; j++)
806 sum += fabs(Asp[i]->values[j]);
807 spvect->values[0] = 1.0 / (sum);
809 case G_MATH_DIAGONAL_PRECONDITION:
811 for (j = 0; j < Asp[i]->cols; j++)
812 if (i == Asp[i]->index[j])
813 spvect->values[0] = 1.0 / Asp[i]->values[j];
817 spvect->index[0] = i;
void G_free(void *buf)
Free allocated memory.
void G_math_d_ax_by(double *x, double *y, double *z, double a, double b, int rows)
Scales vectors x and y with the scalars a and b and adds them.
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
void G_math_d_Ax(double **A, double *x, double *y, int rows, int cols)
Compute the matrix - vector product of matrix A and vector x.
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
void G_message(const char *msg,...)
Print a message to stderr.
void G_warning(const char *msg,...)
Print a warning message to stderr.
#define assert(condition)
int G_math_solver_pcg(double **A, double *x, double *b, int rows, int maxit, double err, int prec)
The iterative preconditioned conjugate gradients solver for symmetric positive definite matrices.
int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit, double err)
The iterative conjugate gradients solver for symmetric positive definite matrices.
int G_math_solver_bicgstab(double **A, double *x, double *b, int rows, int maxit, double err)
The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices.
int G_math_solver_sparse_cg(G_math_spvector **Asp, double *x, double *b, int rows, int maxit, double err)
The iterative conjugate gradients solver for sparse symmetric positive definite matrices.
int G_math_solver_cg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit, double err)
The iterative conjugate gradients solver for symmetric positive definite band matrices.
int G_math_solver_pcg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit, double err, int prec)
The iterative preconditioned conjugate gradients solver for symmetric positive definite band matrices...
int G_math_solver_sparse_bicgstab(G_math_spvector **Asp, double *x, double *b, int rows, int maxit, double err)
The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices.
int G_math_solver_sparse_pcg(G_math_spvector **Asp, double *x, double *b, int rows, int maxit, double err, int prec)
The iterative preconditioned conjugate gradients solver for sparse symmetric positive definite matric...
void G_math_Ax_sparse(G_math_spvector **Asp, double *x, double *y, int rows)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
void G_math_free_spmatrix(G_math_spvector **Asp, int rows)
Release the memory of the sparse matrix.
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse matrix at position row.
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
void G_math_Ax_sband(double **A, double *x, double *y, int rows, int bandwidth)
Compute the matrix - vector product of symmetric band matrix A and vector x.