GRASS Programmer's Manual  6.4.2(2012)
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
N_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/N_pde.h"
24 #include "solvers_local_proto.h"
25 
26 /*local protos */
27 static void scalar_product(double *a, double *b, double *scalar, int rows);
28 static void sub_vectors(double *source_a, double *source_b, double *result,
29  int row);
30 static void sub_vectors_scalar(double *source_a, double *source_b,
31  double *result, double scalar_b, int rows);
32 static void add_vectors(double *source_a, double *source_b, double *result,
33  int rows);
34 static void add_vectors_scalar(double *source_a, double *source_b,
35  double *result, double scalar_b, int rows);
36 static void add_vectors_scalar2(double *source_a, double *source_b,
37  double *result, double scalar_a,
38  double scalar_b, int rows);
39 static void scalar_vector_product(double *a, double *result, double scalar,
40  int rows);
41 static void sync_vectors(double *source, double *target, int rows);
42 
43 
44 /* ******************************************************* *
45  * *** preconditioned conjugate gradients **************** *
46  * ******************************************************* */
64 int N_solver_pcg(N_les * L, int maxit, double err, int prec)
65 {
66  double *r, *z;
67  double *p;
68  double *v;
69  double *x, *b;
70  double s = 0.0;
71  double a0 = 0, a1 = 0, mygamma, tmp = 0;
72  int m, rows, i;
73  int finished = 2;
74  int error_break;
75  N_les *M;
76 
77  if (L->quad != 1) {
78  G_warning(_("The linear equation system is not quadratic"));
79  return -1;
80  }
81 
82  /* check for symmetry */
83  if (check_symmetry(L) != 1) {
84  G_warning(_("Matrix is not symmetric!"));
85  }
86 
87  x = L->x;
88  b = L->b;
89  rows = L->rows;
90 
91  r = vectmem(rows);
92  p = vectmem(rows);
93  v = vectmem(rows);
94  z = vectmem(rows);
95 
96  error_break = 0;
97 
98  /*compute the preconditioning matrix */
99  M = N_create_diag_precond_matrix(L, prec);
100 
101  /*
102  * residual calculation
103  */
104 #pragma omp parallel
105  {
106  /* matrix vector multiplication */
107  if (L->type == N_SPARSE_LES)
109  else
110  N_matrix_vector_product(L, x, v);
111 
112  sub_vectors(b, v, r, rows);
113  /*performe the preconditioning */
115 
116  /* scalar product */
117 #pragma omp for schedule (static) private(i) reduction(+:s)
118  for (i = 0; i < rows; i++) {
119  s += p[i] * r[i];
120  }
121  }
122 
123  a0 = s;
124  s = 0.0;
125 
126  /* ******************* */
127  /* start the iteration */
128  /* ******************* */
129  for (m = 0; m < maxit; m++) {
130 #pragma omp parallel default(shared)
131  {
132  /* matrix vector multiplication */
133  if (L->type == N_SPARSE_LES)
135  else
136  N_matrix_vector_product(L, p, v);
137 
138 
139  /* scalar product */
140 #pragma omp for schedule (static) private(i) reduction(+:s)
141  for (i = 0; i < rows; i++) {
142  s += v[i] * p[i];
143  }
144 
145  /* barrier */
146 #pragma omp single
147  {
148  tmp = s;
149  mygamma = a0 / tmp;
150  s = 0.0;
151  }
152 
153  add_vectors_scalar(x, p, x, mygamma, rows);
154 
155  if (m % 50 == 1) {
156  /* matrix vector multiplication */
157  if (L->type == N_SPARSE_LES)
159  else
160  N_matrix_vector_product(L, x, v);
161 
162  sub_vectors(b, v, r, rows);
163 
164  }
165  else {
166  sub_vectors_scalar(r, v, r, mygamma, rows);
167  }
168 
169  /*performe the preconditioning */
171 
172  /* scalar product */
173 #pragma omp for schedule (static) private(i) reduction(+:s)
174  for (i = 0; i < rows; i++) {
175  s += z[i] * r[i];
176  }
177 
178  /* barrier */
179 #pragma omp single
180  {
181  a1 = s;
182  tmp = a1 / a0;
183  a0 = a1;
184  s = 0.0;
185 
186  if (a1 < 0 || a1 == 0 || a1 > 0) {
187  ;
188  }
189  else {
190  G_warning(_("Unable to solve the linear equation system"));
191  error_break = 1;
192  }
193  }
194  add_vectors_scalar(z, p, p, tmp, rows);
195  }
196 
197  if (L->type == N_SPARSE_LES)
198  G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0);
199  else
200  G_message(_("PCG -- iteration %i error %g\n"), m, a0);
201 
202  if (error_break == 1) {
203  finished = -1;
204  break;
205  }
206 
207 
208  if (a0 < err) {
209  finished = 1;
210  break;
211  }
212  }
213 
214  G_free(r);
215  G_free(p);
216  G_free(v);
217  G_free(z);
218 
219  return finished;
220 }
221 
222 
223 /* ******************************************************* *
224  * ****************** conjugate gradients **************** *
225  * ******************************************************* */
242 int N_solver_cg(N_les * L, int maxit, double err)
243 {
244  double *r;
245  double *p;
246  double *v;
247  double *x, *b;
248  double s = 0.0;
249  double a0 = 0, a1 = 0, mygamma, tmp = 0;
250  int m, rows, i;
251  int finished = 2;
252  int error_break;
253 
254  if (L->quad != 1) {
255  G_warning(_("The linear equation system is not quadratic"));
256  return -1;
257  }
258 
259  /* check for symmetry */
260  if (check_symmetry(L) != 1) {
261  G_warning(_("Matrix is not symmetric!"));
262  }
263 
264  x = L->x;
265  b = L->b;
266  rows = L->rows;
267 
268  r = vectmem(rows);
269  p = vectmem(rows);
270  v = vectmem(rows);
271 
272  error_break = 0;
273  /*
274  * residual calculation
275  */
276 #pragma omp parallel
277  {
278  /* matrix vector multiplication */
279  if (L->type == N_SPARSE_LES)
281  else
282  N_matrix_vector_product(L, x, v);
283 
284  sub_vectors(b, v, r, rows);
285  sync_vectors(r, p, rows);
286 
287  /* scalar product */
288 #pragma omp for schedule (static) private(i) reduction(+:s)
289  for (i = 0; i < rows; i++) {
290  s += r[i] * r[i];
291  }
292  }
293 
294  a0 = s;
295  s = 0.0;
296 
297  /* ******************* */
298  /* start the iteration */
299  /* ******************* */
300  for (m = 0; m < maxit; m++) {
301 #pragma omp parallel default(shared)
302  {
303  /* matrix vector multiplication */
304  if (L->type == N_SPARSE_LES)
306  else
307  N_matrix_vector_product(L, p, v);
308 
309 
310  /* scalar product */
311 #pragma omp for schedule (static) private(i) reduction(+:s)
312  for (i = 0; i < rows; i++) {
313  s += v[i] * p[i];
314  }
315 
316  /* barrier */
317 #pragma omp single
318  {
319  tmp = s;
320  mygamma = a0 / tmp;
321  s = 0.0;
322  }
323 
324  add_vectors_scalar(x, p, x, mygamma, rows);
325 
326  if (m % 50 == 1) {
327  /* matrix vector multiplication */
328  if (L->type == N_SPARSE_LES)
330  else
331  N_matrix_vector_product(L, x, v);
332 
333  sub_vectors(b, v, r, rows);
334  }
335  else {
336  sub_vectors_scalar(r, v, r, mygamma, rows);
337  }
338 
339  /* scalar product */
340 #pragma omp for schedule (static) private(i) reduction(+:s)
341  for (i = 0; i < rows; i++) {
342  s += r[i] * r[i];
343  }
344 
345  /* barrier */
346 #pragma omp single
347  {
348  a1 = s;
349  tmp = a1 / a0;
350  a0 = a1;
351  s = 0.0;
352 
353  if (a1 < 0 || a1 == 0 || a1 > 0) {
354  ;
355  }
356  else {
357  G_warning(_("Unable to solve the linear equation system"));
358  error_break = 1;
359  }
360  }
361  add_vectors_scalar(r, p, p, tmp, rows);
362  }
363 
364  if (L->type == N_SPARSE_LES)
365  G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0);
366  else
367  G_message(_("CG -- iteration %i error %g\n"), m, a0);
368 
369  if (error_break == 1) {
370  finished = -1;
371  break;
372  }
373 
374  if (a0 < err) {
375  finished = 1;
376  break;
377  }
378  }
379 
380  G_free(r);
381  G_free(p);
382  G_free(v);
383 
384  return finished;
385 }
386 
387 /* ******************************************************* *
388  * ************ biconjugate gradients ******************** *
389  * ******************************************************* */
407 int N_solver_bicgstab(N_les * L, int maxit, double err)
408 {
409  double *r;
410  double *r0;
411  double *p;
412  double *v;
413  double *s;
414  double *t;
415  double *x, *b;
416  double s1 = 0.0, s2 = 0.0, s3 = 0.0;
417  double alpha = 0, beta = 0, omega, rr0 = 0, error;
418  int m, rows, i;
419  int finished = 2;
420  int error_break;
421 
422  if (L->quad != 1) {
423  G_warning(_("The linear equation system is not quadratic"));
424  return -1;
425  }
426 
427 
428  x = L->x;
429  b = L->b;
430  rows = L->rows;
431  r = vectmem(rows);
432  r0 = vectmem(rows);
433  p = vectmem(rows);
434  v = vectmem(rows);
435  s = vectmem(rows);
436  t = vectmem(rows);
437 
438  error_break = 0;
439 
440 #pragma omp parallel
441  {
442  if (L->type == N_SPARSE_LES)
444  else
445  N_matrix_vector_product(L, x, v);
446  sub_vectors(b, v, r, rows);
447  sync_vectors(r, r0, rows);
448  sync_vectors(r, p, rows);
449  }
450 
451  s1 = s2 = s3 = 0.0;
452 
453  /* ******************* */
454  /* start the iteration */
455  /* ******************* */
456  for (m = 0; m < maxit; m++) {
457 
458 #pragma omp parallel default(shared)
459  {
460  if (L->type == N_SPARSE_LES)
462  else
463  N_matrix_vector_product(L, p, v);
464 
465  /* scalar product */
466 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3)
467  for (i = 0; i < rows; i++) {
468  s1 += r[i] * r[i];
469  s2 += r[i] * r0[i];
470  s3 += v[i] * r0[i];
471  }
472 
473 #pragma omp single
474  {
475  error = s1;
476 
477  if (error < 0 || error == 0 || error > 0) {
478  ;
479  }
480  else {
481  G_warning(_("Unable to solve the linear equation system"));
482  error_break = 1;
483  }
484 
485  rr0 = s2;
486  alpha = rr0 / s3;
487  s1 = s2 = s3 = 0.0;
488  }
489 
490  sub_vectors_scalar(r, v, s, alpha, rows);
491  if (L->type == N_SPARSE_LES)
493  else
494  N_matrix_vector_product(L, s, t);
495 
496  /* scalar product */
497 #pragma omp for schedule (static) private(i) reduction(+:s1, s2)
498  for (i = 0; i < rows; i++) {
499  s1 += t[i] * s[i];
500  s2 += t[i] * t[i];
501  }
502 
503 #pragma omp single
504  {
505  omega = s1 / s2;
506  s1 = s2 = 0.0;
507  }
508 
509  add_vectors_scalar2(p, s, r, alpha, omega, rows);
510  add_vectors(x, r, x, rows);
511  sub_vectors_scalar(s, t, r, omega, rows);
512 
513 #pragma omp for schedule (static) private(i) reduction(+:s1)
514  for (i = 0; i < rows; i++) {
515  s1 += r[i] * r0[i];
516  }
517 
518 #pragma omp single
519  {
520  beta = alpha / omega * s1 / rr0;
521  s1 = s2 = s3 = 0.0;
522  }
523 
524  sub_vectors_scalar(p, v, p, omega, rows);
525  add_vectors_scalar(r, p, p, beta, rows);
526  }
527 
528 
529  if (L->type == N_SPARSE_LES)
530  G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m,
531  error);
532  else
533  G_message(_("BiCGStab -- iteration %i error %g\n"), m, error);
534 
535  if (error_break == 1) {
536  finished = -1;
537  break;
538  }
539 
540  if (error < err) {
541  finished = 1;
542  break;
543  }
544  }
545 
546  G_free(r);
547  G_free(r0);
548  G_free(p);
549  G_free(v);
550  G_free(s);
551  G_free(t);
552 
553  return finished;
554 }
555 
568 void scalar_product(double *a, double *b, double *scalar, int rows)
569 {
570  int i;
571  double s = 0.0;
572 
573 #pragma omp parallel for schedule (static) reduction(+:s)
574  for (i = 0; i < rows; i++) {
575  s += a[i] * b[i];
576  }
577 
578  *scalar = s;
579  return;
580 }
581 
594 void N_matrix_vector_product(N_les * L, double *x, double *result)
595 {
596  int i, j;
597  double tmp;
598 
599 #pragma omp for schedule (static) private(i, j, tmp)
600  for (i = 0; i < L->rows; i++) {
601  tmp = 0;
602  for (j = 0; j < L->cols; j++) {
603  tmp += L->A[i][j] * x[j];
604  }
605  result[i] = tmp;
606  }
607  return;
608 }
609 
622 void N_sparse_matrix_vector_product(N_les * L, double *x, double *result)
623 {
624  int i, j;
625  double tmp;
626 
627 #pragma omp for schedule (static) private(i, j, tmp)
628  for (i = 0; i < L->rows; i++) {
629  tmp = 0;
630  for (j = 0; j < L->Asp[i]->cols; j++) {
631  tmp += L->Asp[i]->values[j] * x[L->Asp[i]->index[j]];
632  }
633  result[i] = tmp;
634  }
635  return;
636 }
637 
652 void
653 add_vectors_scalar2(double *a, double *b, double *result, double scalar_a,
654  double scalar_b, int rows)
655 {
656  int i;
657 
658 #pragma omp for schedule (static)
659  for (i = 0; i < rows; i++) {
660  result[i] = scalar_a * a[i] + scalar_b * b[i];
661  }
662 
663  return;
664 }
665 
679 void
680 add_vectors_scalar(double *a, double *b, double *result, double scalar_b,
681  int rows)
682 {
683  int i;
684 
685 #pragma omp for schedule (static)
686  for (i = 0; i < rows; i++) {
687  result[i] = a[i] + scalar_b * b[i];
688  }
689 
690  return;
691 }
692 
706 void
707 sub_vectors_scalar(double *a, double *b, double *result, double scalar_b,
708  int rows)
709 {
710  int i;
711 
712 #pragma omp for schedule (static)
713  for (i = 0; i < rows; i++) {
714  result[i] = a[i] - scalar_b * b[i];
715  }
716 
717  return;
718 }
719 
732 void add_vectors(double *a, double *b, double *result, int rows)
733 {
734  int i;
735 
736 #pragma omp for schedule (static)
737  for (i = 0; i < rows; i++) {
738  result[i] = a[i] + b[i];
739  }
740 
741  return;
742 }
743 
756 void sub_vectors(double *a, double *b, double *result, int rows)
757 {
758  int i;
759 
760 #pragma omp for schedule (static)
761  for (i = 0; i < rows; i++) {
762  result[i] = a[i] - b[i];
763  }
764 
765  return;
766 }
767 
780 void scalar_vector_product(double *a, double *result, double scalar, int rows)
781 {
782  int i;
783 
784 #pragma omp for schedule (static)
785  for (i = 0; i < rows; i++) {
786  result[i] = scalar * a[i];
787  }
788 
789  return;
790 }
791 
800 void sync_vectors(double *source, double *target, int rows)
801 {
802  int i;
803 
804 #pragma omp for schedule (static)
805  for (i = 0; i < rows; i++) {
806  target[i] = source[i];
807  }
808 
809  return;
810 }
811 
812 /* ******************************************************* *
813  * ***** Check if matrix is symmetric ******************** *
814  * ******************************************************* */
824 {
825  int i, j, k;
826  double value1 = 0;
827  double value2 = 0;
828  int index;
829  int symm = 0;
830 
831  if (L->quad != 1) {
832  G_warning(_("The linear equation system is not quadratic"));
833  return 0;
834  }
835 
836  G_debug(2, "check_symmetry: Check if matrix is symmetric");
837 
838  if (L->type == N_SPARSE_LES) {
839 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
840  for (j = 0; j < L->rows; j++) {
841  for (i = 1; i < L->Asp[j]->cols; i++) {
842  value1 = 0;
843  value2 = 0;
844  index = 0;
845 
846  index = L->Asp[j]->index[i];
847  value1 = L->Asp[j]->values[i];
848 
849  for (k = 1; k < L->Asp[index]->cols; k++) {
850  if (L->Asp[index]->index[k] == j) {
851  value2 = L->Asp[index]->values[k];
852  if ((value1 != value2)) {
853  if ((fabs((fabs(value1) - fabs(value2))) <
854  SYMM_TOLERANCE)) {
855  G_debug(5,
856  "check_symmetry: sparse matrix is unsymmetric, but within tolerance");
857  }
858  else {
859  G_warning
860  ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf \ndifference = %12.18lf\nStop symmetry calculation.\n",
861  j, index, index, L->Asp[index]->index[k],
862  value1, value2,
863  fabs(fabs(value1) - fabs(value2)));
864  symm++;
865  }
866  }
867  }
868  }
869  }
870  }
871  }
872  else {
873 #pragma omp parallel for schedule (static) private(i, j, k, value1, value2, index) reduction(+:symm) shared(L)
874  for (i = 0; i < L->rows; i++) {
875  for (j = i + 1; j < L->rows; j++) {
876  if (L->A[i][j] != L->A[j][i]) {
877  if ((fabs(fabs(L->A[i][j]) - fabs(L->A[j][i])) <
878  SYMM_TOLERANCE)) {
879  G_debug(5,
880  "check_symmetry: matrix is unsymmetric, but within tolerance");
881  }
882  else {
883  G_warning
884  ("Matrix unsymmetric: Position [%i][%i] : [%i][%i] \nError: %12.18lf != %12.18lf\ndifference = %12.18lf\nStop symmetry calculation.\n",
885  i, j, j, i, L->A[i][j], L->A[j][i],
886  fabs(fabs(L->A[i][j]) - fabs(L->A[j][i])));
887  symm++;
888  }
889  }
890  }
891  }
892  }
893 
894  if (symm > 0)
895  return 0;
896 
897  return 1;
898 }
899 
900 
910 {
911  N_les *L_new;
912  int rows = L->rows;
913  int cols = L->cols;
914  int i, j;
915  double sum;
916 
917  L_new = N_alloc_les_A(rows, N_SPARSE_LES);
918 
919  if (L->type == N_NORMAL_LES) {
920 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
921  for (i = 0; i < rows; i++) {
922  N_spvector *spvect = N_alloc_spvector(1);
923 
924  switch (prec) {
926  sum = 0;
927  for (j = 0; j < cols; j++)
928  sum += L->A[i][j] * L->A[i][j];
929  spvect->values[0] = 1.0 / sqrt(sum);
930  break;
932  sum = 0;
933  for (j = 0; j < cols; j++)
934  sum += fabs(L->A[i][j]);
935  spvect->values[0] = 1.0 / (sum);
936  break;
938  spvect->values[0] = 1.0 / L->A[i][i];
939  break;
940  default:
941  spvect->values[0] = 1.0 / L->A[i][i];
942  }
943 
944 
945  spvect->index[0] = i;
946  spvect->cols = 1;;
947  N_add_spvector_to_les(L_new, spvect, i);
948 
949  }
950  }
951  else {
952 #pragma omp parallel for schedule (static) private(i, sum) shared(L_new, L, rows, prec)
953  for (i = 0; i < rows; i++) {
954  N_spvector *spvect = N_alloc_spvector(1);
955 
956  switch (prec) {
958  sum = 0;
959  for (j = 0; j < L->Asp[i]->cols; j++)
960  sum += L->Asp[i]->values[j] * L->Asp[i]->values[j];
961  spvect->values[0] = 1.0 / sqrt(sum);
962  break;
964  sum = 0;
965  for (j = 0; j < L->Asp[i]->cols; j++)
966  sum += fabs(L->Asp[i]->values[j]);
967  spvect->values[0] = 1.0 / (sum);
968  break;
970  spvect->values[0] = 1.0 / L->Asp[i]->values[0];
971  break;
972  default:
973  spvect->values[0] = 1.0 / L->Asp[i]->values[0];
974  }
975 
976  spvect->index[0] = i;
977  spvect->cols = 1;;
978  N_add_spvector_to_les(L_new, spvect, i);
979  }
980  }
981  return L_new;
982 }