source: CIVL/examples/mpi-omp/AMG2013/krylov/gmres.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 30.5 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14
15/******************************************************************************
16 *
17 * GMRES gmres
18 *
19 *****************************************************************************/
20
21#include "krylov.h"
22#include "utilities.h"
23
24/*--------------------------------------------------------------------------
25 * hypre_GMRESFunctionsCreate
26 *--------------------------------------------------------------------------*/
27
28hypre_GMRESFunctions *
29hypre_GMRESFunctionsCreate(
30 char * (*CAlloc) ( int count, int elt_size ),
31 int (*Free) ( char *ptr ),
32 int (*CommInfo) ( void *A, int *my_id, int *num_procs ),
33 void * (*CreateVector) ( void *vector ),
34 void * (*CreateVectorArray) ( int size, void *vectors ),
35 int (*DestroyVector) ( void *vector ),
36 void * (*MatvecCreate) ( void *A, void *x ),
37 int (*Matvec) ( void *matvec_data, double alpha, void *A,
38 void *x, double beta, void *y ),
39 int (*MatvecDestroy) ( void *matvec_data ),
40 double (*InnerProd) ( void *x, void *y ),
41 int (*CopyVector) ( void *x, void *y ),
42 int (*ClearVector) ( void *x ),
43 int (*ScaleVector) ( double alpha, void *x ),
44 int (*Axpy) ( double alpha, void *x, void *y ),
45 int (*PrecondSetup) ( void *vdata, void *A, void *b, void *x ),
46 int (*Precond) ( void *vdata, void *A, void *b, void *x )
47 )
48{
49 hypre_GMRESFunctions * gmres_functions;
50 gmres_functions = (hypre_GMRESFunctions *)
51 CAlloc( 1, sizeof(hypre_GMRESFunctions) );
52
53 gmres_functions->CAlloc = CAlloc;
54 gmres_functions->Free = Free;
55 gmres_functions->CommInfo = CommInfo; /* not in PCGFunctionsCreate */
56 gmres_functions->CreateVector = CreateVector;
57 gmres_functions->CreateVectorArray = CreateVectorArray; /* not in PCGFunctionsCreate */
58 gmres_functions->DestroyVector = DestroyVector;
59 gmres_functions->MatvecCreate = MatvecCreate;
60 gmres_functions->Matvec = Matvec;
61 gmres_functions->MatvecDestroy = MatvecDestroy;
62 gmres_functions->InnerProd = InnerProd;
63 gmres_functions->CopyVector = CopyVector;
64 gmres_functions->ClearVector = ClearVector;
65 gmres_functions->ScaleVector = ScaleVector;
66 gmres_functions->Axpy = Axpy;
67/* default preconditioner must be set here but can be changed later... */
68 gmres_functions->precond_setup = PrecondSetup;
69 gmres_functions->precond = Precond;
70
71 return gmres_functions;
72}
73
74/*--------------------------------------------------------------------------
75 * hypre_GMRESCreate
76 *--------------------------------------------------------------------------*/
77
78void *
79hypre_GMRESCreate( hypre_GMRESFunctions *gmres_functions )
80{
81 hypre_GMRESData *gmres_data;
82
83 gmres_data = hypre_CTAllocF(hypre_GMRESData, 1, gmres_functions);
84
85 gmres_data->functions = gmres_functions;
86
87 /* set defaults */
88 (gmres_data -> k_dim) = 5;
89 (gmres_data -> tol) = 1.0e-06;
90 (gmres_data -> cf_tol) = 0.0;
91 (gmres_data -> min_iter) = 0;
92 (gmres_data -> max_iter) = 1000;
93 (gmres_data -> rel_change) = 0;
94 (gmres_data -> stop_crit) = 0; /* rel. residual norm */
95 (gmres_data -> converged) = 0;
96 (gmres_data -> precond_data) = NULL;
97 (gmres_data -> print_level) = 0;
98 (gmres_data -> logging) = 0;
99 (gmres_data -> p) = NULL;
100 (gmres_data -> r) = NULL;
101 (gmres_data -> w) = NULL;
102 (gmres_data -> matvec_data) = NULL;
103 (gmres_data -> norms) = NULL;
104 (gmres_data -> log_file_name) = NULL;
105
106 return (void *) gmres_data;
107}
108
109/*--------------------------------------------------------------------------
110 * hypre_GMRESDestroy
111 *--------------------------------------------------------------------------*/
112
113int
114hypre_GMRESDestroy( void *gmres_vdata )
115{
116 hypre_GMRESData *gmres_data = gmres_vdata;
117 hypre_GMRESFunctions *gmres_functions = gmres_data->functions;
118 int i, ierr = 0;
119
120 if (gmres_data)
121 {
122 if ( (gmres_data->logging>0) || (gmres_data->print_level) > 0 )
123 {
124 if ( (gmres_data -> norms) != NULL )
125 hypre_TFreeF( gmres_data -> norms, gmres_functions );
126 }
127
128 if ( (gmres_data -> matvec_data) != NULL )
129 (*(gmres_functions->MatvecDestroy))(gmres_data -> matvec_data);
130
131 if ( (gmres_data -> r) != NULL )
132 (*(gmres_functions->DestroyVector))(gmres_data -> r);
133 if ( (gmres_data -> w) != NULL )
134 (*(gmres_functions->DestroyVector))(gmres_data -> w);
135 if ( (gmres_data -> p) != NULL )
136 {
137 for (i = 0; i < (gmres_data -> k_dim+1); i++)
138 {
139 if ( (gmres_data -> p)[i] != NULL )
140 (*(gmres_functions->DestroyVector))( (gmres_data -> p) [i]);
141 }
142 hypre_TFreeF( gmres_data->p, gmres_functions );
143 }
144 hypre_TFreeF( gmres_data, gmres_functions );
145 hypre_TFreeF( gmres_functions, gmres_functions );
146 }
147
148 return(ierr);
149}
150
151/*--------------------------------------------------------------------------
152 * hypre_GMRESGetResidual
153 *--------------------------------------------------------------------------*/
154
155int hypre_GMRESGetResidual( void *gmres_vdata, void **residual )
156{
157 /* returns a pointer to the residual vector */
158 int ierr = 0;
159 hypre_GMRESData *gmres_data = gmres_vdata;
160 *residual = gmres_data->r;
161 return ierr;
162}
163
164/*--------------------------------------------------------------------------
165 * hypre_GMRESSetup
166 *--------------------------------------------------------------------------*/
167
168int
169hypre_GMRESSetup( void *gmres_vdata,
170 void *A,
171 void *b,
172 void *x )
173{
174 hypre_GMRESData *gmres_data = gmres_vdata;
175 hypre_GMRESFunctions *gmres_functions = gmres_data->functions;
176
177 int k_dim = (gmres_data -> k_dim);
178 int max_iter = (gmres_data -> max_iter);
179 int (*precond_setup)(void *vdata, void *A, void *b, void *x) = (gmres_functions->precond_setup);
180 void *precond_data = (gmres_data -> precond_data);
181 int ierr = 0;
182
183 (gmres_data -> A) = A;
184
185 /*--------------------------------------------------
186 * The arguments for NewVector are important to
187 * maintain consistency between the setup and
188 * compute phases of matvec and the preconditioner.
189 *--------------------------------------------------*/
190
191 if ((gmres_data -> p) == NULL)
192 (gmres_data -> p) = (*(gmres_functions->CreateVectorArray))(k_dim+1,x);
193 if ((gmres_data -> r) == NULL)
194 (gmres_data -> r) = (*(gmres_functions->CreateVector))(b);
195 if ((gmres_data -> w) == NULL)
196 (gmres_data -> w) = (*(gmres_functions->CreateVector))(b);
197
198 if ((gmres_data -> matvec_data) == NULL)
199 (gmres_data -> matvec_data) = (*(gmres_functions->MatvecCreate))(A, x);
200
201 ierr = precond_setup(precond_data, A, b, x);
202
203 /*-----------------------------------------------------
204 * Allocate space for log info
205 *-----------------------------------------------------*/
206
207 if ( (gmres_data->logging)>0 || (gmres_data->print_level) > 0 )
208 {
209 if ((gmres_data -> norms) == NULL)
210 (gmres_data -> norms) = hypre_CTAllocF(double, max_iter + 1,gmres_functions);
211 }
212 if ( (gmres_data->print_level) > 0 ) {
213 if ((gmres_data -> log_file_name) == NULL)
214 (gmres_data -> log_file_name) = "gmres.out.log";
215 }
216
217 return ierr;
218}
219
220/*--------------------------------------------------------------------------
221 * hypre_GMRESSolve
222 *-------------------------------------------------------------------------*/
223
224int
225hypre_GMRESSolve(void *gmres_vdata,
226 void *A,
227 void *b,
228 void *x)
229{
230 hypre_GMRESData *gmres_data = gmres_vdata;
231 hypre_GMRESFunctions *gmres_functions = gmres_data->functions;
232 int k_dim = (gmres_data -> k_dim);
233 int min_iter = (gmres_data -> min_iter);
234 int max_iter = (gmres_data -> max_iter);
235 int rel_change = (gmres_data -> rel_change);
236 int stop_crit = (gmres_data -> stop_crit);
237 double accuracy = (gmres_data -> tol);
238 double cf_tol = (gmres_data -> cf_tol);
239 void *matvec_data = (gmres_data -> matvec_data);
240
241 void *r = (gmres_data -> r);
242 void *w = (gmres_data -> w);
243 void **p = (gmres_data -> p);
244
245 int (*precond)(void*, void*, void*, void*) = (gmres_functions -> precond);
246 int *precond_data = (gmres_data -> precond_data);
247
248 int print_level = (gmres_data -> print_level);
249 int logging = (gmres_data -> logging);
250
251 double *norms = (gmres_data -> norms);
252/* not used yet char *log_file_name = (gmres_data -> log_file_name);*/
253/* FILE *fp; */
254
255 int ierr = 0;
256 int break_value = 0;
257 int i, j, k;
258 double *rs, **hh, *c, *s;
259 int iter;
260 int my_id, num_procs;
261 double epsilon, gamma, t, r_norm, b_norm, den_norm, x_norm;
262 double epsmac = 1.e-16;
263 double ieee_check = 0.;
264
265 double guard_zero_residual;
266 double cf_ave_0 = 0.0;
267 double cf_ave_1 = 0.0;
268 double weight;
269 double r_norm_0;
270 double relative_error;
271
272 (gmres_data -> converged) = 0;
273 /*-----------------------------------------------------------------------
274 * With relative change convergence test on, it is possible to attempt
275 * another iteration with a zero residual. This causes the parameter
276 * alpha to go NaN. The guard_zero_residual parameter is to circumvent
277 * this. Perhaps it should be set to something non-zero (but small).
278 *-----------------------------------------------------------------------*/
279 guard_zero_residual = 0.0;
280
281 (*(gmres_functions->CommInfo))(A,&my_id,&num_procs);
282 if ( logging>0 || print_level>0 )
283 {
284 norms = (gmres_data -> norms);
285 /* not used yet log_file_name = (gmres_data -> log_file_name);*/
286 /* fp = fopen(log_file_name,"w"); */
287 }
288
289 /* initialize work arrays */
290 rs = hypre_CTAllocF(double,k_dim+1,gmres_functions);
291 c = hypre_CTAllocF(double,k_dim,gmres_functions);
292 s = hypre_CTAllocF(double,k_dim,gmres_functions);
293
294 hh = hypre_CTAllocF(double*,k_dim+1,gmres_functions);
295 for (i=0; i < k_dim+1; i++)
296 {
297 hh[i] = hypre_CTAllocF(double,k_dim,gmres_functions);
298 }
299
300 (*(gmres_functions->CopyVector))(b,p[0]);
301
302 /* compute initial residual */
303 (*(gmres_functions->Matvec))(matvec_data,-1.0, A, x, 1.0, p[0]);
304
305 b_norm = sqrt((*(gmres_functions->InnerProd))(b,b));
306
307 /* Since it is does not diminish performance, attempt to return an error flag
308 and notify users when they supply bad input. */
309 if (b_norm != 0.) ieee_check = b_norm/b_norm; /* INF -> NaN conversion */
310 if (ieee_check != ieee_check)
311 {
312 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
313 for ieee_check self-equality works on all IEEE-compliant compilers/
314 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
315 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
316 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
317 if (logging > 0 || print_level > 0)
318 {
319 printf("\n\nERROR detected by Hypre ... BEGIN\n");
320 printf("ERROR -- hypre_GMRESSolve: INFs and/or NaNs detected in input.\n");
321 printf("User probably placed non-numerics in supplied b.\n");
322 printf("Returning error flag += 101. Program not terminated.\n");
323 printf("ERROR detected by Hypre ... END\n\n\n");
324 }
325 ierr += 101;
326 return ierr;
327 }
328
329 r_norm = sqrt((*(gmres_functions->InnerProd))(p[0],p[0]));
330 r_norm_0 = r_norm;
331
332 /* Since it is does not diminish performance, attempt to return an error flag
333 and notify users when they supply bad input. */
334 if (r_norm != 0.) ieee_check = r_norm/r_norm; /* INF -> NaN conversion */
335 if (ieee_check != ieee_check)
336 {
337 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
338 for ieee_check self-equality works on all IEEE-compliant compilers/
339 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
340 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
341 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
342 if (logging > 0 || print_level > 0)
343 {
344 printf("\n\nERROR detected by Hypre ... BEGIN\n");
345 printf("ERROR -- hypre_GMRESSolve: INFs and/or NaNs detected in input.\n");
346 printf("User probably placed non-numerics in supplied A or x_0.\n");
347 printf("Returning error flag += 101. Program not terminated.\n");
348 printf("ERROR detected by Hypre ... END\n\n\n");
349 }
350 ierr += 101;
351 return ierr;
352 }
353
354 if ( logging>0 || print_level > 0)
355 {
356 norms[0] = r_norm;
357 if ( print_level>1 && my_id == 0 )
358 {
359 printf("L2 norm of b: %e\n", b_norm);
360 if (b_norm == 0.0)
361 printf("Rel_resid_norm actually contains the residual norm\n");
362 printf("Initial L2 norm of residual: %e\n", r_norm);
363
364 }
365 }
366 iter = 0;
367
368 if (b_norm > 0.0)
369 {
370/* convergence criterion |r_i|/|b| <= accuracy if |b| > 0 */
371 den_norm= b_norm;
372 }
373 else
374 {
375/* convergence criterion |r_i|/|r0| <= accuracy if |b| = 0 */
376 den_norm= r_norm;
377 };
378
379 epsilon= accuracy;
380
381/* convergence criterion |r_i| <= accuracy , absolute residual norm*/
382 if ( stop_crit && !rel_change )
383 epsilon = accuracy;
384
385 if ( print_level>1 && my_id == 0 )
386 {
387 if (b_norm > 0.0)
388 {printf("=============================================\n\n");
389 printf("Iters resid.norm conv.rate rel.res.norm\n");
390 printf("----- ------------ ---------- ------------\n");
391
392 }
393
394 else
395 {printf("=============================================\n\n");
396 printf("Iters resid.norm conv.rate\n");
397 printf("----- ------------ ----------\n");
398
399 };
400 }
401
402 /* set the relative_error to initially bypass the stopping criterion */
403 if (rel_change)
404 {
405 relative_error= epsilon + 1.;
406 }
407
408 while (iter < max_iter)
409 {
410 /* initialize first term of hessenberg system */
411
412 rs[0] = r_norm;
413 if (r_norm == 0.0)
414 {
415 hypre_TFreeF(c,gmres_functions);
416 hypre_TFreeF(s,gmres_functions);
417 hypre_TFreeF(rs,gmres_functions);
418 for (i=0; i < k_dim+1; i++) hypre_TFreeF(hh[i],gmres_functions);
419 hypre_TFreeF(hh,gmres_functions);
420 ierr = 0;
421 return ierr;
422 }
423
424 if (r_norm/den_norm <= epsilon && iter >= min_iter)
425 {
426 if (rel_change)
427 {
428 if (relative_error <= epsilon)
429 {
430 (*(gmres_functions->CopyVector))(b,r);
431 (*(gmres_functions->Matvec))(matvec_data,-1.0,A,x,1.0,r);
432 r_norm = sqrt((*(gmres_functions->InnerProd))(r,r));
433 if (r_norm/den_norm <= epsilon)
434 {
435 if ( print_level>1 && my_id == 0)
436 {
437 printf("\n\n");
438 printf("Final L2 norm of residual: %e\n\n", r_norm);
439 }
440 break;
441 }
442 else
443 if ( print_level>0 && my_id == 0)
444 printf("false convergence 1\n");
445 }
446 }
447 else
448 {
449 (*(gmres_functions->CopyVector))(b,r);
450 (*(gmres_functions->Matvec))(matvec_data,-1.0,A,x,1.0,r);
451 r_norm = sqrt((*(gmres_functions->InnerProd))(r,r));
452 if (r_norm/den_norm <= epsilon)
453 {
454 if ( print_level>1 && my_id == 0)
455 {
456 printf("\n\n");
457 printf("Final L2 norm of residual: %e\n\n", r_norm);
458 }
459 break;
460 }
461 else
462 if ( print_level>0 && my_id == 0)
463 printf("false convergence 1\n");
464 }
465
466 }
467
468 t = 1.0 / r_norm;
469 (*(gmres_functions->ScaleVector))(t,p[0]);
470 i = 0;
471 while (i < k_dim && ( (r_norm/den_norm > epsilon || iter < min_iter)
472 || ((rel_change) && relative_error > epsilon) )
473 && iter < max_iter)
474 {
475 i++;
476 iter++;
477 (*(gmres_functions->ClearVector))(r);
478 precond(precond_data, A, p[i-1], r);
479 (*(gmres_functions->Matvec))(matvec_data, 1.0, A, r, 0.0, p[i]);
480 /* modified Gram_Schmidt */
481 for (j=0; j < i; j++)
482 {
483 hh[j][i-1] = (*(gmres_functions->InnerProd))(p[j],p[i]);
484 (*(gmres_functions->Axpy))(-hh[j][i-1],p[j],p[i]);
485 }
486 t = sqrt((*(gmres_functions->InnerProd))(p[i],p[i]));
487 hh[i][i-1] = t;
488 if (t != 0.0)
489 {
490 t = 1.0/t;
491 (*(gmres_functions->ScaleVector))(t,p[i]);
492 }
493 /* done with modified Gram_schmidt and Arnoldi step.
494 update factorization of hh */
495 for (j = 1; j < i; j++)
496 {
497 t = hh[j-1][i-1];
498 hh[j-1][i-1] = c[j-1]*t + s[j-1]*hh[j][i-1];
499 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
500 }
501 gamma = sqrt(hh[i-1][i-1]*hh[i-1][i-1] + hh[i][i-1]*hh[i][i-1]);
502 if (gamma == 0.0) gamma = epsmac;
503 c[i-1] = hh[i-1][i-1]/gamma;
504 s[i-1] = hh[i][i-1]/gamma;
505 rs[i] = -s[i-1]*rs[i-1];
506 rs[i-1] = c[i-1]*rs[i-1];
507 /* determine residual norm */
508 hh[i-1][i-1] = c[i-1]*hh[i-1][i-1] + s[i-1]*hh[i][i-1];
509 r_norm = fabs(rs[i]);
510 if ( print_level>0 )
511 {
512 norms[iter] = r_norm;
513 if ( print_level>1 && my_id == 0 )
514 {
515 if (b_norm > 0.0)
516 printf("% 5d %e %f %e\n", iter,
517 norms[iter],norms[iter]/norms[iter-1],
518 norms[iter]/b_norm);
519 else
520 printf("% 5d %e %f\n", iter, norms[iter],
521 norms[iter]/norms[iter-1]);
522 }
523 }
524 if (cf_tol > 0.0)
525 {
526 cf_ave_0 = cf_ave_1;
527 cf_ave_1 = pow( r_norm / r_norm_0, 1.0/(2.0*iter));
528
529 weight = fabs(cf_ave_1 - cf_ave_0);
530 weight = weight / hypre_max(cf_ave_1, cf_ave_0);
531 weight = 1.0 - weight;
532#if 0
533 printf("I = %d: cf_new = %e, cf_old = %e, weight = %e\n",
534 i, cf_ave_1, cf_ave_0, weight );
535#endif
536 if (weight * cf_ave_1 > cf_tol)
537 {
538 break_value = 1;
539 break;
540 }
541 }
542
543 }
544 /* now compute solution, first solve upper triangular system */
545
546 if (break_value) break;
547
548 rs[i-1] = rs[i-1]/hh[i-1][i-1];
549 for (k = i-2; k >= 0; k--)
550 {
551 t = rs[k];
552 for (j = k+1; j < i; j++)
553 {
554 t -= hh[k][j]*rs[j];
555 }
556 rs[k] = t/hh[k][k];
557 }
558
559 (*(gmres_functions->CopyVector))(p[0],w);
560 (*(gmres_functions->ScaleVector))(rs[0],w);
561 for (j = 1; j < i; j++)
562 (*(gmres_functions->Axpy))(rs[j], p[j], w);
563
564 (*(gmres_functions->ClearVector))(r);
565 precond(precond_data, A, w, r);
566
567 (*(gmres_functions->Axpy))(1.0,r,x);
568
569/* check for convergence, evaluate actual residual */
570 if (r_norm/den_norm <= epsilon && iter >= min_iter)
571 {
572 if (rel_change)
573 {
574 x_norm = sqrt( (*(gmres_functions->InnerProd))(x,x) );
575 if ( x_norm<=guard_zero_residual ) break; /* don't divide by 0 */
576 r_norm = sqrt( (*(gmres_functions->InnerProd))(r,r) );
577 relative_error= r_norm/x_norm;
578 }
579
580 (*(gmres_functions->CopyVector))(b,r);
581 (*(gmres_functions->Matvec))(matvec_data,-1.0,A,x,1.0,r);
582 r_norm = sqrt( (*(gmres_functions->InnerProd))(r,r) );
583 if (r_norm/den_norm <= epsilon)
584 {
585 if ( print_level>1 && my_id == 0 )
586 {
587 printf("\n\n");
588 printf("Final L2 norm of residual: %e\n\n", r_norm);
589 }
590 if (rel_change && r_norm > guard_zero_residual)
591 /* Also test on relative change of iterates, x_i - x_(i-1) */
592 { /* At this point r = x_i - x_(i-1) */
593 x_norm = sqrt( (*(gmres_functions->InnerProd))(x,x) );
594 if ( x_norm<=guard_zero_residual ) break; /* don't divide by 0 */
595 if ( relative_error < epsilon )
596 {
597 (gmres_data -> converged) = 1;
598 break;
599 }
600 }
601 else
602 {
603 (gmres_data -> converged) = 1;
604 break;
605 }
606 }
607 else
608 {
609 if ( print_level>0 && my_id == 0)
610 printf("false convergence 2\n");
611 (*(gmres_functions->CopyVector))(r,p[0]);
612 i = 0;
613 }
614 }
615
616/* compute residual vector and continue loop */
617
618 for (j=i ; j > 0; j--)
619 {
620 rs[j-1] = -s[j-1]*rs[j];
621 rs[j] = c[j-1]*rs[j];
622 }
623
624 if (i) (*(gmres_functions->Axpy))(rs[0]-1.0,p[0],p[0]);
625 for (j=1; j < i+1; j++)
626 (*(gmres_functions->Axpy))(rs[j],p[j],p[0]);
627 }
628
629 if ( print_level>1 && my_id == 0 )
630 printf("\n\n");
631
632 (gmres_data -> num_iterations) = iter;
633 if (b_norm > 0.0)
634 (gmres_data -> rel_residual_norm) = r_norm/b_norm;
635 if (b_norm == 0.0)
636 (gmres_data -> rel_residual_norm) = r_norm;
637
638 if (iter >= max_iter && r_norm/den_norm > epsilon) ierr = 1;
639
640 hypre_TFreeF(c,gmres_functions);
641 hypre_TFreeF(s,gmres_functions);
642 hypre_TFreeF(rs,gmres_functions);
643
644 for (i=0; i < k_dim+1; i++)
645 {
646 hypre_TFreeF(hh[i],gmres_functions);
647 }
648 hypre_TFreeF(hh,gmres_functions);
649
650 return ierr;
651}
652
653/*--------------------------------------------------------------------------
654 * hypre_GMRESSetKDim, hypre_GMRESGetKDim
655 *--------------------------------------------------------------------------*/
656
657int
658hypre_GMRESSetKDim( void *gmres_vdata,
659 int k_dim )
660{
661 hypre_GMRESData *gmres_data = gmres_vdata;
662 int ierr = 0;
663
664 (gmres_data -> k_dim) = k_dim;
665
666 return ierr;
667}
668
669int
670hypre_GMRESGetKDim( void *gmres_vdata,
671 int * k_dim )
672{
673 hypre_GMRESData *gmres_data = gmres_vdata;
674 int ierr = 0;
675
676 *k_dim = (gmres_data -> k_dim);
677
678 return ierr;
679}
680
681/*--------------------------------------------------------------------------
682 * hypre_GMRESSetTol, hypre_GMRESGetTol
683 *--------------------------------------------------------------------------*/
684
685int
686hypre_GMRESSetTol( void *gmres_vdata,
687 double tol )
688{
689 hypre_GMRESData *gmres_data = gmres_vdata;
690 int ierr = 0;
691
692 (gmres_data -> tol) = tol;
693
694 return ierr;
695}
696
697int
698hypre_GMRESGetTol( void *gmres_vdata,
699 double * tol )
700{
701 hypre_GMRESData *gmres_data = gmres_vdata;
702 int ierr = 0;
703
704 *tol = (gmres_data -> tol);
705
706 return ierr;
707}
708
709/*--------------------------------------------------------------------------
710 * hypre_GMRESSetConvergenceFactorTol, hypre_GMRESGetConvergenceFactorTol
711 *--------------------------------------------------------------------------*/
712
713int
714hypre_GMRESSetConvergenceFactorTol( void *gmres_vdata,
715 double cf_tol )
716{
717 hypre_GMRESData *gmres_data = gmres_vdata;
718 int ierr = 0;
719
720 (gmres_data -> cf_tol) = cf_tol;
721
722 return ierr;
723}
724
725int
726hypre_GMRESGetConvergenceFactorTol( void *gmres_vdata,
727 double * cf_tol )
728{
729 hypre_GMRESData *gmres_data = gmres_vdata;
730 int ierr = 0;
731
732 *cf_tol = (gmres_data -> cf_tol);
733
734 return ierr;
735}
736
737/*--------------------------------------------------------------------------
738 * hypre_GMRESSetMinIter, hypre_GMRESGetMinIter
739 *--------------------------------------------------------------------------*/
740
741int
742hypre_GMRESSetMinIter( void *gmres_vdata,
743 int min_iter )
744{
745 hypre_GMRESData *gmres_data = gmres_vdata;
746 int ierr = 0;
747
748 (gmres_data -> min_iter) = min_iter;
749
750 return ierr;
751}
752
753int
754hypre_GMRESGetMinIter( void *gmres_vdata,
755 int * min_iter )
756{
757 hypre_GMRESData *gmres_data = gmres_vdata;
758 int ierr = 0;
759
760 *min_iter = (gmres_data -> min_iter);
761
762 return ierr;
763}
764
765/*--------------------------------------------------------------------------
766 * hypre_GMRESSetMaxIter, hypre_GMRESGetMaxIter
767 *--------------------------------------------------------------------------*/
768
769int
770hypre_GMRESSetMaxIter( void *gmres_vdata,
771 int max_iter )
772{
773 hypre_GMRESData *gmres_data = gmres_vdata;
774 int ierr = 0;
775
776 (gmres_data -> max_iter) = max_iter;
777
778 return ierr;
779}
780
781int
782hypre_GMRESGetMaxIter( void *gmres_vdata,
783 int * max_iter )
784{
785 hypre_GMRESData *gmres_data = gmres_vdata;
786 int ierr = 0;
787
788 *max_iter = (gmres_data -> max_iter);
789
790 return ierr;
791}
792
793/*--------------------------------------------------------------------------
794 * hypre_GMRESSetRelChange, hypre_GMRESGetRelChange
795 *--------------------------------------------------------------------------*/
796
797int
798hypre_GMRESSetRelChange( void *gmres_vdata,
799 int rel_change )
800{
801 hypre_GMRESData *gmres_data = gmres_vdata;
802 int ierr = 0;
803
804 (gmres_data -> rel_change) = rel_change;
805
806 return ierr;
807}
808
809int
810hypre_GMRESGetRelChange( void *gmres_vdata,
811 int * rel_change )
812{
813 hypre_GMRESData *gmres_data = gmres_vdata;
814 int ierr = 0;
815
816 *rel_change = (gmres_data -> rel_change);
817
818 return ierr;
819}
820
821/*--------------------------------------------------------------------------
822 * hypre_GMRESSetStopCrit, hypre_GMRESGetStopCrit
823 *--------------------------------------------------------------------------*/
824
825int
826hypre_GMRESSetStopCrit( void *gmres_vdata,
827 int stop_crit )
828{
829 hypre_GMRESData *gmres_data = gmres_vdata;
830 int ierr = 0;
831
832 (gmres_data -> stop_crit) = stop_crit;
833
834 return ierr;
835}
836
837int
838hypre_GMRESGetStopCrit( void *gmres_vdata,
839 int * stop_crit )
840{
841 hypre_GMRESData *gmres_data = gmres_vdata;
842 int ierr = 0;
843
844 *stop_crit = (gmres_data -> stop_crit);
845
846 return ierr;
847}
848
849/*--------------------------------------------------------------------------
850 * hypre_GMRESSetPrecond
851 *--------------------------------------------------------------------------*/
852
853int
854hypre_GMRESSetPrecond( void *gmres_vdata,
855 int (*precond)(void *vdata, void *A, void *b, void *x ),
856 int (*precond_setup)(void *vdata, void *A, void *b, void *x ),
857 void *precond_data )
858{
859 hypre_GMRESData *gmres_data = gmres_vdata;
860 hypre_GMRESFunctions *gmres_functions = gmres_data->functions;
861 int ierr = 0;
862
863 (gmres_functions -> precond) = precond;
864 (gmres_functions -> precond_setup) = precond_setup;
865 (gmres_data -> precond_data) = precond_data;
866
867 return ierr;
868}
869
870/*--------------------------------------------------------------------------
871 * hypre_GMRESGetPrecond
872 *--------------------------------------------------------------------------*/
873
874int
875hypre_GMRESGetPrecond( void *gmres_vdata,
876 HYPRE_Solver *precond_data_ptr )
877{
878 hypre_GMRESData *gmres_data = gmres_vdata;
879 int ierr = 0;
880
881 *precond_data_ptr = (HYPRE_Solver)(gmres_data -> precond_data);
882
883 return ierr;
884}
885
886/*--------------------------------------------------------------------------
887 * hypre_GMRESSetPrintLevel, hypre_GMRESGetPrintLevel
888 *--------------------------------------------------------------------------*/
889
890int
891hypre_GMRESSetPrintLevel( void *gmres_vdata,
892 int level)
893{
894 hypre_GMRESData *gmres_data = gmres_vdata;
895 int ierr = 0;
896
897 (gmres_data -> print_level) = level;
898
899 return ierr;
900}
901
902int
903hypre_GMRESGetPrintLevel( void *gmres_vdata,
904 int * level)
905{
906 hypre_GMRESData *gmres_data = gmres_vdata;
907 int ierr = 0;
908
909 *level = (gmres_data -> print_level);
910
911 return ierr;
912}
913
914/*--------------------------------------------------------------------------
915 * hypre_GMRESSetLogging, hypre_GMRESGetLogging
916 *--------------------------------------------------------------------------*/
917
918int
919hypre_GMRESSetLogging( void *gmres_vdata,
920 int level)
921{
922 hypre_GMRESData *gmres_data = gmres_vdata;
923 int ierr = 0;
924
925 (gmres_data -> logging) = level;
926
927 return ierr;
928}
929
930int
931hypre_GMRESGetLogging( void *gmres_vdata,
932 int * level)
933{
934 hypre_GMRESData *gmres_data = gmres_vdata;
935 int ierr = 0;
936
937 *level = (gmres_data -> logging);
938
939 return ierr;
940}
941
942/*--------------------------------------------------------------------------
943 * hypre_GMRESGetNumIterations
944 *--------------------------------------------------------------------------*/
945
946int
947hypre_GMRESGetNumIterations( void *gmres_vdata,
948 int *num_iterations )
949{
950 hypre_GMRESData *gmres_data = gmres_vdata;
951 int ierr = 0;
952
953 *num_iterations = (gmres_data -> num_iterations);
954
955 return ierr;
956}
957
958/*--------------------------------------------------------------------------
959 * hypre_GMRESGetConverged
960 *--------------------------------------------------------------------------*/
961
962int
963hypre_GMRESGetConverged( void *gmres_vdata,
964 int *converged )
965{
966 hypre_GMRESData *gmres_data = gmres_vdata;
967 int ierr = 0;
968
969 *converged = (gmres_data -> converged);
970
971 return ierr;
972}
973
974/*--------------------------------------------------------------------------
975 * hypre_GMRESGetFinalRelativeResidualNorm
976 *--------------------------------------------------------------------------*/
977
978int
979hypre_GMRESGetFinalRelativeResidualNorm( void *gmres_vdata,
980 double *relative_residual_norm )
981{
982 hypre_GMRESData *gmres_data = gmres_vdata;
983 int ierr = 0;
984
985 *relative_residual_norm = (gmres_data -> rel_residual_norm);
986
987 return ierr;
988}
Note: See TracBrowser for help on using the repository browser.