GRASS Programmer's Manual  6.5.svn(2014)-r66266
solvers_krylov.c
Go to the documentation of this file.
1
2 /*****************************************************************************
3 *
4 * MODULE: Grass PDE Numerical Library
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> gmx <dot> de
7 *
8 * PURPOSE: linear equation system solvers
9 * part of the gpde library
10 *
11 * COPYRIGHT: (C) 2000 by the GRASS Development Team
12 *
13 * This program is free software under the GNU General Public
14 * License (>=v2). Read the file COPYING that comes with GRASS
15 * for details.
16 *
17 *****************************************************************************/
18
19 #include <math.h>
20 #include <unistd.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include <grass/gis.h>
24 #include <grass/gmath.h>
25 #include <grass/glocale.h>
26
27 #define G_MATH_ROWSCALE_L2_NORM_PRECONDITION 1
28 #define G_MATH_ROWSCALE_L1_NORM_PRECONDITION 2
29 #define G_MATH_DIAGONAL_PRECONDITION 3
30
31 static G_math_spvector **create_diag_precond_matrix(double **A,
32  G_math_spvector ** Asp,
33  int rows, int prec);
34 static int solver_pcg(double **A, G_math_spvector ** Asp, double *x,
35  double *b, int rows, int maxit, double err, int prec);
36 static int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
37  int rows, int maxit, double err);
38 static int solver_bicgstab(double **A, G_math_spvector ** Asp, double *x,
39  double *b, int rows, int maxit, double err);
40
41
64 int G_math_solver_pcg(double **A, double *x, double *b, int rows, int maxit,
65  double err, int prec)
66 {
67
68  return solver_pcg(A, NULL, x, b, rows, maxit, err, prec);
69 }
70
93 int G_math_solver_sparse_pcg(G_math_spvector ** Asp, double *x, double *b,
94  int rows, int maxit, double err, int prec)
95 {
96
97  return solver_pcg(NULL, Asp, x, b, rows, maxit, err, prec);
98 }
99
100 int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
101  int rows, int maxit, double err, int prec)
102 {
103  double *r, *z;
104
105  double *p;
106
107  double *v;
108
109  double s = 0.0;
110
111  double a0 = 0, a1 = 0, mygamma, tmp = 0;
112
113  int m, i;
114
115  int finished = 2;
116
117  int error_break;
118
119  G_math_spvector **M;
120
121  r = G_alloc_vector(rows);
122  p = G_alloc_vector(rows);
123  v = G_alloc_vector(rows);
124  z = G_alloc_vector(rows);
125
126  error_break = 0;
127
128  /*compute the preconditioning matrix, this is a sparse matrix */
129  M = create_diag_precond_matrix(A, Asp, rows, prec);
130
131  /*
132  * residual calculation
133  */
134 #pragma omp parallel
135  {
136  if (Asp)
137  G_math_Ax_sparse(Asp, x, v, rows);
138  else
139  G_math_d_Ax(A, x, v, rows, rows);
140
141  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
142  /*performe the preconditioning */
143  G_math_Ax_sparse(M, r, p, rows);
144
145  /* scalar product */
146 #pragma omp for schedule (static) private(i) reduction(+:s)
147  for (i = 0; i < rows; i++) {
148  s += p[i] * r[i];
149  }
150  }
151
152  a0 = s;
153  s = 0.0;
154
155  /* ******************* */
156  /* start the iteration */
157  /* ******************* */
158  for (m = 0; m < maxit; m++) {
159 #pragma omp parallel default(shared)
160  {
161  if (Asp)
162  G_math_Ax_sparse(Asp, p, v, rows);
163  else
164  G_math_d_Ax(A, p, v, rows, rows);
165
166
167
168  /* scalar product */
169 #pragma omp for schedule (static) private(i) reduction(+:s)
170  for (i = 0; i < rows; i++) {
171  s += v[i] * p[i];
172  }
173
174  /* barrier */
175 #pragma omp single
176  {
177  tmp = s;
178  mygamma = a0 / tmp;
179  s = 0.0;
180  }
181
182  G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
183
184  if (m % 50 == 1) {
185  if (Asp)
186  G_math_Ax_sparse(Asp, x, v, rows);
187  else
188  G_math_d_Ax(A, x, v, rows, rows);
189
190  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
191  }
192  else {
193  G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
194  }
195
196  /*performe the preconditioning */
197  G_math_Ax_sparse(M, r, z, rows);
198
199
200  /* scalar product */
201 #pragma omp for schedule (static) private(i) reduction(+:s)
202  for (i = 0; i < rows; i++) {
203  s += z[i] * r[i];
204  }
205
206  /* barrier */
207 #pragma omp single
208  {
209  a1 = s;
210  tmp = a1 / a0;
211  a0 = a1;
212  s = 0.0;
213
214  if (a1 < 0 || a1 == 0 || a1 > 0) {
215  ;
216  }
217  else {
218  G_warning(_
219  ("Unable to solve the linear equation system"));
220  error_break = 1;
221  }
222  }
223  G_math_d_ax_by(p, z, p, tmp, 1.0, rows);
224  }
225
226  if (Asp != NULL)
227  G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0);
228  else
229  G_message(_("PCG -- iteration %i error %g\n"), m, a0);
230
231  if (error_break == 1) {
232  finished = -1;
233  break;
234  }
235
236
237  if (a0 < err) {
238  finished = 1;
239  break;
240  }
241  }
242
243  G_free(r);
244  G_free(p);
245  G_free(v);
246  G_free(z);
247  G_math_free_spmatrix(M, rows);
248
249  return finished;
250 }
251
252
274 int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit,
275  double err)
276 {
277  return solver_cg(A, NULL, x, b, rows, maxit, err);
278 }
279
301 int G_math_solver_sparse_cg(G_math_spvector ** Asp, double *x, double *b,
302  int rows, int maxit, double err)
303 {
304  return solver_cg(NULL, Asp, x, b, rows, maxit, err);
305 }
306
307
308 int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
309  int rows, int maxit, double err)
310 {
311  double *r;
312
313  double *p;
314
315  double *v;
316
317  double s = 0.0;
318
319  double a0 = 0, a1 = 0, mygamma, tmp = 0;
320
321  int m, i;
322
323  int finished = 2;
324
325  int error_break;
326
327  r = G_alloc_vector(rows);
328  p = G_alloc_vector(rows);
329  v = G_alloc_vector(rows);
330
331  error_break = 0;
332  /*
333  * residual calculation
334  */
335 #pragma omp parallel
336  {
337  if (Asp)
338  G_math_Ax_sparse(Asp, x, v, rows);
339  else
340  G_math_d_Ax(A, x, v, rows, rows);
341
342  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
343  G_math_d_copy(r, p, rows);
344
345  /* scalar product */
346 #pragma omp for schedule (static) private(i) reduction(+:s)
347  for (i = 0; i < rows; i++) {
348  s += r[i] * r[i];
349  }
350  }
351
352  a0 = s;
353  s = 0.0;
354
355  /* ******************* */
356  /* start the iteration */
357  /* ******************* */
358  for (m = 0; m < maxit; m++) {
359 #pragma omp parallel default(shared)
360  {
361  if (Asp)
362  G_math_Ax_sparse(Asp, p, v, rows);
363  else
364  G_math_d_Ax(A, p, v, rows, rows);
365
366  /* scalar product */
367 #pragma omp for schedule (static) private(i) reduction(+:s)
368  for (i = 0; i < rows; i++) {
369  s += v[i] * p[i];
370  }
371
372  /* barrier */
373 #pragma omp single
374  {
375  tmp = s;
376  mygamma = a0 / tmp;
377  s = 0.0;
378  }
379
380  G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
381
382  if (m % 50 == 1) {
383  if (Asp)
384  G_math_Ax_sparse(Asp, x, v, rows);
385  else
386  G_math_d_Ax(A, x, v, rows, rows);
387
388  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
389  }
390  else {
391  G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
392  }
393
394  /* scalar product */
395 #pragma omp for schedule (static) private(i) reduction(+:s)
396  for (i = 0; i < rows; i++) {
397  s += r[i] * r[i];
398  }
399
400  /* barrier */
401 #pragma omp single
402  {
403  a1 = s;
404  tmp = a1 / a0;
405  a0 = a1;
406  s = 0.0;
407
408  if (a1 < 0 || a1 == 0 || a1 > 0) {
409  ;
410  }
411  else {
412  G_warning(_
413  ("Unable to solve the linear equation system"));
414  error_break = 1;
415  }
416  }
417  G_math_d_ax_by(p, r, p, tmp, 1.0, rows);
418  }
419
420  if (Asp != NULL)
421  G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0);
422  else
423  G_message(_("CG -- iteration %i error %g\n"), m, a0);
424
425  if (error_break == 1) {
426  finished = -1;
427  break;
428  }
429
430  if (a0 < err) {
431  finished = 1;
432  break;
433  }
434  }
435
436  G_free(r);
437  G_free(p);
438  G_free(v);
439
440  return finished;
441 }
442
443
444
466 int G_math_solver_bicgstab(double **A, double *x, double *b, int rows,
467  int maxit, double err)
468 {
469  return solver_bicgstab(A, NULL, x, b, rows, maxit, err);
470 }
471
493 int G_math_solver_sparse_bicgstab(G_math_spvector ** Asp, double *x,
494  double *b, int rows, int maxit, double err)
495 {
496  return solver_bicgstab(NULL, Asp, x, b, rows, maxit, err);
497 }
498
499
500 int solver_bicgstab(double **A, G_math_spvector ** Asp, double *x, double *b,
501  int rows, int maxit, double err)
502 {
503  double *r;
504
505  double *r0;
506
507  double *p;
508
509  double *v;
510
511  double *s;
512
513  double *t;
514
515  double s1 = 0.0, s2 = 0.0, s3 = 0.0;
516
517  double alpha = 0, beta = 0, omega, rr0 = 0, error;
518
519  int m, i;
520
521  int finished = 2;
522
523  int error_break;
524
525  r = G_alloc_vector(rows);
526  r0 = G_alloc_vector(rows);
527  p = G_alloc_vector(rows);
528  v = G_alloc_vector(rows);
529  s = G_alloc_vector(rows);
530  t = G_alloc_vector(rows);
531
532  error_break = 0;
533
534 #pragma omp parallel
535  {
536  if (Asp)
537  G_math_Ax_sparse(Asp, x, v, rows);
538  else
539  G_math_d_Ax(A, x, v, rows, rows);
540
541  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
542  G_math_d_copy(r, r0, rows);
543  G_math_d_copy(r, p, rows);
544  }
545
546  s1 = s2 = s3 = 0.0;
547
548  /* ******************* */
549  /* start the iteration */
550  /* ******************* */
551  for (m = 0; m < maxit; m++) {
552
553 #pragma omp parallel default(shared)
554  {
555  if (Asp)
556  G_math_Ax_sparse(Asp, p, v, rows);
557  else
558  G_math_d_Ax(A, p, v, rows, rows);
559
560  /* scalar product */
561 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3)
562  for (i = 0; i < rows; i++) {
563  s1 += r[i] * r[i];
564  s2 += r[i] * r0[i];
565  s3 += v[i] * r0[i];
566  }
567
568 #pragma omp single
569  {
570  error = s1;
571
572  if (error < 0 || error == 0 || error > 0) {
573  ;
574  }
575  else {
576  G_warning(_
577  ("Unable to solve the linear equation system"));
578  error_break = 1;
579  }
580
581  rr0 = s2;
582  alpha = rr0 / s3;
583  s1 = s2 = s3 = 0.0;
584  }
585
586  G_math_d_ax_by(r, v, s, 1.0, -1.0 * alpha, rows);
587  if (Asp)
588  G_math_Ax_sparse(Asp, s, t, rows);
589  else
590  G_math_d_Ax(A, s, t, rows, rows);
591
592  /* scalar product */
593 #pragma omp for schedule (static) private(i) reduction(+:s1, s2)
594  for (i = 0; i < rows; i++) {
595  s1 += t[i] * s[i];
596  s2 += t[i] * t[i];
597  }
598
599 #pragma omp single
600  {
601  omega = s1 / s2;
602  s1 = s2 = 0.0;
603  }
604
605  G_math_d_ax_by(p, s, r, alpha, omega, rows);
606  G_math_d_ax_by(x, r, x, 1.0, 1.0, rows);
607  G_math_d_ax_by(s, t, r, 1.0, -1.0 * omega, rows);
608
609 #pragma omp for schedule (static) private(i) reduction(+:s1)
610  for (i = 0; i < rows; i++) {
611  s1 += r[i] * r0[i];
612  }
613
614 #pragma omp single
615  {
616  beta = alpha / omega * s1 / rr0;
617  s1 = s2 = s3 = 0.0;
618  }
619
620  G_math_d_ax_by(p, v, p, 1.0, -1.0 * omega, rows);
621  G_math_d_ax_by(p, r, p, beta, 1.0, rows);
622  }
623
624
625  if (Asp != NULL)
626  G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m,
627  error);
628  else
629  G_message(_("BiCGStab -- iteration %i error %g\n"), m, error);
630
631  if (error_break == 1) {
632  finished = -1;
633  break;
634  }
635
636  if (error < err) {
637  finished = 1;
638  break;
639  }
640  }
641
642  G_free(r);
643  G_free(r0);
644  G_free(p);
645  G_free(v);
646  G_free(s);
647  G_free(t);
648
649  return finished;
650 }
651
652
662 G_math_spvector **create_diag_precond_matrix(double **A,
663  G_math_spvector ** Asp, int rows,
664  int prec)
665 {
666  G_math_spvector **Msp;
667
668  int i, j, cols = rows;
669
670  double sum;
671
672  Msp = G_math_alloc_spmatrix(rows);
673
674  if (A != NULL) {
675 #pragma omp parallel for schedule (static) private(i, j, sum) shared(A, Msp, rows, cols, prec)
676  for (i = 0; i < rows; i++) {
677  G_math_spvector *spvect = G_math_alloc_spvector(1);
678
679  switch (prec) {
681  sum = 0;
682  for (j = 0; j < cols; j++)
683  sum += A[i][j] * A[i][j];
684  spvect->values[0] = 1.0 / sqrt(sum);
685  break;
687  sum = 0;
688  for (j = 0; j < cols; j++)
689  sum += fabs(A[i][j]);
690  spvect->values[0] = 1.0 / (sum);
691  break;
693  default:
694  spvect->values[0] = 1.0 / A[i][i];
695  break;
696  }
697
698
699  spvect->index[0] = i;
700  spvect->cols = 1;;
701  G_math_add_spvector(Msp, spvect, i);
702
703  }
704  }
705  else {
706 #pragma omp parallel for schedule (static) private(i, j, sum) shared(Asp, Msp, rows, cols, prec)
707  for (i = 0; i < rows; i++) {
708  G_math_spvector *spvect = G_math_alloc_spvector(1);
709
710  switch (prec) {
712  sum = 0;
713  for (j = 0; j < Asp[i]->cols; j++)
714  sum += Asp[i]->values[j] * Asp[i]->values[j];
715  spvect->values[0] = 1.0 / sqrt(sum);
716  break;
718  sum = 0;
719  for (j = 0; j < Asp[i]->cols; j++)
720  sum += fabs(Asp[i]->values[j]);
721  spvect->values[0] = 1.0 / (sum);
722  break;
724  default:
725  for (j = 0; j < Asp[i]->cols; j++)
726  if (i == Asp[i]->index[j])
727  spvect->values[0] = 1.0 / Asp[i]->values[j];
728  break;
729  }
730
731  spvect->index[0] = i;
732  spvect->cols = 1;;
733  G_math_add_spvector(Msp, spvect, i);
734  }
735  }
736  return Msp;
737 }
void G_free(void *buf)
Free allocated memory.
Definition: gis/alloc.c:142
float b
Definition: named_colr.c:8
#define G_MATH_DIAGONAL_PRECONDITION
G_math_spvector * G_math_alloc_spvector(int cols)
Allocate memory for a sparse vector.
Definition: sparse_matrix.c:76
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:79
float r
Definition: named_colr.c:8
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...
const char * err
Definition: g3dcolor.c:50
def error
Display an error message using g.message -e
Definition: core.py:370
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...
G_math_spvector ** G_math_alloc_spmatrix(int rows)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
void G_message(const char *msg,...)
Print a message to stderr.
Definition: lib/gis/error.c:74
#define G_MATH_ROWSCALE_L2_NORM_PRECONDITION
void G_math_d_copy(double *x, double *y, int rows)
Copy the vector x to y.
Definition: blas_level_1.c:236
void G_math_free_spmatrix(G_math_spvector **spmatrix, int rows)
Release the memory of the sparse matrix.
return NULL
Definition: dbfopen.c:1394
G_warning("category support for [%s] in mapset [%s] %s", name, mapset, type)
tuple cols
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: blas_level_2.c:45
#define G_MATH_ROWSCALE_L1_NORM_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...
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:172
double * G_alloc_vector(size_t n)
Vector matrix memory allocation.
Definition: dalloc.c:41
int G_math_add_spvector(G_math_spvector **Asp, G_math_spvector *spvector, int row)
Adds a sparse vector to a sparse linear equation system at position row.
Definition: sparse_matrix.c:35
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.
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.