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 * LSICG
11 *
12 *****************************************************************************/
13
14 #include "utilities/_hypre_utilities.h"
15 #include "HYPRE.h"
16 #include "parcsr_ls/_hypre_parcsr_ls.h"
17 #include "parcsr_mv/_hypre_parcsr_mv.h"
18 #include "seq_mv/seq_mv.h"
19
20 /*--------------------------------------------------------------------------
21 * hypre_LSICGData
22 *--------------------------------------------------------------------------*/
23
24 typedef struct
25 {
26 int max_iter;
27 int stop_crit;
28 double tol;
29 double rel_residual_norm;
30
31 void *A;
32 void *r;
33 void *ap;
34 void *p;
35 void *z;
36
37 void *matvec_data;
38
39 int (*precond)(void*, void*, void*, void*);
40 int (*precond_setup)(void*, void*, void*, void*);
41 void *precond_data;
42
43 int num_iterations;
44
45 int logging;
46
47 } hypre_LSICGData;
48
49 /*--------------------------------------------------------------------------
50 * hypre_LSICGCreate
51 *--------------------------------------------------------------------------*/
52
hypre_LSICGCreate()53 void *hypre_LSICGCreate( )
54 {
55 hypre_LSICGData *lsicg_data;
56
57 lsicg_data = hypre_CTAlloc(hypre_LSICGData, 1, HYPRE_MEMORY_HOST);
58
59 /* set defaults */
60 (lsicg_data -> tol) = 1.0e-06;
61 (lsicg_data -> max_iter) = 1000;
62 (lsicg_data -> stop_crit) = 0; /* rel. residual norm */
63 (lsicg_data -> precond) = hypre_ParKrylovIdentity;
64 (lsicg_data -> precond_setup) = hypre_ParKrylovIdentitySetup;
65 (lsicg_data -> precond_data) = NULL;
66 (lsicg_data -> logging) = 0;
67 (lsicg_data -> r) = NULL;
68 (lsicg_data -> p) = NULL;
69 (lsicg_data -> ap) = NULL;
70 (lsicg_data -> z) = NULL;
71 (lsicg_data -> matvec_data) = NULL;
72
73 return (void *) lsicg_data;
74 }
75
76 /*--------------------------------------------------------------------------
77 * hypre_LSICGDestroy
78 *--------------------------------------------------------------------------*/
79
hypre_LSICGDestroy(void * lsicg_vdata)80 int hypre_LSICGDestroy( void *lsicg_vdata )
81 {
82 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
83 int ierr = 0;
84
85 if (lsicg_data)
86 {
87 hypre_ParKrylovMatvecDestroy(lsicg_data -> matvec_data);
88 hypre_ParKrylovDestroyVector(lsicg_data -> r);
89 hypre_ParKrylovDestroyVector(lsicg_data -> p);
90 hypre_ParKrylovDestroyVector(lsicg_data -> ap);
91 hypre_ParKrylovDestroyVector(lsicg_data -> z);
92 hypre_TFree(lsicg_data, HYPRE_MEMORY_HOST);
93 }
94 return(ierr);
95 }
96
97 /*--------------------------------------------------------------------------
98 * hypre_LSICGSetup
99 *--------------------------------------------------------------------------*/
100
hypre_LSICGSetup(void * lsicg_vdata,void * A,void * b,void * x)101 int hypre_LSICGSetup( void *lsicg_vdata, void *A, void *b, void *x )
102 {
103 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
104 int (*precond_setup)(void*, void*, void*, void*) = (lsicg_data -> precond_setup);
105 void *precond_data = (lsicg_data -> precond_data);
106 int ierr = 0;
107
108 (lsicg_data -> A) = A;
109
110 /*--------------------------------------------------
111 * The arguments for NewVector are important to
112 * maintain consistency between the setup and
113 * compute phases of matvec and the preconditioner.
114 *--------------------------------------------------*/
115
116 if ((lsicg_data -> r) == NULL)
117 (lsicg_data -> r) = hypre_ParKrylovCreateVector(b);
118 if ((lsicg_data -> p) == NULL)
119 (lsicg_data -> p) = hypre_ParKrylovCreateVector(b);
120 if ((lsicg_data -> z) == NULL)
121 (lsicg_data -> z) = hypre_ParKrylovCreateVector(b);
122 if ((lsicg_data -> ap) == NULL)
123 (lsicg_data -> ap) = hypre_ParKrylovCreateVector(b);
124 if ((lsicg_data -> matvec_data) == NULL)
125 (lsicg_data -> matvec_data) = hypre_ParKrylovMatvecCreate(A, x);
126
127 ierr = precond_setup(precond_data, A, b, x);
128
129 return ierr;
130 }
131
132 /*--------------------------------------------------------------------------
133 * hypre_LSICGSolve
134 *-------------------------------------------------------------------------*/
135
hypre_LSICGSolve(void * lsicg_vdata,void * A,void * b,void * x)136 int hypre_LSICGSolve(void *lsicg_vdata, void *A, void *b, void *x)
137 {
138 int ierr=0, mypid, nprocs, iter, converged=0;
139 double rhom1, rho, r_norm, b_norm, epsilon;
140 double sigma, alpha, beta, dArray[2], dArray2[2];
141 hypre_Vector *r_local, *z_local;
142 MPI_Comm comm;
143
144 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
145 int max_iter = (lsicg_data -> max_iter);
146 int stop_crit = (lsicg_data -> stop_crit);
147 double accuracy = (lsicg_data -> tol);
148 void *matvec_data = (lsicg_data -> matvec_data);
149 void *r = (lsicg_data -> r);
150 void *p = (lsicg_data -> p);
151 void *z = (lsicg_data -> z);
152 void *ap = (lsicg_data -> ap);
153 int (*precond)(void*, void*, void*, void*) = (lsicg_data -> precond);
154 int *precond_data = (int*)(lsicg_data -> precond_data);
155 int logging = (lsicg_data -> logging);
156
157 /* compute initial residual */
158
159 r_local = hypre_ParVectorLocalVector((hypre_ParVector *) r);
160 z_local = hypre_ParVectorLocalVector((hypre_ParVector *) z);
161 comm = hypre_ParCSRMatrixComm((hypre_ParCSRMatrix *) A);
162 hypre_ParKrylovCommInfo(A,&mypid,&nprocs);
163 hypre_ParKrylovCopyVector(b,r);
164 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, r);
165 r_norm = sqrt(hypre_ParKrylovInnerProd(r,r));
166 b_norm = sqrt(hypre_ParKrylovInnerProd(b,b));
167 if (logging > 0)
168 {
169 if (mypid == 0)
170 {
171 printf("LSICG : L2 norm of b = %e\n", b_norm);
172 if (b_norm == 0.0)
173 printf("Rel_resid_norm actually contains the residual norm\n");
174 printf("LSICG : Initial L2 norm of residual = %e\n", r_norm);
175 }
176 }
177
178 /* set convergence criterion */
179
180 if (b_norm > 0.0) epsilon = accuracy * b_norm;
181 else epsilon = accuracy * r_norm;
182 if ( stop_crit ) epsilon = accuracy;
183
184 iter = 0;
185 hypre_ParKrylovClearVector(p);
186
187 while ( converged == 0 )
188 {
189 while ( r_norm > epsilon && iter < max_iter )
190 {
191 iter++;
192 if ( iter == 1 )
193 {
194 precond(precond_data, A, r, z);
195 rhom1 = rho;
196 rho = hypre_ParKrylovInnerProd(r,z);
197 beta = 0.0;
198 }
199 else beta = rho / rhom1;
200 hypre_ParKrylovScaleVector( beta, p );
201 hypre_ParKrylovAxpy(1.0e0, z, p);
202 hypre_ParKrylovMatvec(matvec_data,1.0e0,A,p,0.0,ap);
203 sigma = hypre_ParKrylovInnerProd(p,ap);
204 alpha = rho / sigma;
205 if ( sigma == 0.0 )
206 {
207 printf("HYPRE::LSICG ERROR - sigma = 0.0.\n");
208 ierr = 2;
209 return ierr;
210 }
211 hypre_ParKrylovAxpy(alpha, p, x);
212 hypre_ParKrylovAxpy(-alpha, ap, r);
213 dArray[0] = hypre_SeqVectorInnerProd( r_local, r_local );
214 precond(precond_data, A, r, z);
215 rhom1 = rho;
216 dArray[1] = hypre_SeqVectorInnerProd( r_local, z_local );
217 MPI_Allreduce(dArray, dArray2, 2, MPI_DOUBLE, MPI_SUM, comm);
218 rho = dArray2[1];
219 r_norm = sqrt( dArray2[0] );
220 if ( iter % 1 == 0 && mypid == 0 )
221 printf("LSICG : iteration %d - residual norm = %e (%e)\n",
222 iter, r_norm, epsilon);
223 }
224 hypre_ParKrylovCopyVector(b,r);
225 hypre_ParKrylovMatvec(matvec_data,-1.0, A, x, 1.0, r);
226 r_norm = sqrt(hypre_ParKrylovInnerProd(r,r));
227 if ( logging >= 1 && mypid == 0 )
228 printf("LSICG actual residual norm = %e \n",r_norm);
229 if ( r_norm < epsilon || iter >= max_iter ) converged = 1;
230 }
231 if ( iter >= max_iter ) ierr = 1;
232 lsicg_data->rel_residual_norm = r_norm;
233 lsicg_data->num_iterations = iter;
234 if ( logging >= 1 && mypid == 0 )
235 printf("LSICG : total number of iterations = %d \n",iter);
236
237 return ierr;
238 }
239
240 /*--------------------------------------------------------------------------
241 * hypre_LSICGSetTol
242 *--------------------------------------------------------------------------*/
243
hypre_LSICGSetTol(void * lsicg_vdata,double tol)244 int hypre_LSICGSetTol( void *lsicg_vdata, double tol )
245 {
246 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
247 (lsicg_data -> tol) = tol;
248 return 0;
249 }
250
251 /*--------------------------------------------------------------------------
252 * hypre_LSICGSetMaxIter
253 *--------------------------------------------------------------------------*/
254
hypre_LSICGSetMaxIter(void * lsicg_vdata,int max_iter)255 int hypre_LSICGSetMaxIter( void *lsicg_vdata, int max_iter )
256 {
257 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
258 (lsicg_data -> max_iter) = max_iter;
259 return 0;
260 }
261
262 /*--------------------------------------------------------------------------
263 * hypre_LSICGSetStopCrit
264 *--------------------------------------------------------------------------*/
265
hypre_LSICGSetStopCrit(void * lsicg_vdata,double stop_crit)266 int hypre_LSICGSetStopCrit( void *lsicg_vdata, double stop_crit )
267 {
268 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
269 (lsicg_data -> stop_crit) = stop_crit;
270 return 0;
271 }
272
273 /*--------------------------------------------------------------------------
274 * hypre_LSICGSetPrecond
275 *--------------------------------------------------------------------------*/
276
hypre_LSICGSetPrecond(void * lsicg_vdata,int (* precond)(void *,void *,void *,void *),int (* precond_setup)(void *,void *,void *,void *),void * precond_data)277 int hypre_LSICGSetPrecond( void *lsicg_vdata, int (*precond)(void*,void*,void*,void*),
278 int (*precond_setup)(void*,void*,void*,void*), void *precond_data )
279 {
280 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
281 (lsicg_data -> precond) = precond;
282 (lsicg_data -> precond_setup) = precond_setup;
283 (lsicg_data -> precond_data) = precond_data;
284 return 0;
285 }
286
287 /*--------------------------------------------------------------------------
288 * hypre_LSICGSetLogging
289 *--------------------------------------------------------------------------*/
290
hypre_LSICGSetLogging(void * lsicg_vdata,int logging)291 int hypre_LSICGSetLogging( void *lsicg_vdata, int logging)
292 {
293 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
294 (lsicg_data -> logging) = logging;
295 return 0;
296 }
297
298 /*--------------------------------------------------------------------------
299 * hypre_LSICGGetNumIterations
300 *--------------------------------------------------------------------------*/
301
hypre_LSICGGetNumIterations(void * lsicg_vdata,int * num_iterations)302 int hypre_LSICGGetNumIterations(void *lsicg_vdata,int *num_iterations)
303 {
304 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
305 *num_iterations = (lsicg_data -> num_iterations);
306 return 0;
307 }
308
309 /*--------------------------------------------------------------------------
310 * hypre_LSICGGetFinalRelativeResidualNorm
311 *--------------------------------------------------------------------------*/
312
hypre_LSICGGetFinalRelativeResidualNorm(void * lsicg_vdata,double * relative_residual_norm)313 int hypre_LSICGGetFinalRelativeResidualNorm(void *lsicg_vdata,
314 double *relative_residual_norm)
315 {
316 hypre_LSICGData *lsicg_data = (hypre_LSICGData *) lsicg_vdata;
317 *relative_residual_norm = (lsicg_data -> rel_residual_norm);
318 return 0;
319 }
320
321