GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
solvers_krylov.c
Go to the documentation of this file.
1 
2 /*****************************************************************************
3 *
4 * MODULE: Grass numerical math interface
5 * AUTHOR(S): Soeren Gebbert, Berlin (GER) Dec 2006
6 * soerengebbert <at> googlemail <dot> com
7 *
8 * PURPOSE: linear equation system solvers
9 * part of the gmath library
10 *
11 * COPYRIGHT: (C) 2010 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 static G_math_spvector **create_diag_precond_matrix(double **A,
28  G_math_spvector ** Asp,
29  int rows, int prec);
30 static int solver_pcg(double **A, G_math_spvector ** Asp, double *x,
31  double *b, int rows, int maxit, double err, int prec, int has_band, int bandwidth);
32 static int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
33  int rows, int maxit, double err, int has_band, int bandwidth);
34 static int solver_bicgstab(double **A, G_math_spvector ** Asp, double *x,
35  double *b, int rows, int maxit, double err);
36 
37 
38 /*!
39  * \brief The iterative preconditioned conjugate gradients solver for symmetric positive definite matrices
40  *
41  * This iterative solver works with symmetric positive definite regular quadratic matrices.
42  *
43  * This solver solves the linear equation system:
44  * A x = b
45  *
46  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
47  * solver will abort the calculation and writes the current result into the vector x.
48  * The parameter <i>err</i> defines the error break criteria for the solver.
49  *
50  * \param A (double **) -- the matrix
51  * \param x (double *) -- the value vector
52  * \param b (double *) -- the right hand side
53  * \param rows (int)
54  * \param maxit (int) -- the maximum number of iterations
55  * \param err (double) -- defines the error break criteria
56  * \param prec (int) -- the preconditioner which should be used 1,2 or 3
57  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
58  *
59  * */
60 int G_math_solver_pcg(double **A, double *x, double *b, int rows, int maxit,
61  double err, int prec)
62 {
63 
64  return solver_pcg(A, NULL, x, b, rows, maxit, err, prec, 0, 0);
65 }
66 
67 /*!
68  * \brief The iterative preconditioned conjugate gradients solver for symmetric positive definite band matrices
69  *
70  * WARNING: The preconditioning of symmetric band matrices is not implemented yet
71  *
72  * This iterative solver works with symmetric positive definite band matrices.
73  *
74  * This solver solves the linear equation system:
75  * A x = b
76  *
77  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
78  * solver will abort the calculation and writes the current result into the vector x.
79  * The parameter <i>err</i> defines the error break criteria for the solver.
80  *
81  * \param A (double **) -- the positive definite band matrix
82  * \param x (double *) -- the value vector
83  * \param b (double *) -- the right hand side
84  * \param rows (int)
85  * \param bandwidth (int) -- bandwidth of matrix A
86  * \param maxit (int) -- the maximum number of iterations
87  * \param err (double) -- defines the error break criteria
88  * \param prec (int) -- the preconditioner which should be used 1,2 or 3
89  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
90  *
91  * */
92 int G_math_solver_pcg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit,
93  double err, int prec)
94 {
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);
97 }
98 
99 
100 /*!
101  * \brief The iterative preconditioned conjugate gradients solver for sparse symmetric positive definite matrices
102  *
103  * This iterative solver works with symmetric positive definite sparse matrices.
104  *
105  * This solver solves the linear equation system:
106  * A x = b
107  *
108  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
109  * solver will abort the calculation and writes the current result into the vector x.
110  * The parameter <i>err</i> defines the error break criteria for the solver.
111  *
112  * \param Asp (G_math_spvector **) -- the sparse matrix
113  * \param x (double *) -- the value vector
114  * \param b (double *) -- the right hand side
115  * \param rows (int)
116  * \param maxit (int) -- the maximum number of iterations
117  * \param err (double) -- defines the error break criteria
118  * \param prec (int) -- the preconditioner which should be used 1,2 or 3
119  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
120  *
121  * */
122 int G_math_solver_sparse_pcg(G_math_spvector ** Asp, double *x, double *b,
123  int rows, int maxit, double err, int prec)
124 {
125 
126  return solver_pcg(NULL, Asp, x, b, rows, maxit, err, prec, 0, 0);
127 }
128 
129 int solver_pcg(double **A, G_math_spvector ** Asp, double *x, double *b,
130  int rows, int maxit, double err, int prec, int has_band, int bandwidth)
131 {
132  double *r, *z;
133 
134  double *p;
135 
136  double *v;
137 
138  double s = 0.0;
139 
140  double a0 = 0, a1 = 0, mygamma, tmp = 0;
141 
142  int m, i;
143 
144  int finished = 2;
145 
146  int error_break;
147 
148  G_math_spvector **M;
149 
150  r = G_alloc_vector(rows);
151  p = G_alloc_vector(rows);
152  v = G_alloc_vector(rows);
153  z = G_alloc_vector(rows);
154 
155  error_break = 0;
156 
157  /*compute the preconditioning matrix, this is a sparse matrix */
158  M = create_diag_precond_matrix(A, Asp, rows, prec);
159 
160  /*
161  * residual calculation
162  */
163 #pragma omp parallel
164  {
165  if (Asp)
166  G_math_Ax_sparse(Asp, x, v, rows);
167  else if(has_band)
168  G_math_Ax_sband(A, x, v, rows, bandwidth);
169  else
170  G_math_d_Ax(A, x, v, rows, rows);
171 
172  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
173  /*performe the preconditioning */
174  G_math_Ax_sparse(M, r, p, rows);
175 
176  /* scalar product */
177 #pragma omp for schedule (static) private(i) reduction(+:s)
178  for (i = 0; i < rows; i++) {
179  s += p[i] * r[i];
180  }
181  }
182 
183  a0 = s;
184  s = 0.0;
185 
186  /* ******************* */
187  /* start the iteration */
188  /* ******************* */
189  for (m = 0; m < maxit; m++) {
190 #pragma omp parallel default(shared)
191  {
192  if (Asp)
193  G_math_Ax_sparse(Asp, p, v, rows);
194  else if(has_band)
195  G_math_Ax_sband(A, p, v, rows, bandwidth);
196  else
197  G_math_d_Ax(A, p, v, rows, rows);
198 
199 
200 
201  /* scalar product */
202 #pragma omp for schedule (static) private(i) reduction(+:s)
203  for (i = 0; i < rows; i++) {
204  s += v[i] * p[i];
205  }
206 
207  /* barrier */
208 #pragma omp single
209  {
210  tmp = s;
211  mygamma = a0 / tmp;
212  s = 0.0;
213  }
214 
215  G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
216 
217  if (m % 50 == 1) {
218  if (Asp)
219  G_math_Ax_sparse(Asp, x, v, rows);
220  else if(has_band)
221  G_math_Ax_sband(A, x, v, rows, bandwidth);
222  else
223  G_math_d_Ax(A, x, v, rows, rows);
224 
225  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
226  }
227  else {
228  G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
229  }
230 
231  /*performe the preconditioning */
232  G_math_Ax_sparse(M, r, z, rows);
233 
234 
235  /* scalar product */
236 #pragma omp for schedule (static) private(i) reduction(+:s)
237  for (i = 0; i < rows; i++) {
238  s += z[i] * r[i];
239  }
240 
241  /* barrier */
242 #pragma omp single
243  {
244  a1 = s;
245  tmp = a1 / a0;
246  a0 = a1;
247  s = 0.0;
248 
249  if (a1 < 0 || a1 == 0 || a1 > 0) {
250  ;
251  }
252  else {
253  G_warning(_
254  ("Unable to solve the linear equation system"));
255  error_break = 1;
256  }
257  }
258  G_math_d_ax_by(p, z, p, tmp, 1.0, rows);
259  }
260 
261  if (Asp != NULL)
262  G_message(_("Sparse PCG -- iteration %i error %g\n"), m, a0);
263  else
264  G_message(_("PCG -- iteration %i error %g\n"), m, a0);
265 
266  if (error_break == 1) {
267  finished = -1;
268  break;
269  }
270 
271 
272  if (a0 < err) {
273  finished = 1;
274  break;
275  }
276  }
277 
278  G_free(r);
279  G_free(p);
280  G_free(v);
281  G_free(z);
282  G_math_free_spmatrix(M, rows);
283 
284  return finished;
285 }
286 
287 
288 /*!
289  * \brief The iterative conjugate gradients solver for symmetric positive definite matrices
290  *
291  * This iterative solver works with symmetric positive definite regular quadratic matrices.
292  *
293  * This solver solves the linear equation system:
294  * A x = b
295  *
296  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
297  * solver will abort the calculation and writes the current result into the vector x.
298  * The parameter <i>err</i> defines the error break criteria for the solver.
299  *
300  * \param A (double **) -- the matrix
301  * \param x (double *) -- the value vector
302  * \param b (double *) -- the right hand side
303  * \param rows (int)
304  * \param maxit (int) -- the maximum number of iterations
305  * \param err (double) -- defines the error break criteria
306  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
307  *
308  * */
309 int G_math_solver_cg(double **A, double *x, double *b, int rows, int maxit,
310  double err)
311 {
312  return solver_cg(A, NULL, x, b, rows, maxit, err, 0, 0);
313 }
314 
315 /*!
316  * \brief The iterative conjugate gradients solver for symmetric positive definite band matrices
317  *
318  * This iterative solver works with symmetric positive definite band matrices.
319  *
320  * This solver solves the linear equation system:
321  * A x = b
322  *
323  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
324  * solver will abort the calculation and writes the current result into the vector x.
325  * The parameter <i>err</i> defines the error break criteria for the solver.
326  *
327  * \param A (double **) -- the symmetric positive definite band matrix
328  * \param x (double *) -- the value vector
329  * \param b (double *) -- the right hand side
330  * \param rows (int)
331  * \param bandwidth (int) -- the bandwidth of matrix A
332  * \param maxit (int) -- the maximum number of iterations
333  * \param err (double) -- defines the error break criteria
334  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
335  *
336  * */
337 int G_math_solver_cg_sband(double **A, double *x, double *b, int rows, int bandwidth, int maxit, double err)
338 {
339  return solver_cg(A, NULL, x, b, rows, maxit, err, 1, bandwidth);
340 }
341 
342 
343 /*!
344  * \brief The iterative conjugate gradients solver for sparse symmetric positive definite matrices
345  *
346  * This iterative solver works with symmetric positive definite sparse matrices.
347  *
348  * This solver solves the linear equation system:
349  * A x = b
350  *
351  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
352  * solver will abort the calculation and writes the current result into the vector x.
353  * The parameter <i>err</i> defines the error break criteria for the solver.
354  *
355  * \param Asp (G_math_spvector **) -- the sparse matrix
356  * \param x (double *) -- the value vector
357  * \param b (double *) -- the right hand side
358  * \param rows (int)
359  * \param maxit (int) -- the maximum number of iterations
360  * \param err (double) -- defines the error break criteria
361  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
362  *
363  * */
364 int G_math_solver_sparse_cg(G_math_spvector ** Asp, double *x, double *b,
365  int rows, int maxit, double err)
366 {
367  return solver_cg(NULL, Asp, x, b, rows, maxit, err, 0, 0);
368 }
369 
370 
371 int solver_cg(double **A, G_math_spvector ** Asp, double *x, double *b,
372  int rows, int maxit, double err, int has_band, int bandwidth)
373 {
374  double *r;
375 
376  double *p;
377 
378  double *v;
379 
380  double s = 0.0;
381 
382  double a0 = 0, a1 = 0, mygamma, tmp = 0;
383 
384  int m, i;
385 
386  int finished = 2;
387 
388  int error_break;
389 
390  r = G_alloc_vector(rows);
391  p = G_alloc_vector(rows);
392  v = G_alloc_vector(rows);
393 
394  error_break = 0;
395  /*
396  * residual calculation
397  */
398 #pragma omp parallel
399  {
400  if (Asp)
401  G_math_Ax_sparse(Asp, x, v, rows);
402  else if(has_band)
403  G_math_Ax_sband(A, x, v, rows, bandwidth);
404  else
405  G_math_d_Ax(A, x, v, rows, rows);
406 
407  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
408  G_math_d_copy(r, p, rows);
409 
410  /* scalar product */
411 #pragma omp for schedule (static) private(i) reduction(+:s)
412  for (i = 0; i < rows; i++) {
413  s += r[i] * r[i];
414  }
415  }
416 
417  a0 = s;
418  s = 0.0;
419 
420  /* ******************* */
421  /* start the iteration */
422  /* ******************* */
423  for (m = 0; m < maxit; m++) {
424 #pragma omp parallel default(shared)
425  {
426  if (Asp)
427  G_math_Ax_sparse(Asp, p, v, rows);
428  else if(has_band)
429  G_math_Ax_sband(A, p, v, rows, bandwidth);
430  else
431  G_math_d_Ax(A, p, v, rows, rows);
432 
433  /* scalar product */
434 #pragma omp for schedule (static) private(i) reduction(+:s)
435  for (i = 0; i < rows; i++) {
436  s += v[i] * p[i];
437  }
438 
439  /* barrier */
440 #pragma omp single
441  {
442  tmp = s;
443  mygamma = a0 / tmp;
444  s = 0.0;
445  }
446 
447  G_math_d_ax_by(p, x, x, mygamma, 1.0, rows);
448 
449  if (m % 50 == 1) {
450  if (Asp)
451  G_math_Ax_sparse(Asp, x, v, rows);
452  else if(has_band)
453  G_math_Ax_sband(A, x, v, rows, bandwidth);
454  else
455  G_math_d_Ax(A, x, v, rows, rows);
456 
457  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
458  }
459  else {
460  G_math_d_ax_by(r, v, r, 1.0, -1.0 * mygamma, rows);
461  }
462 
463  /* scalar product */
464 #pragma omp for schedule (static) private(i) reduction(+:s)
465  for (i = 0; i < rows; i++) {
466  s += r[i] * r[i];
467  }
468 
469  /* barrier */
470 #pragma omp single
471  {
472  a1 = s;
473  tmp = a1 / a0;
474  a0 = a1;
475  s = 0.0;
476 
477  if (a1 < 0 || a1 == 0 || a1 > 0) {
478  ;
479  }
480  else {
481  G_warning(_
482  ("Unable to solve the linear equation system"));
483  error_break = 1;
484  }
485  }
486  G_math_d_ax_by(p, r, p, tmp, 1.0, rows);
487  }
488 
489  if (Asp != NULL)
490  G_message(_("Sparse CG -- iteration %i error %g\n"), m, a0);
491  else
492  G_message(_("CG -- iteration %i error %g\n"), m, a0);
493 
494  if (error_break == 1) {
495  finished = -1;
496  break;
497  }
498 
499  if (a0 < err) {
500  finished = 1;
501  break;
502  }
503  }
504 
505  G_free(r);
506  G_free(p);
507  G_free(v);
508 
509  return finished;
510 }
511 
512 
513 
514 /*!
515  * \brief The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices
516  *
517  * This iterative solver works with regular quadratic matrices.
518  *
519  * This solver solves the linear equation system:
520  * A x = b
521  *
522  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
523  * solver will abort the calculation and writes the current result into the vector x.
524  * The parameter <i>err</i> defines the error break criteria for the solver.
525  *
526  * \param A (double **) -- the matrix
527  * \param x (double *) -- the value vector
528  * \param b (double *) -- the right hand side
529  * \param rows (int)
530  * \param maxit (int) -- the maximum number of iterations
531  * \param err (double) -- defines the error break criteria
532  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
533  *
534  * */
535 int G_math_solver_bicgstab(double **A, double *x, double *b, int rows,
536  int maxit, double err)
537 {
538  return solver_bicgstab(A, NULL, x, b, rows, maxit, err);
539 }
540 
541 /*!
542  * \brief The iterative biconjugate gradients solver with stabilization for unsymmetric non-definite matrices
543  *
544  * This iterative solver works with sparse matrices.
545  *
546  * This solver solves the linear equation system:
547  * A x = b
548  *
549  * The parameter <i>maxit</i> specifies the maximum number of iterations. If the maximum is reached, the
550  * solver will abort the calculation and writes the current result into the vector x.
551  * The parameter <i>err</i> defines the error break criteria for the solver.
552  *
553  * \param Asp (G_math_spvector **) -- the sparse matrix
554  * \param x (double *) -- the value vector
555  * \param b (double *) -- the right hand side
556  * \param rows (int)
557  * \param maxit (int) -- the maximum number of iterations
558  * \param err (double) -- defines the error break criteria
559  * \return (int) -- 1 - success, 2 - not finished but success, 0 - matrix singular, -1 - could not solve the les
560  *
561  * */
563  double *b, int rows, int maxit, double err)
564 {
565  return solver_bicgstab(NULL, Asp, x, b, rows, maxit, err);
566 }
567 
568 
569 int solver_bicgstab(double **A, G_math_spvector ** Asp, double *x, double *b,
570  int rows, int maxit, double err)
571 {
572  double *r;
573 
574  double *r0;
575 
576  double *p;
577 
578  double *v;
579 
580  double *s;
581 
582  double *t;
583 
584  double s1 = 0.0, s2 = 0.0, s3 = 0.0;
585 
586  double alpha = 0, beta = 0, omega, rr0 = 0, error;
587 
588  int m, i;
589 
590  int finished = 2;
591 
592  int error_break;
593 
594  r = G_alloc_vector(rows);
595  r0 = G_alloc_vector(rows);
596  p = G_alloc_vector(rows);
597  v = G_alloc_vector(rows);
598  s = G_alloc_vector(rows);
599  t = G_alloc_vector(rows);
600 
601  error_break = 0;
602 
603 #pragma omp parallel
604  {
605  if (Asp)
606  G_math_Ax_sparse(Asp, x, v, rows);
607  else
608  G_math_d_Ax(A, x, v, rows, rows);
609 
610  G_math_d_ax_by(b, v, r, 1.0, -1.0, rows);
611  G_math_d_copy(r, r0, rows);
612  G_math_d_copy(r, p, rows);
613  }
614 
615  s1 = s2 = s3 = 0.0;
616 
617  /* ******************* */
618  /* start the iteration */
619  /* ******************* */
620  for (m = 0; m < maxit; m++) {
621 
622 #pragma omp parallel default(shared)
623  {
624  if (Asp)
625  G_math_Ax_sparse(Asp, p, v, rows);
626  else
627  G_math_d_Ax(A, p, v, rows, rows);
628 
629  /* scalar product */
630 #pragma omp for schedule (static) private(i) reduction(+:s1, s2, s3)
631  for (i = 0; i < rows; i++) {
632  s1 += r[i] * r[i];
633  s2 += r[i] * r0[i];
634  s3 += v[i] * r0[i];
635  }
636 
637 #pragma omp single
638  {
639  error = s1;
640 
641  if (error < 0 || error == 0 || error > 0) {
642  ;
643  }
644  else {
645  G_warning(_
646  ("Unable to solve the linear equation system"));
647  error_break = 1;
648  }
649 
650  rr0 = s2;
651  alpha = rr0 / s3;
652  s1 = s2 = s3 = 0.0;
653  }
654 
655  G_math_d_ax_by(r, v, s, 1.0, -1.0 * alpha, rows);
656  if (Asp)
657  G_math_Ax_sparse(Asp, s, t, rows);
658  else
659  G_math_d_Ax(A, s, t, rows, rows);
660 
661  /* scalar product */
662 #pragma omp for schedule (static) private(i) reduction(+:s1, s2)
663  for (i = 0; i < rows; i++) {
664  s1 += t[i] * s[i];
665  s2 += t[i] * t[i];
666  }
667 
668 #pragma omp single
669  {
670  omega = s1 / s2;
671  s1 = s2 = 0.0;
672  }
673 
674  G_math_d_ax_by(p, s, r, alpha, omega, rows);
675  G_math_d_ax_by(x, r, x, 1.0, 1.0, rows);
676  G_math_d_ax_by(s, t, r, 1.0, -1.0 * omega, rows);
677 
678 #pragma omp for schedule (static) private(i) reduction(+:s1)
679  for (i = 0; i < rows; i++) {
680  s1 += r[i] * r0[i];
681  }
682 
683 #pragma omp single
684  {
685  beta = alpha / omega * s1 / rr0;
686  s1 = s2 = s3 = 0.0;
687  }
688 
689  G_math_d_ax_by(p, v, p, 1.0, -1.0 * omega, rows);
690  G_math_d_ax_by(p, r, p, beta, 1.0, rows);
691  }
692 
693 
694  if (Asp != NULL)
695  G_message(_("Sparse BiCGStab -- iteration %i error %g\n"), m,
696  error);
697  else
698  G_message(_("BiCGStab -- iteration %i error %g\n"), m, error);
699 
700  if (error_break == 1) {
701  finished = -1;
702  break;
703  }
704 
705  if (error < err) {
706  finished = 1;
707  break;
708  }
709  }
710 
711  G_free(r);
712  G_free(r0);
713  G_free(p);
714  G_free(v);
715  G_free(s);
716  G_free(t);
717 
718  return finished;
719 }
720 
721 
722 /*!
723  * \brief Compute a diagonal preconditioning matrix for krylov space solver
724  *
725  * \param A (double **) -- the matrix for which the precondition should be computed (if the sparse matrix is used, set it to NULL)
726  * \param Asp (G_math_spvector **) -- the matrix for which the precondition should be computed
727  * \param rows (int)
728  * \param prec (int) -- which preconditioner should be used 1, 2 or 3
729  *
730  * */
731 G_math_spvector **create_diag_precond_matrix(double **A,
732  G_math_spvector ** Asp, int rows,
733  int prec)
734 {
735  G_math_spvector **Msp;
736 
737  int i, j, cols = rows;
738 
739  double sum;
740 
741  Msp = G_math_alloc_spmatrix(rows);
742 
743  if (A != NULL) {
744 #pragma omp parallel for schedule (static) private(i, j, sum) shared(A, Msp, rows, cols, prec)
745  for (i = 0; i < rows; i++) {
747 
748  switch (prec) {
750  sum = 0;
751  for (j = 0; j < cols; j++)
752  sum += A[i][j] * A[i][j];
753  spvect->values[0] = 1.0 / sqrt(sum);
754  break;
756  sum = 0;
757  for (j = 0; j < cols; j++)
758  sum += fabs(A[i][j]);
759  spvect->values[0] = 1.0 / (sum);
760  break;
762  default:
763  spvect->values[0] = 1.0 / A[i][i];
764  break;
765  }
766 
767 
768  spvect->index[0] = i;
769  spvect->cols = 1;;
770  G_math_add_spvector(Msp, spvect, i);
771 
772  }
773  }
774  else {
775 #pragma omp parallel for schedule (static) private(i, j, sum) shared(Asp, Msp, rows, cols, prec)
776  for (i = 0; i < rows; i++) {
778 
779  switch (prec) {
781  sum = 0;
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);
785  break;
787  sum = 0;
788  for (j = 0; j < Asp[i]->cols; j++)
789  sum += fabs(Asp[i]->values[j]);
790  spvect->values[0] = 1.0 / (sum);
791  break;
793  default:
794  for (j = 0; j < Asp[i]->cols; j++)
795  if (i == Asp[i]->index[j])
796  spvect->values[0] = 1.0 / Asp[i]->values[j];
797  break;
798  }
799 
800  spvect->index[0] = i;
801  spvect->cols = 1;;
802  G_math_add_spvector(Msp, spvect, i);
803  }
804  }
805  return Msp;
806 }
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.
Definition: sparse_matrix.c:76
The row vector of the sparse matrix.
Definition: gmath.h:54
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.
Definition: gis/alloc.c:149
void G_math_d_copy(double *, double *, int)
Copy the vector x to y.
Definition: blas_level_1.c:237
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.
Definition: blas_level_1.c:173
#define NULL
Definition: ccmath.h:32
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...
#define x
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)
Definition: symbol/read.c:220
void G_message(const char *,...) __attribute__((format(printf
unsigned int * index
Definition: gmath.h:58
double t
Definition: r_raster.c:39
double b
Definition: r_raster.c:39
double * G_alloc_vector(size_t)
Vector matrix memory allocation.
Definition: dalloc.c:41
double * values
Definition: gmath.h:56
#define G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION
Definition: gmath.h:47
unsigned int cols
Definition: gmath.h:57
G_math_spvector ** G_math_alloc_spmatrix(int)
Allocate memory for a sparse matrix.
Definition: sparse_matrix.c:58
#define M(row, col)
Definition: georef.c:49
#define G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION
Definition: gmath.h:48
void G_warning(const char *,...) __attribute__((format(printf
#define _(str)
Definition: glocale.h:10
#define G_MATH_DIAGONAL_PRECONDITION
Definition: gmath.h:46
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.
Definition: blas_level_2.c:47
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.
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...
double r
Definition: r_raster.c:39
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.