23 #include "grass/N_pde.h"
26 static int sparse_jacobi_gauss(
N_les * L,
int maxit,
double sor,
double error,
28 static int jacobi(
double **M,
double *
b,
double *x,
int rows,
int maxit,
29 double sor,
double error);
30 static int gauss_seidel(
double **M,
double *
b,
double *x,
int rows,
int maxit,
31 double sor,
double error);
57 G_warning(_(
"The linear equation system is not quadratic"));
64 return sparse_jacobi_gauss(L, maxit, sor, error,
94 G_warning(_(
"The linear equation system is not quadratic"));
99 return gauss_seidel(L->
A, L->
b, L->
x, L->
rows, maxit, sor, error);
101 return sparse_jacobi_gauss(L, maxit, sor, error,
109 sparse_jacobi_gauss(
N_les * L,
int maxit,
double sor,
double error,
112 int i, j, k, rows, finished = 0;
113 double *Enew, *x, *
b;
122 for (k = 0; k < maxit; k++) {
126 for (j = 0; j < rows; j++) {
130 for (i = 0; i < rows; i++) {
133 for (j = 0; j < L->
Asp[i]->
cols; j++) {
138 for (j = 0; j < L->
Asp[i]->
cols; j++) {
142 Enew[i] = x[i] - sor * (E - b[i]) / L->
Asp[i]->
values[0];
144 for (j = 0; j < rows; j++) {
145 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
152 G_message(_(
"sparse Jacobi -- iteration %5i error %g\n"), k,
err);
154 G_message(_(
"sparse SOR -- iteration %5i error %g\n"), k,
err);
170 int jacobi(
double **M,
double *b,
double *x,
int rows,
int maxit,
double sor,
179 for (j = 0; j < rows; j++) {
183 for (k = 0; k < maxit; k++) {
184 for (i = 0; i < rows; i++) {
186 for (j = 0; j < rows; j++) {
189 Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
192 for (j = 0; j < rows; j++) {
193 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
197 G_message(_(
"Jacobi -- iteration %5i error %g\n"), k, err);
208 int gauss_seidel(
double **M,
double *b,
double *x,
int rows,
int maxit,
209 double sor,
double error)
217 for (j = 0; j < rows; j++) {
221 for (k = 0; k < maxit; k++) {
222 for (i = 0; i < rows; i++) {
224 for (j = 0; j < rows; j++) {
225 E += M[i][j] * Enew[j];
227 Enew[i] = x[i] - sor * (E - b[i]) / M[i][i];
230 for (j = 0; j < rows; j++) {
231 err += (x[j] - Enew[j]) * (x[j] - Enew[j]);
235 G_message(_(
"SOR -- iteration %5i error %g\n"), k, err);
void G_free(void *buf)
Free allocated memory.
#define N_SOLVER_ITERATIVE_JACOBI
int jacobi(double a[MX][MX], long n, double d[MX], double v[MX][MX])
int G_warning(const char *msg,...)
Print a warning message to stderr.
SYMBOL * err(FILE *fp, SYMBOL *s, char *msg)
#define N_SOLVER_ITERATIVE_SOR
void G_message(const char *msg,...)
Print a message to stderr.
double * vectmem(int rows)
Allocate vector memory.
int N_solver_jacobi(N_les *L, int maxit, double sor, double error)
The iterative jacobian solver for regular matrices.
int N_solver_SOR(N_les *L, int maxit, double sor, double error)
The iterative overrelaxed gauss seidel solver for regular matrices.
The linear equation system (les) structure.