1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 *
10 * Preconditioned conjugate gradient (Omin) functions
11 *
12 *****************************************************************************/
13
14 /* This was based on the pcg.c formerly in struct_ls, with
15 changes (GetPrecond and stop_crit) for compatibility with the pcg.c
16 in parcsr_ls and elsewhere. Incompatibilities with the
17 parcsr_ls version:
18 - logging is different; no attempt has been made to be the same
19 - treatment of b=0 in Ax=b is different: this returns x=0; the parcsr
20 version iterates with a special stopping criterion
21 */
22
23 #include "krylov.h"
24 #include "_hypre_utilities.h"
25
26 /*--------------------------------------------------------------------------
27 * hypre_PCGFunctionsCreate
28 *--------------------------------------------------------------------------*/
29
30 hypre_PCGFunctions *
hypre_PCGFunctionsCreate(void * (* CAlloc)(size_t count,size_t elt_size,HYPRE_MemoryLocation location),HYPRE_Int (* Free)(void * ptr),HYPRE_Int (* CommInfo)(void * A,HYPRE_Int * my_id,HYPRE_Int * num_procs),void * (* CreateVector)(void * vector),HYPRE_Int (* DestroyVector)(void * vector),void * (* MatvecCreate)(void * A,void * x),HYPRE_Int (* Matvec)(void * matvec_data,HYPRE_Complex alpha,void * A,void * x,HYPRE_Complex beta,void * y),HYPRE_Int (* MatvecDestroy)(void * matvec_data),HYPRE_Real (* InnerProd)(void * x,void * y),HYPRE_Int (* CopyVector)(void * x,void * y),HYPRE_Int (* ClearVector)(void * x),HYPRE_Int (* ScaleVector)(HYPRE_Complex alpha,void * x),HYPRE_Int (* Axpy)(HYPRE_Complex alpha,void * x,void * y),HYPRE_Int (* PrecondSetup)(void * vdata,void * A,void * b,void * x),HYPRE_Int (* Precond)(void * vdata,void * A,void * b,void * x))31 hypre_PCGFunctionsCreate(
32 void * (*CAlloc) ( size_t count, size_t elt_size, HYPRE_MemoryLocation location ),
33 HYPRE_Int (*Free) ( void *ptr ),
34 HYPRE_Int (*CommInfo) ( void *A, HYPRE_Int *my_id,
35 HYPRE_Int *num_procs ),
36 void * (*CreateVector) ( void *vector ),
37 HYPRE_Int (*DestroyVector) ( void *vector ),
38 void * (*MatvecCreate) ( void *A, void *x ),
39 HYPRE_Int (*Matvec) ( void *matvec_data, HYPRE_Complex alpha, void *A,
40 void *x, HYPRE_Complex beta, void *y ),
41 HYPRE_Int (*MatvecDestroy) ( void *matvec_data ),
42 HYPRE_Real (*InnerProd) ( void *x, void *y ),
43 HYPRE_Int (*CopyVector) ( void *x, void *y ),
44 HYPRE_Int (*ClearVector) ( void *x ),
45 HYPRE_Int (*ScaleVector) ( HYPRE_Complex alpha, void *x ),
46 HYPRE_Int (*Axpy) ( HYPRE_Complex alpha, void *x, void *y ),
47 HYPRE_Int (*PrecondSetup) ( void *vdata, void *A, void *b, void *x ),
48 HYPRE_Int (*Precond) ( void *vdata, void *A, void *b, void *x )
49 )
50 {
51 hypre_PCGFunctions * pcg_functions;
52 pcg_functions = (hypre_PCGFunctions *)
53 CAlloc( 1, sizeof(hypre_PCGFunctions), HYPRE_MEMORY_HOST );
54
55 pcg_functions->CAlloc = CAlloc;
56 pcg_functions->Free = Free;
57 pcg_functions->CommInfo = CommInfo;
58 pcg_functions->CreateVector = CreateVector;
59 pcg_functions->DestroyVector = DestroyVector;
60 pcg_functions->MatvecCreate = MatvecCreate;
61 pcg_functions->Matvec = Matvec;
62 pcg_functions->MatvecDestroy = MatvecDestroy;
63 pcg_functions->InnerProd = InnerProd;
64 pcg_functions->CopyVector = CopyVector;
65 pcg_functions->ClearVector = ClearVector;
66 pcg_functions->ScaleVector = ScaleVector;
67 pcg_functions->Axpy = Axpy;
68 /* default preconditioner must be set here but can be changed later... */
69 pcg_functions->precond_setup = PrecondSetup;
70 pcg_functions->precond = Precond;
71
72 return pcg_functions;
73 }
74
75 /*--------------------------------------------------------------------------
76 * hypre_PCGCreate
77 *--------------------------------------------------------------------------*/
78
79 void *
hypre_PCGCreate(hypre_PCGFunctions * pcg_functions)80 hypre_PCGCreate( hypre_PCGFunctions *pcg_functions )
81 {
82 hypre_PCGData *pcg_data;
83
84 HYPRE_ANNOTATE_FUNC_BEGIN;
85
86 pcg_data = hypre_CTAllocF(hypre_PCGData, 1, pcg_functions, HYPRE_MEMORY_HOST);
87
88 pcg_data -> functions = pcg_functions;
89
90 /* set defaults */
91 (pcg_data -> tol) = 1.0e-06;
92 (pcg_data -> atolf) = 0.0;
93 (pcg_data -> cf_tol) = 0.0;
94 (pcg_data -> a_tol) = 0.0;
95 (pcg_data -> rtol) = 0.0;
96 (pcg_data -> max_iter) = 1000;
97 (pcg_data -> two_norm) = 0;
98 (pcg_data -> rel_change) = 0;
99 (pcg_data -> recompute_residual) = 0;
100 (pcg_data -> recompute_residual_p) = 0;
101 (pcg_data -> stop_crit) = 0;
102 (pcg_data -> converged) = 0;
103 (pcg_data -> hybrid) = 0;
104 (pcg_data -> owns_matvec_data ) = 1;
105 (pcg_data -> matvec_data) = NULL;
106 (pcg_data -> precond_data) = NULL;
107 (pcg_data -> print_level) = 0;
108 (pcg_data -> logging) = 0;
109 (pcg_data -> norms) = NULL;
110 (pcg_data -> rel_norms) = NULL;
111 (pcg_data -> p) = NULL;
112 (pcg_data -> s) = NULL;
113 (pcg_data -> r) = NULL;
114
115 HYPRE_ANNOTATE_FUNC_END;
116
117 return (void *) pcg_data;
118 }
119
120 /*--------------------------------------------------------------------------
121 * hypre_PCGDestroy
122 *--------------------------------------------------------------------------*/
123
124 HYPRE_Int
hypre_PCGDestroy(void * pcg_vdata)125 hypre_PCGDestroy( void *pcg_vdata )
126 {
127 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
128
129 HYPRE_ANNOTATE_FUNC_BEGIN;
130
131 if (pcg_data)
132 {
133 hypre_PCGFunctions *pcg_functions = pcg_data->functions;
134 if ( (pcg_data -> norms) != NULL )
135 {
136 hypre_TFreeF( pcg_data -> norms, pcg_functions );
137 pcg_data -> norms = NULL;
138 }
139 if ( (pcg_data -> rel_norms) != NULL )
140 {
141 hypre_TFreeF( pcg_data -> rel_norms, pcg_functions );
142 pcg_data -> rel_norms = NULL;
143 }
144 if ( pcg_data -> matvec_data != NULL && pcg_data->owns_matvec_data )
145 {
146 (*(pcg_functions->MatvecDestroy))(pcg_data -> matvec_data);
147 pcg_data -> matvec_data = NULL;
148 }
149 if ( pcg_data -> p != NULL )
150 {
151 (*(pcg_functions->DestroyVector))(pcg_data -> p);
152 pcg_data -> p = NULL;
153 }
154 if ( pcg_data -> s != NULL )
155 {
156 (*(pcg_functions->DestroyVector))(pcg_data -> s);
157 pcg_data -> s = NULL;
158 }
159 if ( pcg_data -> r != NULL )
160 {
161 (*(pcg_functions->DestroyVector))(pcg_data -> r);
162 pcg_data -> r = NULL;
163 }
164 hypre_TFreeF( pcg_data, pcg_functions );
165 hypre_TFreeF( pcg_functions, pcg_functions );
166 }
167
168 HYPRE_ANNOTATE_FUNC_END;
169
170 return(hypre_error_flag);
171 }
172
173 /*--------------------------------------------------------------------------
174 * hypre_PCGGetResidual
175 *--------------------------------------------------------------------------*/
176
hypre_PCGGetResidual(void * pcg_vdata,void ** residual)177 HYPRE_Int hypre_PCGGetResidual( void *pcg_vdata, void **residual )
178 {
179 /* returns a pointer to the residual vector */
180
181 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
182 *residual = pcg_data->r;
183 return hypre_error_flag;
184 }
185
186 /*--------------------------------------------------------------------------
187 * hypre_PCGSetup
188 *--------------------------------------------------------------------------*/
189
190 HYPRE_Int
hypre_PCGSetup(void * pcg_vdata,void * A,void * b,void * x)191 hypre_PCGSetup( void *pcg_vdata,
192 void *A,
193 void *b,
194 void *x )
195 {
196 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
197 hypre_PCGFunctions *pcg_functions = pcg_data->functions;
198 HYPRE_Int max_iter = (pcg_data -> max_iter);
199 HYPRE_Int (*precond_setup)(void*,void*,void*,void*) = (pcg_functions -> precond_setup);
200 void *precond_data = (pcg_data -> precond_data);
201
202 HYPRE_ANNOTATE_FUNC_BEGIN;
203
204 (pcg_data -> A) = A;
205
206 /*--------------------------------------------------
207 * The arguments for CreateVector are important to
208 * maintain consistency between the setup and
209 * compute phases of matvec and the preconditioner.
210 *--------------------------------------------------*/
211
212 if ( pcg_data -> p != NULL )
213 (*(pcg_functions->DestroyVector))(pcg_data -> p);
214 (pcg_data -> p) = (*(pcg_functions->CreateVector))(x);
215
216 if ( pcg_data -> s != NULL )
217 (*(pcg_functions->DestroyVector))(pcg_data -> s);
218 (pcg_data -> s) = (*(pcg_functions->CreateVector))(x);
219
220 if ( pcg_data -> r != NULL )
221 (*(pcg_functions->DestroyVector))(pcg_data -> r);
222 (pcg_data -> r) = (*(pcg_functions->CreateVector))(b);
223
224 if ( pcg_data -> matvec_data != NULL && pcg_data->owns_matvec_data )
225 (*(pcg_functions->MatvecDestroy))(pcg_data -> matvec_data);
226 (pcg_data -> matvec_data) = (*(pcg_functions->MatvecCreate))(A, x);
227
228 precond_setup(precond_data, A, b, x);
229
230 /*-----------------------------------------------------
231 * Allocate space for log info
232 *-----------------------------------------------------*/
233
234 if ( (pcg_data->logging)>0 || (pcg_data->print_level)>0 )
235 {
236 if ( (pcg_data -> norms) != NULL )
237 hypre_TFreeF( pcg_data -> norms, pcg_functions );
238 (pcg_data -> norms) = hypre_CTAllocF( HYPRE_Real, max_iter + 1,
239 pcg_functions, HYPRE_MEMORY_HOST);
240
241 if ( (pcg_data -> rel_norms) != NULL )
242 hypre_TFreeF( pcg_data -> rel_norms, pcg_functions );
243 (pcg_data -> rel_norms) = hypre_CTAllocF( HYPRE_Real, max_iter + 1,
244 pcg_functions, HYPRE_MEMORY_HOST );
245 }
246
247 HYPRE_ANNOTATE_FUNC_END;
248
249 return hypre_error_flag;
250 }
251
252 /*--------------------------------------------------------------------------
253 * hypre_PCGSolve
254 *--------------------------------------------------------------------------
255 *
256 * We use the following convergence test as the default (see Ashby, Holst,
257 * Manteuffel, and Saylor):
258 *
259 * ||e||_A ||r||_C
260 * ------- <= [kappa_A(C*A)]^(1/2) ------- < tol
261 * ||x||_A ||b||_C
262 *
263 * where we let (for the time being) kappa_A(CA) = 1.
264 * We implement the test as:
265 *
266 * gamma = <C*r,r>/<C*b,b> < (tol^2) = eps
267 *
268 *--------------------------------------------------------------------------*/
269
270 HYPRE_Int
hypre_PCGSolve(void * pcg_vdata,void * A,void * b,void * x)271 hypre_PCGSolve( void *pcg_vdata,
272 void *A,
273 void *b,
274 void *x )
275 {
276 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
277 hypre_PCGFunctions *pcg_functions = pcg_data->functions;
278
279 HYPRE_Real r_tol = (pcg_data -> tol);
280 HYPRE_Real a_tol = (pcg_data -> a_tol);
281 HYPRE_Real atolf = (pcg_data -> atolf);
282 HYPRE_Real cf_tol = (pcg_data -> cf_tol);
283 HYPRE_Real rtol = (pcg_data -> rtol);
284 HYPRE_Int max_iter = (pcg_data -> max_iter);
285 HYPRE_Int two_norm = (pcg_data -> two_norm);
286 HYPRE_Int rel_change = (pcg_data -> rel_change);
287 HYPRE_Int recompute_residual = (pcg_data -> recompute_residual);
288 HYPRE_Int recompute_residual_p = (pcg_data -> recompute_residual_p);
289 HYPRE_Int stop_crit = (pcg_data -> stop_crit);
290 HYPRE_Int hybrid = (pcg_data -> hybrid);
291 /*
292 HYPRE_Int converged = (pcg_data -> converged);
293 */
294 void *p = (pcg_data -> p);
295 void *s = (pcg_data -> s);
296 void *r = (pcg_data -> r);
297 void *matvec_data = (pcg_data -> matvec_data);
298 HYPRE_Int (*precond)(void*,void*,void*,void*) = (pcg_functions -> precond);
299 void *precond_data = (pcg_data -> precond_data);
300 HYPRE_Int print_level = (pcg_data -> print_level);
301 HYPRE_Int logging = (pcg_data -> logging);
302 HYPRE_Real *norms = (pcg_data -> norms);
303 HYPRE_Real *rel_norms = (pcg_data -> rel_norms);
304
305 HYPRE_Real alpha, beta;
306 HYPRE_Real gamma, gamma_old;
307 HYPRE_Real bi_prod, eps;
308 HYPRE_Real pi_prod, xi_prod;
309 HYPRE_Real ieee_check = 0.;
310
311 HYPRE_Real i_prod = 0.0;
312 HYPRE_Real i_prod_0 = 0.0;
313 HYPRE_Real cf_ave_0 = 0.0;
314 HYPRE_Real cf_ave_1 = 0.0;
315 HYPRE_Real weight;
316 HYPRE_Real ratio;
317
318 HYPRE_Real guard_zero_residual, sdotp;
319 HYPRE_Int tentatively_converged = 0;
320 HYPRE_Int recompute_true_residual = 0;
321
322 HYPRE_Int i = 0;
323 HYPRE_Int my_id, num_procs;
324
325 HYPRE_ANNOTATE_FUNC_BEGIN;
326
327 (pcg_data -> converged) = 0;
328
329 (*(pcg_functions->CommInfo))(A,&my_id,&num_procs);
330
331 /*-----------------------------------------------------------------------
332 * With relative change convergence test on, it is possible to attempt
333 * another iteration with a zero residual. This causes the parameter
334 * alpha to go NaN. The guard_zero_residual parameter is to circumvent
335 * this. Perhaps it should be set to something non-zero (but small).
336 *-----------------------------------------------------------------------*/
337
338 guard_zero_residual = 0.0;
339
340 /*-----------------------------------------------------------------------
341 * Start pcg solve
342 *-----------------------------------------------------------------------*/
343
344 /* compute eps */
345 if (two_norm)
346 {
347 /* bi_prod = <b,b> */
348 bi_prod = (*(pcg_functions->InnerProd))(b, b);
349 if (print_level > 1 && my_id == 0)
350 hypre_printf("<b,b>: %e\n",bi_prod);
351 }
352 else
353 {
354 /* bi_prod = <C*b,b> */
355 (*(pcg_functions->ClearVector))(p);
356 precond(precond_data, A, b, p);
357 bi_prod = (*(pcg_functions->InnerProd))(p, b);
358 if (print_level > 1 && my_id == 0)
359 hypre_printf("<C*b,b>: %e\n",bi_prod);
360 };
361
362 /* Since it is does not diminish performance, attempt to return an error flag
363 and notify users when they supply bad input. */
364 if (bi_prod != 0.) ieee_check = bi_prod/bi_prod; /* INF -> NaN conversion */
365 if (ieee_check != ieee_check)
366 {
367 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
368 for ieee_check self-equality works on all IEEE-compliant compilers/
369 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
370 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
371 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
372 if (print_level > 0 || logging > 0)
373 {
374 hypre_printf("\n\nERROR detected by Hypre ... BEGIN\n");
375 hypre_printf("ERROR -- hypre_PCGSolve: INFs and/or NaNs detected in input.\n");
376 hypre_printf("User probably placed non-numerics in supplied b.\n");
377 hypre_printf("Returning error flag += 101. Program not terminated.\n");
378 hypre_printf("ERROR detected by Hypre ... END\n\n\n");
379 }
380 hypre_error(HYPRE_ERROR_GENERIC);
381 HYPRE_ANNOTATE_FUNC_END;
382
383 return hypre_error_flag;
384 }
385
386 eps = r_tol*r_tol; /* note: this may be re-assigned below */
387 if ( bi_prod > 0.0 )
388 {
389 if ( stop_crit && !rel_change && atolf<=0 ) /* pure absolute tolerance */
390 {
391 eps = eps / bi_prod;
392 /* Note: this section is obsolete. Aside from backwards comatability
393 concerns, we could delete the stop_crit parameter and related code,
394 using tol & atolf instead. */
395 }
396 else if ( atolf>0 ) /* mixed relative and absolute tolerance */
397 {
398 bi_prod += atolf;
399 }
400 else /* DEFAULT (stop_crit and atolf exist for backwards compatibilty
401 and are not in the reference manual) */
402 {
403 /* convergence criteria: <C*r,r> <= max( a_tol^2, r_tol^2 * <C*b,b> )
404 note: default for a_tol is 0.0, so relative residual criteria is used unless
405 user specifies a_tol, or sets r_tol = 0.0, which means absolute
406 tol only is checked */
407 eps = hypre_max(r_tol*r_tol, a_tol*a_tol/bi_prod);
408
409 }
410 }
411 else /* bi_prod==0.0: the rhs vector b is zero */
412 {
413 /* Set x equal to zero and return */
414 (*(pcg_functions->CopyVector))(b, x);
415 if (logging>0 || print_level>0)
416 {
417 norms[0] = 0.0;
418 rel_norms[i] = 0.0;
419 }
420 HYPRE_ANNOTATE_FUNC_END;
421
422 return hypre_error_flag;
423 /* In this case, for the original parcsr pcg, the code would take special
424 action to force iterations even though the exact value was known. */
425 };
426
427 /* r = b - Ax */
428 (*(pcg_functions->CopyVector))(b, r);
429
430 (*(pcg_functions->Matvec))(matvec_data, -1.0, A, x, 1.0, r);
431
432 //hypre_ParVectorUpdateHost(r);
433 /* p = C*r */
434 (*(pcg_functions->ClearVector))(p);
435 precond(precond_data, A, r, p);
436
437 /* gamma = <r,p> */
438 gamma = (*(pcg_functions->InnerProd))(r,p);
439
440 /* Since it is does not diminish performance, attempt to return an error flag
441 and notify users when they supply bad input. */
442 if (gamma != 0.) ieee_check = gamma/gamma; /* INF -> NaN conversion */
443 if (ieee_check != ieee_check)
444 {
445 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
446 for ieee_check self-equality works on all IEEE-compliant compilers/
447 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
448 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
449 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
450 if (print_level > 0 || logging > 0)
451 {
452 hypre_printf("\n\nERROR detected by Hypre ... BEGIN\n");
453 hypre_printf("ERROR -- hypre_PCGSolve: INFs and/or NaNs detected in input.\n");
454 hypre_printf("User probably placed non-numerics in supplied A or x_0.\n");
455 hypre_printf("Returning error flag += 101. Program not terminated.\n");
456 hypre_printf("ERROR detected by Hypre ... END\n\n\n");
457 }
458 hypre_error(HYPRE_ERROR_GENERIC);
459 HYPRE_ANNOTATE_FUNC_END;
460
461 return hypre_error_flag;
462 }
463
464 /* Set initial residual norm */
465 if ( logging>0 || print_level > 0 || cf_tol > 0.0 )
466 {
467 if (two_norm)
468 i_prod_0 = (*(pcg_functions->InnerProd))(r,r);
469 else
470 i_prod_0 = gamma;
471
472 if ( logging>0 || print_level>0 ) norms[0] = sqrt(i_prod_0);
473 }
474 if ( print_level > 1 && my_id==0 )
475 {
476 hypre_printf("\n\n");
477 if (two_norm)
478 {
479 if ( stop_crit && !rel_change && atolf==0 ) /* pure absolute tolerance */
480 {
481 hypre_printf("Iters ||r||_2 conv.rate\n");
482 hypre_printf("----- ------------ ---------\n");
483 }
484 else
485 {
486 hypre_printf("Iters ||r||_2 conv.rate ||r||_2/||b||_2\n");
487 hypre_printf("----- ------------ --------- ------------ \n");
488 }
489 }
490 else /* !two_norm */
491 {
492 hypre_printf("Iters ||r||_C conv.rate ||r||_C/||b||_C\n");
493 hypre_printf("----- ------------ --------- ------------ \n");
494 }
495 /* hypre_printf("% 5d %e\n", i, norms[i]); */
496 }
497
498 while ((i+1) <= max_iter)
499 {
500 /*--------------------------------------------------------------------
501 * the core CG calculations...
502 *--------------------------------------------------------------------*/
503 i++;
504
505 /* At user request, periodically recompute the residual from the formula
506 r = b - A x (instead of using the recursive definition). Note that this
507 is potentially expensive and can lead to degraded convergence (since it
508 essentially a "restarted CG"). */
509 recompute_true_residual = recompute_residual_p && !(i%recompute_residual_p);
510
511 /* s = A*p */
512 (*(pcg_functions->Matvec))(matvec_data, 1.0, A, p, 0.0, s);
513
514 /* alpha = gamma / <s,p> */
515 sdotp = (*(pcg_functions->InnerProd))(s, p);
516 if ( sdotp==0.0 )
517 {
518 hypre_error_w_msg(HYPRE_ERROR_CONV, "Zero sdotp value in PCG");
519 if (i==1) i_prod=i_prod_0;
520 break;
521 }
522 alpha = gamma / sdotp;
523 if (! (alpha > HYPRE_REAL_MIN) )
524 {
525 hypre_error_w_msg(HYPRE_ERROR_CONV, "Subnormal alpha value in PCG");
526 if (i==1) i_prod=i_prod_0;
527 break;
528 }
529
530 gamma_old = gamma;
531
532 /* x = x + alpha*p */
533 (*(pcg_functions->Axpy))(alpha, p, x);
534
535 /* r = r - alpha*s */
536 if ( !recompute_true_residual )
537 {
538 (*(pcg_functions->Axpy))(-alpha, s, r);
539 }
540 else
541 {
542 if (print_level > 1 && my_id == 0)
543 {
544 hypre_printf("Recomputing the residual...\n");
545 }
546 (*(pcg_functions->CopyVector))(b, r);
547 (*(pcg_functions->Matvec))(matvec_data, -1.0, A, x, 1.0, r);
548 }
549
550 /* residual-based stopping criteria: ||r_new-r_old|| < rtol ||b|| */
551 if (rtol && two_norm)
552 {
553 /* use that r_new-r_old = alpha * s */
554 HYPRE_Real drob2 = alpha*alpha*(*(pcg_functions->InnerProd))(s,s)/bi_prod;
555 if ( drob2 < rtol*rtol )
556 {
557 if (print_level > 1 && my_id == 0)
558 {
559 hypre_printf("\n\n||r_old-r_new||/||b||: %e\n", sqrt(drob2));
560 }
561 break;
562 }
563 }
564
565 /* s = C*r */
566 (*(pcg_functions->ClearVector))(s);
567 precond(precond_data, A, r, s);
568
569 /* gamma = <r,s> */
570 gamma = (*(pcg_functions->InnerProd))(r, s);
571
572 /* residual-based stopping criteria: ||r_new-r_old||_C < rtol ||b||_C */
573 if (rtol && !two_norm)
574 {
575 /* use that ||r_new-r_old||_C^2 = (r_new ,C r_new) + (r_old, C r_old) */
576 HYPRE_Real r2ob2 = (gamma + gamma_old)/bi_prod;
577 if ( r2ob2 < rtol*rtol)
578 {
579 if (print_level > 1 && my_id == 0)
580 {
581 hypre_printf("\n\n||r_old-r_new||_C/||b||_C: %e\n", sqrt(r2ob2));
582 }
583 break;
584 }
585 }
586
587 /* set i_prod for convergence test */
588 if (two_norm)
589 i_prod = (*(pcg_functions->InnerProd))(r,r);
590 else
591 i_prod = gamma;
592
593 /*--------------------------------------------------------------------
594 * optional output
595 *--------------------------------------------------------------------*/
596 #if 0
597 if (two_norm)
598 hypre_printf("Iter (%d): ||r||_2 = %e, ||r||_2/||b||_2 = %e\n",
599 i, sqrt(i_prod), (bi_prod ? sqrt(i_prod/bi_prod) : 0));
600 else
601 hypre_printf("Iter (%d): ||r||_C = %e, ||r||_C/||b||_C = %e\n",
602 i, sqrt(i_prod), (bi_prod ? sqrt(i_prod/bi_prod) : 0));
603 #endif
604
605 /* print norm info */
606 if ( logging>0 || print_level>0 )
607 {
608 norms[i] = sqrt(i_prod);
609 rel_norms[i] = bi_prod ? sqrt(i_prod/bi_prod) : 0;
610 }
611 if ( print_level > 1 && my_id==0 )
612 {
613 if (two_norm)
614 {
615 if ( stop_crit && !rel_change && atolf==0 ) { /* pure absolute tolerance */
616 hypre_printf("% 5d %e %f\n", i, norms[i],
617 norms[i]/norms[i-1] );
618 }
619 else
620 {
621 hypre_printf("% 5d %e %f %e\n", i, norms[i],
622 norms[i]/norms[i-1], rel_norms[i] );
623 }
624 }
625 else
626 {
627 hypre_printf("% 5d %e %f %e\n", i, norms[i],
628 norms[i]/norms[i-1], rel_norms[i] );
629 }
630 }
631
632
633 /*--------------------------------------------------------------------
634 * check for convergence
635 *--------------------------------------------------------------------*/
636 if (i_prod / bi_prod < eps) /* the basic convergence test */
637 tentatively_converged = 1;
638 if ( tentatively_converged && recompute_residual )
639 /* At user request, don't trust the convergence test until we've recomputed
640 the residual from scratch. This is expensive in the usual case where an
641 the norm is the energy norm.
642 This calculation is coded on the assumption that r's accuracy is only a
643 concern for problems where CG takes many iterations. */
644 {
645 /* r = b - Ax */
646 (*(pcg_functions->CopyVector))(b, r);
647 (*(pcg_functions->Matvec))(matvec_data, -1.0, A, x, 1.0, r);
648
649 /* set i_prod for convergence test */
650 if (two_norm)
651 {
652 i_prod = (*(pcg_functions->InnerProd))(r,r);
653 }
654 else
655 {
656 /* s = C*r */
657 (*(pcg_functions->ClearVector))(s);
658 precond(precond_data, A, r, s);
659 /* iprod = gamma = <r,s> */
660 i_prod = (*(pcg_functions->InnerProd))(r, s);
661 }
662 if (i_prod / bi_prod >= eps) tentatively_converged = 0;
663 }
664 if ( tentatively_converged && rel_change && (i_prod > guard_zero_residual ))
665 /* At user request, don't treat this as converged unless x didn't change
666 much in the last iteration. */
667 {
668 pi_prod = (*(pcg_functions->InnerProd))(p,p);
669 xi_prod = (*(pcg_functions->InnerProd))(x,x);
670 ratio = alpha*alpha*pi_prod/xi_prod;
671 if (ratio >= eps) tentatively_converged = 0;
672 }
673 if ( tentatively_converged )
674 /* we've passed all the convergence tests, it's for real */
675 {
676 (pcg_data -> converged) = 1;
677 break;
678 }
679
680 if (! (gamma > HYPRE_REAL_MIN) )
681 {
682 hypre_error_w_msg(HYPRE_ERROR_CONV, "Subnormal gamma value in PCG");
683
684 break;
685 }
686 /* ... gamma should be >=0. IEEE subnormal numbers are < 2**(-1022)=2.2e-308
687 (and >= 2**(-1074)=4.9e-324). So a gamma this small means we're getting
688 dangerously close to subnormal or zero numbers (usually if gamma is small,
689 so will be other variables). Thus further calculations risk a crash.
690 Such small gamma generally means no hope of progress anyway. */
691
692 /*--------------------------------------------------------------------
693 * Optional test to see if adequate progress is being made.
694 * The average convergence factor is recorded and compared
695 * against the tolerance 'cf_tol'. The weighting factor is
696 * intended to pay more attention to the test when an accurate
697 * estimate for average convergence factor is available.
698 *--------------------------------------------------------------------*/
699
700 if (cf_tol > 0.0)
701 {
702 cf_ave_0 = cf_ave_1;
703 if (! (i_prod_0 > HYPRE_REAL_MIN) )
704 {
705 /* i_prod_0 is zero, or (almost) subnormal, yet i_prod wasn't small
706 enough to pass the convergence test. Therefore initial guess was good,
707 and we're just calculating garbage - time to bail out before the
708 next step, which will be a divide by zero (or close to it). */
709 hypre_error_w_msg(HYPRE_ERROR_CONV, "Subnormal i_prod value in PCG");
710
711 break;
712 }
713 cf_ave_1 = pow( i_prod / i_prod_0, 1.0/(2.0*i) );
714
715 weight = fabs(cf_ave_1 - cf_ave_0);
716 weight = weight / hypre_max(cf_ave_1, cf_ave_0);
717 weight = 1.0 - weight;
718 #if 0
719 hypre_printf("I = %d: cf_new = %e, cf_old = %e, weight = %e\n",
720 i, cf_ave_1, cf_ave_0, weight );
721 #endif
722 if (weight * cf_ave_1 > cf_tol) break;
723 }
724
725 /*--------------------------------------------------------------------
726 * back to the core CG calculations
727 *--------------------------------------------------------------------*/
728
729 /* beta = gamma / gamma_old */
730 beta = gamma / gamma_old;
731
732 /* p = s + beta p */
733 if ( !recompute_true_residual )
734 {
735 (*(pcg_functions->ScaleVector))(beta, p);
736 (*(pcg_functions->Axpy))(1.0, s, p);
737 }
738 else
739 (*(pcg_functions->CopyVector))(s, p);
740 }
741
742 /*--------------------------------------------------------------------
743 * Finish up with some outputs.
744 *--------------------------------------------------------------------*/
745
746 if ( print_level > 1 && my_id==0 )
747 hypre_printf("\n\n");
748
749 if (i >= max_iter && (i_prod/bi_prod) >= eps && eps > 0 && hybrid != -1)
750 {
751 hypre_error_w_msg(HYPRE_ERROR_CONV, "Reached max iterations in PCG before convergence");
752 }
753
754 (pcg_data -> num_iterations) = i;
755 if (bi_prod > 0.0)
756 (pcg_data -> rel_residual_norm) = sqrt(i_prod/bi_prod);
757 else /* actually, we'll never get here... */
758 (pcg_data -> rel_residual_norm) = 0.0;
759
760 HYPRE_ANNOTATE_FUNC_END;
761
762 return hypre_error_flag;
763 }
764
765 /*--------------------------------------------------------------------------
766 * hypre_PCGSetTol, hypre_PCGGetTol
767 *--------------------------------------------------------------------------*/
768
769 HYPRE_Int
hypre_PCGSetTol(void * pcg_vdata,HYPRE_Real tol)770 hypre_PCGSetTol( void *pcg_vdata,
771 HYPRE_Real tol )
772 {
773 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
774
775 (pcg_data -> tol) = tol;
776
777 return hypre_error_flag;
778 }
779
780 HYPRE_Int
hypre_PCGGetTol(void * pcg_vdata,HYPRE_Real * tol)781 hypre_PCGGetTol( void *pcg_vdata,
782 HYPRE_Real * tol )
783 {
784 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
785
786 *tol = (pcg_data -> tol);
787
788 return hypre_error_flag;
789 }
790 /*--------------------------------------------------------------------------
791 * hypre_PCGSetAbsoluteTol, hypre_PCGGetAbsoluteTol
792 *--------------------------------------------------------------------------*/
793
794 HYPRE_Int
hypre_PCGSetAbsoluteTol(void * pcg_vdata,HYPRE_Real a_tol)795 hypre_PCGSetAbsoluteTol( void *pcg_vdata,
796 HYPRE_Real a_tol )
797 {
798 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
799
800 (pcg_data -> a_tol) = a_tol;
801
802 return hypre_error_flag;
803 }
804
805 HYPRE_Int
hypre_PCGGetAbsoluteTol(void * pcg_vdata,HYPRE_Real * a_tol)806 hypre_PCGGetAbsoluteTol( void *pcg_vdata,
807 HYPRE_Real * a_tol )
808 {
809 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
810
811 *a_tol = (pcg_data -> a_tol);
812
813 return hypre_error_flag;
814 }
815
816 /*--------------------------------------------------------------------------
817 * hypre_PCGSetAbsoluteTolFactor, hypre_PCGGetAbsoluteTolFactor
818 *--------------------------------------------------------------------------*/
819
820 HYPRE_Int
hypre_PCGSetAbsoluteTolFactor(void * pcg_vdata,HYPRE_Real atolf)821 hypre_PCGSetAbsoluteTolFactor( void *pcg_vdata,
822 HYPRE_Real atolf )
823 {
824 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
825
826 (pcg_data -> atolf) = atolf;
827
828 return hypre_error_flag;
829 }
830
831 HYPRE_Int
hypre_PCGGetAbsoluteTolFactor(void * pcg_vdata,HYPRE_Real * atolf)832 hypre_PCGGetAbsoluteTolFactor( void *pcg_vdata,
833 HYPRE_Real * atolf )
834 {
835 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
836
837 *atolf = (pcg_data -> atolf);
838
839 return hypre_error_flag;
840 }
841
842 /*--------------------------------------------------------------------------
843 * hypre_PCGSetResidualTol, hypre_PCGGetResidualTol
844 *--------------------------------------------------------------------------*/
845
846 HYPRE_Int
hypre_PCGSetResidualTol(void * pcg_vdata,HYPRE_Real rtol)847 hypre_PCGSetResidualTol( void *pcg_vdata,
848 HYPRE_Real rtol )
849 {
850 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
851
852 (pcg_data -> rtol) = rtol;
853
854 return hypre_error_flag;
855 }
856
857 HYPRE_Int
hypre_PCGGetResidualTol(void * pcg_vdata,HYPRE_Real * rtol)858 hypre_PCGGetResidualTol( void *pcg_vdata,
859 HYPRE_Real * rtol )
860 {
861 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
862
863 *rtol = (pcg_data -> rtol);
864
865 return hypre_error_flag;
866 }
867
868 /*--------------------------------------------------------------------------
869 * hypre_PCGSetConvergenceFactorTol, hypre_PCGGetConvergenceFactorTol
870 *--------------------------------------------------------------------------*/
871
872 HYPRE_Int
hypre_PCGSetConvergenceFactorTol(void * pcg_vdata,HYPRE_Real cf_tol)873 hypre_PCGSetConvergenceFactorTol( void *pcg_vdata,
874 HYPRE_Real cf_tol )
875 {
876 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
877
878 (pcg_data -> cf_tol) = cf_tol;
879
880 return hypre_error_flag;
881 }
882
883 HYPRE_Int
hypre_PCGGetConvergenceFactorTol(void * pcg_vdata,HYPRE_Real * cf_tol)884 hypre_PCGGetConvergenceFactorTol( void *pcg_vdata,
885 HYPRE_Real * cf_tol )
886 {
887 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
888
889 *cf_tol = (pcg_data -> cf_tol);
890
891 return hypre_error_flag;
892 }
893
894 /*--------------------------------------------------------------------------
895 * hypre_PCGSetMaxIter, hypre_PCGGetMaxIter
896 *--------------------------------------------------------------------------*/
897
898 HYPRE_Int
hypre_PCGSetMaxIter(void * pcg_vdata,HYPRE_Int max_iter)899 hypre_PCGSetMaxIter( void *pcg_vdata,
900 HYPRE_Int max_iter )
901 {
902 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
903
904 (pcg_data -> max_iter) = max_iter;
905
906 return hypre_error_flag;
907 }
908
909 HYPRE_Int
hypre_PCGGetMaxIter(void * pcg_vdata,HYPRE_Int * max_iter)910 hypre_PCGGetMaxIter( void *pcg_vdata,
911 HYPRE_Int * max_iter )
912 {
913 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
914
915
916 *max_iter = (pcg_data -> max_iter);
917
918 return hypre_error_flag;
919 }
920
921 /*--------------------------------------------------------------------------
922 * hypre_PCGSetTwoNorm, hypre_PCGGetTwoNorm
923 *--------------------------------------------------------------------------*/
924
925 HYPRE_Int
hypre_PCGSetTwoNorm(void * pcg_vdata,HYPRE_Int two_norm)926 hypre_PCGSetTwoNorm( void *pcg_vdata,
927 HYPRE_Int two_norm )
928 {
929 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
930
931
932 (pcg_data -> two_norm) = two_norm;
933
934 return hypre_error_flag;
935 }
936
937 HYPRE_Int
hypre_PCGGetTwoNorm(void * pcg_vdata,HYPRE_Int * two_norm)938 hypre_PCGGetTwoNorm( void *pcg_vdata,
939 HYPRE_Int * two_norm )
940 {
941 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
942
943
944 *two_norm = (pcg_data -> two_norm);
945
946 return hypre_error_flag;
947 }
948
949 /*--------------------------------------------------------------------------
950 * hypre_PCGSetRelChange, hypre_PCGGetRelChange
951 *--------------------------------------------------------------------------*/
952
953 HYPRE_Int
hypre_PCGSetRelChange(void * pcg_vdata,HYPRE_Int rel_change)954 hypre_PCGSetRelChange( void *pcg_vdata,
955 HYPRE_Int rel_change )
956 {
957 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
958
959
960 (pcg_data -> rel_change) = rel_change;
961
962 return hypre_error_flag;
963 }
964
965 HYPRE_Int
hypre_PCGGetRelChange(void * pcg_vdata,HYPRE_Int * rel_change)966 hypre_PCGGetRelChange( void *pcg_vdata,
967 HYPRE_Int * rel_change )
968 {
969 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
970
971
972 *rel_change = (pcg_data -> rel_change);
973
974 return hypre_error_flag;
975 }
976
977 /*--------------------------------------------------------------------------
978 * hypre_PCGSetRecomputeResidual, hypre_PCGGetRecomputeResidual
979 *--------------------------------------------------------------------------*/
980
981 HYPRE_Int
hypre_PCGSetRecomputeResidual(void * pcg_vdata,HYPRE_Int recompute_residual)982 hypre_PCGSetRecomputeResidual( void *pcg_vdata,
983 HYPRE_Int recompute_residual )
984 {
985 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
986
987
988 (pcg_data -> recompute_residual) = recompute_residual;
989
990 return hypre_error_flag;
991 }
992
993 HYPRE_Int
hypre_PCGGetRecomputeResidual(void * pcg_vdata,HYPRE_Int * recompute_residual)994 hypre_PCGGetRecomputeResidual( void *pcg_vdata,
995 HYPRE_Int * recompute_residual )
996 {
997 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
998
999
1000 *recompute_residual = (pcg_data -> recompute_residual);
1001
1002 return hypre_error_flag;
1003 }
1004
1005 /*--------------------------------------------------------------------------
1006 * hypre_PCGSetRecomputeResidualP, hypre_PCGGetRecomputeResidualP
1007 *--------------------------------------------------------------------------*/
1008
1009 HYPRE_Int
hypre_PCGSetRecomputeResidualP(void * pcg_vdata,HYPRE_Int recompute_residual_p)1010 hypre_PCGSetRecomputeResidualP( void *pcg_vdata,
1011 HYPRE_Int recompute_residual_p )
1012 {
1013 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1014
1015 (pcg_data -> recompute_residual_p) = recompute_residual_p;
1016
1017 return hypre_error_flag;
1018 }
1019
1020 HYPRE_Int
hypre_PCGGetRecomputeResidualP(void * pcg_vdata,HYPRE_Int * recompute_residual_p)1021 hypre_PCGGetRecomputeResidualP( void *pcg_vdata,
1022 HYPRE_Int * recompute_residual_p )
1023 {
1024 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1025
1026 *recompute_residual_p = (pcg_data -> recompute_residual_p);
1027
1028 return hypre_error_flag;
1029 }
1030
1031 /*--------------------------------------------------------------------------
1032 * hypre_PCGSetStopCrit, hypre_PCGGetStopCrit
1033 *--------------------------------------------------------------------------*/
1034
1035 HYPRE_Int
hypre_PCGSetStopCrit(void * pcg_vdata,HYPRE_Int stop_crit)1036 hypre_PCGSetStopCrit( void *pcg_vdata,
1037 HYPRE_Int stop_crit )
1038 {
1039 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1040
1041
1042 (pcg_data -> stop_crit) = stop_crit;
1043
1044 return hypre_error_flag;
1045 }
1046
1047 HYPRE_Int
hypre_PCGGetStopCrit(void * pcg_vdata,HYPRE_Int * stop_crit)1048 hypre_PCGGetStopCrit( void *pcg_vdata,
1049 HYPRE_Int * stop_crit )
1050 {
1051 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1052
1053
1054 *stop_crit = (pcg_data -> stop_crit);
1055
1056 return hypre_error_flag;
1057 }
1058
1059 /*--------------------------------------------------------------------------
1060 * hypre_PCGGetPrecond
1061 *--------------------------------------------------------------------------*/
1062
1063 HYPRE_Int
hypre_PCGGetPrecond(void * pcg_vdata,HYPRE_Solver * precond_data_ptr)1064 hypre_PCGGetPrecond( void *pcg_vdata,
1065 HYPRE_Solver *precond_data_ptr )
1066 {
1067 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1068
1069
1070 *precond_data_ptr = (HYPRE_Solver)(pcg_data -> precond_data);
1071
1072 return hypre_error_flag;
1073 }
1074
1075 /*--------------------------------------------------------------------------
1076 * hypre_PCGSetPrecond
1077 *--------------------------------------------------------------------------*/
1078
1079 HYPRE_Int
hypre_PCGSetPrecond(void * pcg_vdata,HYPRE_Int (* precond)(void *,void *,void *,void *),HYPRE_Int (* precond_setup)(void *,void *,void *,void *),void * precond_data)1080 hypre_PCGSetPrecond( void *pcg_vdata,
1081 HYPRE_Int (*precond)(void*,void*,void*,void*),
1082 HYPRE_Int (*precond_setup)(void*,void*,void*,void*),
1083 void *precond_data )
1084 {
1085 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1086 hypre_PCGFunctions *pcg_functions = pcg_data->functions;
1087
1088
1089 (pcg_functions -> precond) = precond;
1090 (pcg_functions -> precond_setup) = precond_setup;
1091 (pcg_data -> precond_data) = precond_data;
1092
1093 return hypre_error_flag;
1094 }
1095
1096 /*--------------------------------------------------------------------------
1097 * hypre_PCGSetPrintLevel, hypre_PCGGetPrintLevel
1098 *--------------------------------------------------------------------------*/
1099
1100 HYPRE_Int
hypre_PCGSetPrintLevel(void * pcg_vdata,HYPRE_Int level)1101 hypre_PCGSetPrintLevel( void *pcg_vdata,
1102 HYPRE_Int level)
1103 {
1104 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1105
1106
1107 (pcg_data -> print_level) = level;
1108
1109 return hypre_error_flag;
1110 }
1111
1112 HYPRE_Int
hypre_PCGGetPrintLevel(void * pcg_vdata,HYPRE_Int * level)1113 hypre_PCGGetPrintLevel( void *pcg_vdata,
1114 HYPRE_Int * level)
1115 {
1116 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1117
1118
1119 *level = (pcg_data -> print_level);
1120
1121 return hypre_error_flag;
1122 }
1123
1124 /*--------------------------------------------------------------------------
1125 * hypre_PCGSetLogging, hypre_PCGGetLogging
1126 *--------------------------------------------------------------------------*/
1127
1128 HYPRE_Int
hypre_PCGSetLogging(void * pcg_vdata,HYPRE_Int level)1129 hypre_PCGSetLogging( void *pcg_vdata,
1130 HYPRE_Int level)
1131 {
1132 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1133
1134 (pcg_data -> logging) = level;
1135
1136 return hypre_error_flag;
1137 }
1138
1139 HYPRE_Int
hypre_PCGGetLogging(void * pcg_vdata,HYPRE_Int * level)1140 hypre_PCGGetLogging( void *pcg_vdata,
1141 HYPRE_Int * level)
1142 {
1143 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1144
1145 *level = (pcg_data -> logging);
1146
1147 return hypre_error_flag;
1148 }
1149
1150 HYPRE_Int
hypre_PCGSetHybrid(void * pcg_vdata,HYPRE_Int level)1151 hypre_PCGSetHybrid( void *pcg_vdata,
1152 HYPRE_Int level)
1153 {
1154 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1155
1156 (pcg_data -> hybrid) = level;
1157
1158 return hypre_error_flag;
1159 }
1160
1161 /*--------------------------------------------------------------------------
1162 * hypre_PCGGetNumIterations
1163 *--------------------------------------------------------------------------*/
1164
1165 HYPRE_Int
hypre_PCGGetNumIterations(void * pcg_vdata,HYPRE_Int * num_iterations)1166 hypre_PCGGetNumIterations( void *pcg_vdata,
1167 HYPRE_Int *num_iterations )
1168 {
1169 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1170
1171 *num_iterations = (pcg_data -> num_iterations);
1172
1173 return hypre_error_flag;
1174 }
1175
1176 /*--------------------------------------------------------------------------
1177 * hypre_PCGGetConverged
1178 *--------------------------------------------------------------------------*/
1179
1180 HYPRE_Int
hypre_PCGGetConverged(void * pcg_vdata,HYPRE_Int * converged)1181 hypre_PCGGetConverged( void *pcg_vdata,
1182 HYPRE_Int *converged)
1183 {
1184 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1185
1186 *converged = (pcg_data -> converged);
1187
1188 return hypre_error_flag;
1189 }
1190
1191 /*--------------------------------------------------------------------------
1192 * hypre_PCGPrintLogging
1193 *--------------------------------------------------------------------------*/
1194
1195 HYPRE_Int
hypre_PCGPrintLogging(void * pcg_vdata,HYPRE_Int myid)1196 hypre_PCGPrintLogging( void *pcg_vdata,
1197 HYPRE_Int myid)
1198 {
1199 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1200
1201 HYPRE_Int num_iterations = (pcg_data -> num_iterations);
1202 HYPRE_Int print_level = (pcg_data -> print_level);
1203 HYPRE_Real *norms = (pcg_data -> norms);
1204 HYPRE_Real *rel_norms = (pcg_data -> rel_norms);
1205
1206 HYPRE_Int i;
1207
1208 if (myid == 0)
1209 {
1210 if (print_level > 0)
1211 {
1212 for (i = 0; i < num_iterations; i++)
1213 {
1214 hypre_printf("Residual norm[%d] = %e ", i, norms[i]);
1215 hypre_printf("Relative residual norm[%d] = %e\n", i, rel_norms[i]);
1216 }
1217 }
1218 }
1219
1220 return hypre_error_flag;
1221 }
1222
1223 /*--------------------------------------------------------------------------
1224 * hypre_PCGGetFinalRelativeResidualNorm
1225 *--------------------------------------------------------------------------*/
1226
1227 HYPRE_Int
hypre_PCGGetFinalRelativeResidualNorm(void * pcg_vdata,HYPRE_Real * relative_residual_norm)1228 hypre_PCGGetFinalRelativeResidualNorm( void *pcg_vdata,
1229 HYPRE_Real *relative_residual_norm )
1230 {
1231 hypre_PCGData *pcg_data = (hypre_PCGData *)pcg_vdata;
1232
1233 HYPRE_Real rel_residual_norm = (pcg_data -> rel_residual_norm);
1234
1235 *relative_residual_norm = rel_residual_norm;
1236
1237 return hypre_error_flag;
1238 }
1239