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 #include "_hypre_FEI.h"
9
10 /******************************************************************************
11 *
12 * FGMRES - flexible gmres
13 *
14 *****************************************************************************/
15
16 #include "utilities/_hypre_utilities.h"
17 #include "HYPRE.h"
18 #include "IJ_mv/HYPRE_IJ_mv.h"
19 #include "parcsr_mv/HYPRE_parcsr_mv.h"
20 #include "parcsr_mv/_hypre_parcsr_mv.h"
21 #include "parcsr_ls/_hypre_parcsr_ls.h"
22 #include "parcsr_ls/HYPRE_parcsr_ls.h"
23
24 /*--------------------------------------------------------------------------
25 * hypre_FGMRESData
26 *--------------------------------------------------------------------------*/
27
28 typedef struct
29 {
30 int max_iter;
31 int stop_crit;
32 int k_dim;
33 double tol;
34 double rel_residual_norm;
35 void *A;
36 void *w;
37 void **p;
38 void **z;
39 void *r;
40 void *matvec_data;
41 int (*precond)(void*, void*, void*, void*);
42 int (*precond_setup)(void*, void*, void*, void*);
43 void *precond_data;
44 int num_iterations;
45 int logging;
46 double *norms;
47 char *log_file_name;
48 int precond_tol_update;
49 int (*precond_update_tol)(int*,double);
50
51 } hypre_FGMRESData;
52
53 /*--------------------------------------------------------------------------
54 * hypre_FGMRESCreate
55 *--------------------------------------------------------------------------*/
56
hypre_FGMRESCreate()57 void *hypre_FGMRESCreate()
58 {
59 hypre_FGMRESData *fgmres_data;
60
61 fgmres_data = hypre_CTAlloc(hypre_FGMRESData, 1, HYPRE_MEMORY_HOST);
62
63 /* set defaults */
64
65 (fgmres_data -> k_dim) = 5;
66 (fgmres_data -> tol) = 1.0e-06;
67 (fgmres_data -> max_iter) = 1000;
68 (fgmres_data -> stop_crit) = 0; /* rel. residual norm */
69 (fgmres_data -> precond) = hypre_ParKrylovIdentity;
70 (fgmres_data -> precond_setup) = hypre_ParKrylovIdentitySetup;
71 (fgmres_data -> precond_data) = NULL;
72 (fgmres_data -> logging) = 0;
73 (fgmres_data -> p) = NULL;
74 (fgmres_data -> z) = NULL;
75 (fgmres_data -> r) = NULL;
76 (fgmres_data -> w) = NULL;
77 (fgmres_data -> matvec_data) = NULL;
78 (fgmres_data -> norms) = NULL;
79 (fgmres_data -> log_file_name) = NULL;
80 (fgmres_data -> logging) = 0;
81 (fgmres_data -> precond_tol_update) = 0;
82 (fgmres_data -> precond_update_tol) = NULL;
83 return (void *) fgmres_data;
84 }
85
86 /*--------------------------------------------------------------------------
87 * hypre_FGMRESDestroy
88 *--------------------------------------------------------------------------*/
89
hypre_FGMRESDestroy(void * fgmres_vdata)90 int hypre_FGMRESDestroy( void *fgmres_vdata )
91 {
92 int i, ierr=0;
93 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
94
95 if (fgmres_data)
96 {
97 if ( (fgmres_data->logging) > 0 && (fgmres_data->norms != NULL) )
98 hypre_TFree( fgmres_data -> norms , HYPRE_MEMORY_HOST);
99 if ( (fgmres_data->matvec_data) != NULL )
100 hypre_ParKrylovMatvecDestroy(fgmres_data -> matvec_data);
101 if ( (fgmres_data-> r) != NULL )
102 hypre_ParKrylovDestroyVector(fgmres_data -> r);
103 if ( (fgmres_data-> w) != NULL )
104 hypre_ParKrylovDestroyVector(fgmres_data -> w);
105 if ( (fgmres_data-> p) != NULL )
106 {
107 for (i = 0; i < (fgmres_data -> k_dim+1); i++)
108 hypre_ParKrylovDestroyVector((fgmres_data -> p)[i]);
109 hypre_TFree( fgmres_data -> p , HYPRE_MEMORY_HOST);
110 }
111 if ( (fgmres_data-> z) != NULL )
112 {
113 for (i = 0; i < (fgmres_data -> k_dim+1); i++)
114 hypre_ParKrylovDestroyVector((fgmres_data -> z)[i]);
115 hypre_TFree( fgmres_data -> z , HYPRE_MEMORY_HOST);
116 }
117 hypre_TFree( fgmres_data , HYPRE_MEMORY_HOST);
118 }
119 return(ierr);
120 }
121
122 /*--------------------------------------------------------------------------
123 * hypre_FGMRESSetup
124 *--------------------------------------------------------------------------*/
125
hypre_FGMRESSetup(void * fgmres_vdata,void * A,void * b,void * x)126 int hypre_FGMRESSetup( void *fgmres_vdata, void *A, void *b, void *x )
127 {
128 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
129 int k_dim = (fgmres_data -> k_dim);
130 int max_iter = (fgmres_data -> max_iter);
131 int (*precond_setup)(void*, void*, void*, void*) = (fgmres_data -> precond_setup);
132 void *precond_data = (fgmres_data -> precond_data);
133 int ierr = 0;
134
135 (fgmres_data -> A) = A;
136
137 if ((fgmres_data -> r) == NULL)
138 (fgmres_data -> r) = hypre_ParKrylovCreateVector(b);
139 if ((fgmres_data -> w) == NULL)
140 (fgmres_data -> w) = hypre_ParKrylovCreateVector(b);
141 if ((fgmres_data -> p) == NULL)
142 (fgmres_data -> p) = (void**) hypre_ParKrylovCreateVectorArray(k_dim+1,b);
143 if ((fgmres_data -> z) == NULL)
144 (fgmres_data -> z) = (void**) hypre_ParKrylovCreateVectorArray(k_dim+1,b);
145
146 if ((fgmres_data -> matvec_data) == NULL)
147 (fgmres_data -> matvec_data) = hypre_ParKrylovMatvecCreate(A, x);
148
149 ierr = precond_setup(precond_data, A, b, x);
150
151 if ((fgmres_data -> logging) > 0)
152 {
153 if ((fgmres_data -> norms) == NULL)
154 (fgmres_data -> norms) = hypre_CTAlloc(double, max_iter + 1, HYPRE_MEMORY_HOST);
155 if ((fgmres_data -> log_file_name) == NULL)
156 (fgmres_data -> log_file_name) = (char*) "fgmres.out.log";
157 }
158 return ierr;
159 }
160
161 /*--------------------------------------------------------------------------
162 * hypre_FGMRESSolve
163 *-------------------------------------------------------------------------*/
164
hypre_FGMRESSolve(void * fgmres_vdata,void * A,void * b,void * x)165 int hypre_FGMRESSolve(void *fgmres_vdata, void *A, void *b, void *x)
166 {
167 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
168 int k_dim = (fgmres_data -> k_dim);
169 int max_iter = (fgmres_data -> max_iter);
170 int stop_crit = (fgmres_data -> stop_crit);
171 double accuracy = (fgmres_data -> tol);
172 void *matvec_data = (fgmres_data -> matvec_data);
173
174 void *r = (fgmres_data -> r);
175 void **p = (fgmres_data -> p);
176 void **z = (fgmres_data -> z);
177
178 int (*precond)(void*, void*, void*, void*) = (fgmres_data -> precond);
179 int *precond_data = (int*) (fgmres_data -> precond_data);
180
181 int logging = (fgmres_data -> logging);
182 double *norms = (fgmres_data -> norms);
183
184 int tol_update = (fgmres_data -> precond_tol_update);
185 int (*update_tol)(int*,double)= (fgmres_data -> precond_update_tol);
186
187 int i, j, k, ierr = 0, iter, my_id, num_procs;
188 double *rs, **hh, *c, *s, t;
189 double epsilon, gamma, r_norm, b_norm, epsmac = 1.e-16;
190
191 hypre_ParKrylovCommInfo(A,&my_id,&num_procs);
192
193 /* initialize work arrays */
194
195 if (logging > 0) norms = (fgmres_data -> norms);
196 rs = hypre_CTAlloc(double, k_dim+1, HYPRE_MEMORY_HOST);
197 c = hypre_CTAlloc(double, k_dim, HYPRE_MEMORY_HOST);
198 s = hypre_CTAlloc(double, k_dim, HYPRE_MEMORY_HOST);
199 hh = hypre_CTAlloc(double*, k_dim+1, HYPRE_MEMORY_HOST);
200 for (i=0; i < k_dim+1; i++) hh[i] = hypre_CTAlloc(double, k_dim, HYPRE_MEMORY_HOST);
201 hypre_ParKrylovCopyVector(b,p[0]);
202
203 /* compute initial residual */
204
205 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, p[0]);
206 r_norm = sqrt(hypre_ParKrylovInnerProd(p[0],p[0]));
207 b_norm = sqrt(hypre_ParKrylovInnerProd(b,b));
208 if (logging > 0)
209 {
210 norms[0] = r_norm;
211 if (my_id == 0)
212 {
213 printf("FGMRES : L2 norm of b: %e\n", b_norm);
214 if (b_norm == 0.0)
215 printf("Rel_resid_norm actually contains the residual norm\n");
216 printf("FGMRES : Initial L2 norm of residual: %e\n", r_norm);
217 }
218
219 }
220 iter = 0;
221
222 if (b_norm > 0.0)
223 {
224 /* convergence criterion |r_i| <= accuracy*|b| if |b| > 0 */
225 epsilon = accuracy * b_norm;
226 }
227 else
228 {
229 /* convergence criterion |r_i| <= accuracy*|r0| if |b| = 0 */
230 epsilon = accuracy * r_norm;
231 };
232
233 /* convergence criterion |r_i| <= accuracy , absolute residual norm*/
234
235 if ( stop_crit ) epsilon = accuracy;
236
237 while (iter < max_iter)
238 {
239 /* initialize first term of hessenberg system */
240
241 rs[0] = r_norm;
242 if (r_norm == 0.0)
243 {
244 ierr = 0;
245 return ierr;
246 }
247
248 if (r_norm <= epsilon && iter > 0)
249 {
250 hypre_ParKrylovCopyVector(b,r);
251 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, r);
252 r_norm = sqrt(hypre_ParKrylovInnerProd(r,r));
253 if (r_norm <= epsilon)
254 {
255 if (logging > 0 && my_id == 0)
256 printf("Final L2 norm of residual: %e\n\n", r_norm);
257 break;
258 }
259 }
260
261 t = 1.0 / r_norm;
262 hypre_ParKrylovScaleVector(t,p[0]);
263 i = 0;
264 while (i < k_dim && r_norm > epsilon && iter < max_iter)
265 {
266 i++;
267 iter++;
268 hypre_ParKrylovClearVector(z[i-1]);
269
270 if ( tol_update != 0 && update_tol != NULL )
271 update_tol(precond_data,r_norm/b_norm);
272
273 precond(precond_data, A, p[i-1], z[i-1]);
274 hypre_ParKrylovMatvec(matvec_data, 1.0, A, z[i-1], 0.0, p[i]);
275
276 /* modified Gram_Schmidt */
277
278 for (j=0; j < i; j++)
279 {
280 hh[j][i-1] = hypre_ParKrylovInnerProd(p[j],p[i]);
281 hypre_ParKrylovAxpy(-hh[j][i-1],p[j],p[i]);
282 }
283 t = sqrt(hypre_ParKrylovInnerProd(p[i],p[i]));
284 hh[i][i-1] = t;
285 if (t != 0.0)
286 {
287 t = 1.0/t;
288 hypre_ParKrylovScaleVector(t, p[i]);
289 }
290
291 /* done with modified Gram_schmidt. update factorization of hh */
292
293 for (j = 1; j < i; j++)
294 {
295 t = hh[j-1][i-1];
296 hh[j-1][i-1] = c[j-1]*t + s[j-1]*hh[j][i-1];
297 hh[j][i-1] = -s[j-1]*t + c[j-1]*hh[j][i-1];
298 }
299 gamma = sqrt(hh[i-1][i-1]*hh[i-1][i-1] + hh[i][i-1]*hh[i][i-1]);
300 if (gamma == 0.0) gamma = epsmac;
301 c[i-1] = hh[i-1][i-1]/gamma;
302 s[i-1] = hh[i][i-1]/gamma;
303 rs[i] = -s[i-1]*rs[i-1];
304 rs[i-1] = c[i-1]*rs[i-1];
305
306 /* determine residual norm */
307
308 hh[i-1][i-1] = c[i-1]*hh[i-1][i-1] + s[i-1]*hh[i][i-1];
309 r_norm = fabs(rs[i]);
310 if (logging > 0)
311 {
312 norms[iter] = r_norm;
313 if (my_id == 0)
314 printf("FGMRES : iteration = %6d, norm of r = %e\n", iter,
315 r_norm);
316 }
317 }
318
319 /* now compute solution, first solve upper triangular system */
320
321 rs[i-1] = rs[i-1]/hh[i-1][i-1];
322 for (k = i-2; k >= 0; k--)
323 {
324 t = rs[k];
325 for (j = k+1; j < i; j++) t -= hh[k][j]*rs[j];
326 rs[k] = t/hh[k][k];
327 }
328
329
330 for (j = 0; j < i; j++) hypre_ParKrylovAxpy(rs[j], z[j], x);
331
332 /* check for convergence, evaluate actual residual */
333
334 hypre_ParKrylovCopyVector(b,p[0]);
335 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, p[0]);
336 r_norm = sqrt(hypre_ParKrylovInnerProd(p[0],p[0]));
337 if (r_norm <= epsilon)
338 {
339 if (logging > 0 && my_id == 0)
340 printf("FGMRES Final L2 norm of residual: %e\n\n", r_norm);
341 break;
342 }
343 }
344
345 (fgmres_data -> num_iterations) = iter;
346 if (b_norm > 0.0)
347 (fgmres_data -> rel_residual_norm) = r_norm/b_norm;
348 if (b_norm == 0.0)
349 (fgmres_data -> rel_residual_norm) = r_norm;
350
351 if (iter >= max_iter && r_norm > epsilon) ierr = 1;
352
353 hypre_TFree(c, HYPRE_MEMORY_HOST);
354 hypre_TFree(s, HYPRE_MEMORY_HOST);
355 hypre_TFree(rs, HYPRE_MEMORY_HOST);
356
357 for (i=0; i < k_dim+1; i++) hypre_TFree(hh[i], HYPRE_MEMORY_HOST);
358 hypre_TFree(hh, HYPRE_MEMORY_HOST);
359
360 return ierr;
361 }
362
363 /*--------------------------------------------------------------------------
364 * hypre_FGMRESSetKDim
365 *--------------------------------------------------------------------------*/
366
hypre_FGMRESSetKDim(void * fgmres_vdata,int k_dim)367 int hypre_FGMRESSetKDim( void *fgmres_vdata, int k_dim )
368 {
369 int ierr = 0;
370 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
371
372 (fgmres_data -> k_dim) = k_dim;
373
374 return ierr;
375 }
376
377 /*--------------------------------------------------------------------------
378 * hypre_FGMRESSetTol
379 *--------------------------------------------------------------------------*/
380
hypre_FGMRESSetTol(void * fgmres_vdata,double tol)381 int hypre_FGMRESSetTol( void *fgmres_vdata, double tol )
382 {
383 int ierr = 0;
384 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
385
386 (fgmres_data -> tol) = tol;
387
388 return ierr;
389 }
390
391 /*--------------------------------------------------------------------------
392 * hypre_FGMRESSetMaxIter
393 *--------------------------------------------------------------------------*/
394
hypre_FGMRESSetMaxIter(void * fgmres_vdata,int max_iter)395 int hypre_FGMRESSetMaxIter( void *fgmres_vdata, int max_iter )
396 {
397 int ierr = 0;
398 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
399
400 (fgmres_data -> max_iter) = max_iter;
401
402 return ierr;
403 }
404
405 /*--------------------------------------------------------------------------
406 * hypre_FGMRESSetStopCrit
407 *--------------------------------------------------------------------------*/
408
hypre_FGMRESSetStopCrit(void * fgmres_vdata,double stop_crit)409 int hypre_FGMRESSetStopCrit( void *fgmres_vdata, double stop_crit )
410 {
411 int ierr = 0;
412 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
413
414 (fgmres_data -> stop_crit) = stop_crit;
415
416 return ierr;
417 }
418
419 /*--------------------------------------------------------------------------
420 * hypre_FGMRESSetPrecond
421 *--------------------------------------------------------------------------*/
422
hypre_FGMRESSetPrecond(void * fgmres_vdata,int (* precond)(void *,void *,void *,void *),int (* precond_setup)(void *,void *,void *,void *),void * precond_data)423 int hypre_FGMRESSetPrecond( void *fgmres_vdata, int (*precond)(void*,void*,void*,void*),
424 int (*precond_setup)(void*,void*,void*,void*), void *precond_data )
425 {
426 int ierr = 0;
427 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
428
429 (fgmres_data -> precond) = precond;
430 (fgmres_data -> precond_setup) = precond_setup;
431 (fgmres_data -> precond_data) = precond_data;
432
433 return ierr;
434 }
435
436 /*--------------------------------------------------------------------------
437 * hypre_FGMRESGetPrecond
438 *--------------------------------------------------------------------------*/
439
hypre_FGMRESGetPrecond(void * fgmres_vdata,HYPRE_Solver * precond_data_ptr)440 int hypre_FGMRESGetPrecond(void *fgmres_vdata, HYPRE_Solver *precond_data_ptr)
441 {
442 int ierr = 0;
443 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
444
445 *precond_data_ptr = (HYPRE_Solver)(fgmres_data -> precond_data);
446
447 return ierr;
448 }
449
450 /*--------------------------------------------------------------------------
451 * hypre_FGMRESSetLogging
452 *--------------------------------------------------------------------------*/
453
hypre_FGMRESSetLogging(void * fgmres_vdata,int logging)454 int hypre_FGMRESSetLogging( void *fgmres_vdata, int logging )
455 {
456 int ierr = 0;
457 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
458
459 (fgmres_data -> logging) = logging;
460
461 return ierr;
462 }
463
464 /*--------------------------------------------------------------------------
465 * hypre_FGMRESGetNumIterations
466 *--------------------------------------------------------------------------*/
467
hypre_FGMRESGetNumIterations(void * fgmres_vdata,int * num_iterations)468 int hypre_FGMRESGetNumIterations( void *fgmres_vdata, int *num_iterations )
469 {
470 int ierr = 0;
471 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
472
473 *num_iterations = (fgmres_data -> num_iterations);
474
475 return ierr;
476 }
477
478 /*--------------------------------------------------------------------------
479 * hypre_FGMRESGetFinalRelativeResidualNorm
480 *--------------------------------------------------------------------------*/
481
hypre_FGMRESGetFinalRelativeResidualNorm(void * fgmres_vdata,double * relative_residual_norm)482 int hypre_FGMRESGetFinalRelativeResidualNorm(void *fgmres_vdata,
483 double *relative_residual_norm )
484 {
485 int ierr = 0;
486 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
487
488 *relative_residual_norm = (fgmres_data -> rel_residual_norm);
489
490 return ierr;
491 }
492
493 /*--------------------------------------------------------------------------
494 * hypre_FGMRESUpdatePrecondTolerance
495 *--------------------------------------------------------------------------*/
496
hypre_FGMRESUpdatePrecondTolerance(void * fgmres_vdata,int (* update_tol)(int *,double))497 int hypre_FGMRESUpdatePrecondTolerance(void *fgmres_vdata, int (*update_tol)(int*, double))
498 {
499 int ierr = 0;
500 hypre_FGMRESData *fgmres_data = (hypre_FGMRESData *) fgmres_vdata;
501
502 (fgmres_data -> precond_tol_update) = 1;
503 (fgmres_data -> precond_update_tol) = update_tol;
504 return ierr;
505 }
506
507