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  * COGMRES cogmres
11  *
12  *****************************************************************************/
13 
14 #include "krylov.h"
15 #include "_hypre_utilities.h"
16 
17 /*--------------------------------------------------------------------------
18  * hypre_COGMRESFunctionsCreate
19  *--------------------------------------------------------------------------*/
20 
21 hypre_COGMRESFunctions *
hypre_COGMRESFunctionsCreate(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),void * (* CreateVectorArray)(HYPRE_Int size,void * vectors),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 (* MassInnerProd)(void * x,void ** y,HYPRE_Int k,HYPRE_Int unroll,void * result),HYPRE_Int (* MassDotpTwo)(void * x,void * y,void ** z,HYPRE_Int k,HYPRE_Int unroll,void * result_x,void * result_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 (* MassAxpy)(HYPRE_Complex * alpha,void ** x,void * y,HYPRE_Int k,HYPRE_Int unroll),HYPRE_Int (* PrecondSetup)(void * vdata,void * A,void * b,void * x),HYPRE_Int (* Precond)(void * vdata,void * A,void * b,void * x))22 hypre_COGMRESFunctionsCreate(
23    void *       (*CAlloc)        ( size_t count, size_t elt_size, HYPRE_MemoryLocation location ),
24    HYPRE_Int    (*Free)          ( void *ptr ),
25    HYPRE_Int    (*CommInfo)      ( void  *A, HYPRE_Int   *my_id,
26                                    HYPRE_Int   *num_procs ),
27    void *       (*CreateVector)  ( void *vector ),
28    void *       (*CreateVectorArray)  ( HYPRE_Int size, void *vectors ),
29    HYPRE_Int    (*DestroyVector) ( void *vector ),
30    void *       (*MatvecCreate)  ( void *A, void *x ),
31    HYPRE_Int    (*Matvec)        ( void *matvec_data, HYPRE_Complex alpha, void *A,
32                                    void *x, HYPRE_Complex beta, void *y ),
33    HYPRE_Int    (*MatvecDestroy) ( void *matvec_data ),
34    HYPRE_Real   (*InnerProd)     ( void *x, void *y ),
35    HYPRE_Int    (*MassInnerProd) (void *x, void **y, HYPRE_Int k, HYPRE_Int unroll, void *result),
36    HYPRE_Int    (*MassDotpTwo)   (void *x, void *y, void **z, HYPRE_Int k, HYPRE_Int unroll, void *result_x, void *result_y),
37    HYPRE_Int    (*CopyVector)    ( void *x, void *y ),
38    HYPRE_Int    (*ClearVector)   ( void *x ),
39    HYPRE_Int    (*ScaleVector)   ( HYPRE_Complex alpha, void *x ),
40    HYPRE_Int    (*Axpy)          ( HYPRE_Complex alpha, void *x, void *y ),
41    HYPRE_Int    (*MassAxpy)      ( HYPRE_Complex *alpha, void **x, void *y, HYPRE_Int k, HYPRE_Int unroll),
42    HYPRE_Int    (*PrecondSetup)  ( void *vdata, void *A, void *b, void *x ),
43    HYPRE_Int    (*Precond)       ( void *vdata, void *A, void *b, void *x )
44    )
45 {
46    hypre_COGMRESFunctions * cogmres_functions;
47    cogmres_functions = (hypre_COGMRESFunctions *)
48     CAlloc( 1, sizeof(hypre_COGMRESFunctions), HYPRE_MEMORY_HOST );
49 
50    cogmres_functions->CAlloc            = CAlloc;
51    cogmres_functions->Free              = Free;
52    cogmres_functions->CommInfo          = CommInfo; /* not in PCGFunctionsCreate */
53    cogmres_functions->CreateVector      = CreateVector;
54    cogmres_functions->CreateVectorArray = CreateVectorArray; /* not in PCGFunctionsCreate */
55    cogmres_functions->DestroyVector     = DestroyVector;
56    cogmres_functions->MatvecCreate      = MatvecCreate;
57    cogmres_functions->Matvec            = Matvec;
58    cogmres_functions->MatvecDestroy     = MatvecDestroy;
59    cogmres_functions->InnerProd         = InnerProd;
60    cogmres_functions->MassInnerProd     = MassInnerProd;
61    cogmres_functions->MassDotpTwo       = MassDotpTwo;
62    cogmres_functions->CopyVector        = CopyVector;
63    cogmres_functions->ClearVector       = ClearVector;
64    cogmres_functions->ScaleVector       = ScaleVector;
65    cogmres_functions->Axpy              = Axpy;
66    cogmres_functions->MassAxpy          = MassAxpy;
67    /* default preconditioner must be set here but can be changed later... */
68    cogmres_functions->precond_setup     = PrecondSetup;
69    cogmres_functions->precond           = Precond;
70 
71    return cogmres_functions;
72 }
73 
74 /*--------------------------------------------------------------------------
75  * hypre_COGMRESCreate
76  *--------------------------------------------------------------------------*/
77 
78 void *
hypre_COGMRESCreate(hypre_COGMRESFunctions * cogmres_functions)79 hypre_COGMRESCreate( hypre_COGMRESFunctions *cogmres_functions )
80 {
81    hypre_COGMRESData *cogmres_data;
82 
83    HYPRE_ANNOTATE_FUNC_BEGIN;
84 
85    cogmres_data = hypre_CTAllocF(hypre_COGMRESData, 1, cogmres_functions, HYPRE_MEMORY_HOST);
86    cogmres_data->functions = cogmres_functions;
87 
88    /* set defaults */
89    (cogmres_data -> k_dim)          = 5;
90    (cogmres_data -> cgs)            = 1; /* if 2 performs reorthogonalization */
91    (cogmres_data -> tol)            = 1.0e-06; /* relative residual tol */
92    (cogmres_data -> cf_tol)         = 0.0;
93    (cogmres_data -> a_tol)          = 0.0; /* abs. residual tol */
94    (cogmres_data -> min_iter)       = 0;
95    (cogmres_data -> max_iter)       = 1000;
96    (cogmres_data -> rel_change)     = 0;
97    (cogmres_data -> skip_real_r_check) = 0;
98    (cogmres_data -> converged)      = 0;
99    (cogmres_data -> precond_data)   = NULL;
100    (cogmres_data -> print_level)    = 0;
101    (cogmres_data -> logging)        = 0;
102    (cogmres_data -> p)              = NULL;
103    (cogmres_data -> r)              = NULL;
104    (cogmres_data -> w)              = NULL;
105    (cogmres_data -> w_2)            = NULL;
106    (cogmres_data -> matvec_data)    = NULL;
107    (cogmres_data -> norms)          = NULL;
108    (cogmres_data -> log_file_name)  = NULL;
109    (cogmres_data -> unroll)         = 0;
110 
111    HYPRE_ANNOTATE_FUNC_END;
112 
113    return (void *) cogmres_data;
114 }
115 
116 /*--------------------------------------------------------------------------
117  * hypre_COGMRESDestroy
118  *--------------------------------------------------------------------------*/
119 
120 HYPRE_Int
hypre_COGMRESDestroy(void * cogmres_vdata)121 hypre_COGMRESDestroy( void *cogmres_vdata )
122 {
123    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
124    HYPRE_Int i;
125 
126    HYPRE_ANNOTATE_FUNC_BEGIN;
127    if (cogmres_data)
128    {
129       hypre_COGMRESFunctions *cogmres_functions = cogmres_data->functions;
130       if ( (cogmres_data->logging>0) || (cogmres_data->print_level) > 0 )
131       {
132          if ( (cogmres_data -> norms) != NULL )
133          hypre_TFreeF( cogmres_data -> norms, cogmres_functions );
134       }
135 
136       if ( (cogmres_data -> matvec_data) != NULL )
137          (*(cogmres_functions->MatvecDestroy))(cogmres_data -> matvec_data);
138 
139       if ( (cogmres_data -> r) != NULL )
140          (*(cogmres_functions->DestroyVector))(cogmres_data -> r);
141       if ( (cogmres_data -> w) != NULL )
142          (*(cogmres_functions->DestroyVector))(cogmres_data -> w);
143       if ( (cogmres_data -> w_2) != NULL )
144          (*(cogmres_functions->DestroyVector))(cogmres_data -> w_2);
145 
146 
147       if ( (cogmres_data -> p) != NULL )
148       {
149          for (i = 0; i < (cogmres_data -> k_dim+1); i++)
150          {
151             if ( (cogmres_data -> p)[i] != NULL )
152             (*(cogmres_functions->DestroyVector))( (cogmres_data -> p) [i]);
153          }
154          hypre_TFreeF( cogmres_data->p, cogmres_functions );
155       }
156       hypre_TFreeF( cogmres_data, cogmres_functions );
157       hypre_TFreeF( cogmres_functions, cogmres_functions );
158    }
159 
160    HYPRE_ANNOTATE_FUNC_END;
161 
162    return hypre_error_flag;
163 }
164 
165 /*--------------------------------------------------------------------------
166  * hypre_COGMRESGetResidual
167  *--------------------------------------------------------------------------*/
168 
hypre_COGMRESGetResidual(void * cogmres_vdata,void ** residual)169 HYPRE_Int hypre_COGMRESGetResidual( void *cogmres_vdata, void **residual )
170 {
171    hypre_COGMRESData  *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
172    *residual = cogmres_data->r;
173 
174    return hypre_error_flag;
175 }
176 
177 /*--------------------------------------------------------------------------
178  * hypre_COGMRESSetup
179  *--------------------------------------------------------------------------*/
180 
181 HYPRE_Int
hypre_COGMRESSetup(void * cogmres_vdata,void * A,void * b,void * x)182 hypre_COGMRESSetup( void *cogmres_vdata,
183                     void *A,
184                     void *b,
185                     void *x         )
186 {
187    hypre_COGMRESData *cogmres_data     = (hypre_COGMRESData *)cogmres_vdata;
188    hypre_COGMRESFunctions *cogmres_functions = cogmres_data->functions;
189 
190    HYPRE_Int k_dim            = (cogmres_data -> k_dim);
191    HYPRE_Int max_iter         = (cogmres_data -> max_iter);
192    HYPRE_Int (*precond_setup)(void*,void*,void*,void*) = (cogmres_functions->precond_setup);
193    void       *precond_data   = (cogmres_data -> precond_data);
194    HYPRE_Int rel_change       = (cogmres_data -> rel_change);
195 
196    HYPRE_ANNOTATE_FUNC_BEGIN;
197 
198    (cogmres_data -> A) = A;
199 
200    /*--------------------------------------------------
201     * The arguments for NewVector are important to
202     * maintain consistency between the setup and
203     * compute phases of matvec and the preconditioner.
204     *--------------------------------------------------*/
205 
206    if ((cogmres_data -> p) == NULL)
207       (cogmres_data -> p) = (void**)(*(cogmres_functions->CreateVectorArray))(k_dim+1,x);
208    if ((cogmres_data -> r) == NULL)
209       (cogmres_data -> r) = (*(cogmres_functions->CreateVector))(b);
210    if ((cogmres_data -> w) == NULL)
211       (cogmres_data -> w) = (*(cogmres_functions->CreateVector))(b);
212 
213    if (rel_change)
214    {
215       if ((cogmres_data -> w_2) == NULL)
216          (cogmres_data -> w_2) = (*(cogmres_functions->CreateVector))(b);
217    }
218 
219 
220    if ((cogmres_data -> matvec_data) == NULL)
221       (cogmres_data -> matvec_data) = (*(cogmres_functions->MatvecCreate))(A, x);
222 
223    precond_setup(precond_data, A, b, x);
224 
225    /*-----------------------------------------------------
226     * Allocate space for log info
227     *-----------------------------------------------------*/
228 
229    if ( (cogmres_data->logging)>0 || (cogmres_data->print_level) > 0 )
230    {
231       if ((cogmres_data -> norms) == NULL)
232          (cogmres_data -> norms) = hypre_CTAllocF(HYPRE_Real, max_iter + 1,cogmres_functions, HYPRE_MEMORY_HOST);
233    }
234    if ( (cogmres_data->print_level) > 0 )
235    {
236       if ((cogmres_data -> log_file_name) == NULL)
237          (cogmres_data -> log_file_name) = (char*)"cogmres.out.log";
238    }
239 
240    HYPRE_ANNOTATE_FUNC_END;
241 
242    return hypre_error_flag;
243 }
244 /*--------------------------------------------------------------------------
245  * hypre_COGMRESSolve
246  *-------------------------------------------------------------------------*/
247 
248 HYPRE_Int
hypre_COGMRESSolve(void * cogmres_vdata,void * A,void * b,void * x)249 hypre_COGMRESSolve(void  *cogmres_vdata,
250                    void  *A,
251                    void  *b,
252                    void  *x)
253 {
254 
255    hypre_COGMRESData      *cogmres_data      = (hypre_COGMRESData *)cogmres_vdata;
256    hypre_COGMRESFunctions *cogmres_functions = cogmres_data->functions;
257    HYPRE_Int     k_dim             = (cogmres_data -> k_dim);
258    HYPRE_Int     unroll            = (cogmres_data -> unroll);
259    HYPRE_Int     cgs               = (cogmres_data -> cgs);
260    HYPRE_Int     min_iter          = (cogmres_data -> min_iter);
261    HYPRE_Int     max_iter          = (cogmres_data -> max_iter);
262    HYPRE_Int     rel_change        = (cogmres_data -> rel_change);
263    HYPRE_Int     skip_real_r_check = (cogmres_data -> skip_real_r_check);
264    HYPRE_Real    r_tol             = (cogmres_data -> tol);
265    HYPRE_Real    cf_tol            = (cogmres_data -> cf_tol);
266    HYPRE_Real    a_tol             = (cogmres_data -> a_tol);
267    void         *matvec_data       = (cogmres_data -> matvec_data);
268 
269    void         *r                 = (cogmres_data -> r);
270    void         *w                 = (cogmres_data -> w);
271    /* note: w_2 is only allocated if rel_change = 1 */
272    void         *w_2               = (cogmres_data -> w_2);
273 
274    void        **p                 = (cogmres_data -> p);
275 
276    HYPRE_Int (*precond)(void*,void*,void*,void*) = (cogmres_functions -> precond);
277    HYPRE_Int  *precond_data       = (HYPRE_Int*)(cogmres_data -> precond_data);
278 
279    HYPRE_Int print_level = (cogmres_data -> print_level);
280    HYPRE_Int logging     = (cogmres_data -> logging);
281 
282    HYPRE_Real     *norms          = (cogmres_data -> norms);
283   /* not used yet   char           *log_file_name  = (cogmres_data -> log_file_name);*/
284   /*   FILE           *fp; */
285 
286    HYPRE_Int  break_value = 0;
287    HYPRE_Int  i, j, k;
288   /*KS: rv is the norm history */
289    HYPRE_Real *rs, *hh, *uu, *c, *s, *rs_2, *rv;
290   //, *tmp;
291    HYPRE_Int  iter;
292    HYPRE_Int  my_id, num_procs;
293    HYPRE_Real epsilon, gamma, t, r_norm, b_norm, den_norm, x_norm;
294    HYPRE_Real w_norm;
295 
296    HYPRE_Real epsmac = 1.e-16;
297    HYPRE_Real ieee_check = 0.;
298 
299    HYPRE_Real guard_zero_residual;
300    HYPRE_Real cf_ave_0 = 0.0;
301    HYPRE_Real cf_ave_1 = 0.0;
302    HYPRE_Real weight;
303    HYPRE_Real r_norm_0;
304    HYPRE_Real relative_error = 1.0;
305 
306    HYPRE_Int        rel_change_passed = 0, num_rel_change_check = 0;
307    HYPRE_Int    itmp = 0;
308 
309    HYPRE_Real real_r_norm_old, real_r_norm_new;
310 
311    HYPRE_ANNOTATE_FUNC_BEGIN;
312 
313    (cogmres_data -> converged) = 0;
314    /*-----------------------------------------------------------------------
315     * With relative change convergence test on, it is possible to attempt
316     * another iteration with a zero residual. This causes the parameter
317     * alpha to go NaN. The guard_zero_residual parameter is to circumvent
318     * this. Perhaps it should be set to something non-zero (but small).
319     *-----------------------------------------------------------------------*/
320    guard_zero_residual = 0.0;
321 
322    (*(cogmres_functions->CommInfo))(A,&my_id,&num_procs);
323    if ( logging>0 || print_level>0 )
324    {
325       norms = (cogmres_data -> norms);
326    }
327 
328    /* initialize work arrays */
329    rs = hypre_CTAllocF(HYPRE_Real,k_dim+1,cogmres_functions, HYPRE_MEMORY_HOST);
330    c  = hypre_CTAllocF(HYPRE_Real,k_dim,cogmres_functions, HYPRE_MEMORY_HOST);
331    s  = hypre_CTAllocF(HYPRE_Real,k_dim,cogmres_functions, HYPRE_MEMORY_HOST);
332    if (rel_change) rs_2 = hypre_CTAllocF(HYPRE_Real,k_dim+1,cogmres_functions, HYPRE_MEMORY_HOST);
333 
334    rv = hypre_CTAllocF(HYPRE_Real, k_dim+1, cogmres_functions, HYPRE_MEMORY_HOST);
335 
336    hh = hypre_CTAllocF(HYPRE_Real, (k_dim+1)*k_dim, cogmres_functions, HYPRE_MEMORY_HOST);
337    uu = hypre_CTAllocF(HYPRE_Real, (k_dim+1)*k_dim, cogmres_functions, HYPRE_MEMORY_HOST);
338 
339    (*(cogmres_functions->CopyVector))(b,p[0]);
340 
341    /* compute initial residual */
342    (*(cogmres_functions->Matvec))(matvec_data,-1.0, A, x, 1.0, p[0]);
343 
344    b_norm = sqrt((*(cogmres_functions->InnerProd))(b,b));
345    real_r_norm_old = b_norm;
346 
347    /* Since it is does not diminish performance, attempt to return an error flag
348       and notify users when they supply bad input. */
349    if (b_norm != 0.) ieee_check = b_norm/b_norm; /* INF -> NaN conversion */
350    if (ieee_check != ieee_check)
351    {
352       /* ...INFs or NaNs in input can make ieee_check a NaN.  This test
353          for ieee_check self-equality works on all IEEE-compliant compilers/
354          machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
355          by W. Kahan, May 31, 1996.  Currently (July 2002) this paper may be
356          found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
357       if (logging > 0 || print_level > 0)
358       {
359          hypre_printf("\n\nERROR detected by Hypre ... BEGIN\n");
360          hypre_printf("ERROR -- hypre_COGMRESSolve: INFs and/or NaNs detected in input.\n");
361          hypre_printf("User probably placed non-numerics in supplied b.\n");
362          hypre_printf("Returning error flag += 101.  Program not terminated.\n");
363          hypre_printf("ERROR detected by Hypre ... END\n\n\n");
364       }
365       hypre_error(HYPRE_ERROR_GENERIC);
366       HYPRE_ANNOTATE_FUNC_END;
367 
368       return hypre_error_flag;
369    }
370 
371    r_norm   = sqrt((*(cogmres_functions->InnerProd))(p[0],p[0]));
372    r_norm_0 = r_norm;
373 
374    /* Since it is does not diminish performance, attempt to return an error flag
375       and notify users when they supply bad input. */
376    if (r_norm != 0.) ieee_check = r_norm/r_norm; /* INF -> NaN conversion */
377    if (ieee_check != ieee_check)
378    {
379       /* ...INFs or NaNs in input can make ieee_check a NaN.  This test
380          for ieee_check self-equality works on all IEEE-compliant compilers/
381          machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
382          by W. Kahan, May 31, 1996.  Currently (July 2002) this paper may be
383          found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
384       if (logging > 0 || print_level > 0)
385       {
386          hypre_printf("\n\nERROR detected by Hypre ... BEGIN\n");
387          hypre_printf("ERROR -- hypre_COGMRESSolve: INFs and/or NaNs detected in input.\n");
388          hypre_printf("User probably placed non-numerics in supplied A or x_0.\n");
389          hypre_printf("Returning error flag += 101.  Program not terminated.\n");
390          hypre_printf("ERROR detected by Hypre ... END\n\n\n");
391       }
392       hypre_error(HYPRE_ERROR_GENERIC);
393       HYPRE_ANNOTATE_FUNC_END;
394 
395       return hypre_error_flag;
396    }
397 
398    if ( logging>0 || print_level > 0)
399    {
400       norms[0] = r_norm;
401       if ( print_level>1 && my_id == 0 )
402       {
403          hypre_printf("L2 norm of b: %e\n", b_norm);
404          if (b_norm == 0.0)
405             hypre_printf("Rel_resid_norm actually contains the residual norm\n");
406          hypre_printf("Initial L2 norm of residual: %e\n", r_norm);
407       }
408    }
409    iter = 0;
410 
411    if (b_norm > 0.0)
412    {
413       /* convergence criterion |r_i|/|b| <= accuracy if |b| > 0 */
414       den_norm = b_norm;
415    }
416    else
417    {
418       /* convergence criterion |r_i|/|r0| <= accuracy if |b| = 0 */
419       den_norm = r_norm;
420    };
421 
422    /* convergence criteria: |r_i| <= max( a_tol, r_tol * den_norm)
423       den_norm = |r_0| or |b|
424       note: default for a_tol is 0.0, so relative residual criteria is used unless
425       user specifies a_tol, or sets r_tol = 0.0, which means absolute
426       tol only is checked  */
427 
428    epsilon = hypre_max(a_tol,r_tol*den_norm);
429 
430    /* so now our stop criteria is |r_i| <= epsilon */
431 
432    if ( print_level>1 && my_id == 0 )
433    {
434       if (b_norm > 0.0)
435       {
436          hypre_printf("=============================================\n\n");
437          hypre_printf("Iters     resid.norm     conv.rate  rel.res.norm\n");
438          hypre_printf("-----    ------------    ---------- ------------\n");
439 
440       }
441       else
442       {
443          hypre_printf("=============================================\n\n");
444          hypre_printf("Iters     resid.norm     conv.rate\n");
445          hypre_printf("-----    ------------    ----------\n");
446       };
447    }
448 
449 
450    /* once the rel. change check has passed, we do not want to check it again */
451    rel_change_passed = 0;
452 
453    while (iter < max_iter)
454    {
455       /* initialize first term of hessenberg system */
456       rs[0] = r_norm;
457       if (r_norm == 0.0)
458       {
459          hypre_TFreeF(c,cogmres_functions);
460          hypre_TFreeF(s,cogmres_functions);
461          hypre_TFreeF(rs,cogmres_functions);
462          hypre_TFreeF(rv,cogmres_functions);
463          if (rel_change)  hypre_TFreeF(rs_2,cogmres_functions);
464          hypre_TFreeF(hh,cogmres_functions);
465          hypre_TFreeF(uu,cogmres_functions);
466          HYPRE_ANNOTATE_FUNC_END;
467 
468          return hypre_error_flag;
469       }
470 
471       /* see if we are already converged and
472          should print the final norm and exit */
473 
474       if (r_norm  <= epsilon && iter >= min_iter)
475       {
476          if (!rel_change) /* shouldn't exit after no iterations if
477                            * relative change is on*/
478          {
479             (*(cogmres_functions->CopyVector))(b,r);
480             (*(cogmres_functions->Matvec))(matvec_data,-1.0,A,x,1.0,r);
481             r_norm = sqrt((*(cogmres_functions->InnerProd))(r,r));
482             if (r_norm  <= epsilon)
483             {
484                if ( print_level>1 && my_id == 0)
485                {
486                   hypre_printf("\n\n");
487                   hypre_printf("Final L2 norm of residual: %e\n\n", r_norm);
488                }
489                break;
490             }
491             else if ( print_level>0 && my_id == 0)
492                hypre_printf("false convergence 1\n");
493          }
494       }
495 
496 
497 
498       t = 1.0 / r_norm;
499       (*(cogmres_functions->ScaleVector))(t,p[0]);
500       i = 0;
501       /***RESTART CYCLE (right-preconditioning) ***/
502       while (i < k_dim && iter < max_iter)
503       {
504          i++;
505          iter++;
506          itmp = (i-1)*(k_dim+1);
507 
508          (*(cogmres_functions->ClearVector))(r);
509 
510          precond(precond_data, A, p[i-1], r);
511          (*(cogmres_functions->Matvec))(matvec_data, 1.0, A, r, 0.0, p[i]);
512          for (j=0; j<i; j++)
513             rv[j]  = 0;
514 
515          if (cgs > 1)
516          {
517             (*(cogmres_functions->MassDotpTwo))((void *) p[i], p[i-1], p, i, unroll, &hh[itmp], &uu[itmp]);
518             for (j=0; j<i-1; j++) uu[j*(k_dim+1)+i-1] = uu[itmp+j];
519             for (j=0; j<i; j++) rv[j] = hh[itmp+j];
520             for (k=0; k < i; k++)
521             {
522                for (j=0; j < i; j++)
523                {
524                   hh[itmp+j] -= (uu[k*(k_dim+1)+j]*rv[j]);
525                }
526             }
527             for (j=0; j<i; j++)
528                hh[itmp+j]  = -rv[j]-hh[itmp+j];
529          }
530          else
531          {
532             (*(cogmres_functions->MassInnerProd))((void *) p[i], p, i, unroll, &hh[itmp]);
533             for (j=0; j<i; j++)
534                hh[itmp+j]  = -hh[itmp+j];
535          }
536 
537          (*(cogmres_functions->MassAxpy))(&hh[itmp],p,p[i], i, unroll);
538          for (j=0; j<i; j++)
539             hh[itmp+j]  = -hh[itmp+j];
540          t = sqrt( (*(cogmres_functions->InnerProd))(p[i],p[i]) );
541          hh[itmp+i] = t;
542 
543          if (hh[itmp+i] != 0.0)
544          {
545             t = 1.0/t;
546             (*(cogmres_functions->ScaleVector))(t,p[i]);
547          }
548          for (j = 1; j < i; j++)
549          {
550             t = hh[itmp+j-1];
551             hh[itmp+j-1] = s[j-1]*hh[itmp+j] + c[j-1]*t;
552             hh[itmp+j] = -s[j-1]*t + c[j-1]*hh[itmp+j];
553          }
554          t= hh[itmp+i]*hh[itmp+i];
555          t+= hh[itmp+i-1]*hh[itmp+i-1];
556          gamma = sqrt(t);
557          if (gamma == 0.0) gamma = epsmac;
558          c[i-1] = hh[itmp+i-1]/gamma;
559          s[i-1] = hh[itmp+i]/gamma;
560          rs[i] = -hh[itmp+i]*rs[i-1];
561          rs[i] /=  gamma;
562          rs[i-1] = c[i-1]*rs[i-1];
563          // determine residual norm
564          hh[itmp+i-1] = s[i-1]*hh[itmp+i] + c[i-1]*hh[itmp+i-1];
565          r_norm = fabs(rs[i]);
566          if ( print_level>0 )
567          {
568             norms[iter] = r_norm;
569             if ( print_level>1 && my_id == 0 )
570             {
571                if (b_norm > 0.0)
572                   hypre_printf("% 5d    %e    %f   %e\n", iter,
573                      norms[iter],norms[iter]/norms[iter-1],
574                      norms[iter]/b_norm);
575                else
576                   hypre_printf("% 5d    %e    %f\n", iter, norms[iter],
577                      norms[iter]/norms[iter-1]);
578             }
579          }
580          /*convergence factor tolerance */
581          if (cf_tol > 0.0)
582          {
583             cf_ave_0 = cf_ave_1;
584             cf_ave_1 = pow( r_norm / r_norm_0, 1.0/(2.0*iter));
585 
586             weight = fabs(cf_ave_1 - cf_ave_0);
587             weight = weight / hypre_max(cf_ave_1, cf_ave_0);
588 
589             weight = 1.0 - weight;
590 #if 0
591            hypre_printf("I = %d: cf_new = %e, cf_old = %e, weight = %e\n",
592               i, cf_ave_1, cf_ave_0, weight );
593 #endif
594             if (weight * cf_ave_1 > cf_tol)
595             {
596                break_value = 1;
597                break;
598             }
599          }
600          /* should we exit the restart cycle? (conv. check) */
601          if (r_norm <= epsilon && iter >= min_iter)
602          {
603             if (rel_change && !rel_change_passed)
604             {
605                /* To decide whether to break here: to actually
606                   determine the relative change requires the approx
607                   solution (so a triangular solve) and a
608                   precond. solve - so if we have to do this many
609                   times, it will be expensive...(unlike cg where is
610                   is relatively straightforward)
611                   previously, the intent (there was a bug), was to
612                   exit the restart cycle based on the residual norm
613                   and check the relative change outside the cycle.
614                   Here we will check the relative here as we don't
615                   want to exit the restart cycle prematurely */
616                for (k=0; k<i; k++) /* extra copy of rs so we don't need
617                                    to change the later solve */
618                   rs_2[k] = rs[k];
619 
620                /* solve tri. system*/
621                rs_2[i-1] = rs_2[i-1]/hh[itmp+i-1];
622                for (k = i-2; k >= 0; k--)
623                {
624                   t = 0.0;
625                   for (j = k+1; j < i; j++)
626                   {
627                      t -= hh[j*(k_dim+1)+k]*rs_2[j];
628                   }
629                   t+= rs_2[k];
630                   rs_2[k] = t/hh[k*(k_dim+1)+k];
631                }
632                (*(cogmres_functions->CopyVector))(p[i-1],w);
633                (*(cogmres_functions->ScaleVector))(rs_2[i-1],w);
634                for (j = i-2; j >=0; j--)
635                   (*(cogmres_functions->Axpy))(rs_2[j], p[j], w);
636 
637                (*(cogmres_functions->ClearVector))(r);
638                /* find correction (in r) */
639                precond(precond_data, A, w, r);
640                /* copy current solution (x) to w (don't want to over-write x)*/
641                (*(cogmres_functions->CopyVector))(x,w);
642 
643                /* add the correction */
644                (*(cogmres_functions->Axpy))(1.0,r,w);
645 
646                /* now w is the approx solution  - get the norm*/
647                x_norm = sqrt( (*(cogmres_functions->InnerProd))(w,w) );
648 
649                if ( !(x_norm <= guard_zero_residual ))
650                   /* don't divide by zero */
651                {  /* now get  x_i - x_i-1 */
652                   if (num_rel_change_check)
653                   {
654                      /* have already checked once so we can avoid another precond.
655                         solve */
656                      (*(cogmres_functions->CopyVector))(w, r);
657                      (*(cogmres_functions->Axpy))(-1.0, w_2, r);
658                      /* now r contains x_i - x_i-1*/
659 
660                      /* save current soln w in w_2 for next time */
661                      (*(cogmres_functions->CopyVector))(w, w_2);
662                   }
663                   else
664                   {
665                      /* first time to check rel change*/
666                      /* first save current soln w in w_2 for next time */
667                      (*(cogmres_functions->CopyVector))(w, w_2);
668 
669                      (*(cogmres_functions->ClearVector))(w);
670                      (*(cogmres_functions->Axpy))(rs_2[i-1], p[i-1], w);
671                      (*(cogmres_functions->ClearVector))(r);
672                      /* apply the preconditioner */
673                      precond(precond_data, A, w, r);
674                      /* now r contains x_i - x_i-1 */
675                   }
676                   /* find the norm of x_i - x_i-1 */
677                   w_norm = sqrt( (*(cogmres_functions->InnerProd))(r,r) );
678                   relative_error = w_norm/x_norm;
679                   if (relative_error <= r_tol)
680                   {
681                      rel_change_passed = 1;
682                      break;
683                   }
684                }
685                else
686                {
687                   rel_change_passed = 1;
688                   break;
689                }
690                num_rel_change_check++;
691             }
692             else /* no relative change */
693             {
694                break;
695             }
696          }
697       } /*** end of restart cycle ***/
698 
699       /* now compute solution, first solve upper triangular system */
700       if (break_value) break;
701 
702       rs[i-1] = rs[i-1]/hh[itmp+i-1];
703       for (k = i-2; k >= 0; k--)
704       {
705          t = 0.0;
706          for (j = k+1; j < i; j++)
707          {
708             t -= hh[j*(k_dim+1)+k]*rs[j];
709          }
710          t+= rs[k];
711          rs[k] = t/hh[k*(k_dim+1)+k];
712       }
713 
714       (*(cogmres_functions->CopyVector))(p[i-1],w);
715       (*(cogmres_functions->ScaleVector))(rs[i-1],w);
716       for (j = i-2; j >=0; j--)
717          (*(cogmres_functions->Axpy))(rs[j], p[j], w);
718 
719       (*(cogmres_functions->ClearVector))(r);
720       /* find correction (in r) */
721       precond(precond_data, A, w, r);
722 
723       /* update current solution x (in x) */
724       (*(cogmres_functions->Axpy))(1.0,r,x);
725 
726 
727       /* check for convergence by evaluating the actual residual */
728       if (r_norm  <= epsilon && iter >= min_iter)
729       {
730          if (skip_real_r_check)
731          {
732             (cogmres_data -> converged) = 1;
733             break;
734          }
735 
736          /* calculate actual residual norm*/
737          (*(cogmres_functions->CopyVector))(b,r);
738          (*(cogmres_functions->Matvec))(matvec_data,-1.0,A,x,1.0,r);
739          real_r_norm_new = r_norm = sqrt( (*(cogmres_functions->InnerProd))(r,r) );
740 
741          if (r_norm <= epsilon)
742          {
743             if (rel_change && !rel_change_passed) /* calculate the relative change */
744             {
745                /* calculate the norm of the solution */
746                x_norm = sqrt( (*(cogmres_functions->InnerProd))(x,x) );
747 
748                if ( !(x_norm <= guard_zero_residual ))
749                /* don't divide by zero */
750                {
751                   (*(cogmres_functions->ClearVector))(w);
752                   (*(cogmres_functions->Axpy))(rs[i-1], p[i-1], w);
753                   (*(cogmres_functions->ClearVector))(r);
754                   /* apply the preconditioner */
755                   precond(precond_data, A, w, r);
756                   /* find the norm of x_i - x_i-1 */
757                   w_norm = sqrt( (*(cogmres_functions->InnerProd))(r,r) );
758                   relative_error= w_norm/x_norm;
759                   if ( relative_error < r_tol )
760                   {
761                      (cogmres_data -> converged) = 1;
762                      if ( print_level>1 && my_id == 0 )
763                      {
764                         hypre_printf("\n\n");
765                         hypre_printf("Final L2 norm of residual: %e\n\n", r_norm);
766                      }
767                      break;
768                   }
769                }
770                else
771                {
772                   (cogmres_data -> converged) = 1;
773                   if ( print_level>1 && my_id == 0 )
774                   {
775                      hypre_printf("\n\n");
776                      hypre_printf("Final L2 norm of residual: %e\n\n", r_norm);
777                   }
778                   break;
779                }
780             }
781             else /* don't need to check rel. change */
782             {
783                if ( print_level>1 && my_id == 0 )
784                {
785                   hypre_printf("\n\n");
786                   hypre_printf("Final L2 norm of residual: %e\n\n", r_norm);
787                }
788                (cogmres_data -> converged) = 1;
789                break;
790             }
791          }
792          else /* conv. has not occurred, according to true residual */
793          {
794             /* exit if the real residual norm has not decreased */
795             if (real_r_norm_new >= real_r_norm_old)
796             {
797                if (print_level > 1 && my_id == 0)
798                {
799                   hypre_printf("\n\n");
800                   hypre_printf("Final L2 norm of residual: %e\n\n", r_norm);
801                }
802                (cogmres_data -> converged) = 1;
803                break;
804             }
805             /* report discrepancy between real/COGMRES residuals and restart */
806             if ( print_level>0 && my_id == 0)
807                hypre_printf("false convergence 2, L2 norm of residual: %e\n", r_norm);
808             (*(cogmres_functions->CopyVector))(r,p[0]);
809             i = 0;
810             real_r_norm_old = real_r_norm_new;
811          }
812       } /* end of convergence check */
813 
814       /* compute residual vector and continue loop */
815       for (j=i ; j > 0; j--)
816       {
817          rs[j-1] = -s[j-1]*rs[j];
818          rs[j] = c[j-1]*rs[j];
819       }
820 
821       if (i) (*(cogmres_functions->Axpy))(rs[i]-1.0,p[i],p[i]);
822       for (j=i-1 ; j > 0; j--)
823          (*(cogmres_functions->Axpy))(rs[j],p[j],p[i]);
824 
825       if (i)
826       {
827          (*(cogmres_functions->Axpy))(rs[0]-1.0,p[0],p[0]);
828          (*(cogmres_functions->Axpy))(1.0,p[i],p[0]);
829       }
830 
831    } /* END of iteration while loop */
832 
833 
834    (cogmres_data -> num_iterations) = iter;
835    if (b_norm > 0.0)
836       (cogmres_data -> rel_residual_norm) = r_norm/b_norm;
837    if (b_norm == 0.0)
838       (cogmres_data -> rel_residual_norm) = r_norm;
839 
840    if (iter >= max_iter && r_norm > epsilon && epsilon > 0) hypre_error(HYPRE_ERROR_CONV);
841 
842    hypre_TFreeF(c,cogmres_functions);
843    hypre_TFreeF(s,cogmres_functions);
844    hypre_TFreeF(rs,cogmres_functions);
845    hypre_TFreeF(rv,cogmres_functions);
846    if (rel_change)  hypre_TFreeF(rs_2,cogmres_functions);
847 
848    /*for (i=0; i < k_dim+1; i++)
849    {
850       hypre_TFreeF(hh[i],cogmres_functions);
851       hypre_TFreeF(uu[i],cogmres_functions);
852    }*/
853    hypre_TFreeF(hh,cogmres_functions);
854    hypre_TFreeF(uu,cogmres_functions);
855 
856    HYPRE_ANNOTATE_FUNC_END;
857 
858    return hypre_error_flag;
859 }
860 
861 /*--------------------------------------------------------------------------
862  * hypre_COGMRESSetKDim, hypre_COGMRESGetKDim
863  *--------------------------------------------------------------------------*/
864 
865 HYPRE_Int
hypre_COGMRESSetKDim(void * cogmres_vdata,HYPRE_Int k_dim)866 hypre_COGMRESSetKDim( void   *cogmres_vdata,
867         HYPRE_Int   k_dim )
868 {
869    hypre_COGMRESData *cogmres_data =(hypre_COGMRESData *) cogmres_vdata;
870    (cogmres_data -> k_dim) = k_dim;
871    return hypre_error_flag;
872 }
873 
874 HYPRE_Int
hypre_COGMRESGetKDim(void * cogmres_vdata,HYPRE_Int * k_dim)875 hypre_COGMRESGetKDim( void   *cogmres_vdata,
876         HYPRE_Int * k_dim )
877 {
878    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
879    *k_dim = (cogmres_data -> k_dim);
880    return hypre_error_flag;
881 }
882 
883 /*--------------------------------------------------------------------------
884  * hypre_COGMRESSetUnroll, hypre_COGMRESGetUnroll
885  *--------------------------------------------------------------------------*/
886 
887 HYPRE_Int
hypre_COGMRESSetUnroll(void * cogmres_vdata,HYPRE_Int unroll)888 hypre_COGMRESSetUnroll( void   *cogmres_vdata,
889         HYPRE_Int   unroll )
890 {
891    hypre_COGMRESData *cogmres_data =(hypre_COGMRESData *) cogmres_vdata;
892    (cogmres_data -> unroll) = unroll;
893    return hypre_error_flag;
894 }
895 
896 HYPRE_Int
hypre_COGMRESGetUnroll(void * cogmres_vdata,HYPRE_Int * unroll)897 hypre_COGMRESGetUnroll( void   *cogmres_vdata,
898         HYPRE_Int * unroll )
899 {
900    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
901    *unroll = (cogmres_data -> unroll);
902    return hypre_error_flag;
903 }
904 
905 /*--------------------------------------------------------------------------
906  * hypre_COGMRESSetCGS, hypre_COGMRESGetCGS
907  *--------------------------------------------------------------------------*/
908 
909 HYPRE_Int
hypre_COGMRESSetCGS(void * cogmres_vdata,HYPRE_Int cgs)910 hypre_COGMRESSetCGS( void   *cogmres_vdata,
911         HYPRE_Int   cgs )
912 {
913    hypre_COGMRESData *cogmres_data =(hypre_COGMRESData *) cogmres_vdata;
914    (cogmres_data -> cgs) = cgs;
915    return hypre_error_flag;
916 }
917 
918 HYPRE_Int
hypre_COGMRESGetCGS(void * cogmres_vdata,HYPRE_Int * cgs)919 hypre_COGMRESGetCGS( void   *cogmres_vdata,
920         HYPRE_Int * cgs )
921 {
922    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
923    *cgs = (cogmres_data -> cgs);
924    return hypre_error_flag;
925 }
926 
927 /*--------------------------------------------------------------------------
928  * hypre_COGMRESSetTol, hypre_COGMRESGetTol
929  *--------------------------------------------------------------------------*/
930 
931 HYPRE_Int
hypre_COGMRESSetTol(void * cogmres_vdata,HYPRE_Real tol)932 hypre_COGMRESSetTol( void   *cogmres_vdata,
933         HYPRE_Real  tol       )
934 {
935    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
936    (cogmres_data -> tol) = tol;
937    return hypre_error_flag;
938 }
939 
940 HYPRE_Int
hypre_COGMRESGetTol(void * cogmres_vdata,HYPRE_Real * tol)941 hypre_COGMRESGetTol( void   *cogmres_vdata,
942         HYPRE_Real  * tol      )
943 {
944    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
945    *tol = (cogmres_data -> tol);
946    return hypre_error_flag;
947 }
948 /*--------------------------------------------------------------------------
949  * hypre_COGMRESSetAbsoluteTol, hypre_COGMRESGetAbsoluteTol
950  *--------------------------------------------------------------------------*/
951 
952 HYPRE_Int
hypre_COGMRESSetAbsoluteTol(void * cogmres_vdata,HYPRE_Real a_tol)953 hypre_COGMRESSetAbsoluteTol( void   *cogmres_vdata,
954         HYPRE_Real  a_tol       )
955 {
956    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
957    (cogmres_data -> a_tol) = a_tol;
958    return hypre_error_flag;
959 }
960 
961 HYPRE_Int
hypre_COGMRESGetAbsoluteTol(void * cogmres_vdata,HYPRE_Real * a_tol)962 hypre_COGMRESGetAbsoluteTol( void   *cogmres_vdata,
963         HYPRE_Real  * a_tol      )
964 {
965    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
966    *a_tol = (cogmres_data -> a_tol);
967    return hypre_error_flag;
968 }
969 /*--------------------------------------------------------------------------
970  * hypre_COGMRESSetConvergenceFactorTol, hypre_COGMRESGetConvergenceFactorTol
971  *--------------------------------------------------------------------------*/
972 
973 HYPRE_Int
hypre_COGMRESSetConvergenceFactorTol(void * cogmres_vdata,HYPRE_Real cf_tol)974 hypre_COGMRESSetConvergenceFactorTol( void   *cogmres_vdata,
975         HYPRE_Real  cf_tol       )
976 {
977    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
978    (cogmres_data -> cf_tol) = cf_tol;
979    return hypre_error_flag;
980 }
981 
982 HYPRE_Int
hypre_COGMRESGetConvergenceFactorTol(void * cogmres_vdata,HYPRE_Real * cf_tol)983 hypre_COGMRESGetConvergenceFactorTol( void   *cogmres_vdata,
984         HYPRE_Real * cf_tol       )
985 {
986    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
987    *cf_tol = (cogmres_data -> cf_tol);
988    return hypre_error_flag;
989 }
990 
991 /*--------------------------------------------------------------------------
992  * hypre_COGMRESSetMinIter, hypre_COGMRESGetMinIter
993  *--------------------------------------------------------------------------*/
994 
995 HYPRE_Int
hypre_COGMRESSetMinIter(void * cogmres_vdata,HYPRE_Int min_iter)996 hypre_COGMRESSetMinIter( void *cogmres_vdata,
997         HYPRE_Int   min_iter  )
998 {
999    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1000    (cogmres_data -> min_iter) = min_iter;
1001    return hypre_error_flag;
1002 }
1003 
1004 HYPRE_Int
hypre_COGMRESGetMinIter(void * cogmres_vdata,HYPRE_Int * min_iter)1005 hypre_COGMRESGetMinIter( void *cogmres_vdata,
1006         HYPRE_Int * min_iter  )
1007 {
1008    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1009    *min_iter = (cogmres_data -> min_iter);
1010    return hypre_error_flag;
1011 }
1012 
1013 /*--------------------------------------------------------------------------
1014  * hypre_COGMRESSetMaxIter, hypre_COGMRESGetMaxIter
1015  *--------------------------------------------------------------------------*/
1016 
1017 HYPRE_Int
hypre_COGMRESSetMaxIter(void * cogmres_vdata,HYPRE_Int max_iter)1018 hypre_COGMRESSetMaxIter( void *cogmres_vdata,
1019         HYPRE_Int   max_iter  )
1020 {
1021    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1022    (cogmres_data -> max_iter) = max_iter;
1023    return hypre_error_flag;
1024 }
1025 
1026 HYPRE_Int
hypre_COGMRESGetMaxIter(void * cogmres_vdata,HYPRE_Int * max_iter)1027 hypre_COGMRESGetMaxIter( void *cogmres_vdata,
1028         HYPRE_Int * max_iter  )
1029 {
1030    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1031    *max_iter = (cogmres_data -> max_iter);
1032    return hypre_error_flag;
1033 }
1034 
1035 /*--------------------------------------------------------------------------
1036  * hypre_COGMRESSetRelChange, hypre_COGMRESGetRelChange
1037  *--------------------------------------------------------------------------*/
1038 
1039 HYPRE_Int
hypre_COGMRESSetRelChange(void * cogmres_vdata,HYPRE_Int rel_change)1040 hypre_COGMRESSetRelChange( void *cogmres_vdata,
1041         HYPRE_Int   rel_change  )
1042 {
1043    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1044    (cogmres_data -> rel_change) = rel_change;
1045    return hypre_error_flag;
1046 }
1047 
1048 HYPRE_Int
hypre_COGMRESGetRelChange(void * cogmres_vdata,HYPRE_Int * rel_change)1049 hypre_COGMRESGetRelChange( void *cogmres_vdata,
1050         HYPRE_Int * rel_change  )
1051 {
1052    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1053    *rel_change = (cogmres_data -> rel_change);
1054    return hypre_error_flag;
1055 }
1056 
1057 /*--------------------------------------------------------------------------
1058  * hypre_COGMRESSetSkipRealResidualCheck, hypre_COGMRESGetSkipRealResidualCheck
1059  *--------------------------------------------------------------------------*/
1060 
1061 HYPRE_Int
hypre_COGMRESSetSkipRealResidualCheck(void * cogmres_vdata,HYPRE_Int skip_real_r_check)1062 hypre_COGMRESSetSkipRealResidualCheck( void *cogmres_vdata,
1063         HYPRE_Int skip_real_r_check )
1064 {
1065    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1066    (cogmres_data -> skip_real_r_check) = skip_real_r_check;
1067    return hypre_error_flag;
1068 }
1069 
1070 HYPRE_Int
hypre_COGMRESGetSkipRealResidualCheck(void * cogmres_vdata,HYPRE_Int * skip_real_r_check)1071 hypre_COGMRESGetSkipRealResidualCheck( void *cogmres_vdata,
1072         HYPRE_Int *skip_real_r_check)
1073 {
1074    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1075    *skip_real_r_check = (cogmres_data -> skip_real_r_check);
1076    return hypre_error_flag;
1077 }
1078 
1079 /*--------------------------------------------------------------------------
1080  * hypre_COGMRESSetPrecond
1081  *--------------------------------------------------------------------------*/
1082 
1083 HYPRE_Int
hypre_COGMRESSetPrecond(void * cogmres_vdata,HYPRE_Int (* precond)(void *,void *,void *,void *),HYPRE_Int (* precond_setup)(void *,void *,void *,void *),void * precond_data)1084 hypre_COGMRESSetPrecond( void  *cogmres_vdata,
1085         HYPRE_Int  (*precond)(void*,void*,void*,void*),
1086         HYPRE_Int  (*precond_setup)(void*,void*,void*,void*),
1087         void  *precond_data )
1088 {
1089    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1090    hypre_COGMRESFunctions *cogmres_functions = cogmres_data->functions;
1091    (cogmres_functions -> precond)        = precond;
1092    (cogmres_functions -> precond_setup)  = precond_setup;
1093    (cogmres_data -> precond_data)   = precond_data;
1094    return hypre_error_flag;
1095 }
1096 
1097 /*--------------------------------------------------------------------------
1098  * hypre_COGMRESGetPrecond
1099  *--------------------------------------------------------------------------*/
1100 
1101 HYPRE_Int
hypre_COGMRESGetPrecond(void * cogmres_vdata,HYPRE_Solver * precond_data_ptr)1102 hypre_COGMRESGetPrecond( void         *cogmres_vdata,
1103         HYPRE_Solver *precond_data_ptr )
1104 {
1105    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1106    *precond_data_ptr = (HYPRE_Solver)(cogmres_data -> precond_data);
1107    return hypre_error_flag;
1108 }
1109 
1110 /*--------------------------------------------------------------------------
1111  * hypre_COGMRESSetPrintLevel, hypre_COGMRESGetPrintLevel
1112  *--------------------------------------------------------------------------*/
1113 
1114 HYPRE_Int
hypre_COGMRESSetPrintLevel(void * cogmres_vdata,HYPRE_Int level)1115 hypre_COGMRESSetPrintLevel( void *cogmres_vdata,
1116         HYPRE_Int   level)
1117 {
1118    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1119    (cogmres_data -> print_level) = level;
1120    return hypre_error_flag;
1121 }
1122 
1123 HYPRE_Int
hypre_COGMRESGetPrintLevel(void * cogmres_vdata,HYPRE_Int * level)1124 hypre_COGMRESGetPrintLevel( void *cogmres_vdata,
1125         HYPRE_Int * level)
1126 {
1127    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1128    *level = (cogmres_data -> print_level);
1129    return hypre_error_flag;
1130 }
1131 
1132 /*--------------------------------------------------------------------------
1133  * hypre_COGMRESSetLogging, hypre_COGMRESGetLogging
1134  *--------------------------------------------------------------------------*/
1135 
1136 HYPRE_Int
hypre_COGMRESSetLogging(void * cogmres_vdata,HYPRE_Int level)1137 hypre_COGMRESSetLogging( void *cogmres_vdata,
1138         HYPRE_Int   level)
1139 {
1140    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1141    (cogmres_data -> logging) = level;
1142    return hypre_error_flag;
1143 }
1144 
1145 HYPRE_Int
hypre_COGMRESGetLogging(void * cogmres_vdata,HYPRE_Int * level)1146 hypre_COGMRESGetLogging( void *cogmres_vdata,
1147         HYPRE_Int * level)
1148 {
1149    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1150    *level = (cogmres_data -> logging);
1151    return hypre_error_flag;
1152 }
1153 
1154 /*--------------------------------------------------------------------------
1155  * hypre_COGMRESGetNumIterations
1156  *--------------------------------------------------------------------------*/
1157 
1158 HYPRE_Int
hypre_COGMRESGetNumIterations(void * cogmres_vdata,HYPRE_Int * num_iterations)1159 hypre_COGMRESGetNumIterations( void *cogmres_vdata,
1160         HYPRE_Int  *num_iterations )
1161 {
1162    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1163    *num_iterations = (cogmres_data -> num_iterations);
1164    return hypre_error_flag;
1165 }
1166 
1167 /*--------------------------------------------------------------------------
1168  * hypre_COGMRESGetConverged
1169  *--------------------------------------------------------------------------*/
1170 
1171 HYPRE_Int
hypre_COGMRESGetConverged(void * cogmres_vdata,HYPRE_Int * converged)1172 hypre_COGMRESGetConverged( void *cogmres_vdata,
1173         HYPRE_Int  *converged )
1174 {
1175    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1176    *converged = (cogmres_data -> converged);
1177    return hypre_error_flag;
1178 }
1179 
1180 /*--------------------------------------------------------------------------
1181  * hypre_COGMRESGetFinalRelativeResidualNorm
1182  *--------------------------------------------------------------------------*/
1183 
1184 HYPRE_Int
hypre_COGMRESGetFinalRelativeResidualNorm(void * cogmres_vdata,HYPRE_Real * relative_residual_norm)1185 hypre_COGMRESGetFinalRelativeResidualNorm( void   *cogmres_vdata,
1186         HYPRE_Real *relative_residual_norm )
1187 {
1188    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1189    *relative_residual_norm = (cogmres_data -> rel_residual_norm);
1190    return hypre_error_flag;
1191 }
1192 
1193 
1194 HYPRE_Int
hypre_COGMRESSetModifyPC(void * cogmres_vdata,HYPRE_Int (* modify_pc)(void * precond_data,HYPRE_Int iteration,HYPRE_Real rel_residual_norm))1195 hypre_COGMRESSetModifyPC(void *cogmres_vdata,
1196       HYPRE_Int (*modify_pc)(void *precond_data, HYPRE_Int iteration, HYPRE_Real rel_residual_norm))
1197 {
1198    hypre_COGMRESData *cogmres_data = (hypre_COGMRESData *)cogmres_vdata;
1199    hypre_COGMRESFunctions *cogmres_functions = cogmres_data->functions;
1200    (cogmres_functions -> modify_pc)        = modify_pc;
1201    return hypre_error_flag;
1202 }
1203