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 * AMG solve routine
11 *
12 *****************************************************************************/
13
14 #include "_hypre_parcsr_ls.h"
15 #include "par_amg.h"
16
17 /*--------------------------------------------------------------------
18 * hypre_BoomerAMGSolve
19 *--------------------------------------------------------------------*/
20
21 HYPRE_Int
hypre_BoomerAMGSolve(void * amg_vdata,hypre_ParCSRMatrix * A,hypre_ParVector * f,hypre_ParVector * u)22 hypre_BoomerAMGSolve( void *amg_vdata,
23 hypre_ParCSRMatrix *A,
24 hypre_ParVector *f,
25 hypre_ParVector *u )
26 {
27 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
28
29 hypre_ParAMGData *amg_data = (hypre_ParAMGData*) amg_vdata;
30
31 /* Data Structure variables */
32
33 HYPRE_Int amg_print_level;
34 HYPRE_Int amg_logging;
35 HYPRE_Int cycle_count;
36 HYPRE_Int num_levels;
37 /* HYPRE_Int num_unknowns; */
38 HYPRE_Int converge_type;
39 HYPRE_Real tol;
40
41 HYPRE_Int block_mode;
42
43
44 hypre_ParCSRMatrix **A_array;
45 hypre_ParVector **F_array;
46 hypre_ParVector **U_array;
47
48 hypre_ParCSRBlockMatrix **A_block_array;
49
50
51 /* Local variables */
52
53 HYPRE_Int j;
54 HYPRE_Int Solve_err_flag;
55 HYPRE_Int min_iter;
56 HYPRE_Int max_iter;
57 HYPRE_Int num_procs, my_id;
58 HYPRE_Int additive;
59 HYPRE_Int mult_additive;
60 HYPRE_Int simple;
61
62 HYPRE_Real alpha = 1.0;
63 HYPRE_Real beta = -1.0;
64 HYPRE_Real cycle_op_count;
65 HYPRE_Real total_coeffs;
66 HYPRE_Real total_variables;
67 HYPRE_Real *num_coeffs;
68 HYPRE_Real *num_variables;
69 HYPRE_Real cycle_cmplxty = 0.0;
70 HYPRE_Real operat_cmplxty;
71 HYPRE_Real grid_cmplxty;
72 HYPRE_Real conv_factor = 0.0;
73 HYPRE_Real resid_nrm = 1.0;
74 HYPRE_Real resid_nrm_init = 0.0;
75 HYPRE_Real relative_resid;
76 HYPRE_Real rhs_norm = 0.0;
77 HYPRE_Real old_resid;
78 HYPRE_Real ieee_check = 0.;
79
80 hypre_ParVector *Vtemp;
81 hypre_ParVector *Residual;
82
83 HYPRE_ANNOTATE_FUNC_BEGIN;
84 hypre_MPI_Comm_size(comm, &num_procs);
85 hypre_MPI_Comm_rank(comm,&my_id);
86
87 amg_print_level = hypre_ParAMGDataPrintLevel(amg_data);
88 amg_logging = hypre_ParAMGDataLogging(amg_data);
89 if ( amg_logging > 1 )
90 Residual = hypre_ParAMGDataResidual(amg_data);
91 /* num_unknowns = hypre_ParAMGDataNumUnknowns(amg_data); */
92 num_levels = hypre_ParAMGDataNumLevels(amg_data);
93 A_array = hypre_ParAMGDataAArray(amg_data);
94 F_array = hypre_ParAMGDataFArray(amg_data);
95 U_array = hypre_ParAMGDataUArray(amg_data);
96
97 converge_type = hypre_ParAMGDataConvergeType(amg_data);
98 tol = hypre_ParAMGDataTol(amg_data);
99 min_iter = hypre_ParAMGDataMinIter(amg_data);
100 max_iter = hypre_ParAMGDataMaxIter(amg_data);
101 additive = hypre_ParAMGDataAdditive(amg_data);
102 simple = hypre_ParAMGDataSimple(amg_data);
103 mult_additive = hypre_ParAMGDataMultAdditive(amg_data);
104
105 A_array[0] = A;
106 F_array[0] = f;
107 U_array[0] = u;
108
109 block_mode = hypre_ParAMGDataBlockMode(amg_data);
110
111 A_block_array = hypre_ParAMGDataABlockArray(amg_data);
112
113
114 /* Vtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
115 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
116 hypre_ParCSRMatrixRowStarts(A_array[0]));
117 hypre_ParVectorInitialize(Vtemp);
118 hypre_ParAMGDataVtemp(amg_data) = Vtemp;
119 */
120 Vtemp = hypre_ParAMGDataVtemp(amg_data);
121
122 /*-----------------------------------------------------------------------
123 * Write the solver parameters
124 *-----------------------------------------------------------------------*/
125
126 if (my_id == 0 && amg_print_level > 1)
127 {
128 hypre_BoomerAMGWriteSolverParams(amg_data);
129 }
130
131 /*-----------------------------------------------------------------------
132 * Initialize the solver error flag and assorted bookkeeping variables
133 *-----------------------------------------------------------------------*/
134
135 Solve_err_flag = 0;
136
137 total_coeffs = 0;
138 total_variables = 0;
139 cycle_count = 0;
140 operat_cmplxty = 0;
141 grid_cmplxty = 0;
142
143 /*-----------------------------------------------------------------------
144 * write some initial info
145 *-----------------------------------------------------------------------*/
146
147 if (my_id == 0 && amg_print_level > 1 && tol > 0.)
148 hypre_printf("\n\nAMG SOLUTION INFO:\n");
149
150
151 /*-----------------------------------------------------------------------
152 * Compute initial fine-grid residual and print
153 *-----------------------------------------------------------------------*/
154
155 if (amg_print_level > 1 || amg_logging > 1 || tol > 0.)
156 {
157 if ( amg_logging > 1 )
158 {
159 hypre_ParVectorCopy(F_array[0], Residual );
160 if (tol > 0)
161 {
162 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Residual );
163 }
164 resid_nrm = sqrt(hypre_ParVectorInnerProd( Residual, Residual ));
165 }
166 else
167 {
168 hypre_ParVectorCopy(F_array[0], Vtemp);
169 if (tol > 0)
170 {
171 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Vtemp);
172 }
173 resid_nrm = sqrt(hypre_ParVectorInnerProd(Vtemp, Vtemp));
174 }
175
176 /* Since it is does not diminish performance, attempt to return an error flag
177 and notify users when they supply bad input. */
178 if (resid_nrm != 0.)
179 {
180 ieee_check = resid_nrm/resid_nrm; /* INF -> NaN conversion */
181 }
182
183 if (ieee_check != ieee_check)
184 {
185 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
186 for ieee_check self-equality works on all IEEE-compliant compilers/
187 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
188 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
189 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
190 if (amg_print_level > 0)
191 {
192 hypre_printf("\n\nERROR detected by Hypre ... BEGIN\n");
193 hypre_printf("ERROR -- hypre_BoomerAMGSolve: INFs and/or NaNs detected in input.\n");
194 hypre_printf("User probably placed non-numerics in supplied A, x_0, or b.\n");
195 hypre_printf("ERROR detected by Hypre ... END\n\n\n");
196 }
197 hypre_error(HYPRE_ERROR_GENERIC);
198 HYPRE_ANNOTATE_FUNC_END;
199
200 return hypre_error_flag;
201 }
202
203 /* r0 */
204 resid_nrm_init = resid_nrm;
205
206 if (0 == converge_type)
207 {
208 rhs_norm = sqrt(hypre_ParVectorInnerProd(f, f));
209 if (rhs_norm)
210 {
211 relative_resid = resid_nrm_init / rhs_norm;
212 }
213 else
214 {
215 relative_resid = resid_nrm_init;
216 }
217 }
218 else
219 {
220 /* converge_type != 0, test convergence with ||r|| / ||r0|| */
221 relative_resid = 1.0;
222 }
223 }
224 else
225 {
226 relative_resid = 1.;
227 }
228
229 if (my_id == 0 && amg_print_level > 1)
230 {
231 hypre_printf(" relative\n");
232 hypre_printf(" residual factor residual\n");
233 hypre_printf(" -------- ------ --------\n");
234 hypre_printf(" Initial %e %e\n", resid_nrm_init,
235 relative_resid);
236 }
237
238 /*-----------------------------------------------------------------------
239 * Main V-cycle loop
240 *-----------------------------------------------------------------------*/
241
242 while ( (relative_resid >= tol || cycle_count < min_iter) && cycle_count < max_iter )
243 {
244 hypre_ParAMGDataCycleOpCount(amg_data) = 0;
245 /* Op count only needed for one cycle */
246 if ( (additive < 0 || additive >= num_levels) &&
247 (mult_additive < 0 || mult_additive >= num_levels) &&
248 (simple < 0 || simple >= num_levels) )
249 {
250 hypre_BoomerAMGCycle(amg_data, F_array, U_array);
251 }
252 else
253 {
254 hypre_BoomerAMGAdditiveCycle(amg_data);
255 }
256
257 /*---------------------------------------------------------------
258 * Compute fine-grid residual and residual norm
259 *----------------------------------------------------------------*/
260
261 if (amg_print_level > 1 || amg_logging > 1 || tol > 0.)
262 {
263 old_resid = resid_nrm;
264
265 if ( amg_logging > 1 )
266 {
267 hypre_ParCSRMatrixMatvecOutOfPlace(alpha, A_array[0], U_array[0], beta, F_array[0], Residual );
268 resid_nrm = sqrt(hypre_ParVectorInnerProd( Residual, Residual ));
269 }
270 else
271 {
272 hypre_ParCSRMatrixMatvecOutOfPlace(alpha, A_array[0], U_array[0], beta, F_array[0], Vtemp);
273 resid_nrm = sqrt(hypre_ParVectorInnerProd(Vtemp, Vtemp));
274 }
275
276 if (old_resid)
277 {
278 conv_factor = resid_nrm / old_resid;
279 }
280 else
281 {
282 conv_factor = resid_nrm;
283 }
284
285 if (0 == converge_type)
286 {
287 if (rhs_norm)
288 {
289 relative_resid = resid_nrm / rhs_norm;
290 }
291 else
292 {
293 relative_resid = resid_nrm;
294 }
295 }
296 else
297 {
298 relative_resid = resid_nrm / resid_nrm_init;
299 }
300
301 hypre_ParAMGDataRelativeResidualNorm(amg_data) = relative_resid;
302 }
303
304 ++cycle_count;
305
306 hypre_ParAMGDataNumIterations(amg_data) = cycle_count;
307 #ifdef CUMNUMIT
308 ++hypre_ParAMGDataCumNumIterations(amg_data);
309 #endif
310
311 if (my_id == 0 && amg_print_level > 1)
312 {
313 hypre_printf(" Cycle %2d %e %f %e \n", cycle_count,
314 resid_nrm, conv_factor, relative_resid);
315 }
316 }
317
318 if (cycle_count == max_iter && tol > 0.)
319 {
320 Solve_err_flag = 1;
321 hypre_error(HYPRE_ERROR_CONV);
322 }
323
324 /*-----------------------------------------------------------------------
325 * Compute closing statistics
326 *-----------------------------------------------------------------------*/
327
328 if (cycle_count > 0 && resid_nrm_init)
329 conv_factor = pow((resid_nrm/resid_nrm_init),(1.0/(HYPRE_Real) cycle_count));
330 else
331 conv_factor = 1.;
332
333 if (amg_print_level > 1)
334 {
335 num_coeffs = hypre_CTAlloc(HYPRE_Real, num_levels, HYPRE_MEMORY_HOST);
336 num_variables = hypre_CTAlloc(HYPRE_Real, num_levels, HYPRE_MEMORY_HOST);
337 num_coeffs[0] = hypre_ParCSRMatrixDNumNonzeros(A);
338 num_variables[0] = hypre_ParCSRMatrixGlobalNumRows(A);
339
340 if (block_mode)
341 {
342 for (j = 1; j < num_levels; j++)
343 {
344 num_coeffs[j] = (HYPRE_Real) hypre_ParCSRBlockMatrixNumNonzeros(A_block_array[j]);
345 num_variables[j] = (HYPRE_Real) hypre_ParCSRBlockMatrixGlobalNumRows(A_block_array[j]);
346 }
347 num_coeffs[0] = hypre_ParCSRBlockMatrixDNumNonzeros(A_block_array[0]);
348 num_variables[0] = hypre_ParCSRBlockMatrixGlobalNumRows(A_block_array[0]);
349
350 }
351 else
352 {
353 for (j = 1; j < num_levels; j++)
354 {
355 num_coeffs[j] = (HYPRE_Real) hypre_ParCSRMatrixNumNonzeros(A_array[j]);
356 num_variables[j] = (HYPRE_Real) hypre_ParCSRMatrixGlobalNumRows(A_array[j]);
357 }
358 }
359
360
361 for (j=0;j<hypre_ParAMGDataNumLevels(amg_data);j++)
362 {
363 total_coeffs += num_coeffs[j];
364 total_variables += num_variables[j];
365 }
366
367 cycle_op_count = hypre_ParAMGDataCycleOpCount(amg_data);
368
369 if (num_variables[0])
370 grid_cmplxty = total_variables / num_variables[0];
371 if (num_coeffs[0])
372 {
373 operat_cmplxty = total_coeffs / num_coeffs[0];
374 cycle_cmplxty = cycle_op_count / num_coeffs[0];
375 }
376
377 if (my_id == 0)
378 {
379 if (Solve_err_flag == 1)
380 {
381 hypre_printf("\n\n==============================================");
382 hypre_printf("\n NOTE: Convergence tolerance was not achieved\n");
383 hypre_printf(" within the allowed %d V-cycles\n",max_iter);
384 hypre_printf("==============================================");
385 }
386 hypre_printf("\n\n Average Convergence Factor = %f",conv_factor);
387 hypre_printf("\n\n Complexity: grid = %f\n",grid_cmplxty);
388 hypre_printf(" operator = %f\n",operat_cmplxty);
389 hypre_printf(" cycle = %f\n\n\n\n",cycle_cmplxty);
390 }
391
392 hypre_TFree(num_coeffs, HYPRE_MEMORY_HOST);
393 hypre_TFree(num_variables, HYPRE_MEMORY_HOST);
394 }
395 HYPRE_ANNOTATE_FUNC_END;
396
397 return hypre_error_flag;
398 }
399