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 * */
G_math_solver_pcg(double ** A,double * x,double * b,int rows,int maxit,double err,int prec)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 * */
G_math_solver_pcg_sband(double ** A,double * x,double * b,int rows,int bandwidth,int maxit,double err,int prec)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 * */
G_math_solver_sparse_pcg(G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err,int prec)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
solver_pcg(double ** A,G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err,int prec,int has_band,int bandwidth)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 * */
G_math_solver_cg(double ** A,double * x,double * b,int rows,int maxit,double err)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 * */
G_math_solver_cg_sband(double ** A,double * x,double * b,int rows,int bandwidth,int maxit,double err)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 * */
G_math_solver_sparse_cg(G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err)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
solver_cg(double ** A,G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err,int has_band,int bandwidth)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 * */
G_math_solver_bicgstab(double ** A,double * x,double * b,int rows,int maxit,double err)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 * */
G_math_solver_sparse_bicgstab(G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err)562 int G_math_solver_sparse_bicgstab(G_math_spvector ** Asp, double *x,
563 double *b, int rows, int maxit, double err)
564 {
565 return solver_bicgstab(NULL, Asp, x, b, rows, maxit, err);
566 }
567
568
solver_bicgstab(double ** A,G_math_spvector ** Asp,double * x,double * b,int rows,int maxit,double err)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 * */
create_diag_precond_matrix(double ** A,G_math_spvector ** Asp,int rows,int prec)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++) {
746 G_math_spvector *spvect = G_math_alloc_spvector(1);
747
748 switch (prec) {
749 case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
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;
755 case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
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;
761 case G_MATH_DIAGONAL_PRECONDITION:
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++) {
777 G_math_spvector *spvect = G_math_alloc_spvector(1);
778
779 switch (prec) {
780 case G_MATH_ROWSCALE_EUKLIDNORM_PRECONDITION:
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;
786 case G_MATH_ROWSCALE_ABSSUMNORM_PRECONDITION:
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;
792 case G_MATH_DIAGONAL_PRECONDITION:
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 }
807