source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_amg_solve.c

main
Last change on this file was ea777aa, checked in by Alex Wilton <awilton@…>, 3 years ago

Moved examples, include, build_default.properties, common.xml, and README out from dev.civl.com into the root of the repo.

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 10.4 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13#define MPIP_SOLVE_ON 0
14
15
16
17/******************************************************************************
18 *
19 * AMG solve routine
20 *
21 *****************************************************************************/
22
23#include "headers.h"
24#include "par_amg.h"
25
26/*--------------------------------------------------------------------
27 * hypre_BoomerAMGSolve
28 *--------------------------------------------------------------------*/
29
30int
31hypre_BoomerAMGSolve( void *amg_vdata,
32 hypre_ParCSRMatrix *A,
33 hypre_ParVector *f,
34 hypre_ParVector *u )
35{
36
37 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
38
39 hypre_ParAMGData *amg_data = amg_vdata;
40
41 /* Data Structure variables */
42
43 int amg_print_level;
44 int amg_logging;
45 int cycle_count;
46 int num_levels;
47 /* int num_unknowns; */
48 double tol;
49
50
51 hypre_ParCSRMatrix **A_array;
52 hypre_ParVector **F_array;
53 hypre_ParVector **U_array;
54
55 /*hypre_ParCSRBlockMatrix **A_block_array;*/
56
57
58 /* Local variables */
59
60 int j;
61 int Solve_err_flag;
62 int min_iter;
63 int max_iter;
64 int num_procs, my_id;
65
66 double alpha = 1.0;
67 double beta = -1.0;
68 double cycle_op_count;
69 double total_coeffs;
70 double total_variables;
71 double *num_coeffs;
72 double *num_variables;
73 double cycle_cmplxty;
74 double operat_cmplxty;
75 double grid_cmplxty;
76 double conv_factor;
77 double resid_nrm;
78 double resid_nrm_init;
79 double relative_resid;
80 double rhs_norm;
81 double old_resid;
82 double ieee_check = 0.;
83
84 hypre_ParVector *Vtemp;
85 hypre_ParVector *Residual;
86
87 MPI_Comm_size(comm, &num_procs);
88 MPI_Comm_rank(comm,&my_id);
89
90 amg_print_level = hypre_ParAMGDataPrintLevel(amg_data);
91 amg_logging = hypre_ParAMGDataLogging(amg_data);
92 if ( amg_logging > 1 )
93 Residual = hypre_ParAMGDataResidual(amg_data);
94 /* num_unknowns = hypre_ParAMGDataNumUnknowns(amg_data); */
95 num_levels = hypre_ParAMGDataNumLevels(amg_data);
96 A_array = hypre_ParAMGDataAArray(amg_data);
97 F_array = hypre_ParAMGDataFArray(amg_data);
98 U_array = hypre_ParAMGDataUArray(amg_data);
99
100 tol = hypre_ParAMGDataTol(amg_data);
101 min_iter = hypre_ParAMGDataMinIter(amg_data);
102 max_iter = hypre_ParAMGDataMaxIter(amg_data);
103
104 num_coeffs = hypre_CTAlloc(double, num_levels);
105 num_variables = hypre_CTAlloc(double, num_levels);
106 num_coeffs[0] = hypre_ParCSRMatrixDNumNonzeros(A);
107 num_variables[0] = hypre_ParCSRMatrixGlobalNumRows(A);
108
109 A_array[0] = A;
110 F_array[0] = f;
111 U_array[0] = u;
112
113 Vtemp = hypre_ParAMGDataVtemp(amg_data);
114
115
116 for (j = 1; j < num_levels; j++)
117 {
118 num_coeffs[j] = (double) hypre_ParCSRMatrixNumNonzeros(A_array[j]);
119 num_variables[j] = (double) hypre_ParCSRMatrixGlobalNumRows(A_array[j]);
120 }
121
122
123 /*-----------------------------------------------------------------------
124 * Write the solver parameters
125 *-----------------------------------------------------------------------*/
126
127
128 if (my_id == 0 && amg_print_level > 1)
129 hypre_BoomerAMGWriteSolverParams(amg_data);
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 printf("\n\nAMG SOLUTION INFO:\n");
149
150
151 /*-----------------------------------------------------------------------
152 * Compute initial fine-grid residual and print
153 *-----------------------------------------------------------------------*/
154
155 if (tol >= 0.)
156 {
157 if ( amg_logging > 1 ) {
158 hypre_ParVectorCopy(F_array[0], Residual );
159 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Residual );
160 resid_nrm = sqrt(hypre_ParVectorInnerProd( Residual, Residual ));
161 }
162 else {
163 hypre_ParVectorCopy(F_array[0], Vtemp);
164 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Vtemp);
165 resid_nrm = sqrt(hypre_ParVectorInnerProd(Vtemp, Vtemp));
166 }
167
168 /* Since it is does not diminish performance, attempt to return an error flag
169 and notify users when they supply bad input. */
170 if (resid_nrm != 0.) ieee_check = resid_nrm/resid_nrm; /* INF -> NaN conversion */
171 if (ieee_check != ieee_check)
172 {
173 /* ...INFs or NaNs in input can make ieee_check a NaN. This test
174 for ieee_check self-equality works on all IEEE-compliant compilers/
175 machines, c.f. page 8 of "Lecture Notes on the Status of IEEE 754"
176 by W. Kahan, May 31, 1996. Currently (July 2002) this paper may be
177 found at http://HTTP.CS.Berkeley.EDU/~wkahan/ieee754status/IEEE754.PDF */
178 if (amg_print_level > 0)
179 {
180 printf("\n\nERROR detected by Hypre ... BEGIN\n");
181 printf("ERROR -- hypre_BoomerAMGSolve: INFs and/or NaNs detected in input.\n");
182 printf("User probably placed non-numerics in supplied A, x_0, or b.\n");
183 printf("ERROR detected by Hypre ... END\n\n\n");
184 }
185 hypre_error(HYPRE_ERROR_GENERIC);
186 return hypre_error_flag;
187 }
188
189 resid_nrm_init = resid_nrm;
190 rhs_norm = sqrt(hypre_ParVectorInnerProd(f, f));
191 if (rhs_norm)
192 {
193 relative_resid = resid_nrm_init / rhs_norm;
194 }
195 else
196 {
197 relative_resid = resid_nrm_init;
198 }
199 }
200 else
201 {
202 relative_resid = 1.;
203 }
204
205 if (my_id == 0 && amg_print_level > 1 && tol >= 0.)
206 {
207 printf(" relative\n");
208 printf(" residual factor residual\n");
209 printf(" -------- ------ --------\n");
210 printf(" Initial %e %e\n",resid_nrm_init,
211 relative_resid);
212 }
213
214 /*-----------------------------------------------------------------------
215 * Main V-cycle loop
216 *-----------------------------------------------------------------------*/
217
218 while ((relative_resid >= tol || cycle_count < min_iter)
219 && cycle_count < max_iter)
220 {
221 hypre_ParAMGDataCycleOpCount(amg_data) = 0;
222 /* Op count only needed for one cycle */
223
224#if MPIP_SOLVE_ON
225 MPI_Pcontrol(2);
226 MPI_Pcontrol(1);
227#endif
228
229
230 hypre_BoomerAMGCycle(amg_data, F_array, U_array);
231
232#if MPIP_SOLVE_ON
233 MPI_Pcontrol(3);
234 MPI_Pcontrol(0);
235#endif
236
237 /*---------------------------------------------------------------
238 * Compute fine-grid residual and residual norm
239 *----------------------------------------------------------------*/
240
241 if (tol >= 0.)
242 {
243 old_resid = resid_nrm;
244
245 if ( amg_logging > 1 ) {
246 hypre_ParVectorCopy(F_array[0], Residual);
247 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Residual );
248 resid_nrm = sqrt(hypre_ParVectorInnerProd( Residual, Residual ));
249 }
250 else {
251 hypre_ParVectorCopy(F_array[0], Vtemp);
252 hypre_ParCSRMatrixMatvec(alpha, A_array[0], U_array[0], beta, Vtemp);
253 resid_nrm = sqrt(hypre_ParVectorInnerProd(Vtemp, Vtemp));
254 }
255
256 conv_factor = resid_nrm / old_resid;
257 if (rhs_norm)
258 {
259 relative_resid = resid_nrm / rhs_norm;
260 }
261 else
262 {
263 relative_resid = resid_nrm;
264 }
265
266 hypre_ParAMGDataRelativeResidualNorm(amg_data) = relative_resid;
267 }
268
269 ++cycle_count;
270
271 hypre_ParAMGDataNumIterations(amg_data) = cycle_count;
272#ifdef CUMNUMIT
273 ++hypre_ParAMGDataCumNumIterations(amg_data);
274#endif
275
276 if (my_id == 0 && amg_print_level > 1 && tol >= 0.)
277 {
278 printf(" Cycle %2d %e %f %e \n", cycle_count,
279 resid_nrm, conv_factor, relative_resid);
280 }
281 }
282
283 if (cycle_count == max_iter && tol > 0.)
284 {
285 Solve_err_flag = 1;
286 hypre_error(HYPRE_ERROR_CONV);
287 }
288
289 /*-----------------------------------------------------------------------
290 * Compute closing statistics
291 *-----------------------------------------------------------------------*/
292
293 if (cycle_count > 0 && tol >= 0.)
294 conv_factor = pow((resid_nrm/resid_nrm_init),(1.0/(double) cycle_count));
295 else
296 conv_factor = 1.;
297
298
299 for (j=0;j<hypre_ParAMGDataNumLevels(amg_data);j++)
300 {
301 total_coeffs += num_coeffs[j];
302 total_variables += num_variables[j];
303 }
304
305 cycle_op_count = hypre_ParAMGDataCycleOpCount(amg_data);
306
307 if (num_variables[0])
308 grid_cmplxty = total_variables / num_variables[0];
309 if (num_coeffs[0])
310 {
311 operat_cmplxty = total_coeffs / num_coeffs[0];
312 cycle_cmplxty = cycle_op_count / num_coeffs[0];
313 }
314
315 if (my_id == 0 && amg_print_level > 1)
316 {
317 if (Solve_err_flag == 1)
318 {
319 printf("\n\n==============================================");
320 printf("\n NOTE: Convergence tolerance was not achieved\n");
321 printf(" within the allowed %d V-cycles\n",max_iter);
322 printf("==============================================");
323 }
324 if (tol >= 0.)
325 printf("\n\n Average Convergence Factor = %f",conv_factor);
326 printf("\n\n Complexity: grid = %f\n",grid_cmplxty);
327 printf(" operator = %f\n",operat_cmplxty);
328 printf(" cycle = %f\n\n\n\n",cycle_cmplxty);
329 }
330
331 hypre_TFree(num_coeffs);
332 hypre_TFree(num_variables);
333
334 return hypre_error_flag;
335}
336
Note: See TracBrowser for help on using the repository browser.