Loading...
Searching...
No Matches
solvers_krylov.c
Go to the documentation of this file.
/*****************************************************************************
*
* MODULE: Grass numerical math interface
* AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
* soerengebbert <at> googlemail <dot> com
*
* PURPOSE: linear equation system solvers
* part of the gmath library
*
* COPYRIGHT: (C) 2010 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
* for details.
*
*****************************************************************************/
#include <assert.h>
#include <math.h>
#include <unistd.h>
#include <stdio.h>
#include <string.h>
#include <grass/gis.h>
#include <grass/gmath.h>
#include <grass/glocale.h>
static G_math_spvector **create_diag_precond_matrix(double **A,
G_math_spvector **Asp,
int rows, int prec);
static int solver_pcg(double **A, G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err, int prec, int has_band,
int bandwidth);
static int solver_cg(double **A, G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err, int has_band,
int bandwidth);
static int solver_bicgstab(double **A, G_math_spvector **Asp, double *x,
double *b, int rows, int maxit, double err);
/*!
* \brief The iterative preconditioned conjugate gradients solver for symmetric
* positive definite matrices
*
* This iterative solver works with symmetric positive definite regular
* quadratic matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param A (double **) -- the matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \param prec (int) -- the preconditioner which should be used 1,2 or 3
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_pcg(double **A, double *x, double *b, int rows, int maxit,
double err, int prec)
{
return solver_pcg(A, NULL, x, b, rows, maxit, err, prec, 0, 0);
}
/*!
* \brief The iterative preconditioned conjugate gradients solver for symmetric
* positive definite band matrices
*
* WARNING: The preconditioning of symmetric band matrices is not implemented
* yet
*
* This iterative solver works with symmetric positive definite band matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param A (double **) -- the positive definite band matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param bandwidth (int) -- bandwidth of matrix A
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \param prec (int) -- the preconditioner which should be used 1,2 or 3
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_pcg_sband(double **A, double *x, double *b, int rows,
int bandwidth, int maxit, double err, int prec)
{
G_fatal_error("Preconditioning of band matrices is not implemented yet");
return solver_pcg(A, NULL, x, b, rows, maxit, err, prec, 1, bandwidth);
}
/*!
* \brief The iterative preconditioned conjugate gradients solver for sparse
* symmetric positive definite matrices
*
* This iterative solver works with symmetric positive definite sparse matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param Asp (G_math_spvector **) -- the sparse matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \param prec (int) -- the preconditioner which should be used 1,2 or 3
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_sparse_pcg(G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err, int prec)
{
return solver_pcg(NULL, Asp, x, b, rows, maxit, err, prec, 0, 0);
}
int solver_pcg(double **A, G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err, int prec, int has_band,
int bandwidth)
{
double *r, *z;
double *p;
double *v;
double s = 0.0;
double a0 = 0, a1 = 0, mygamma, tmp = 0;
int m, i;
int finished = 2;
int error_break;
G_math_spvector **M;
r = G_alloc_vector(rows);
p = G_alloc_vector(rows);
v = G_alloc_vector(rows);
z = G_alloc_vector(rows);
error_break = 0;
/*compute the preconditioning matrix, this is a sparse matrix */
M = create_diag_precond_matrix(A, Asp, rows, prec);
/*
* residual calculation
*/
#pragma omp parallel
{
if (Asp)
G_math_Ax_sparse(Asp, x, v, rows);
else if (has_band)
G_math_Ax_sband(A, x, v, rows, bandwidth);
else
G_math_d_Ax(A, x, v, rows, rows);
G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
/*perform the preconditioning */
G_math_Ax_sparse(M, r, p, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += p[i] * r[i];
}
}
a0 = s;
s = 0.0;
/* ******************* */
/* start the iteration */
/* ******************* */
for (m = 0; m < maxit; m++) {
#pragma omp parallel default(shared)
{
if (Asp)
G_math_Ax_sparse(Asp, p, v, rows);
else if (has_band)
G_math_Ax_sband(A, p, v, rows, bandwidth);
else
G_math_d_Ax(A, p, v, rows, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += v[i] * p[i];
}
/* barrier */
#pragma omp single
{
tmp = s;
mygamma = a0 / tmp;
s = 0.0;
}
G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
if (m % 50 == 1) {
if (Asp)
G_math_Ax_sparse(Asp, x, v, rows);
else if (has_band)
G_math_Ax_sband(A, x, v, rows, bandwidth);
else
G_math_d_Ax(A, x, v, rows, rows);
G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
}
else {
G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
}
/*perform the preconditioning */
G_math_Ax_sparse(M, r, z, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += z[i] * r[i];
}
/* barrier */
#pragma omp single
{
a1 = s;
tmp = a1 / a0;
a0 = a1;
s = 0.0;
if (a1 < 0 || a1 == 0 || a1 > 0) {
;
}
else {
G_warning(_("Unable to solve the linear equation system"));
error_break = 1;
}
}
G_math_d_ax_by(p, z, p, tmp, 1.0, rows);
}
if (Asp != NULL)
G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0);
else
G_message(_("PCG -- iteration %i error %g\n"), m, a0);
if (error_break == 1) {
finished = -1;
break;
}
if (a0 < err) {
finished = 1;
break;
}
}
G_free(r);
G_free(p);
G_free(v);
G_free(z);
G_math_free_spmatrix(M, rows);
return finished;
}
/*!
* \brief The iterative conjugate gradients solver for symmetric positive
* definite matrices
*
* This iterative solver works with symmetric positive definite regular
* quadratic matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param A (double **) -- the matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit,
double err)
{
return solver_cg(A, NULL, x, b, rows, maxit, err, 0, 0);
}
/*!
* \brief The iterative conjugate gradients solver for symmetric positive
* definite band matrices
*
* This iterative solver works with symmetric positive definite band matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param A (double **) -- the symmetric positive definite band matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param bandwidth (int) -- the bandwidth of matrix A
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_cg_sband(double **A, double *x, double *b, int rows,
int bandwidth, int maxit, double err)
{
return solver_cg(A, NULL, x, b, rows, maxit, err, 1, bandwidth);
}
/*!
* \brief The iterative conjugate gradients solver for sparse symmetric positive
* definite matrices
*
* This iterative solver works with symmetric positive definite sparse matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param Asp (G_math_spvector **) -- the sparse matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_sparse_cg(G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err)
{
return solver_cg(NULL, Asp, x, b, rows, maxit, err, 0, 0);
}
int solver_cg(double **A, G_math_spvector **Asp, double *x, double *b, int rows,
int maxit, double err, int has_band, int bandwidth)
{
double *r;
double *p;
double *v;
double s = 0.0;
double a0 = 0, a1 = 0, mygamma, tmp = 0;
int m, i;
int finished = 2;
int error_break;
r = G_alloc_vector(rows);
p = G_alloc_vector(rows);
v = G_alloc_vector(rows);
error_break = 0;
/*
* residual calculation
*/
#pragma omp parallel
{
if (Asp)
G_math_Ax_sparse(Asp, x, v, rows);
else if (has_band)
G_math_Ax_sband(A, x, v, rows, bandwidth);
else
G_math_d_Ax(A, x, v, rows, rows);
G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
G_math_d_copy(r, p, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += r[i] * r[i];
}
}
a0 = s;
s = 0.0;
/* ******************* */
/* start the iteration */
/* ******************* */
for (m = 0; m < maxit; m++) {
#pragma omp parallel default(shared)
{
if (Asp)
G_math_Ax_sparse(Asp, p, v, rows);
else if (has_band)
G_math_Ax_sband(A, p, v, rows, bandwidth);
else
G_math_d_Ax(A, p, v, rows, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += v[i] * p[i];
}
/* barrier */
#pragma omp single
{
tmp = s;
mygamma = a0 / tmp;
s = 0.0;
}
G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
if (m % 50 == 1) {
if (Asp)
G_math_Ax_sparse(Asp, x, v, rows);
else if (has_band)
G_math_Ax_sband(A, x, v, rows, bandwidth);
else
G_math_d_Ax(A, x, v, rows, rows);
G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
}
else {
G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
}
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s)
for (i = 0; i < rows; i++) {
s += r[i] * r[i];
}
/* barrier */
#pragma omp single
{
a1 = s;
tmp = a1 / a0;
a0 = a1;
s = 0.0;
if (a1 < 0 || a1 == 0 || a1 > 0) {
;
}
else {
G_warning(_("Unable to solve the linear equation system"));
error_break = 1;
}
}
G_math_d_ax_by(p, r, p, tmp, 1.0, rows);
}
if (Asp != NULL)
G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0);
else
G_message(_("CG -- iteration %i error %g\n"), m, a0);
if (error_break == 1) {
finished = -1;
break;
}
if (a0 < err) {
finished = 1;
break;
}
}
G_free(r);
G_free(p);
G_free(v);
return finished;
}
/*!
* \brief The iterative biconjugate gradients solver with stabilization for
* unsymmetric non-definite matrices
*
* This iterative solver works with regular quadratic matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param A (double **) -- the matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_bicgstab(double **A, double *x, double *b, int rows,
int maxit, double err)
{
return solver_bicgstab(A, NULL, x, b, rows, maxit, err);
}
/*!
* \brief The iterative biconjugate gradients solver with stabilization for
* unsymmetric non-definite matrices
*
* This iterative solver works with sparse matrices.
*
* This solver solves the linear equation system:
* A x = b
*
* The parameter <i>maxit</i> specifies the maximum number of iterations. If the
* maximum is reached, the solver will abort the calculation and writes the
* current result into the vector x. The parameter <i>err</i> defines the error
* break criteria for the solver.
*
* \param Asp (G_math_spvector **) -- the sparse matrix
* \param x (double *) -- the value vector
* \param b (double *) -- the right hand side
* \param rows (int)
* \param maxit (int) -- the maximum number of iterations
* \param err (double) -- defines the error break criteria
* \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix
* singular, -1 - could not solve the les
*
* */
int G_math_solver_sparse_bicgstab(G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err)
{
return solver_bicgstab(NULL, Asp, x, b, rows, maxit, err);
}
int solver_bicgstab(double **A, G_math_spvector **Asp, double *x, double *b,
int rows, int maxit, double err)
{
double *r;
double *r0;
double *p;
double *v;
double *s;
double *t;
double s1 = 0.0, s2 = 0.0, s3 = 0.0;
double alpha = 0, beta = 0, omega, rr0 = 0, error;
int m, i;
int finished = 2;
int error_break;
r = G_alloc_vector(rows);
r0 = G_alloc_vector(rows);
p = G_alloc_vector(rows);
v = G_alloc_vector(rows);
s = G_alloc_vector(rows);
t = G_alloc_vector(rows);
error_break = 0;
#pragma omp parallel
{
if (Asp)
G_math_Ax_sparse(Asp, x, v, rows);
else
G_math_d_Ax(A, x, v, rows, rows);
G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
G_math_d_copy(r, r0, rows);
G_math_d_copy(r, p, rows);
}
s1 = s2 = s3 = 0.0;
/* ******************* */
/* start the iteration */
/* ******************* */
for (m = 0; m < maxit; m++) {
#pragma omp parallel default(shared)
{
if (Asp)
G_math_Ax_sparse(Asp, p, v, rows);
else
G_math_d_Ax(A, p, v, rows, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s1, s2, s3)
for (i = 0; i < rows; i++) {
s1 += r[i] * r[i];
s2 += r[i] * r0[i];
s3 += v[i] * r0[i];
}
#pragma omp single
{
error = s1;
if (error < 0 || error == 0 || error > 0) {
;
}
else {
G_warning(_("Unable to solve the linear equation system"));
error_break = 1;
}
rr0 = s2;
alpha = rr0 / s3;
s1 = s2 = s3 = 0.0;
}
G_math_d_ax_by(r, v, s, 1.0, -1.0 * alpha, rows);
if (Asp)
G_math_Ax_sparse(Asp, s, t, rows);
else
G_math_d_Ax(A, s, t, rows, rows);
/* scalar product */
#pragma omp for schedule(static) private(i) reduction(+ : s1, s2)
for (i = 0; i < rows; i++) {
s1 += t[i] * s[i];
s2 += t[i] * t[i];
}
#pragma omp single
{
omega = s1 / s2;
s1 = s2 = 0.0;
}
G_math_d_ax_by(p, s, r, alpha, omega, rows);
G_math_d_ax_by(x, r, x, 1.0, 1.0, rows);
G_math_d_ax_by(s, t, r, 1.0, -1.0 * omega, rows);
#pragma omp for schedule(static) private(i) reduction(+ : s1)
for (i = 0; i < rows; i++) {
s1 += r[i] * r0[i];
}
#pragma omp single
{
beta = alpha / omega * s1 / rr0;
s1 = s2 = s3 = 0.0;
}
G_math_d_ax_by(p, v, p, 1.0, -1.0 * omega, rows);
G_math_d_ax_by(p, r, p, beta, 1.0, rows);
}
if (Asp != NULL)
G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m,
error);
else
G_message(_("BiCGStab -- iteration %i error %g\n"), m, error);
if (error_break == 1) {
finished = -1;
break;
}
if (error < err) {
finished = 1;
break;
}
}
G_free(r);
G_free(r0);
G_free(p);
G_free(v);
G_free(s);
G_free(t);
return finished;
}
/*!
* \brief Compute a diagonal preconditioning matrix for krylov space solver
*
* \param A (double **) -- the matrix for which the precondition should be
* computed (if the sparse matrix is used, set it to NULL) \param Asp
* (G_math_spvector **) -- the matrix for which the precondition should be
* computed \param rows (int) \param prec (int) -- which preconditioner should
* be used 1, 2 or 3
*
* */
G_math_spvector **create_diag_precond_matrix(double **A, G_math_spvector **Asp,
int rows, int prec)
{
G_math_spvector **Msp;
unsigned int i, j, cols = (unsigned int)rows;
double sum;
assert(rows >= 0);
Msp = G_math_alloc_spmatrix(rows);
if (A != NULL) {
#pragma omp parallel for schedule(static) private(i, j, sum) \
shared(A, Msp, rows, cols, prec)
for (i = 0; i < (unsigned int)rows; i++) {
G_math_spvector *spvect = G_math_alloc_spvector(1);
switch (prec) {
case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
sum = 0;
for (j = 0; j < cols; j++)
sum += A[i][j] * A[i][j];
spvect->values[0] = 1.0 / sqrt(sum);
break;
case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
sum = 0;
for (j = 0; j < cols; j++)
sum += fabs(A[i][j]);
spvect->values[0] = 1.0 / (sum);
break;
case G_MATH_DIAGONAL_PRECONDITION:
default:
spvect->values[0] = 1.0 / A[i][i];
break;
}
spvect->index[0] = i;
spvect->cols = 1;
;
G_math_add_spvector(Msp, spvect, i);
}
}
else {
#pragma omp parallel for schedule(static) private(i, j, sum) \
shared(Asp, Msp, rows, cols, prec)
for (i = 0; i < (unsigned int)rows; i++) {
G_math_spvector *spvect = G_math_alloc_spvector(1);
switch (prec) {
case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
sum = 0;
for (j = 0; j < Asp[i]->cols; j++)
sum += Asp[i]->values[j] * Asp[i]->values[j];
spvect->values[0] = 1.0 / sqrt(sum);
break;
case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
sum = 0;
for (j = 0; j < Asp[i]->cols; j++)
sum += fabs(Asp[i]->values[j]);
spvect->values[0] = 1.0 / (sum);
break;
case G_MATH_DIAGONAL_PRECONDITION:
default:
for (j = 0; j < Asp[i]->cols; j++)
if (i == Asp[i]->index[j])
spvect->values[0] = 1.0 / Asp[i]->values[j];
break;
}
spvect->index[0] = i;
spvect->cols = 1;
;
G_math_add_spvector(Msp, spvect, i);
}
}
return Msp;
}
G_free
void G_free(void *buf)
Free allocated memory.
Definition alloc.c:147
G_math_d_ax_by
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.
Definition blas_level_1.c:173
G_math_d_copy
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
Definition blas_level_1.c:237
G_math_d_Ax
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.
Definition blas_level_2.c:46
NULL
#define NULL
Definition ccmath.h:32
G_alloc_vector
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition dalloc.c:38
b
double b
Definition driver/set_window.c:5
t
double t
Definition driver/set_window.c:5
r
double r
Definition driver/set_window.c:5
G_fatal_error
void G_fatal_error(const char *msg,...)
Print a fatal error message to stderr.
Definition gis/error.c:159
G_message
void G_message(const char *msg,...)
Print a message to stderr.
Definition gis/error.c:89
G_warning
void G_warning(const char *msg,...)
Print a warning message to stderr.
Definition gis/error.c:203
assert
#define assert(condition)
Definition lz4.c:291
G_math_solver_pcg
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.
Definition solvers_krylov.c:65
G_math_solver_cg
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.
Definition solvers_krylov.c:319
G_math_solver_bicgstab
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.
Definition solvers_krylov.c:550
G_math_solver_sparse_cg
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.
Definition solvers_krylov.c:380
G_math_solver_cg_sband
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.
Definition solvers_krylov.c:350
G_math_solver_pcg_sband
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...
Definition solvers_krylov.c:101
G_math_solver_sparse_bicgstab
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.
Definition solvers_krylov.c:580
G_math_solver_sparse_pcg
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...
Definition solvers_krylov.c:133
G_math_Ax_sparse
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.
Definition sparse_matrix.c:369
G_math_alloc_spvector
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
Definition sparse_matrix.c:77
G_math_free_spmatrix
void G_math_free_spmatrix(G_math_spvector **Asp, int rows)
Release the memory of the sparse matrix.
Definition sparse_matrix.c:122
G_math_add_spvector
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.
Definition sparse_matrix.c:35
G_math_alloc_spmatrix
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition sparse_matrix.c:59
err
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
Definition symbol/read.c:216
G_math_Ax_sband
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.
Definition symmetric_band_matrix.c:123
x
#define x