31 double *
b,
int rows,
int maxit,
double err,
int prec,
int has_band,
int bandwidth);
33 int rows,
int maxit,
double err,
int has_band,
int bandwidth);
35 double *
b,
int rows,
int maxit,
double err);
64 return solver_pcg(A,
NULL, x, b, rows, maxit, err, prec, 0, 0);
95 G_fatal_error(
"Preconditioning of band matrics is not implemented yet");
96 return solver_pcg(A,
NULL, x, b, rows, maxit, err, prec, 1, bandwidth);
123 int rows,
int maxit,
double err,
int prec)
126 return solver_pcg(
NULL, Asp, x, b, rows, maxit, err, prec, 0, 0);
130 int rows,
int maxit,
double err,
int prec,
int has_band,
int bandwidth)
140 double a0 = 0, a1 = 0, mygamma, tmp = 0;
158 M = create_diag_precond_matrix(A, Asp, rows, prec);
177 #pragma omp for schedule (static) private(i) reduction(+:s) 178 for (i = 0; i < rows; i++) {
189 for (m = 0; m < maxit; m++) {
190 #pragma omp parallel default(shared) 202 #pragma omp for schedule (static) private(i) reduction(+:s) 203 for (i = 0; i < rows; i++) {
236 #pragma omp for schedule (static) private(i) reduction(+:s) 237 for (i = 0; i < rows; i++) {
249 if (a1 < 0 || a1 == 0 || a1 > 0) {
254 (
"Unable to solve the linear equation system"));
262 G_message(
_(
"Sparse PCG -- iteration %i error %g\n"), m, a0);
264 G_message(
_(
"PCG -- iteration %i error %g\n"), m, a0);
266 if (error_break == 1) {
312 return solver_cg(A,
NULL, x, b, rows, maxit, err, 0, 0);
339 return solver_cg(A,
NULL, x, b, rows, maxit, err, 1, bandwidth);
365 int rows,
int maxit,
double err)
367 return solver_cg(
NULL, Asp, x, b, rows, maxit, err, 0, 0);
372 int rows,
int maxit,
double err,
int has_band,
int bandwidth)
382 double a0 = 0, a1 = 0, mygamma, tmp = 0;
411 #pragma omp for schedule (static) private(i) reduction(+:s) 412 for (i = 0; i < rows; i++) {
423 for (m = 0; m < maxit; m++) {
424 #pragma omp parallel default(shared) 434 #pragma omp for schedule (static) private(i) reduction(+:s) 435 for (i = 0; i < rows; i++) {
464 #pragma omp for schedule (static) private(i) reduction(+:s) 465 for (i = 0; i < rows; i++) {
477 if (a1 < 0 || a1 == 0 || a1 > 0) {
482 (
"Unable to solve the linear equation system"));
490 G_message(
_(
"Sparse CG -- iteration %i error %g\n"), m, a0);
492 G_message(
_(
"CG -- iteration %i error %g\n"), m, a0);
494 if (error_break == 1) {
536 int maxit,
double err)
538 return solver_bicgstab(A,
NULL, x, b, rows, maxit, err);
563 double *
b,
int rows,
int maxit,
double err)
565 return solver_bicgstab(
NULL, Asp, x, b, rows, maxit, err);
570 int rows,
int maxit,
double err)
584 double s1 = 0.0, s2 = 0.0, s3 = 0.0;
586 double alpha = 0, beta = 0, omega, rr0 = 0, error;
620 for (m = 0; m < maxit; m++) {
622 #pragma omp parallel default(shared) 630 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3) 631 for (i = 0; i < rows; i++) {
641 if (error < 0 || error == 0 || error > 0) {
646 (
"Unable to solve the linear equation system"));
662 #pragma omp for schedule (static) private(i) reduction(+:s1, s2) 663 for (i = 0; i < rows; i++) {
678 #pragma omp for schedule (static) private(i) reduction(+:s1) 679 for (i = 0; i < rows; i++) {
685 beta = alpha / omega * s1 / rr0;
695 G_message(
_(
"Sparse BiCGStab -- iteration %i error %g\n"), m,
698 G_message(
_(
"BiCGStab -- iteration %i error %g\n"), m, error);
700 if (error_break == 1) {
737 int i, j, cols = rows;
744 #pragma omp parallel for schedule (static) private(i, j, sum) shared(A, Msp, rows, cols, prec) 745 for (i = 0; i < rows; i++) {
751 for (j = 0; j < cols; j++)
752 sum += A[i][j] * A[i][j];
753 spvect->
values[0] = 1.0 / sqrt(sum);
757 for (j = 0; j < cols; j++)
758 sum += fabs(A[i][j]);
759 spvect->
values[0] = 1.0 / (sum);
763 spvect->
values[0] = 1.0 / A[i][i];
768 spvect->
index[0] = i;
775 #pragma omp parallel for schedule (static) private(i, j, sum) shared(Asp, Msp, rows, cols, prec) 776 for (i = 0; i < rows; i++) {
782 for (j = 0; j < Asp[i]->
cols; j++)
783 sum += Asp[i]->values[j] * Asp[i]->values[j];
784 spvect->
values[0] = 1.0 / sqrt(sum);
788 for (j = 0; j < Asp[i]->
cols; j++)
789 sum += fabs(Asp[i]->values[j]);
790 spvect->
values[0] = 1.0 / (sum);
794 for (j = 0; j < Asp[i]->
cols; j++)
795 if (i == Asp[i]->index[j])
800 spvect->
index[0] = i;
void void void void G_fatal_error(const char *,...) __attribute__((format(printf
G_math_spvector * G_math_alloc_spvector(int)
Allocate memory for a sparse vector.
The row vector of the sparse matrix.
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_free_spmatrix(G_math_spvector **, int)
Release the memory of the sparse matrix.
void G_free(void *)
Free allocated memory.
void G_math_d_copy(double *, double *, int)
Copy the vector x to y.
void G_math_d_ax_by(double *, double *, double *, double, double, int)
Scales vectors x and y with the scalars a and b and adds them.
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...
void G_math_Ax_sparse(G_math_spvector **, double *, double *, int)
Compute the matrix - vector product of sparse matrix **Asp and vector x.
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
void G_message(const char *,...) __attribute__((format(printf
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
#define G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION
G_math_spvector ** G_math_alloc_spmatrix(int)
Allocate memory for a sparse matrix.
#define G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION
void G_warning(const char *,...) __attribute__((format(printf
#define G_MATH_DIAGONAL_PRECONDITION
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_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...
void G_math_d_Ax(double **, double *, double *, int, int)
Compute the matrix - vector product of matrix A and vector x.
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_add_spvector(G_math_spvector **, G_math_spvector *, int)
Adds a sparse vector to a sparse matrix at position row.
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_cg(double **A, double *x, double *b, int rows, int maxit, double err)
The iterative conjugate gradients solver for symmetric positive definite matrices.
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.
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.