source: CIVL/examples/mpi-omp/AMG2013/parcsr_ls/par_cycle.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: 17.3 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_CYCLE_ON 0
14
15
16
17/******************************************************************************
18 *
19 * ParAMG cycling routine
20 *
21 *****************************************************************************/
22
23#include "headers.h"
24#include "par_amg.h"
25/*--------------------------------------------------------------------------
26 * hypre_BoomerAMGCycle
27 *--------------------------------------------------------------------------*/
28
29int
30hypre_BoomerAMGCycle( void *amg_vdata,
31 hypre_ParVector **F_array,
32 hypre_ParVector **U_array )
33{
34 hypre_ParAMGData *amg_data = amg_vdata;
35
36 hypre_SeqAMGData *seq_data = hypre_ParAMGDataSeqData(amg_data);
37
38 MPI_Comm comm;
39 HYPRE_Solver *smoother;
40 /* Data Structure variables */
41
42 hypre_ParCSRMatrix **A_array;
43 hypre_ParCSRMatrix **P_array;
44 hypre_ParCSRMatrix **R_array;
45 /*hypre_ParVector *Utemp;*/
46 hypre_ParVector *Vtemp;
47 hypre_ParVector *Rtemp;
48 hypre_ParVector *Ptemp;
49 hypre_ParVector *Ztemp;
50 hypre_ParVector *Aux_U;
51 hypre_ParVector *Aux_F;
52
53 int **CF_marker_array;
54 /* int **unknown_map_array;
55 int **point_map_array;
56 int **v_at_point_array; */
57
58 double cycle_op_count;
59 int cycle_type;
60 int num_levels;
61 int max_levels;
62
63 double *num_coeffs;
64 int *num_grid_sweeps;
65 int *grid_relax_type;
66 int **grid_relax_points;
67
68 /*int block_mode;*/
69
70 double *max_eig_est;
71 double *min_eig_est;
72 int cheby_order;
73 double cheby_eig_ratio;
74
75
76 /* Local variables */
77 int *lev_counter;
78 int Solve_err_flag;
79 int k;
80 int j, jj;
81 int level;
82 int cycle_param;
83 int coarse_grid;
84 int fine_grid;
85 int Not_Finished;
86 int num_sweep;
87 int cg_num_sweep = 1;
88 int relax_type;
89 int relax_points;
90 int relax_order;
91 int relax_local;
92 int old_version = 0;
93 double *relax_weight;
94 double *omega;
95 double beta;
96 int local_size;
97
98 double alpha;
99 double **l1_norms = NULL;
100 double *l1_norms_level;
101
102#if 0
103 double *D_mat;
104 double *S_vec;
105#endif
106
107
108 int myid;
109 int seq_cg = 0;
110
111 if (seq_data)
112 seq_cg = 1;
113
114
115 MPI_Comm_rank(MPI_COMM_WORLD, &myid);
116
117 /* Acquire data and allocate storage */
118
119 A_array = hypre_ParAMGDataAArray(amg_data);
120 P_array = hypre_ParAMGDataPArray(amg_data);
121 R_array = hypre_ParAMGDataRArray(amg_data);
122 CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
123 Vtemp = hypre_ParAMGDataVtemp(amg_data);
124 Rtemp = hypre_ParAMGDataRtemp(amg_data);
125 Ptemp = hypre_ParAMGDataPtemp(amg_data);
126 Ztemp = hypre_ParAMGDataZtemp(amg_data);
127 num_levels = hypre_ParAMGDataNumLevels(amg_data);
128 max_levels = hypre_ParAMGDataMaxLevels(amg_data);
129 cycle_type = hypre_ParAMGDataCycleType(amg_data);
130
131 num_grid_sweeps = hypre_ParAMGDataNumGridSweeps(amg_data);
132 grid_relax_type = hypre_ParAMGDataGridRelaxType(amg_data);
133 grid_relax_points = hypre_ParAMGDataGridRelaxPoints(amg_data);
134 relax_order = hypre_ParAMGDataRelaxOrder(amg_data);
135 relax_weight = hypre_ParAMGDataRelaxWeight(amg_data);
136 omega = hypre_ParAMGDataOmega(amg_data);
137 l1_norms = hypre_ParAMGDataL1Norms(amg_data);
138
139 max_eig_est = hypre_ParAMGDataMaxEigEst(amg_data);
140 min_eig_est = hypre_ParAMGDataMinEigEst(amg_data);
141 cheby_order = hypre_ParAMGDataChebyOrder(amg_data);
142 cheby_eig_ratio = hypre_ParAMGDataChebyEigRatio(amg_data);
143
144 cycle_op_count = hypre_ParAMGDataCycleOpCount(amg_data);
145
146 lev_counter = hypre_CTAlloc(int, num_levels);
147
148 /* Initialize */
149
150 Solve_err_flag = 0;
151
152 if (grid_relax_points) old_version = 1;
153
154 num_coeffs = hypre_CTAlloc(double, num_levels);
155 num_coeffs[0] = hypre_ParCSRMatrixDNumNonzeros(A_array[0]);
156 comm = hypre_ParCSRMatrixComm(A_array[0]);
157
158 for (j = 1; j < num_levels; j++)
159 num_coeffs[j] = hypre_ParCSRMatrixDNumNonzeros(A_array[j]);
160
161
162 /*---------------------------------------------------------------------
163 * Initialize cycling control counter
164 *
165 * Cycling is controlled using a level counter: lev_counter[k]
166 *
167 * Each time relaxation is performed on level k, the
168 * counter is decremented by 1. If the counter is then
169 * negative, we go to the next finer level. If non-
170 * negative, we go to the next coarser level. The
171 * following actions control cycling:
172 *
173 * a. lev_counter[0] is initialized to 1.
174 * b. lev_counter[k] is initialized to cycle_type for k>0.
175 *
176 * c. During cycling, when going down to level k, lev_counter[k]
177 * is set to the max of (lev_counter[k],cycle_type)
178 *---------------------------------------------------------------------*/
179
180 Not_Finished = 1;
181
182 smoother = hypre_ParAMGDataSmoother(amg_data);
183
184 lev_counter[0] = 1;
185 for (k = 1; k < num_levels; ++k)
186 {
187 lev_counter[k] = cycle_type;
188 }
189
190 level = 0;
191 cycle_param = 1;
192
193 /*---------------------------------------------------------------------
194 * Main loop of cycling
195 *--------------------------------------------------------------------*/
196
197
198 while (Not_Finished)
199 {
200
201 if (num_levels > 1)
202 {
203 local_size
204 = hypre_VectorSize(hypre_ParVectorLocalVector(F_array[level]));
205 hypre_VectorSize(hypre_ParVectorLocalVector(Vtemp)) = local_size;
206 cg_num_sweep = 1;
207 num_sweep = num_grid_sweeps[cycle_param];
208 Aux_U = U_array[level];
209 Aux_F = F_array[level];
210 relax_type = grid_relax_type[cycle_param];
211 }
212 else /* do this for max levels = 1 also */
213 {
214 /* If no coarsening occurred, apply a simple smoother once */
215 Aux_U = U_array[level];
216 Aux_F = F_array[level];
217 num_sweep = 1;
218 relax_type = 0;
219 }
220
221 if (l1_norms != NULL)
222 l1_norms_level = l1_norms[level];
223 else
224 l1_norms_level = NULL;
225
226
227 if (cycle_param == 3 && seq_cg)
228 {
229 //do a seq amg solve
230 hypre_seqAMGCycle(amg_data,
231 level,
232 F_array,
233 U_array);
234
235
236
237 }
238 else
239 {
240 /*------------------------------------------------------------------
241 * Do the relaxation num_sweep times
242 *-----------------------------------------------------------------*/
243 for (jj = 0; jj < cg_num_sweep; jj++)
244 {
245
246 for (j = 0; j < num_sweep; j++)
247 {
248 if (num_levels == 1 && max_levels > 1)
249 {
250 relax_points = 0;
251 relax_local = 0;
252 }
253 else
254 {
255 if (old_version)
256 relax_points = grid_relax_points[cycle_param][j];
257 relax_local = relax_order;
258 }
259
260 /*-----------------------------------------------
261 * VERY sloppy approximation to cycle complexity
262 *-----------------------------------------------*/
263
264 if (old_version && level < num_levels -1)
265 {
266 switch (relax_points)
267 {
268 case 1:
269 cycle_op_count += num_coeffs[level+1];
270 break;
271
272 case -1:
273 cycle_op_count += (num_coeffs[level]-num_coeffs[level+1]);
274 break;
275 }
276 }
277 else
278 {
279 cycle_op_count += num_coeffs[level];
280 }
281 if (relax_type == 11 || relax_type == 15 || relax_type == 16
282 || relax_type == 17 )
283
284 { /* Chebyshev */
285
286 int scale = 0;
287 int variant = 0;
288
289 if (relax_type == 15 || relax_type == 17) /*modified Cheby */
290 {
291 variant = 1;
292
293 }
294
295 if (relax_type == 16 || relax_type == 17 ) /* scaled Cheby */
296 {
297 scale = 1;
298 }
299
300 hypre_ParCSRRelax_Cheby(A_array[level],
301 Aux_F,
302 max_eig_est[level],
303 min_eig_est[level],
304 cheby_eig_ratio, cheby_order, scale,
305 variant, Aux_U, Vtemp, Ztemp );
306 }
307
308 else if (relax_type ==12)
309 {
310 hypre_BoomerAMGRelax_FCFJacobi(A_array[level],
311 Aux_F,
312 CF_marker_array[level],
313 relax_weight[level],
314 Aux_U,
315 Vtemp);
316 }
317 else if (relax_type == 13 || relax_type == 14)
318 {
319 if (j ==0) /* do num sweep iterations of CG */
320 hypre_ParCSRRelax_CG( smoother[level],
321 A_array[level],
322 Aux_F,
323 Aux_U,
324 num_sweep);
325 }
326 else if (relax_type == 8)
327 {
328 hypre_ParCSRRelax_L1(A_array[level],
329 Aux_F,
330 relax_weight[level],
331 omega[level],
332 l1_norms_level,
333 Aux_U,
334 Vtemp,
335 Ztemp);
336 }
337 else if (relax_type == 19)
338 {
339 hypre_ParCSRRelax_L1_GS(A_array[level],
340 Aux_F,
341 relax_weight[level],
342 omega[level],
343 l1_norms_level,
344 Aux_U,
345 Vtemp,
346 Ztemp);
347 }
348 else if (relax_type == 18)
349 {
350 if (relax_order == 1 && cycle_type < 3)
351 {
352 int i;
353 int loc_relax_points[2];
354
355 if (cycle_type < 2)
356 {
357 loc_relax_points[0] = 1;
358 loc_relax_points[1] = -1;
359 }
360 else
361 {
362 loc_relax_points[0] = -1;
363 loc_relax_points[1] = 1;
364 }
365 for (i=0; i < 2; i++)
366 hypre_ParCSRRelax_L1_Jacobi(A_array[level],
367 Aux_F,
368 CF_marker_array[level],
369 loc_relax_points[i],
370 relax_weight[level],
371 l1_norms_level,
372 Aux_U,
373 Vtemp);
374 }
375 else
376 {
377 hypre_ParCSRRelax_L1_Jacobi(A_array[level],
378 Aux_F,
379 CF_marker_array[level],
380 0,
381 relax_weight[level],
382 l1_norms_level,
383 Aux_U,
384 Vtemp);
385 }
386 }
387 else
388 {
389
390 if (old_version)
391 {
392 Solve_err_flag = hypre_BoomerAMGRelax(A_array[level],
393 Aux_F,
394 CF_marker_array[level],
395 relax_type,
396 relax_points,
397 relax_weight[level],
398 omega[level],
399 l1_norms_level,
400 Aux_U,
401 Vtemp, Ztemp);
402 }
403 else
404 {
405 Solve_err_flag = hypre_BoomerAMGRelaxIF(A_array[level],
406 Aux_F,
407 CF_marker_array[level],
408 relax_type,
409 relax_local,
410 cycle_param,
411 relax_weight[level],
412 omega[level],
413 l1_norms_level,
414 Aux_U,
415 Vtemp, Ztemp);
416
417 }
418 }
419
420 if (Solve_err_flag != 0)
421 return(Solve_err_flag);
422 }
423 }
424 }//end of relaxation
425
426
427 /*------------------------------------------------------------------
428 * Decrement the control counter and determine which grid to visit next
429 *-----------------------------------------------------------------*/
430
431 --lev_counter[level];
432
433 if (lev_counter[level] >= 0 && level != num_levels-1)
434 {
435
436 /*---------------------------------------------------------------
437 * Visit coarser level next.
438 * Compute residual using hypre_ParCSRMatrixMatvec.
439 * Perform restriction using hypre_ParCSRMatrixMatvecT.
440 * Reset counters and cycling parameters for coarse level
441 *--------------------------------------------------------------*/
442
443 fine_grid = level;
444 coarse_grid = level + 1;
445
446 hypre_ParVectorSetConstantValues(U_array[coarse_grid], 0.0);
447
448 hypre_ParVectorCopy(F_array[fine_grid],Vtemp);
449 alpha = -1.0;
450 beta = 1.0;
451
452 hypre_ParCSRMatrixMatvec(alpha, A_array[fine_grid], U_array[fine_grid],
453 beta, Vtemp);
454
455 alpha = 1.0;
456 beta = 0.0;
457
458 hypre_ParCSRMatrixMatvecT(alpha,R_array[fine_grid],Vtemp,
459 beta,F_array[coarse_grid]);
460
461 ++level;
462 lev_counter[level] = hypre_max(lev_counter[level],cycle_type);
463 cycle_param = 1;
464 if (level == num_levels-1) cycle_param = 3;
465 }
466
467 else if (level != 0)
468 {
469 /*---------------------------------------------------------------
470 * Visit finer level next.
471 * Interpolate and add correction using hypre_ParCSRMatrixMatvec.
472 * Reset counters and cycling parameters for finer level.
473 *--------------------------------------------------------------*/
474
475 fine_grid = level - 1;
476 coarse_grid = level;
477 alpha = 1.0;
478 beta = 1.0;
479 hypre_ParCSRMatrixMatvec(alpha, P_array[fine_grid],
480 U_array[coarse_grid],
481 beta, U_array[fine_grid]);
482
483 --level;
484 cycle_param = 2;
485
486 }
487 else
488 {
489 Not_Finished = 0;
490 }
491
492 }
493
494
495 hypre_ParAMGDataCycleOpCount(amg_data) = cycle_op_count;
496
497 hypre_TFree(lev_counter);
498 hypre_TFree(num_coeffs);
499
500 return(Solve_err_flag);
501}
Note: See TracBrowser for help on using the repository browser.