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