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 * Symmetric QMR
11 *
12 *****************************************************************************/
13
14 #include "utilities/_hypre_utilities.h"
15 #include "HYPRE.h"
16 #include "IJ_mv/HYPRE_IJ_mv.h"
17 #include "parcsr_mv/HYPRE_parcsr_mv.h"
18 #include "parcsr_mv/_hypre_parcsr_mv.h"
19 #include "parcsr_ls/_hypre_parcsr_ls.h"
20 #include "parcsr_ls/HYPRE_parcsr_ls.h"
21
22 /*--------------------------------------------------------------------------
23 * hypre_SymQMRData
24 *--------------------------------------------------------------------------*/
25
26 typedef struct
27 {
28 int max_iter;
29 int stop_crit;
30 double tol;
31 double rel_residual_norm;
32
33 void *A;
34 void *r;
35 void *q;
36 void *u;
37 void *d;
38 void *t;
39 void *rq;
40
41 void *matvec_data;
42
43 int (*precond)(void*,void*,void*,void*);
44 int (*precond_setup)(void*,void*,void*,void*);
45 void *precond_data;
46
47 /* log info (always logged) */
48 int num_iterations;
49
50 /* additional log info (logged when `logging' > 0) */
51 int logging;
52 double *norms;
53 char *log_file_name;
54
55 } hypre_SymQMRData;
56
57 /*--------------------------------------------------------------------------
58 * hypre_SymQMRCreate
59 *--------------------------------------------------------------------------*/
60
hypre_SymQMRCreate()61 void * hypre_SymQMRCreate( )
62 {
63 hypre_SymQMRData *symqmr_data;
64
65 symqmr_data = hypre_CTAlloc(hypre_SymQMRData, 1, HYPRE_MEMORY_HOST);
66
67 /* set defaults */
68 (symqmr_data -> tol) = 1.0e-06;
69 (symqmr_data -> max_iter) = 1000;
70 (symqmr_data -> stop_crit) = 0; /* rel. residual norm */
71 (symqmr_data -> precond) = hypre_ParKrylovIdentity;
72 (symqmr_data -> precond_setup) = hypre_ParKrylovIdentitySetup;
73 (symqmr_data -> precond_data) = NULL;
74 (symqmr_data -> logging) = 0;
75 (symqmr_data -> r) = NULL;
76 (symqmr_data -> q) = NULL;
77 (symqmr_data -> u) = NULL;
78 (symqmr_data -> d) = NULL;
79 (symqmr_data -> t) = NULL;
80 (symqmr_data -> rq) = NULL;
81 (symqmr_data -> matvec_data) = NULL;
82 (symqmr_data -> norms) = NULL;
83 (symqmr_data -> log_file_name) = NULL;
84
85 return (void *) symqmr_data;
86 }
87
88 /*--------------------------------------------------------------------------
89 * hypre_SymQMRDestroy
90 *--------------------------------------------------------------------------*/
91
hypre_SymQMRDestroy(void * symqmr_vdata)92 int hypre_SymQMRDestroy( void *symqmr_vdata )
93 {
94 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
95 int ierr = 0;
96
97 if (symqmr_data)
98 {
99 if ((symqmr_data -> logging) > 0)
100 {
101 hypre_TFree(symqmr_data -> norms, HYPRE_MEMORY_HOST);
102 }
103
104 hypre_ParKrylovMatvecDestroy(symqmr_data -> matvec_data);
105
106 hypre_ParKrylovDestroyVector(symqmr_data -> r);
107 hypre_ParKrylovDestroyVector(symqmr_data -> q);
108 hypre_ParKrylovDestroyVector(symqmr_data -> u);
109 hypre_ParKrylovDestroyVector(symqmr_data -> d);
110 hypre_ParKrylovDestroyVector(symqmr_data -> t);
111 hypre_ParKrylovDestroyVector(symqmr_data -> rq);
112
113 hypre_TFree(symqmr_data, HYPRE_MEMORY_HOST);
114 }
115
116 return(ierr);
117 }
118
119 /*--------------------------------------------------------------------------
120 * hypre_SymQMRSetup
121 *--------------------------------------------------------------------------*/
122
hypre_SymQMRSetup(void * symqmr_vdata,void * A,void * b,void * x)123 int hypre_SymQMRSetup( void *symqmr_vdata, void *A, void *b, void *x )
124 {
125 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
126 int max_iter = (symqmr_data -> max_iter);
127 int (*precond_setup)(void*, void*, void*, void*) = (symqmr_data -> precond_setup);
128 void *precond_data = (symqmr_data -> precond_data);
129 int ierr = 0;
130
131 (symqmr_data -> A) = A;
132
133 /*--------------------------------------------------
134 * The arguments for NewVector are important to
135 * maintain consistency between the setup and
136 * compute phases of matvec and the preconditioner.
137 *--------------------------------------------------*/
138
139 if ((symqmr_data -> r) == NULL)
140 (symqmr_data -> r) = hypre_ParKrylovCreateVector(b);
141 if ((symqmr_data -> q) == NULL)
142 (symqmr_data -> q) = hypre_ParKrylovCreateVector(b);
143 if ((symqmr_data -> u) == NULL)
144 (symqmr_data -> u) = hypre_ParKrylovCreateVector(b);
145 if ((symqmr_data -> d) == NULL)
146 (symqmr_data -> d) = hypre_ParKrylovCreateVector(b);
147 if ((symqmr_data -> t) == NULL)
148 (symqmr_data -> t) = hypre_ParKrylovCreateVector(b);
149 if ((symqmr_data -> rq) == NULL)
150 (symqmr_data -> rq) = hypre_ParKrylovCreateVector(b);
151 if ((symqmr_data -> matvec_data) == NULL)
152 (symqmr_data -> matvec_data) = hypre_ParKrylovMatvecCreate(A, x);
153
154 ierr = precond_setup(precond_data, A, b, x);
155
156 /*-----------------------------------------------------
157 * Allocate space for log info
158 *-----------------------------------------------------*/
159
160 if ((symqmr_data -> logging) > 0)
161 {
162 if ((symqmr_data -> norms) == NULL)
163 (symqmr_data -> norms) = hypre_CTAlloc(double, max_iter + 1, HYPRE_MEMORY_HOST);
164 if ((symqmr_data -> log_file_name) == NULL)
165 (symqmr_data -> log_file_name) = (char*)"symqmr.out.log";
166 }
167
168 return ierr;
169 }
170
171 /*--------------------------------------------------------------------------
172 * hypre_SymQMRSolve
173 *-------------------------------------------------------------------------*/
174
hypre_SymQMRSolve(void * symqmr_vdata,void * A,void * b,void * x)175 int hypre_SymQMRSolve(void *symqmr_vdata, void *A, void *b, void *x)
176 {
177 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
178 int max_iter = (symqmr_data -> max_iter);
179 int stop_crit = (symqmr_data -> stop_crit);
180 double accuracy = (symqmr_data -> tol);
181 void *matvec_data = (symqmr_data -> matvec_data);
182
183 void *r = (symqmr_data -> r);
184 void *q = (symqmr_data -> q);
185 void *u = (symqmr_data -> u);
186 void *d = (symqmr_data -> d);
187 void *t = (symqmr_data -> t);
188 void *rq = (symqmr_data -> rq);
189 int (*precond)(void*, void*, void*, void*) = (symqmr_data -> precond);
190 int *precond_data = (int*)(symqmr_data -> precond_data);
191
192 /* logging variables */
193 int logging = (symqmr_data -> logging);
194 double *norms = (symqmr_data -> norms);
195
196 int ierr=0, my_id, num_procs, iter;
197 double theta, tau, rhom1, rho, dtmp, r_norm;
198 double thetam1, c, epsilon;
199 double sigma, alpha, beta;
200
201 hypre_ParKrylovCommInfo(A,&my_id,&num_procs);
202 if (logging > 0)
203 {
204 norms = (symqmr_data -> norms);
205 }
206
207 /* initialize work arrays */
208
209 hypre_ParKrylovCopyVector(b,r);
210
211 /* compute initial residual */
212
213 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, r);
214 r_norm = sqrt(hypre_ParKrylovInnerProd(r,r));
215 if (logging > 0)
216 {
217 norms[0] = r_norm;
218 if (my_id == 0)
219 printf("SymQMR : Initial L2 norm of residual = %e\n", r_norm);
220 }
221 iter = 0;
222 epsilon = accuracy * r_norm;
223
224 /* convergence criterion |r_i| <= accuracy , absolute residual norm*/
225 if (stop_crit) epsilon = accuracy;
226
227 while ( iter < max_iter && r_norm > epsilon )
228 {
229 if ( my_id == 0 && iter > 0 && logging ) printf("SymQMR restart... \n");
230
231 tau = r_norm;
232 precond(precond_data, A, r, q);
233 rho = hypre_ParKrylovInnerProd(r,q);
234 theta = 0.0;
235 hypre_ParKrylovClearVector(d);
236 hypre_ParKrylovCopyVector(r,rq);
237
238 while ( iter < max_iter && r_norm > epsilon )
239 {
240 iter++;
241
242 hypre_ParKrylovMatvec(matvec_data,1.0,A,q,0.0,t);
243 sigma = hypre_ParKrylovInnerProd(q,t);
244 if ( sigma == 0.0 )
245 {
246 printf("SymQMR ERROR : sigma = 0.0\n");
247 exit(1);
248 }
249 alpha = rho / sigma;
250 dtmp = - alpha;
251 hypre_ParKrylovAxpy(dtmp,t,r);
252 thetam1 = theta;
253 theta = sqrt(hypre_ParKrylovInnerProd(r,r)) / tau;
254 c = 1.0 / sqrt(1.0 + theta * theta );
255 tau = tau * theta * c;
256 dtmp = c * c * thetam1 * thetam1;
257 hypre_ParKrylovScaleVector(dtmp,d);
258 dtmp = c * c * alpha;
259 hypre_ParKrylovAxpy(dtmp,q,d);
260 dtmp = 1.0;
261 hypre_ParKrylovAxpy(dtmp,d,x);
262
263 precond(precond_data, A, r, u);
264 rhom1 = rho;
265 rho = hypre_ParKrylovInnerProd(r,u);
266 beta = rho / rhom1;
267 hypre_ParKrylovScaleVector(beta,q);
268 dtmp = 1.0;
269 hypre_ParKrylovAxpy(dtmp,u,q);
270
271 dtmp = 1.0 - c * c;
272 hypre_ParKrylovScaleVector(dtmp,rq);
273 dtmp = c * c;
274 hypre_ParKrylovAxpy(dtmp,r,rq);
275 r_norm = sqrt(hypre_ParKrylovInnerProd(rq,rq));
276 norms[iter] = r_norm;
277
278 if ( my_id == 0 && logging )
279 printf(" SymQMR : iteration %4d - residual norm = %e \n",
280 iter, r_norm);
281 }
282
283 /* compute true residual */
284
285 hypre_ParKrylovCopyVector(b,r);
286 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, r);
287 r_norm = sqrt(hypre_ParKrylovInnerProd(r,r));
288 }
289
290 (symqmr_data -> num_iterations) = iter;
291 (symqmr_data -> rel_residual_norm) = r_norm;
292
293 if (iter >= max_iter && r_norm > epsilon) ierr = 1;
294
295 return ierr;
296 }
297
298 /*--------------------------------------------------------------------------
299 * hypre_SymQMRSetTol
300 *--------------------------------------------------------------------------*/
301
hypre_SymQMRSetTol(void * symqmr_vdata,double tol)302 int hypre_SymQMRSetTol( void *symqmr_vdata, double tol )
303 {
304 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
305 int ierr = 0;
306
307 (symqmr_data -> tol) = tol;
308
309 return ierr;
310 }
311
312 /*--------------------------------------------------------------------------
313 * hypre_SymQMRSetMaxIter
314 *--------------------------------------------------------------------------*/
315
hypre_SymQMRSetMaxIter(void * symqmr_vdata,int max_iter)316 int hypre_SymQMRSetMaxIter( void *symqmr_vdata, int max_iter )
317 {
318 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
319 int ierr = 0;
320
321 (symqmr_data -> max_iter) = max_iter;
322
323 return ierr;
324 }
325
326 /*--------------------------------------------------------------------------
327 * hypre_SymQMRSetStopCrit
328 *--------------------------------------------------------------------------*/
329
hypre_SymQMRSetStopCrit(void * symqmr_vdata,double stop_crit)330 int hypre_SymQMRSetStopCrit( void *symqmr_vdata, double stop_crit )
331 {
332 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
333 int ierr = 0;
334
335 (symqmr_data -> stop_crit) = stop_crit;
336
337 return ierr;
338 }
339
340 /*--------------------------------------------------------------------------
341 * hypre_SymQMRSetPrecond
342 *--------------------------------------------------------------------------*/
343
hypre_SymQMRSetPrecond(void * symqmr_vdata,int (* precond)(void *,void *,void *,void *),int (* precond_setup)(void *,void *,void *,void *),void * precond_data)344 int hypre_SymQMRSetPrecond( void *symqmr_vdata, int (*precond)(void*,void*,void*,void*),
345 int (*precond_setup)(void*,void*,void*,void*), void *precond_data )
346 {
347 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
348 int ierr = 0;
349
350 (symqmr_data -> precond) = precond;
351 (symqmr_data -> precond_setup) = precond_setup;
352 (symqmr_data -> precond_data) = precond_data;
353
354 return ierr;
355 }
356
357 /*--------------------------------------------------------------------------
358 * hypre_SymQMRSetLogging
359 *--------------------------------------------------------------------------*/
360
hypre_SymQMRSetLogging(void * symqmr_vdata,int logging)361 int hypre_SymQMRSetLogging( void *symqmr_vdata, int logging)
362 {
363 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
364 int ierr = 0;
365
366 (symqmr_data -> logging) = logging;
367
368 return ierr;
369 }
370
371 /*--------------------------------------------------------------------------
372 * hypre_SymQMRGetNumIterations
373 *--------------------------------------------------------------------------*/
374
hypre_SymQMRGetNumIterations(void * symqmr_vdata,int * num_iterations)375 int hypre_SymQMRGetNumIterations(void *symqmr_vdata,int *num_iterations)
376 {
377 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
378 int ierr = 0;
379
380 *num_iterations = (symqmr_data -> num_iterations);
381
382 return ierr;
383 }
384
385 /*--------------------------------------------------------------------------
386 * hypre_SymQMRGetFinalRelativeResidualNorm
387 *--------------------------------------------------------------------------*/
388
hypre_SymQMRGetFinalRelativeResidualNorm(void * symqmr_vdata,double * relative_residual_norm)389 int hypre_SymQMRGetFinalRelativeResidualNorm( void *symqmr_vdata,
390 double *relative_residual_norm )
391 {
392 hypre_SymQMRData *symqmr_data = (hypre_SymQMRData*) symqmr_vdata;
393 int ierr = 0;
394
395 *relative_residual_norm = (symqmr_data -> rel_residual_norm);
396
397 return ierr;
398 }
399
400