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 #include "_hypre_parcsr_ls.h"
9 #include "par_amg.h"
10 #include "../parcsr_block_mv/par_csr_block_matrix.h"
11
12 #define DEBUG 0
13 #define PRINT_CF 0
14
15 #define DEBUG_SAVE_ALL_OPS 0
16 /*****************************************************************************
17 *
18 * Routine for driving the setup phase of AMG
19 *
20 *****************************************************************************/
21
22 /*****************************************************************************
23 * hypre_BoomerAMGSetup
24 *****************************************************************************/
25
26 HYPRE_Int
hypre_BoomerAMGSetup(void * amg_vdata,hypre_ParCSRMatrix * A,hypre_ParVector * f,hypre_ParVector * u)27 hypre_BoomerAMGSetup( void *amg_vdata,
28 hypre_ParCSRMatrix *A,
29 hypre_ParVector *f,
30 hypre_ParVector *u )
31 {
32 MPI_Comm comm = hypre_ParCSRMatrixComm(A);
33 hypre_ParAMGData *amg_data = (hypre_ParAMGData*) amg_vdata;
34
35 /* Data Structure variables */
36 hypre_ParCSRMatrix **A_array;
37 hypre_ParVector **F_array;
38 hypre_ParVector **U_array;
39 hypre_ParVector *Vtemp = NULL;
40 hypre_ParVector *Rtemp = NULL;
41 hypre_ParVector *Ptemp = NULL;
42 hypre_ParVector *Ztemp = NULL;
43 hypre_ParCSRMatrix **P_array;
44 hypre_ParCSRMatrix **R_array;
45 hypre_ParVector *Residual_array;
46 hypre_IntArray **CF_marker_array;
47 hypre_IntArray **dof_func_array;
48 hypre_IntArray *dof_func;
49 HYPRE_Int *dof_func_data;
50 HYPRE_Real *relax_weight;
51 HYPRE_Real *omega;
52 HYPRE_Real schwarz_relax_wt = 1;
53 HYPRE_Real strong_threshold;
54 HYPRE_Int coarsen_cut_factor;
55 HYPRE_Int useSabs;
56 HYPRE_Real CR_strong_th;
57 HYPRE_Real max_row_sum;
58 HYPRE_Real trunc_factor, jacobi_trunc_threshold;
59 HYPRE_Real agg_trunc_factor, agg_P12_trunc_factor;
60 HYPRE_Real CR_rate;
61 HYPRE_Int relax_order;
62 HYPRE_Int max_levels;
63 HYPRE_Int amg_logging;
64 HYPRE_Int amg_print_level;
65 HYPRE_Int debug_flag;
66 HYPRE_Int dbg_flg;
67 HYPRE_Int local_num_vars;
68 HYPRE_Int P_max_elmts;
69 HYPRE_Int agg_P_max_elmts;
70 HYPRE_Int agg_P12_max_elmts;
71 HYPRE_Int IS_type;
72 HYPRE_Int num_CR_relax_steps;
73 HYPRE_Int CR_use_CG;
74 HYPRE_Int cgc_its; /* BM Aug 25, 2006 */
75 HYPRE_Int mult_additive = hypre_ParAMGDataMultAdditive(amg_data);
76 HYPRE_Int additive = hypre_ParAMGDataAdditive(amg_data);
77 HYPRE_Int simple = hypre_ParAMGDataSimple(amg_data);
78 HYPRE_Int add_last_lvl = hypre_ParAMGDataAddLastLvl(amg_data);
79 HYPRE_Int add_P_max_elmts = hypre_ParAMGDataMultAddPMaxElmts(amg_data);
80 HYPRE_Int keep_same_sign = hypre_ParAMGDataKeepSameSign(amg_data);
81 HYPRE_Real add_trunc_factor = hypre_ParAMGDataMultAddTruncFactor(amg_data);
82 HYPRE_Int add_rlx = hypre_ParAMGDataAddRelaxType(amg_data);
83 HYPRE_Real add_rlx_wt = hypre_ParAMGDataAddRelaxWt(amg_data);
84
85 hypre_ParCSRBlockMatrix **A_block_array, **P_block_array, **R_block_array;
86
87 HYPRE_MemoryLocation memory_location = hypre_ParCSRMatrixMemoryLocation(A);
88
89 /* Local variables */
90 HYPRE_Int *CF_marker;
91 hypre_IntArray *CF_marker_host;
92 hypre_IntArray *CFN_marker = NULL;
93 hypre_IntArray *CF2_marker = NULL;
94 hypre_IntArray *CF3_marker = NULL;
95 hypre_ParCSRMatrix *S = NULL, *Sabs = NULL;
96 hypre_ParCSRMatrix *S2;
97 hypre_ParCSRMatrix *SN = NULL;
98 hypre_ParCSRMatrix *SCR;
99 hypre_ParCSRMatrix *P = NULL;
100 hypre_ParCSRMatrix *R = NULL;
101 hypre_ParCSRMatrix *A_H;
102 hypre_ParCSRMatrix *AN = NULL;
103 hypre_ParCSRMatrix *P1;
104 hypre_ParCSRMatrix *P2;
105 hypre_ParCSRMatrix *Pnew = NULL;
106 HYPRE_Real *SmoothVecs = NULL;
107 hypre_Vector **l1_norms = NULL;
108 hypre_Vector **cheby_ds = NULL;
109 HYPRE_Real **cheby_coefs = NULL;
110
111 HYPRE_Int old_num_levels, num_levels;
112 HYPRE_Int level;
113 HYPRE_Int local_size, i, row;
114 HYPRE_BigInt first_local_row;
115 HYPRE_BigInt coarse_size;
116 HYPRE_Int coarsen_type;
117 HYPRE_Int measure_type;
118 HYPRE_Int setup_type;
119 HYPRE_BigInt fine_size;
120 HYPRE_Int offset;
121 HYPRE_Real size;
122 HYPRE_Int not_finished_coarsening = 1;
123 HYPRE_Int coarse_threshold = hypre_ParAMGDataMaxCoarseSize(amg_data);
124 HYPRE_Int min_coarse_size = hypre_ParAMGDataMinCoarseSize(amg_data);
125 HYPRE_Int seq_threshold = hypre_ParAMGDataSeqThreshold(amg_data);
126 HYPRE_Int j, k;
127 HYPRE_Int num_procs,my_id;
128 #if !defined(HYPRE_USING_CUDA) && !defined(HYPRE_USING_HIP)
129 HYPRE_Int num_threads = hypre_NumThreads();
130 #endif
131 HYPRE_Int *grid_relax_type = hypre_ParAMGDataGridRelaxType(amg_data);
132 HYPRE_Int num_functions = hypre_ParAMGDataNumFunctions(amg_data);
133 HYPRE_Int nodal = hypre_ParAMGDataNodal(amg_data);
134 HYPRE_Int nodal_levels = hypre_ParAMGDataNodalLevels(amg_data);
135 HYPRE_Int nodal_diag = hypre_ParAMGDataNodalDiag(amg_data);
136 HYPRE_Int num_paths = hypre_ParAMGDataNumPaths(amg_data);
137 HYPRE_Int agg_num_levels = hypre_ParAMGDataAggNumLevels(amg_data);
138 HYPRE_Int agg_interp_type = hypre_ParAMGDataAggInterpType(amg_data);
139 HYPRE_Int sep_weight = hypre_ParAMGDataSepWeight(amg_data);
140 hypre_IntArray *coarse_dof_func = NULL;
141 HYPRE_BigInt *coarse_pnts_global = NULL;
142 HYPRE_BigInt *coarse_pnts_global1 = NULL;
143 HYPRE_Int num_cg_sweeps;
144
145 HYPRE_Real *max_eig_est = NULL;
146 HYPRE_Real *min_eig_est = NULL;
147
148 HYPRE_Solver *smoother = NULL;
149 HYPRE_Int smooth_type = hypre_ParAMGDataSmoothType(amg_data);
150 HYPRE_Int smooth_num_levels = hypre_ParAMGDataSmoothNumLevels(amg_data);
151 HYPRE_Int sym;
152 HYPRE_Int nlevel;
153 HYPRE_Real thresh;
154 HYPRE_Real filter;
155 HYPRE_Real drop_tol;
156 HYPRE_Int max_nz_per_row;
157 char *euclidfile;
158 HYPRE_Int eu_level;
159 HYPRE_Int eu_bj;
160 HYPRE_Real eu_sparse_A;
161 HYPRE_Int ilu_type;
162 HYPRE_Int ilu_lfil;
163 HYPRE_Int ilu_max_row_nnz;
164 HYPRE_Int ilu_max_iter;
165 HYPRE_Real ilu_droptol;
166 HYPRE_Int ilu_reordering_type;
167 HYPRE_Int needZ = 0;
168
169 HYPRE_Int interp_type, restri_type;
170 HYPRE_Int post_interp_type; /* what to do after computing the interpolation matrix
171 0 for nothing, 1 for a Jacobi step */
172
173 /*for fittting interp vectors */
174 /*HYPRE_Int smooth_interp_vectors= hypre_ParAMGSmoothInterpVectors(amg_data); */
175 HYPRE_Real abs_q_trunc= hypre_ParAMGInterpVecAbsQTrunc(amg_data);
176 HYPRE_Int q_max = hypre_ParAMGInterpVecQMax(amg_data);
177 HYPRE_Int num_interp_vectors= hypre_ParAMGNumInterpVectors(amg_data);
178 HYPRE_Int num_levels_interp_vectors = hypre_ParAMGNumLevelsInterpVectors(amg_data);
179 hypre_ParVector **interp_vectors = hypre_ParAMGInterpVectors(amg_data);
180 hypre_ParVector ***interp_vectors_array= hypre_ParAMGInterpVectorsArray(amg_data);
181 HYPRE_Int interp_vec_variant= hypre_ParAMGInterpVecVariant(amg_data);
182 HYPRE_Int interp_refine= hypre_ParAMGInterpRefine(amg_data);
183 HYPRE_Int interp_vec_first_level= hypre_ParAMGInterpVecFirstLevel(amg_data);
184 HYPRE_Real *expandp_weights = hypre_ParAMGDataExpandPWeights(amg_data);
185
186 /* parameters for non-Galerkin stuff */
187 HYPRE_Int nongalerk_num_tol = hypre_ParAMGDataNonGalerkNumTol (amg_data);
188 HYPRE_Real *nongalerk_tol = hypre_ParAMGDataNonGalerkTol (amg_data);
189 HYPRE_Real nongalerk_tol_l = 0.0;
190 HYPRE_Real *nongal_tol_array = hypre_ParAMGDataNonGalTolArray (amg_data);
191
192 hypre_ParCSRBlockMatrix *A_H_block;
193
194 HYPRE_Int block_mode = 0;
195
196 HYPRE_Int mult_addlvl = hypre_max(mult_additive, simple);
197 HYPRE_Int addlvl = hypre_max(mult_addlvl, additive);
198 HYPRE_Int rap2 = hypre_ParAMGDataRAP2(amg_data);
199 HYPRE_Int keepTranspose = hypre_ParAMGDataKeepTranspose(amg_data);
200
201 HYPRE_Int local_coarse_size;
202 HYPRE_Int num_C_points_coarse = hypre_ParAMGDataNumCPoints(amg_data);
203 HYPRE_Int *C_points_local_marker = hypre_ParAMGDataCPointsLocalMarker(amg_data);
204 HYPRE_BigInt *C_points_marker = hypre_ParAMGDataCPointsMarker(amg_data);
205 HYPRE_Int num_F_points = hypre_ParAMGDataNumFPoints(amg_data);
206 HYPRE_BigInt *F_points_marker = hypre_ParAMGDataFPointsMarker(amg_data);
207 HYPRE_Int num_isolated_F_points = hypre_ParAMGDataNumIsolatedFPoints(amg_data);
208 HYPRE_BigInt *isolated_F_points_marker = hypre_ParAMGDataIsolatedFPointsMarker(amg_data);
209
210 HYPRE_Int *num_grid_sweeps = hypre_ParAMGDataNumGridSweeps(amg_data);
211 HYPRE_Int ns = num_grid_sweeps[1];
212 HYPRE_Real wall_time; /* for debugging instrumentation */
213 HYPRE_Int add_end;
214
215 #ifdef HYPRE_USING_DSUPERLU
216 HYPRE_Int dslu_threshold = hypre_ParAMGDataDSLUThreshold(amg_data);
217 #endif
218
219 hypre_MPI_Comm_size(comm, &num_procs);
220 hypre_MPI_Comm_rank(comm,&my_id);
221
222 /*A_new = hypre_CSRMatrixDeleteZeros(hypre_ParCSRMatrixDiag(A), 1.e-16);
223 hypre_CSRMatrixPrint(A_new, "Atestnew"); */
224 old_num_levels = hypre_ParAMGDataNumLevels(amg_data);
225 max_levels = hypre_ParAMGDataMaxLevels(amg_data);
226 add_end = hypre_min(add_last_lvl, max_levels-1);
227 if (add_end == -1) add_end = max_levels-1;
228 amg_logging = hypre_ParAMGDataLogging(amg_data);
229 amg_print_level = hypre_ParAMGDataPrintLevel(amg_data);
230 coarsen_type = hypre_ParAMGDataCoarsenType(amg_data);
231 measure_type = hypre_ParAMGDataMeasureType(amg_data);
232 setup_type = hypre_ParAMGDataSetupType(amg_data);
233 debug_flag = hypre_ParAMGDataDebugFlag(amg_data);
234 relax_weight = hypre_ParAMGDataRelaxWeight(amg_data);
235 omega = hypre_ParAMGDataOmega(amg_data);
236 sym = hypre_ParAMGDataSym(amg_data);
237 nlevel = hypre_ParAMGDataLevel(amg_data);
238 filter = hypre_ParAMGDataFilter(amg_data);
239 thresh = hypre_ParAMGDataThreshold(amg_data);
240 drop_tol = hypre_ParAMGDataDropTol(amg_data);
241 max_nz_per_row = hypre_ParAMGDataMaxNzPerRow(amg_data);
242 euclidfile = hypre_ParAMGDataEuclidFile(amg_data);
243 eu_level = hypre_ParAMGDataEuLevel(amg_data);
244 eu_sparse_A = hypre_ParAMGDataEuSparseA(amg_data);
245 eu_bj = hypre_ParAMGDataEuBJ(amg_data);
246 ilu_type = hypre_ParAMGDataILUType(amg_data);
247 ilu_lfil = hypre_ParAMGDataILULevel(amg_data);
248 ilu_max_row_nnz = hypre_ParAMGDataILUMaxRowNnz(amg_data);
249 ilu_droptol = hypre_ParAMGDataILUDroptol(amg_data);
250 ilu_max_iter = hypre_ParAMGDataILUMaxIter(amg_data);
251 ilu_reordering_type = hypre_ParAMGDataILULocalReordering(amg_data);
252 interp_type = hypre_ParAMGDataInterpType(amg_data);
253 restri_type = hypre_ParAMGDataRestriction(amg_data); /* RL */
254 post_interp_type = hypre_ParAMGDataPostInterpType(amg_data);
255 IS_type = hypre_ParAMGDataISType(amg_data);
256 num_CR_relax_steps = hypre_ParAMGDataNumCRRelaxSteps(amg_data);
257 CR_rate = hypre_ParAMGDataCRRate(amg_data);
258 CR_use_CG = hypre_ParAMGDataCRUseCG(amg_data);
259 cgc_its = hypre_ParAMGDataCGCIts(amg_data);
260
261 relax_order = hypre_ParAMGDataRelaxOrder(amg_data);
262
263 hypre_ParCSRMatrixSetNumNonzeros(A);
264 hypre_ParCSRMatrixSetDNumNonzeros(A);
265 hypre_ParAMGDataNumVariables(amg_data) = hypre_ParCSRMatrixNumRows(A);
266
267 if (num_procs == 1) seq_threshold = 0;
268 if (setup_type == 0) return hypre_error_flag;
269
270 S = NULL;
271
272 A_array = hypre_ParAMGDataAArray(amg_data);
273 P_array = hypre_ParAMGDataPArray(amg_data);
274 R_array = hypre_ParAMGDataRArray(amg_data);
275 CF_marker_array = hypre_ParAMGDataCFMarkerArray(amg_data);
276 dof_func_array = hypre_ParAMGDataDofFuncArray(amg_data);
277 dof_func = hypre_ParAMGDataDofFunc(amg_data);
278 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
279 first_local_row = hypre_ParCSRMatrixFirstRowIndex(A);
280
281 /* set size of dof_func hypre_IntArray if necessary */
282 if (dof_func && hypre_IntArraySize(dof_func) < 0)
283 {
284 hypre_IntArraySize(dof_func) = local_size;
285 }
286
287 A_block_array = hypre_ParAMGDataABlockArray(amg_data);
288 P_block_array = hypre_ParAMGDataPBlockArray(amg_data);
289 R_block_array = hypre_ParAMGDataRBlockArray(amg_data);
290
291 grid_relax_type[3] = hypre_ParAMGDataUserCoarseRelaxType(amg_data);
292
293 HYPRE_ANNOTATE_FUNC_BEGIN;
294
295 /* change in definition of standard and multipass interpolation, by
296 eliminating interp_type 9 and 5 and setting sep_weight instead
297 when using separation of weights option */
298 if (interp_type == 9)
299 {
300 interp_type = 8;
301 sep_weight = 1;
302 }
303 else if (interp_type == 5)
304 {
305 interp_type = 4;
306 sep_weight = 1;
307 }
308
309 /* Verify that if the user has selected the interp_vec_variant > 0
310 (so GM or LN interpolation) then they have nodal coarsening
311 selected also */
312 if (interp_vec_variant > 0 && nodal < 1)
313 {
314 nodal = 1;
315 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"WARNING: Changing to node-based coarsening because LN of GM interpolation has been specified via HYPRE_BoomerAMGSetInterpVecVariant.\n");
316 }
317
318 /* Verify that settings are correct for solving systmes */
319 /* If the user has specified either a block interpolation or a block relaxation then
320 we need to make sure the other has been choosen as well - so we can be
321 in "block mode" - storing only block matrices on the coarse levels*/
322 /* Furthermore, if we are using systems and nodal = 0, then
323 we will change nodal to 1 */
324 /* probably should disable stuff like smooth num levels at some point */
325
326
327 if (grid_relax_type[0] >= 20) /* block relaxation choosen */
328 {
329
330 if (!((interp_type >= 20 && interp_type != 100) || interp_type == 11 || interp_type == 10 ) )
331 {
332 hypre_ParAMGDataInterpType(amg_data) = 20;
333 interp_type = hypre_ParAMGDataInterpType(amg_data) ;
334 }
335
336 for (i=1; i < 3; i++)
337 {
338 if (grid_relax_type[i] < 20)
339 {
340 grid_relax_type[i] = 23;
341 }
342
343 }
344 if (grid_relax_type[3] < 20) grid_relax_type[3] = 29; /* GE */
345
346 block_mode = 1;
347 }
348
349 if ((interp_type >= 20 && interp_type != 100) || interp_type == 11 || interp_type == 10 ) /* block interp choosen */
350 {
351 if (!(nodal))
352 {
353 hypre_ParAMGDataNodal(amg_data) = 1;
354 nodal = hypre_ParAMGDataNodal(amg_data);
355 }
356 for (i=0; i < 3; i++)
357 {
358 if (grid_relax_type[i] < 20)
359 grid_relax_type[i] = 23;
360 }
361
362 if (grid_relax_type[3] < 20) grid_relax_type[3] = 29; /* GE */
363
364 block_mode = 1;
365
366 }
367
368 hypre_ParAMGDataBlockMode(amg_data) = block_mode;
369
370
371 /* end of systems checks */
372
373 /* free up storage in case of new setup without previous destroy */
374
375 if (A_array || A_block_array || P_array || P_block_array || CF_marker_array ||
376 dof_func_array || R_array || R_block_array)
377 {
378 for (j = 1; j < old_num_levels; j++)
379 {
380 if (A_array[j])
381 {
382 hypre_ParCSRMatrixDestroy(A_array[j]);
383 A_array[j] = NULL;
384 }
385
386 if (A_block_array[j])
387 {
388 hypre_ParCSRBlockMatrixDestroy(A_block_array[j]);
389 A_block_array[j] = NULL;
390 }
391
392 hypre_IntArrayDestroy(dof_func_array[j]);
393 dof_func_array[j] = NULL;
394 }
395
396 for (j = 0; j < old_num_levels-1; j++)
397 {
398 if (P_array[j])
399 {
400 hypre_ParCSRMatrixDestroy(P_array[j]);
401 P_array[j] = NULL;
402 }
403
404 if (P_block_array[j])
405 {
406 hypre_ParCSRBlockMatrixDestroy(P_block_array[j]);
407 P_block_array[j] = NULL;
408 }
409 /* RL */
410 if (R_array[j])
411 {
412 hypre_ParCSRMatrixDestroy(R_array[j]);
413 R_array[j] = NULL;
414 }
415
416 if (R_block_array[j])
417 {
418 hypre_ParCSRBlockMatrixDestroy(R_block_array[j]);
419 R_block_array[j] = NULL;
420 }
421 }
422
423 /* Special case use of CF_marker_array when old_num_levels == 1
424 requires us to attempt this deallocation every time */
425 hypre_IntArrayDestroy(CF_marker_array[0]);
426 CF_marker_array[0] = NULL;
427
428 for (j = 1; j < old_num_levels-1; j++)
429 {
430 hypre_IntArrayDestroy(CF_marker_array[j]);
431 CF_marker_array[j] = NULL;
432 }
433 }
434
435 {
436 MPI_Comm new_comm = hypre_ParAMGDataNewComm(amg_data);
437 void *amg = hypre_ParAMGDataCoarseSolver(amg_data);
438 if (hypre_ParAMGDataRtemp(amg_data))
439 {
440 hypre_ParVectorDestroy(hypre_ParAMGDataRtemp(amg_data));
441 hypre_ParAMGDataRtemp(amg_data) = NULL;
442 }
443 if (hypre_ParAMGDataPtemp(amg_data))
444 {
445 hypre_ParVectorDestroy(hypre_ParAMGDataPtemp(amg_data));
446 hypre_ParAMGDataPtemp(amg_data) = NULL;
447 }
448 if (hypre_ParAMGDataZtemp(amg_data))
449 {
450 hypre_ParVectorDestroy(hypre_ParAMGDataZtemp(amg_data));
451 hypre_ParAMGDataZtemp(amg_data) = NULL;
452 }
453
454 if (hypre_ParAMGDataACoarse(amg_data))
455 {
456 hypre_ParCSRMatrixDestroy(hypre_ParAMGDataACoarse(amg_data));
457 hypre_ParAMGDataACoarse(amg_data) = NULL;
458 }
459
460 if (hypre_ParAMGDataUCoarse(amg_data))
461 {
462 hypre_ParVectorDestroy(hypre_ParAMGDataUCoarse(amg_data));
463 hypre_ParAMGDataUCoarse(amg_data) = NULL;
464 }
465
466 if (hypre_ParAMGDataFCoarse(amg_data))
467 {
468 hypre_ParVectorDestroy(hypre_ParAMGDataFCoarse(amg_data));
469 hypre_ParAMGDataFCoarse(amg_data) = NULL;
470 }
471
472 hypre_TFree(hypre_ParAMGDataAMat(amg_data), HYPRE_MEMORY_HOST);
473 hypre_TFree(hypre_ParAMGDataAInv(amg_data), HYPRE_MEMORY_HOST);
474 hypre_TFree(hypre_ParAMGDataBVec(amg_data), HYPRE_MEMORY_HOST);
475 hypre_TFree(hypre_ParAMGDataCommInfo(amg_data), HYPRE_MEMORY_HOST);
476
477 if (new_comm != hypre_MPI_COMM_NULL)
478 {
479 hypre_MPI_Comm_free (&new_comm);
480 hypre_ParAMGDataNewComm(amg_data) = hypre_MPI_COMM_NULL;
481 }
482
483 if (amg)
484 {
485 hypre_BoomerAMGDestroy (amg);
486 hypre_ParAMGDataCoarseSolver(amg_data) = NULL;
487 }
488
489 hypre_TFree(hypre_ParAMGDataMaxEigEst(amg_data), HYPRE_MEMORY_HOST);
490 hypre_TFree(hypre_ParAMGDataMinEigEst(amg_data), HYPRE_MEMORY_HOST);
491
492 if (hypre_ParAMGDataL1Norms(amg_data))
493 {
494 for (i = 0; i < old_num_levels; i++)
495 {
496 hypre_SeqVectorDestroy(hypre_ParAMGDataL1Norms(amg_data)[i]);
497 }
498 hypre_TFree(hypre_ParAMGDataL1Norms(amg_data), HYPRE_MEMORY_HOST);
499 }
500 if (smooth_num_levels && smoother)
501 {
502 if (smooth_num_levels > old_num_levels-1)
503 {
504 smooth_num_levels = old_num_levels -1;
505 }
506 if (hypre_ParAMGDataSmoothType(amg_data) == 7)
507 {
508 for (i=0; i < smooth_num_levels; i++)
509 {
510 if (smoother[i])
511 {
512 HYPRE_ParCSRPilutDestroy(smoother[i]);
513 smoother[i] = NULL;
514 }
515 }
516 }
517 else if (hypre_ParAMGDataSmoothType(amg_data) == 8)
518 {
519 for (i=0; i < smooth_num_levels; i++)
520 {
521 if (smoother[i])
522 {
523 HYPRE_ParCSRParaSailsDestroy(smoother[i]);
524 smoother[i] = NULL;
525 }
526 }
527 }
528 else if (hypre_ParAMGDataSmoothType(amg_data) == 9)
529 {
530 for (i=0; i < smooth_num_levels; i++)
531 {
532 if (smoother[i])
533 {
534 HYPRE_EuclidDestroy(smoother[i]);
535 smoother[i] = NULL;
536 }
537 }
538 }
539 else if (hypre_ParAMGDataSmoothType(amg_data) == 5)
540 {
541 for (i=0; i < smooth_num_levels; i++)
542 {
543 if (smoother[i])
544 {
545 HYPRE_ILUDestroy(smoother[i]);
546 smoother[i] = NULL;
547 }
548 }
549 }
550 else if (hypre_ParAMGDataSmoothType(amg_data) == 6)
551 {
552 for (i=0; i < smooth_num_levels; i++)
553 {
554 if (smoother[i])
555 {
556 HYPRE_SchwarzDestroy(smoother[i]);
557 smoother[i] = NULL;
558 }
559 }
560 }
561 hypre_TFree(hypre_ParAMGDataSmoother(amg_data), HYPRE_MEMORY_HOST);
562 }
563 if ( hypre_ParAMGDataResidual(amg_data) )
564 {
565 hypre_ParVectorDestroy( hypre_ParAMGDataResidual(amg_data) );
566 hypre_ParAMGDataResidual(amg_data) = NULL;
567 }
568 }
569
570 if (A_array == NULL)
571 {
572 A_array = hypre_CTAlloc(hypre_ParCSRMatrix*, max_levels, HYPRE_MEMORY_HOST);
573 }
574 if (A_block_array == NULL)
575 {
576 A_block_array = hypre_CTAlloc(hypre_ParCSRBlockMatrix*, max_levels, HYPRE_MEMORY_HOST);
577 }
578
579 if (P_array == NULL && max_levels > 1)
580 {
581 P_array = hypre_CTAlloc(hypre_ParCSRMatrix*, max_levels-1, HYPRE_MEMORY_HOST);
582 }
583 if (P_block_array == NULL && max_levels > 1)
584 {
585 P_block_array = hypre_CTAlloc(hypre_ParCSRBlockMatrix*, max_levels-1, HYPRE_MEMORY_HOST);
586 }
587
588 /* RL: if retri_type != 0, R != P^T, allocate R matrices */
589 if (restri_type)
590 {
591 if (R_array == NULL && max_levels > 1)
592 {
593 R_array = hypre_CTAlloc(hypre_ParCSRMatrix*, max_levels-1, HYPRE_MEMORY_HOST);
594 }
595 if (R_block_array == NULL && max_levels > 1)
596 {
597 R_block_array = hypre_CTAlloc(hypre_ParCSRBlockMatrix*, max_levels-1, HYPRE_MEMORY_HOST);
598 }
599 }
600
601 if (CF_marker_array == NULL)
602 {
603 CF_marker_array = hypre_CTAlloc(hypre_IntArray*, max_levels, HYPRE_MEMORY_HOST);
604 }
605
606 if (num_C_points_coarse > 0)
607 {
608 k = 0;
609 for (j = 0; j < num_C_points_coarse; j++)
610 {
611 row = (HYPRE_Int) (C_points_marker[j] - first_local_row);
612 if ((row >= 0) && (row < local_size))
613 {
614 C_points_local_marker[k++] = row;
615 }
616 }
617 num_C_points_coarse = k;
618 }
619
620 if (dof_func_array == NULL)
621 {
622 dof_func_array = hypre_CTAlloc(hypre_IntArray*, max_levels, HYPRE_MEMORY_HOST);
623 }
624
625 if (num_functions > 1 && dof_func == NULL)
626 {
627 dof_func = hypre_IntArrayCreate(local_size);
628 hypre_IntArrayInitialize(dof_func);
629
630 offset = (HYPRE_Int) ( first_local_row % ((HYPRE_BigInt) num_functions) );
631
632 #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
633 hypre_BoomerAMGInitDofFuncDevice(hypre_IntArrayData(dof_func), local_size, offset, num_functions);
634 #else
635 for (i = 0; i < local_size; i++)
636 {
637 hypre_IntArrayData(dof_func)[i] = (i + offset) % num_functions;
638 }
639 #endif /* defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP) */
640 }
641
642 A_array[0] = A;
643
644 /* interp vectors setup */
645 if (interp_vec_variant == 1)
646 {
647 num_levels_interp_vectors = interp_vec_first_level + 1;
648 hypre_ParAMGNumLevelsInterpVectors(amg_data) = num_levels_interp_vectors;
649 }
650 if ( interp_vec_variant > 0 && num_interp_vectors > 0)
651 {
652 interp_vectors_array = hypre_CTAlloc(hypre_ParVector**, num_levels_interp_vectors, HYPRE_MEMORY_HOST);
653 interp_vectors_array[0] = interp_vectors;
654 hypre_ParAMGInterpVectorsArray(amg_data)= interp_vectors_array;
655 }
656
657 if (block_mode)
658 {
659 A_block_array[0] = hypre_ParCSRBlockMatrixConvertFromParCSRMatrix(A_array[0],
660 num_functions);
661 hypre_ParCSRBlockMatrixSetNumNonzeros(A_block_array[0]);
662 hypre_ParCSRBlockMatrixSetDNumNonzeros(A_block_array[0]);
663 }
664
665 dof_func_array[0] = dof_func;
666 hypre_ParAMGDataCFMarkerArray(amg_data) = CF_marker_array;
667 hypre_ParAMGDataNumCPoints(amg_data) = num_C_points_coarse;
668 hypre_ParAMGDataDofFunc(amg_data) = dof_func;
669 hypre_ParAMGDataDofFuncArray(amg_data) = dof_func_array;
670 hypre_ParAMGDataAArray(amg_data) = A_array;
671 hypre_ParAMGDataPArray(amg_data) = P_array;
672
673 /* RL: if R != P^T */
674 if (restri_type)
675 {
676 hypre_ParAMGDataRArray(amg_data) = R_array;
677 }
678 else
679 {
680 hypre_ParAMGDataRArray(amg_data) = P_array;
681 }
682
683 hypre_ParAMGDataABlockArray(amg_data) = A_block_array;
684 hypre_ParAMGDataPBlockArray(amg_data) = P_block_array;
685
686 /* RL: if R != P^T */
687 if (restri_type)
688 {
689 hypre_ParAMGDataRBlockArray(amg_data) = R_block_array;
690 }
691 else
692 {
693 hypre_ParAMGDataRBlockArray(amg_data) = P_block_array;
694 }
695
696 Vtemp = hypre_ParAMGDataVtemp(amg_data);
697
698 if (Vtemp != NULL)
699 {
700 hypre_ParVectorDestroy(Vtemp);
701 Vtemp = NULL;
702 }
703
704 Vtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
705 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
706 hypre_ParCSRMatrixRowStarts(A_array[0]));
707 hypre_ParVectorInitialize_v2(Vtemp, memory_location);
708 hypre_ParAMGDataVtemp(amg_data) = Vtemp;
709
710 /* If we are doing Cheby relaxation, we also need up two more temp vectors.
711 * If cheby_scale is false, only need one, otherwise need two */
712 if ((smooth_num_levels > 0 && smooth_type > 9) || relax_weight[0] < 0 || omega[0] < 0 ||
713 hypre_ParAMGDataSchwarzRlxWeight(amg_data) < 0 ||
714 (grid_relax_type[0] == 16 || grid_relax_type[1] == 16 || grid_relax_type[2] == 16 || grid_relax_type[3] == 16))
715 {
716 Ptemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
717 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
718 hypre_ParCSRMatrixRowStarts(A_array[0]));
719 hypre_ParVectorInitialize_v2(Ptemp, memory_location);
720 hypre_ParAMGDataPtemp(amg_data) = Ptemp;
721
722 /* If not doing chebyshev relaxation, or (doing chebyshev relaxation and scaling) */
723 if (!(grid_relax_type[0] == 16 || grid_relax_type[1] == 16 || grid_relax_type[2] == 16 ||
724 grid_relax_type[3] == 16) ||
725 (hypre_ParAMGDataChebyScale(amg_data)))
726 {
727 Rtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
728 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
729 hypre_ParCSRMatrixRowStarts(A_array[0]));
730 hypre_ParVectorInitialize_v2(Rtemp, memory_location);
731 hypre_ParAMGDataRtemp(amg_data) = Rtemp;
732 }
733 }
734
735 /* See if we need the Ztemp vector */
736 if ( (smooth_num_levels > 0 && smooth_type > 6) || relax_weight[0] < 0 || omega[0] < 0 || hypre_ParAMGDataSchwarzRlxWeight(amg_data) < 0 )
737 {
738 needZ = hypre_max(needZ, 1);
739 }
740
741 if ( grid_relax_type[0] == 16 || grid_relax_type[1] == 16 || grid_relax_type[2] == 16 || grid_relax_type[3] == 16 )
742 {
743 /* Chebyshev */
744 needZ = hypre_max(needZ, 1);
745 }
746
747 #if !defined(HYPRE_USING_CUDA) && !defined(HYPRE_USING_HIP)
748 /* GPU impl. needs Z */
749 if (num_threads > 1)
750 #endif
751 {
752 /* we need the temp Z vector for relaxation 3 and 6 now if we are using threading */
753 for (j = 0; j < 4; j++)
754 {
755 if (grid_relax_type[j] == 3 || grid_relax_type[j] == 4 || grid_relax_type[j] == 6 ||
756 grid_relax_type[j] == 8 || grid_relax_type[j] == 13 || grid_relax_type[j] == 14 ||
757 grid_relax_type[j] == 11 || grid_relax_type[j] == 12)
758 {
759 needZ = hypre_max(needZ, 1);
760 break;
761 }
762 }
763 }
764
765 if (needZ)
766 {
767 Ztemp = hypre_ParMultiVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
768 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
769 hypre_ParCSRMatrixRowStarts(A_array[0]),
770 needZ);
771 hypre_ParVectorInitialize_v2(Ztemp, memory_location);
772 hypre_ParAMGDataZtemp(amg_data) = Ztemp;
773 }
774
775 F_array = hypre_ParAMGDataFArray(amg_data);
776 U_array = hypre_ParAMGDataUArray(amg_data);
777
778 if (F_array != NULL || U_array != NULL)
779 {
780 for (j = 1; j < old_num_levels; j++)
781 {
782 if (F_array[j] != NULL)
783 {
784 hypre_ParVectorDestroy(F_array[j]);
785 F_array[j] = NULL;
786 }
787 if (U_array[j] != NULL)
788 {
789 hypre_ParVectorDestroy(U_array[j]);
790 U_array[j] = NULL;
791 }
792 }
793 }
794
795 if (F_array == NULL)
796 {
797 F_array = hypre_CTAlloc(hypre_ParVector*, max_levels, HYPRE_MEMORY_HOST);
798 }
799 if (U_array == NULL)
800 {
801 U_array = hypre_CTAlloc(hypre_ParVector*, max_levels, HYPRE_MEMORY_HOST);
802 }
803
804 F_array[0] = f;
805 U_array[0] = u;
806
807 hypre_ParAMGDataFArray(amg_data) = F_array;
808 hypre_ParAMGDataUArray(amg_data) = U_array;
809
810 /*----------------------------------------------------------
811 * Initialize hypre_ParAMGData
812 *----------------------------------------------------------*/
813
814 not_finished_coarsening = 1;
815 level = 0;
816 HYPRE_ANNOTATE_MGLEVEL_BEGIN(level);
817
818 strong_threshold = hypre_ParAMGDataStrongThreshold(amg_data);
819 coarsen_cut_factor = hypre_ParAMGDataCoarsenCutFactor(amg_data);
820 useSabs = hypre_ParAMGDataSabs(amg_data);
821 CR_strong_th = hypre_ParAMGDataCRStrongTh(amg_data);
822 max_row_sum = hypre_ParAMGDataMaxRowSum(amg_data);
823 trunc_factor = hypre_ParAMGDataTruncFactor(amg_data);
824 agg_trunc_factor = hypre_ParAMGDataAggTruncFactor(amg_data);
825 agg_P12_trunc_factor = hypre_ParAMGDataAggP12TruncFactor(amg_data);
826 P_max_elmts = hypre_ParAMGDataPMaxElmts(amg_data);
827 agg_P_max_elmts = hypre_ParAMGDataAggPMaxElmts(amg_data);
828 agg_P12_max_elmts = hypre_ParAMGDataAggP12MaxElmts(amg_data);
829 jacobi_trunc_threshold = hypre_ParAMGDataJacobiTruncThreshold(amg_data);
830 if (smooth_num_levels > level)
831 {
832 smoother = hypre_CTAlloc(HYPRE_Solver, smooth_num_levels, HYPRE_MEMORY_HOST);
833 hypre_ParAMGDataSmoother(amg_data) = smoother;
834 }
835
836 /*-----------------------------------------------------
837 * Enter Coarsening Loop
838 *-----------------------------------------------------*/
839
840 while (not_finished_coarsening)
841 {
842 /* only do nodal coarsening on a fixed number of levels */
843 if (level >= nodal_levels)
844 {
845 nodal = 0;
846 }
847
848 if (block_mode)
849 {
850 fine_size = hypre_ParCSRBlockMatrixGlobalNumRows(A_block_array[level]);
851 }
852 else
853 {
854 fine_size = hypre_ParCSRMatrixGlobalNumRows(A_array[level]);
855 }
856
857 if (level > 0)
858 {
859
860 if (block_mode)
861 {
862 F_array[level] =
863 hypre_ParVectorCreateFromBlock(hypre_ParCSRBlockMatrixComm(A_block_array[level]),
864 hypre_ParCSRMatrixGlobalNumRows(A_block_array[level]),
865 hypre_ParCSRBlockMatrixRowStarts(A_block_array[level]),
866 hypre_ParCSRBlockMatrixBlockSize(A_block_array[level]));
867 hypre_ParVectorInitialize(F_array[level]);
868
869 U_array[level] =
870 hypre_ParVectorCreateFromBlock(hypre_ParCSRBlockMatrixComm(A_block_array[level]),
871 hypre_ParCSRMatrixGlobalNumRows(A_block_array[level]),
872 hypre_ParCSRBlockMatrixRowStarts(A_block_array[level]),
873 hypre_ParCSRBlockMatrixBlockSize(A_block_array[level]));
874
875 hypre_ParVectorInitialize(U_array[level]);
876 }
877 else
878 {
879 F_array[level] =
880 hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
881 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
882 hypre_ParCSRMatrixRowStarts(A_array[level]));
883 hypre_ParVectorInitialize_v2(F_array[level], memory_location);
884
885 U_array[level] =
886 hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
887 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
888 hypre_ParCSRMatrixRowStarts(A_array[level]));
889 hypre_ParVectorInitialize_v2(U_array[level], memory_location);
890 }
891 }
892
893 /*-------------------------------------------------------------
894 * Select coarse-grid points on 'level' : returns CF_marker
895 * for the level. Returns strength matrix, S
896 *--------------------------------------------------------------*/
897
898 dof_func_data = NULL;
899 if (dof_func_array[level] != NULL)
900 {
901 dof_func_data = hypre_IntArrayData(dof_func_array[level]);
902 }
903
904 if (debug_flag==1) wall_time = time_getWallclockSeconds();
905 if (debug_flag==3)
906 {
907 hypre_printf("\n ===== Proc = %d Level = %d =====\n",
908 my_id, level);
909 fflush(NULL);
910 }
911
912 if (max_levels == 1)
913 {
914 S = NULL;
915 coarse_pnts_global = NULL;
916 CF_marker_array[level] = hypre_IntArrayCreate(local_size);
917 hypre_IntArrayInitialize(CF_marker_array[level]);
918 hypre_IntArraySetConstantValues(CF_marker_array[level], 1);
919 coarse_size = fine_size;
920 }
921 else /* max_levels > 1 */
922 {
923 if (block_mode)
924 {
925 local_num_vars =
926 hypre_CSRBlockMatrixNumRows(hypre_ParCSRBlockMatrixDiag(A_block_array[level]));
927 }
928 else
929 {
930 local_num_vars =
931 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level]));
932 }
933 if (hypre_ParAMGDataGSMG(amg_data) ||
934 hypre_ParAMGDataInterpType(amg_data) == 1)
935 {
936 hypre_BoomerAMGCreateSmoothVecs(amg_data, A_array[level],
937 hypre_ParAMGDataNumGridSweeps(amg_data)[1],
938 level, &SmoothVecs);
939 }
940
941 /**** Get the Strength Matrix ****/
942 if (hypre_ParAMGDataGSMG(amg_data) == 0)
943 {
944 if (nodal) /* if we are solving systems and
945 not using the unknown approach then we need to
946 convert A to a nodal matrix - values that represent the
947 blocks - before getting the strength matrix*/
948 {
949
950 if (block_mode)
951 {
952 hypre_BoomerAMGBlockCreateNodalA(A_block_array[level], hypre_abs(nodal), nodal_diag, &AN);
953 }
954 else
955 {
956 hypre_BoomerAMGCreateNodalA(A_array[level],num_functions,
957 dof_func_data, hypre_abs(nodal), nodal_diag, &AN);
958 }
959
960 /* dof array not needed for creating S because we pass in that
961 the number of functions is 1 */
962 /* creat s two different ways - depending on if any entries in AN are negative: */
963
964 /* first: positive and negative entries */
965 if (nodal == 3 || nodal == 6 || nodal_diag > 0)
966 {
967 hypre_BoomerAMGCreateS(AN, strong_threshold, max_row_sum,
968 1, NULL, &SN);
969 }
970 else /* all entries are positive */
971 {
972 hypre_BoomerAMGCreateSabs(AN, strong_threshold, max_row_sum,
973 1, NULL, &SN);
974 }
975 }
976 else /* standard AMG or unknown approach */
977 {
978 if (!useSabs)
979 {
980 hypre_BoomerAMGCreateS(A_array[level], strong_threshold, max_row_sum,
981 num_functions, dof_func_data, &S);
982 }
983 else
984 {
985 /*
986 hypre_BoomerAMGCreateSabs(A_array[level], strong_threshold, max_row_sum,
987 num_functions, dof_func_array[level], &S);
988 */
989 hypre_BoomerAMGCreateSabs(A_array[level], strong_threshold, 1.0,
990 1, NULL, &S);
991 }
992 }
993
994 /* for AIR, need absolute value SOC: use a different threshold */
995 if (restri_type == 1 || restri_type == 2 || restri_type == 15)
996 {
997 HYPRE_Real strong_thresholdR;
998 strong_thresholdR = hypre_ParAMGDataStrongThresholdR(amg_data);
999 hypre_BoomerAMGCreateSabs(A_array[level], strong_thresholdR, 1.0,
1000 1, NULL, &Sabs);
1001 }
1002 }
1003 else
1004 {
1005 hypre_BoomerAMGCreateSmoothDirs(amg_data, A_array[level],
1006 SmoothVecs, strong_threshold,
1007 num_functions, dof_func_data, &S);
1008 }
1009
1010 /* Allocate CF_marker for the current level */
1011 CF_marker_array[level] = hypre_IntArrayCreate(local_num_vars);
1012 hypre_IntArrayInitialize(CF_marker_array[level]);
1013 CF_marker = hypre_IntArrayData(CF_marker_array[level]);
1014
1015 /* Set isolated fine points (SF_PT) given by the user */
1016 if ((num_isolated_F_points > 0) && (level == 0))
1017 {
1018 if (block_mode)
1019 {
1020 first_local_row = hypre_ParCSRBlockMatrixFirstRowIndex(A_block_array[level]);
1021 }
1022 else
1023 {
1024 first_local_row = hypre_ParCSRMatrixFirstRowIndex(A_array[level]);
1025 }
1026
1027 /* copy CF_marker to the host if needed */
1028 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1029 {
1030 CF_marker_host = hypre_IntArrayCloneDeep_v2(CF_marker_array[level], HYPRE_MEMORY_HOST);
1031 }
1032 else
1033 {
1034 CF_marker_host = CF_marker_array[level];
1035 }
1036
1037 for (j = 0; j < num_isolated_F_points; j++)
1038 {
1039 row = (HYPRE_Int) (isolated_F_points_marker[j] - first_local_row);
1040 if ((row >= 0) && (row < local_size))
1041 {
1042 hypre_IntArrayData(CF_marker_host)[row] = -3; // Assumes SF_PT == -3
1043 }
1044 }
1045
1046 /* copy back to device and destroy host copy */
1047 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1048 {
1049 hypre_IntArrayCopy(CF_marker_host, CF_marker_array[level]);
1050 hypre_IntArrayDestroy(CF_marker_host);
1051 }
1052
1053
1054 }
1055
1056 /**** Do the appropriate coarsening ****/
1057 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Coarsening");
1058
1059 if (nodal == 0) /* no nodal coarsening */
1060 {
1061 if (coarsen_type == 6)
1062 hypre_BoomerAMGCoarsenFalgout(S, A_array[level], measure_type,
1063 coarsen_cut_factor, debug_flag, &(CF_marker_array[level]));
1064 else if (coarsen_type == 7)
1065 hypre_BoomerAMGCoarsen(S, A_array[level], 2,
1066 debug_flag, &(CF_marker_array[level]));
1067 else if (coarsen_type == 8)
1068 hypre_BoomerAMGCoarsenPMIS(S, A_array[level], 0,
1069 debug_flag, &(CF_marker_array[level]));
1070 else if (coarsen_type == 9)
1071 hypre_BoomerAMGCoarsenPMIS(S, A_array[level], 2,
1072 debug_flag, &(CF_marker_array[level]));
1073 else if (coarsen_type == 10)
1074 hypre_BoomerAMGCoarsenHMIS(S, A_array[level], measure_type,
1075 coarsen_cut_factor, debug_flag, &(CF_marker_array[level]));
1076 else if (coarsen_type == 21 || coarsen_type == 22)
1077 {
1078 #ifdef HYPRE_MIXEDINT
1079 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"CGC coarsening is not available in mixedint mode!");
1080 return hypre_error_flag;
1081 #endif
1082 hypre_BoomerAMGCoarsenCGCb(S, A_array[level], measure_type, coarsen_type,
1083 cgc_its, debug_flag, &(CF_marker_array[level]));
1084 }
1085 else if (coarsen_type == 98)
1086 hypre_BoomerAMGCoarsenCR1(A_array[level], &(CF_marker_array[level]),
1087 &coarse_size, num_CR_relax_steps, IS_type, 0);
1088 else if (coarsen_type == 99)
1089 {
1090 hypre_BoomerAMGCreateS(A_array[level],
1091 CR_strong_th, 1,
1092 num_functions, dof_func_data,&SCR);
1093 hypre_BoomerAMGCoarsenCR(A_array[level], &(CF_marker_array[level]),
1094 &coarse_size,
1095 num_CR_relax_steps, IS_type, 1, grid_relax_type[0],
1096 relax_weight[level], omega[level], CR_rate,
1097 NULL,NULL,CR_use_CG,SCR);
1098 hypre_ParCSRMatrixDestroy(SCR);
1099 }
1100 #if DEBUG
1101 else if (coarsen_type == 999)
1102 {
1103 /* RL_DEBUG: read C/F splitting from files */
1104 /* read from file */
1105 HYPRE_Int my_id;
1106 MPI_Comm comm = hypre_ParCSRMatrixComm(A_array[level]);
1107 hypre_MPI_Comm_rank(comm, &my_id);
1108 HYPRE_Int first_local_row = hypre_ParCSRMatrixFirstRowIndex(A_array[level]);
1109 HYPRE_Int local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level]));
1110 char CFfile[256], line[1024];
1111 hypre_sprintf(CFfile, "CF_%d.txt", level);
1112 hypre_printf("myid %d: level %d, read C/F from file %s, first_row %d, local_size %d\n",
1113 my_id, level, CFfile, first_local_row, local_size);
1114 FILE *fp;
1115 if ((fp = fopen(CFfile, "r")) == NULL)
1116 {
1117 hypre_printf("cannot open file %s\n", CFfile);
1118 exit(0);
1119 }
1120 HYPRE_Int i;
1121 for (i=0; i<first_local_row; i++)
1122 {
1123 if (fgets(line, 1024, fp) == NULL)
1124 {
1125 exit(-1);
1126 }
1127 /*HYPRE_Real tmp; fscanf(fp, "%le\n", &tmp);*/
1128 }
1129
1130 /* copy CF_marker to the host if needed */
1131 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1132 {
1133 CF_marker_host = hypre_IntArrayCloneDeep_v2(CF_marker_array[level], HYPRE_MEMORY_HOST);
1134 }
1135 else
1136 {
1137 CF_marker_host = CF_marker_array[level];
1138 }
1139
1140 for (i=0; i<local_size; i++)
1141 {
1142 HYPRE_Real dj;
1143 HYPRE_Int j;
1144 if (fgets(line, 1024, fp) == NULL)
1145 {
1146 hypre_printf("CF file read error\n");
1147 exit(0);
1148 }
1149 dj = atof(line);
1150 j = (HYPRE_Int) dj;
1151 /* 1: C, 0: F*/
1152 if (j == 1)
1153 {
1154 hypre_IntArrayData(CF_marker_host)[i] = 1;
1155 } else if (j == 0)
1156 {
1157 hypre_IntArrayData(CF_marker_host)[i] = -1;
1158 }
1159 else
1160 {
1161 hypre_printf("CF Error: %d\n", j);
1162 exit(0);
1163 }
1164 }
1165
1166 /* copy back to device and destroy host copy */
1167 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1168 {
1169 hypre_IntArrayCopy(CF_marker_host, CF_marker_array[level]);
1170 hypre_IntArrayDestroy(CF_marker_host);
1171 }
1172
1173 fclose(fp);
1174 }
1175 #endif
1176 else if (coarsen_type)
1177 {
1178 hypre_BoomerAMGCoarsenRuge(S, A_array[level], measure_type, coarsen_type,
1179 coarsen_cut_factor, debug_flag, &(CF_marker_array[level]));
1180 /* DEBUG: SAVE CF the splitting
1181 HYPRE_Int my_id;
1182 MPI_Comm comm = hypre_ParCSRMatrixComm(A_array[level]);
1183 hypre_MPI_Comm_rank(comm, &my_id);
1184 char CFfile[256];
1185 hypre_sprintf(CFfile, "hypreCF_%d.txt.%d", level, my_id);
1186 FILE *fp = fopen(CFfile, "w");
1187 for (i=0; i<local_size; i++)
1188 {
1189 HYPRE_Int k = CF_marker[i];
1190 HYPRE_Real j;
1191 if (k == 1) {
1192 j = 1.0;
1193 } else if (k == -1) {
1194 j = 0.0;
1195 } else {
1196 if (k < 0) {
1197 CF_marker[i] = -1;
1198 }
1199 j = (HYPRE_Real) k;
1200 }
1201 hypre_fprintf(fp, "%.18e\n", j);
1202 }
1203 fclose(fp);
1204 */
1205 }
1206 else
1207 {
1208 hypre_BoomerAMGCoarsen(S, A_array[level], 0,
1209 debug_flag, &(CF_marker_array[level]));
1210 }
1211
1212 if (level < agg_num_levels)
1213 {
1214 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1215 1, dof_func_array[level], CF_marker_array[level],
1216 &coarse_dof_func, &coarse_pnts_global1);
1217 hypre_BoomerAMGCreate2ndS(S, CF_marker, num_paths,
1218 coarse_pnts_global1, &S2);
1219 if (coarsen_type == 10)
1220 {
1221 hypre_BoomerAMGCoarsenHMIS(S2, S2, measure_type+3, coarsen_cut_factor,
1222 debug_flag, &CFN_marker);
1223 }
1224 else if (coarsen_type == 8)
1225 {
1226 hypre_BoomerAMGCoarsenPMIS(S2, S2, 3,
1227 debug_flag, &CFN_marker);
1228 }
1229 else if (coarsen_type == 9)
1230 {
1231 hypre_BoomerAMGCoarsenPMIS(S2, S2, 4,
1232 debug_flag, &CFN_marker);
1233 }
1234 else if (coarsen_type == 6)
1235 {
1236 hypre_BoomerAMGCoarsenFalgout(S2, S2, measure_type, coarsen_cut_factor,
1237 debug_flag, &CFN_marker);
1238 }
1239 else if (coarsen_type == 21 || coarsen_type == 22)
1240 {
1241 hypre_BoomerAMGCoarsenCGCb(S2, S2, measure_type,
1242 coarsen_type, cgc_its, debug_flag, &CFN_marker);
1243 }
1244 else if (coarsen_type == 7)
1245 {
1246 hypre_BoomerAMGCoarsen(S2, S2, 2, debug_flag, &CFN_marker);
1247 }
1248 else if (coarsen_type)
1249 {
1250 hypre_BoomerAMGCoarsenRuge(S2, S2, measure_type, coarsen_type,
1251 coarsen_cut_factor, debug_flag, &CFN_marker);
1252 }
1253 else
1254 {
1255 hypre_BoomerAMGCoarsen(S2, S2, 0, debug_flag, &CFN_marker);
1256 }
1257
1258 hypre_ParCSRMatrixDestroy(S2);
1259 }
1260 }
1261 else if (block_mode)
1262 {
1263 if (coarsen_type == 6)
1264 hypre_BoomerAMGCoarsenFalgout(SN, SN, measure_type, coarsen_cut_factor,
1265 debug_flag, &(CF_marker_array[level]));
1266 else if (coarsen_type == 7)
1267 hypre_BoomerAMGCoarsen(SN, SN, 2,
1268 debug_flag, &(CF_marker_array[level]));
1269 else if (coarsen_type == 8)
1270 hypre_BoomerAMGCoarsenPMIS(SN, SN, 0,
1271 debug_flag, &(CF_marker_array[level]));
1272 else if (coarsen_type == 9)
1273 hypre_BoomerAMGCoarsenPMIS(SN, SN, 2,
1274 debug_flag, &(CF_marker_array[level]));
1275 else if (coarsen_type == 10)
1276 hypre_BoomerAMGCoarsenHMIS(SN, SN, measure_type, coarsen_cut_factor,
1277 debug_flag, &(CF_marker_array[level]));
1278 else if (coarsen_type == 21 || coarsen_type == 22)
1279 hypre_BoomerAMGCoarsenCGCb(SN, SN, measure_type,
1280 coarsen_type, cgc_its, debug_flag, &(CF_marker_array[level]));
1281 else if (coarsen_type)
1282 hypre_BoomerAMGCoarsenRuge(SN, SN, measure_type, coarsen_type,
1283 coarsen_cut_factor, debug_flag, &(CF_marker_array[level]));
1284 else
1285 hypre_BoomerAMGCoarsen(SN, SN, 0, debug_flag, &(CF_marker_array[level]));
1286 }
1287 else if (nodal > 0)
1288 {
1289 if (coarsen_type == 6)
1290 hypre_BoomerAMGCoarsenFalgout(SN, SN, measure_type, coarsen_cut_factor,
1291 debug_flag, &CFN_marker);
1292 else if (coarsen_type == 7)
1293 hypre_BoomerAMGCoarsen(SN, SN, 2, debug_flag, &CFN_marker);
1294 else if (coarsen_type == 8)
1295 hypre_BoomerAMGCoarsenPMIS(SN, SN, 0, debug_flag, &CFN_marker);
1296 else if (coarsen_type == 9)
1297 hypre_BoomerAMGCoarsenPMIS(SN, SN, 2, debug_flag, &CFN_marker);
1298 else if (coarsen_type == 10)
1299 hypre_BoomerAMGCoarsenHMIS(SN, SN, measure_type, coarsen_cut_factor,
1300 debug_flag, &CFN_marker);
1301 else if (coarsen_type == 21 || coarsen_type == 22)
1302 hypre_BoomerAMGCoarsenCGCb(SN, SN, measure_type,
1303 coarsen_type, cgc_its, debug_flag, &CFN_marker);
1304 else if (coarsen_type)
1305 hypre_BoomerAMGCoarsenRuge(SN, SN, measure_type, coarsen_type,
1306 coarsen_cut_factor, debug_flag, &CFN_marker);
1307 else
1308 hypre_BoomerAMGCoarsen(SN, SN, 0,
1309 debug_flag, &CFN_marker);
1310 if (level < agg_num_levels)
1311 {
1312 hypre_BoomerAMGCoarseParms(comm, local_num_vars/num_functions,
1313 1, dof_func_array[level], CFN_marker,
1314 &coarse_dof_func, &coarse_pnts_global1);
1315 hypre_BoomerAMGCreate2ndS(SN, hypre_IntArrayData(CFN_marker), num_paths,
1316 coarse_pnts_global1, &S2);
1317 if (coarsen_type == 10)
1318 {
1319 hypre_BoomerAMGCoarsenHMIS(S2, S2, measure_type+3, coarsen_cut_factor,
1320 debug_flag, &CF2_marker);
1321 }
1322 else if (coarsen_type == 8)
1323 {
1324 hypre_BoomerAMGCoarsenPMIS(S2, S2, 3,
1325 debug_flag, &CF2_marker);
1326 }
1327 else if (coarsen_type == 9)
1328 {
1329 hypre_BoomerAMGCoarsenPMIS(S2, S2, 4,
1330 debug_flag, &CF2_marker);
1331 }
1332 else if (coarsen_type == 6)
1333 {
1334 hypre_BoomerAMGCoarsenFalgout(S2, S2, measure_type, coarsen_cut_factor,
1335 debug_flag, &CF2_marker);
1336 }
1337 else if (coarsen_type == 21 || coarsen_type == 22)
1338 {
1339 hypre_BoomerAMGCoarsenCGCb(S2, S2, measure_type,
1340 coarsen_type, cgc_its, debug_flag, &CF2_marker);
1341 }
1342 else if (coarsen_type == 7)
1343 {
1344 hypre_BoomerAMGCoarsen(S2, S2, 2, debug_flag, &CF2_marker);
1345 }
1346 else if (coarsen_type)
1347 {
1348 hypre_BoomerAMGCoarsenRuge(S2, S2, measure_type, coarsen_type,
1349 coarsen_cut_factor, debug_flag, &CF2_marker);
1350 }
1351 else
1352 {
1353 hypre_BoomerAMGCoarsen(S2, S2, 0, debug_flag, &CF2_marker);
1354 }
1355
1356 hypre_ParCSRMatrixDestroy(S2);
1357 S2 = NULL;
1358 }
1359 else
1360 {
1361 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1362 num_functions, nodal, keep_same_sign,
1363 &dof_func, &(CF_marker_array[level]),
1364 &S);
1365 hypre_IntArrayDestroy(CFN_marker);
1366 CFN_marker = NULL;
1367 hypre_ParCSRMatrixDestroy(SN);
1368 SN = NULL;
1369 hypre_ParCSRMatrixDestroy(AN);
1370 AN = NULL;
1371 }
1372 }
1373
1374 /**************************************************/
1375 /*********Set the fixed index to CF_marker*********/
1376 /* copy CF_marker to the host if needed */
1377 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1378 {
1379 CF_marker_host = hypre_IntArrayCloneDeep_v2(CF_marker_array[level], HYPRE_MEMORY_HOST);
1380 }
1381 else
1382 {
1383 CF_marker_host = CF_marker_array[level];
1384 }
1385
1386 /* Set fine points (F_PT) given by the user */
1387 if ((num_F_points > 0) && (level == 0))
1388 {
1389 for (j = 0; j < num_F_points; j++)
1390 {
1391 row = (HYPRE_Int) (F_points_marker[j] - first_local_row);
1392 if ((row >= 0) && (row < local_size))
1393 {
1394 hypre_IntArrayData(CF_marker_host)[row] = -1; // Assumes F_PT == -1
1395 }
1396 }
1397 }
1398
1399
1400 if (num_C_points_coarse > 0)
1401 {
1402 if (block_mode)
1403 {
1404 hypre_error_w_msg(HYPRE_ERROR_GENERIC, "Keeping coarse nodes in block mode is not implemented\n");
1405 }
1406 else if (level < hypre_ParAMGDataCPointsLevel(amg_data))
1407 {
1408 for (j = 0; j < num_C_points_coarse; j++)
1409 {
1410 hypre_IntArrayData(CF_marker_host)[C_points_local_marker[j]] = 2;
1411 }
1412
1413 local_coarse_size = 0;
1414 k = 0;
1415 for (j = 0; j < local_num_vars; j ++)
1416 {
1417 if (hypre_IntArrayData(CF_marker_host)[j] == 1)
1418 {
1419 local_coarse_size++;
1420 }
1421 else if (hypre_IntArrayData(CF_marker_host)[j] == 2)
1422 {
1423 if ((level + 1) < hypre_ParAMGDataCPointsLevel(amg_data))
1424 {
1425 C_points_local_marker[k++] = local_coarse_size;
1426 }
1427 local_coarse_size++;
1428 hypre_IntArrayData(CF_marker_host)[j] = 1;
1429 }
1430 }
1431 }
1432 }
1433
1434 /* copy back to device and destroy host copy */
1435 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
1436 {
1437 hypre_IntArrayCopy(CF_marker_host, CF_marker_array[level]);
1438 hypre_IntArrayDestroy(CF_marker_host);
1439 }
1440
1441 /*****xxxxxxxxxxxxx changes for min_coarse_size */
1442 /* here we will determine the coarse grid size to be able to
1443 determine if it is not smaller than requested minimal size */
1444
1445 if (level >= agg_num_levels)
1446 {
1447 if (block_mode)
1448 {
1449 hypre_BoomerAMGCoarseParms(comm,
1450 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(AN)),
1451 1, NULL, CF_marker_array[level], NULL, &coarse_pnts_global);
1452 }
1453 else
1454 {
1455 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1456 num_functions, dof_func_array[level], CF_marker_array[level],
1457 &coarse_dof_func, &coarse_pnts_global);
1458 }
1459 if (my_id == (num_procs -1)) coarse_size = coarse_pnts_global[1];
1460 hypre_MPI_Bcast(&coarse_size, 1, HYPRE_MPI_BIG_INT, num_procs-1, comm);
1461 /* if no coarse-grid, stop coarsening, and set the
1462 * coarsest solve to be a single sweep of default smoother or smoother set by user */
1463 if ((coarse_size == 0) || (coarse_size == fine_size))
1464 {
1465 HYPRE_Int *num_grid_sweeps = hypre_ParAMGDataNumGridSweeps(amg_data);
1466 HYPRE_Int **grid_relax_points = hypre_ParAMGDataGridRelaxPoints(amg_data);
1467 if (grid_relax_type[3] == 9 || grid_relax_type[3] == 99 ||
1468 grid_relax_type[3] == 19 || grid_relax_type[3] == 98)
1469 {
1470 grid_relax_type[3] = grid_relax_type[0];
1471 num_grid_sweeps[3] = 1;
1472 if (grid_relax_points) grid_relax_points[3][0] = 0;
1473 }
1474 if (S) hypre_ParCSRMatrixDestroy(S);
1475 if (SN) hypre_ParCSRMatrixDestroy(SN);
1476 if (AN) hypre_ParCSRMatrixDestroy(AN);
1477 //hypre_TFree(CF_marker, HYPRE_MEMORY_HOST);
1478 if (level > 0)
1479 {
1480 /* note special case treatment of CF_marker is necessary
1481 * to do CF relaxation correctly when num_levels = 1 */
1482 hypre_IntArrayDestroy(CF_marker_array[level]);
1483 CF_marker_array[level] = NULL;
1484 hypre_ParVectorDestroy(F_array[level]);
1485 hypre_ParVectorDestroy(U_array[level]);
1486 }
1487 coarse_size = fine_size;
1488
1489 if (Sabs)
1490 {
1491 hypre_ParCSRMatrixDestroy(Sabs);
1492 Sabs = NULL;
1493 }
1494
1495 HYPRE_ANNOTATE_REGION_END("%s", "Coarsening");
1496 break;
1497 }
1498
1499 if (coarse_size < min_coarse_size)
1500 {
1501 if (S) hypre_ParCSRMatrixDestroy(S);
1502 if (SN) hypre_ParCSRMatrixDestroy(SN);
1503 if (AN) hypre_ParCSRMatrixDestroy(AN);
1504 if (num_functions > 1) hypre_IntArrayDestroy(coarse_dof_func);
1505 hypre_IntArrayDestroy(CF_marker_array[level]);
1506 CF_marker_array[level] = NULL;
1507 hypre_TFree(coarse_pnts_global, HYPRE_MEMORY_HOST);
1508 if (level > 0)
1509 {
1510 hypre_ParVectorDestroy(F_array[level]);
1511 hypre_ParVectorDestroy(U_array[level]);
1512 }
1513 coarse_size = fine_size;
1514
1515 if (Sabs)
1516 {
1517 hypre_ParCSRMatrixDestroy(Sabs);
1518 Sabs = NULL;
1519 }
1520
1521 HYPRE_ANNOTATE_REGION_END("%s", "Coarsening");
1522 break;
1523 }
1524 }
1525
1526 /*****xxxxxxxxxxxxx changes for min_coarse_size end */
1527 HYPRE_ANNOTATE_REGION_END("%s", "Coarsening");
1528 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Interpolation");
1529
1530 if (level < agg_num_levels)
1531 {
1532 if (nodal == 0)
1533 {
1534 if (agg_interp_type == 1)
1535 {
1536 hypre_BoomerAMGBuildExtPIInterp(A_array[level],
1537 CF_marker, S, coarse_pnts_global1,
1538 num_functions, dof_func_data, debug_flag,
1539 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1540 }
1541 else if (agg_interp_type == 2)
1542 {
1543 hypre_BoomerAMGBuildStdInterp(A_array[level],
1544 CF_marker, S, coarse_pnts_global1,
1545 num_functions, dof_func_data, debug_flag,
1546 agg_P12_trunc_factor, agg_P12_max_elmts, 0, &P1);
1547 }
1548 else if (agg_interp_type == 3)
1549 {
1550 hypre_BoomerAMGBuildExtInterp(A_array[level],
1551 CF_marker, S, coarse_pnts_global1,
1552 num_functions, dof_func_data, debug_flag,
1553 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1554 }
1555 else if (agg_interp_type == 5)
1556 {
1557 hypre_BoomerAMGBuildModExtInterp(A_array[level],
1558 CF_marker, S, coarse_pnts_global1,
1559 num_functions, dof_func_data,
1560 debug_flag,
1561 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1562 }
1563 else if (agg_interp_type == 6)
1564 {
1565 hypre_BoomerAMGBuildModExtPIInterp(A_array[level],
1566 CF_marker, S, coarse_pnts_global1,
1567 num_functions, dof_func_data,
1568 debug_flag,
1569 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1570 }
1571 else if (agg_interp_type == 7)
1572 {
1573 hypre_BoomerAMGBuildModExtPEInterp(A_array[level],
1574 CF_marker, S, coarse_pnts_global1,
1575 num_functions, dof_func_data,
1576 debug_flag,
1577 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1578 }
1579
1580 if (agg_interp_type == 4)
1581 {
1582 hypre_BoomerAMGCorrectCFMarker(CF_marker_array[level], CFN_marker);
1583 hypre_IntArrayDestroy(CFN_marker);
1584 CFN_marker = NULL;
1585 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1586 num_functions, dof_func_array[level], CF_marker_array[level],
1587 &coarse_dof_func,&coarse_pnts_global);
1588 hypre_BoomerAMGBuildMultipass(A_array[level],
1589 CF_marker, S, coarse_pnts_global,
1590 num_functions, dof_func_data, debug_flag,
1591 agg_trunc_factor, agg_P_max_elmts, sep_weight,
1592 &P);
1593 }
1594 else if (agg_interp_type == 8)
1595 {
1596 hypre_BoomerAMGCorrectCFMarker(CF_marker_array[level], CFN_marker);
1597 hypre_IntArrayDestroy(CFN_marker);
1598 CFN_marker = NULL;
1599 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1600 num_functions, dof_func_array[level], CF_marker_array[level],
1601 &coarse_dof_func,&coarse_pnts_global);
1602 hypre_BoomerAMGBuildModMultipass(A_array[level],
1603 CF_marker, S, coarse_pnts_global,
1604 agg_trunc_factor, agg_P_max_elmts, 8,
1605 num_functions, dof_func_data, &P);
1606 }
1607 else if (agg_interp_type == 9)
1608 {
1609 hypre_BoomerAMGCorrectCFMarker(CF_marker_array[level], CFN_marker);
1610 hypre_IntArrayDestroy(CFN_marker);
1611 CFN_marker = NULL;
1612 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1613 num_functions, dof_func_array[level], CF_marker_array[level],
1614 &coarse_dof_func,&coarse_pnts_global);
1615 hypre_BoomerAMGBuildModMultipass(A_array[level],
1616 CF_marker, S, coarse_pnts_global,
1617 agg_trunc_factor, agg_P_max_elmts, 9,
1618 num_functions, dof_func_data, &P);
1619 }
1620 else
1621 {
1622 hypre_BoomerAMGCorrectCFMarker2 (CF_marker_array[level], (CFN_marker));
1623 hypre_IntArrayDestroy(CFN_marker);
1624 CFN_marker = NULL;
1625 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1626 num_functions, dof_func_array[level], CF_marker_array[level],
1627 &coarse_dof_func,&coarse_pnts_global);
1628 if (agg_interp_type == 1 || agg_interp_type == 6 )
1629 {
1630 hypre_BoomerAMGBuildPartialExtPIInterp(A_array[level],
1631 CF_marker, S, coarse_pnts_global,
1632 coarse_pnts_global1, num_functions,
1633 dof_func_data, debug_flag, agg_P12_trunc_factor,
1634 agg_P12_max_elmts, &P2);
1635 }
1636 else if (agg_interp_type == 2)
1637 {
1638 hypre_BoomerAMGBuildPartialStdInterp(A_array[level],
1639 CF_marker, S, coarse_pnts_global,
1640 coarse_pnts_global1, num_functions,
1641 dof_func_data, debug_flag, agg_P12_trunc_factor,
1642 agg_P12_max_elmts, sep_weight, &P2);
1643 }
1644 else if (agg_interp_type == 3)
1645 {
1646 hypre_BoomerAMGBuildPartialExtInterp(A_array[level],
1647 CF_marker, S, coarse_pnts_global,
1648 coarse_pnts_global1, num_functions,
1649 dof_func_data, debug_flag, agg_P12_trunc_factor,
1650 agg_P12_max_elmts, &P2);
1651 }
1652 else if (agg_interp_type == 5)
1653 {
1654 hypre_BoomerAMGBuildModPartialExtInterp(A_array[level],
1655 CF_marker, S, coarse_pnts_global,
1656 coarse_pnts_global1,
1657 num_functions, dof_func_data,
1658 debug_flag,
1659 agg_P12_trunc_factor, agg_P12_max_elmts, &P2);
1660 }
1661 else if (agg_interp_type == 7)
1662 {
1663 hypre_BoomerAMGBuildModPartialExtPEInterp(A_array[level],
1664 CF_marker, S, coarse_pnts_global,
1665 coarse_pnts_global1,
1666 num_functions, dof_func_data,
1667 debug_flag,
1668 agg_P12_trunc_factor, agg_P12_max_elmts, &P2);
1669 }
1670
1671 if (hypre_ParAMGDataModularizedMatMat(amg_data))
1672 {
1673 P = hypre_ParCSRMatMat(P1, P2);
1674 }
1675 else
1676 {
1677 P = hypre_ParMatmul(P1, P2);
1678 }
1679
1680 hypre_BoomerAMGInterpTruncation(P, agg_trunc_factor, agg_P_max_elmts);
1681
1682 hypre_MatvecCommPkgCreate(P);
1683 hypre_ParCSRMatrixDestroy(P1);
1684 hypre_ParCSRMatrixDestroy(P2);
1685 }
1686 }
1687 else if (nodal > 0)
1688 {
1689 if (agg_interp_type == 4)
1690 {
1691 hypre_BoomerAMGCorrectCFMarker(CFN_marker, CF2_marker);
1692 hypre_IntArrayDestroy(CF2_marker);
1693 CF2_marker = NULL;
1694
1695 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1696 num_functions, nodal, keep_same_sign,
1697 &dof_func, &(CF_marker_array[level]), &S);
1698 hypre_IntArrayDestroy(CFN_marker);
1699 CFN_marker = NULL;
1700 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1701 num_functions, dof_func_array[level], CF_marker_array[level],
1702 &coarse_dof_func, &coarse_pnts_global);
1703 hypre_BoomerAMGBuildMultipass(A_array[level],
1704 CF_marker, S, coarse_pnts_global,
1705 num_functions, dof_func_data, debug_flag,
1706 agg_trunc_factor, agg_P_max_elmts, sep_weight,
1707 &P);
1708 }
1709 else if (agg_interp_type == 8)
1710 {
1711 hypre_BoomerAMGCorrectCFMarker(CFN_marker, CF2_marker);
1712 hypre_IntArrayDestroy(CF2_marker);
1713 CF2_marker = NULL;
1714
1715 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1716 num_functions, nodal, keep_same_sign,
1717 &dof_func, &(CF_marker_array[level]), &S);
1718 hypre_IntArrayDestroy(CFN_marker);
1719 CFN_marker = NULL;
1720 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1721 num_functions, dof_func_array[level], CF_marker_array[level],
1722 &coarse_dof_func, &coarse_pnts_global);
1723 hypre_BoomerAMGBuildModMultipass(A_array[level],
1724 CF_marker, S, coarse_pnts_global,
1725 agg_trunc_factor, agg_P_max_elmts, 8,
1726 num_functions, dof_func_data, &P);
1727 }
1728 else if (agg_interp_type == 9)
1729 {
1730 hypre_BoomerAMGCorrectCFMarker(CFN_marker, CF2_marker);
1731 hypre_IntArrayDestroy(CF2_marker);
1732 CF2_marker = NULL;
1733
1734 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1735 num_functions, nodal, keep_same_sign,
1736 &dof_func, &(CF_marker_array[level]), &S);
1737 hypre_IntArrayDestroy(CFN_marker);
1738 CFN_marker = NULL;
1739 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1740 num_functions, dof_func_array[level], CF_marker_array[level],
1741 &coarse_dof_func,&coarse_pnts_global);
1742 hypre_BoomerAMGBuildModMultipass(A_array[level],
1743 CF_marker, S, coarse_pnts_global,
1744 agg_trunc_factor, agg_P_max_elmts, 9,
1745 num_functions, dof_func_data, &P);
1746 }
1747 else
1748 {
1749 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1750 num_functions, nodal, keep_same_sign,
1751 &dof_func, &CF3_marker, &S);
1752 for (i = 0; i < 2; i++)
1753 {
1754 coarse_pnts_global1[i] *= num_functions;
1755 }
1756 if (agg_interp_type == 1)
1757 {
1758 hypre_BoomerAMGBuildExtPIInterp(A_array[level],
1759 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1760 num_functions, dof_func_data, debug_flag,
1761 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1762 }
1763 else if (agg_interp_type == 2)
1764 {
1765 hypre_BoomerAMGBuildStdInterp(A_array[level],
1766 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1767 num_functions, dof_func_data, debug_flag,
1768 agg_P12_trunc_factor, agg_P12_max_elmts, 0, &P1);
1769 }
1770 else if (agg_interp_type == 3)
1771 {
1772 hypre_BoomerAMGBuildExtInterp(A_array[level],
1773 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1774 num_functions, dof_func_data, debug_flag,
1775 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1776 }
1777 else if (agg_interp_type == 5)
1778 {
1779 hypre_BoomerAMGBuildModExtInterp(A_array[level],
1780 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1781 num_functions, dof_func_data,
1782 debug_flag,
1783 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1784 }
1785 else if (agg_interp_type == 6)
1786 {
1787 hypre_BoomerAMGBuildModExtPIInterp(A_array[level],
1788 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1789 num_functions, dof_func_data,
1790 debug_flag,
1791 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1792 }
1793 else if (agg_interp_type == 7 )
1794 {
1795 hypre_BoomerAMGBuildModExtPEInterp(A_array[level],
1796 hypre_IntArrayData(CF3_marker), S, coarse_pnts_global1,
1797 num_functions, dof_func_data,
1798 debug_flag,
1799 agg_P12_trunc_factor, agg_P12_max_elmts, &P1);
1800 }
1801
1802 hypre_BoomerAMGCorrectCFMarker2 (CFN_marker, CF2_marker);
1803 hypre_IntArrayDestroy(CF2_marker);
1804 CF2_marker = NULL;
1805 hypre_IntArrayDestroy(CF3_marker);
1806 CF3_marker = NULL;
1807 hypre_ParCSRMatrixDestroy(S);
1808 hypre_BoomerAMGCreateScalarCFS(SN, A_array[level], hypre_IntArrayData(CFN_marker),
1809 num_functions, nodal, keep_same_sign,
1810 &dof_func, &(CF_marker_array[level]), &S);
1811
1812 hypre_IntArrayDestroy(CFN_marker);
1813 CFN_marker = NULL;
1814 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1815 num_functions, dof_func_array[level], CF_marker_array[level],
1816 &coarse_dof_func,&coarse_pnts_global);
1817 if (agg_interp_type == 1 || agg_interp_type == 6)
1818 {
1819 hypre_BoomerAMGBuildPartialExtPIInterp(A_array[level],
1820 CF_marker, S, coarse_pnts_global,
1821 coarse_pnts_global1, num_functions,
1822 dof_func_data, debug_flag, agg_P12_trunc_factor,
1823 agg_P12_max_elmts, &P2);
1824 }
1825 else if (agg_interp_type == 2)
1826 {
1827 hypre_BoomerAMGBuildPartialStdInterp(A_array[level],
1828 CF_marker, S, coarse_pnts_global,
1829 coarse_pnts_global1, num_functions,
1830 dof_func_data, debug_flag, agg_P12_trunc_factor,
1831 agg_P12_max_elmts, sep_weight, &P2);
1832 }
1833 else if (agg_interp_type == 3)
1834 {
1835 hypre_BoomerAMGBuildPartialExtInterp(A_array[level],
1836 CF_marker, S, coarse_pnts_global,
1837 coarse_pnts_global1, num_functions,
1838 dof_func_data, debug_flag, agg_P12_trunc_factor,
1839 agg_P12_max_elmts, &P2);
1840 }
1841 else if (agg_interp_type == 5)
1842 {
1843 hypre_BoomerAMGBuildModPartialExtInterp(A_array[level],
1844 CF_marker, S, coarse_pnts_global, coarse_pnts_global1,
1845 num_functions, dof_func_data,
1846 debug_flag,
1847 agg_P12_trunc_factor, agg_P12_max_elmts, &P2);
1848 }
1849 else if (agg_interp_type == 7)
1850 {
1851 hypre_BoomerAMGBuildModPartialExtPEInterp(A_array[level],
1852 CF_marker, S, coarse_pnts_global, coarse_pnts_global1,
1853 num_functions, dof_func_data,
1854 debug_flag,
1855 agg_P12_trunc_factor, agg_P12_max_elmts, &P2);
1856 }
1857
1858 if (hypre_ParAMGDataModularizedMatMat(amg_data))
1859 {
1860 P = hypre_ParCSRMatMat(P1, P2);
1861 }
1862 else
1863 {
1864 P = hypre_ParMatmul(P1, P2);
1865 }
1866 hypre_BoomerAMGInterpTruncation(P, agg_trunc_factor,
1867 agg_P_max_elmts);
1868 hypre_MatvecCommPkgCreate(P);
1869 hypre_ParCSRMatrixDestroy(P1);
1870 hypre_ParCSRMatrixDestroy(P2);
1871 }
1872 if (SN)
1873 {
1874 hypre_ParCSRMatrixDestroy(SN);
1875 }
1876 SN = NULL;
1877 if (AN)
1878 {
1879 hypre_ParCSRMatrixDestroy(AN);
1880 }
1881 AN = NULL;
1882 }
1883 if (my_id == (num_procs -1)) coarse_size = coarse_pnts_global[1];
1884 hypre_MPI_Bcast(&coarse_size, 1, HYPRE_MPI_BIG_INT, num_procs-1, comm);
1885 }
1886 else /* no aggressive coarsening */
1887 {
1888 /**** Get the coarse parameters ****/
1889 /* xxxxxxxxxxxxxxxxxxxxxxxxx change for min_coarse_size
1890 if (block_mode )
1891 {
1892 hypre_BoomerAMGCoarseParms(comm,
1893 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(AN)),
1894 1, NULL, CF_marker, NULL, &coarse_pnts_global);
1895 }
1896 else
1897 {
1898 hypre_BoomerAMGCoarseParms(comm, local_num_vars,
1899 num_functions, dof_func_array[level], CF_marker,
1900 &coarse_dof_func,&coarse_pnts_global);
1901 }
1902 if (my_id == (num_procs -1)) coarse_size = coarse_pnts_global[1];
1903 hypre_MPI_Bcast(&coarse_size, 1, HYPRE_MPI_BIG_INT, num_procs-1, comm);
1904 xxxxxxxxxxxxxxxxxxxxxxxxx change for min_coarse_size */
1905 if (debug_flag==1)
1906 {
1907 wall_time = time_getWallclockSeconds() - wall_time;
1908 hypre_printf("Proc = %d Level = %d Coarsen Time = %f\n",
1909 my_id,level, wall_time);
1910 fflush(NULL);
1911 }
1912
1913 /* RL: build restriction */
1914 if (restri_type)
1915 {
1916 HYPRE_Real filter_thresholdR;
1917 filter_thresholdR = hypre_ParAMGDataFilterThresholdR(amg_data);
1918 HYPRE_Int is_triangular = hypre_ParAMGDataIsTriangular(amg_data);
1919 HYPRE_Int gmres_switch = hypre_ParAMGDataGMRESSwitchR(amg_data);
1920 /* !!! RL: ensure that CF_marker contains -1 or 1 !!! */
1921 #if defined(HYPRE_USING_CUDA) || defined(HYPRE_USING_HIP)
1922 if (hypre_IntArrayMemoryLocation(CF_marker_array[level]) == HYPRE_MEMORY_DEVICE)
1923 {
1924 hypre_BoomerAMGCFMarkerTo1minus1Device(CF_marker, hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level])));
1925 }
1926 else
1927 #endif
1928 {
1929 for (i = 0; i < hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level])); i++)
1930 {
1931 CF_marker[i] = CF_marker[i] > 0 ? 1 : -1;
1932 }
1933 }
1934
1935 if (restri_type == 1) /* distance-1 AIR */
1936 {
1937 hypre_BoomerAMGBuildRestrAIR(A_array[level], CF_marker,
1938 Sabs, coarse_pnts_global, 1, NULL,
1939 filter_thresholdR, debug_flag,
1940 &R,
1941 is_triangular, gmres_switch );
1942 }
1943 else if (restri_type == 2 || restri_type == 15) /* distance-2, 1.5 AIR */
1944 {
1945 hypre_BoomerAMGBuildRestrDist2AIR(A_array[level], CF_marker,
1946 Sabs, coarse_pnts_global, 1, NULL,
1947 filter_thresholdR, debug_flag,
1948 &R, restri_type == 15,
1949 is_triangular, gmres_switch);
1950 }
1951 else
1952 {
1953 HYPRE_Int NeumannAIRDeg = restri_type - 3;
1954 hypre_assert(NeumannAIRDeg >= 0);
1955 HYPRE_Real strong_thresholdR;
1956 strong_thresholdR = hypre_ParAMGDataStrongThresholdR(amg_data);
1957 hypre_BoomerAMGBuildRestrNeumannAIR(A_array[level], CF_marker,
1958 coarse_pnts_global, 1, NULL,
1959 NeumannAIRDeg, strong_thresholdR,
1960 filter_thresholdR, debug_flag,
1961 &R );
1962 }
1963
1964 #if DEBUG_SAVE_ALL_OPS
1965 char file[256];
1966 hypre_sprintf(file, "R_%d.mtx", level);
1967 hypre_ParCSRMatrixPrintIJ(R, 1, 1, file);
1968 #endif
1969 if (Sabs)
1970 {
1971 hypre_ParCSRMatrixDestroy(Sabs);
1972 Sabs = NULL;
1973 }
1974 }
1975
1976 if (debug_flag==1) wall_time = time_getWallclockSeconds();
1977
1978 if (interp_type == 4)
1979 {
1980 hypre_BoomerAMGBuildMultipass(A_array[level], CF_marker,
1981 S, coarse_pnts_global, num_functions, dof_func_data,
1982 debug_flag, trunc_factor, P_max_elmts, sep_weight, &P);
1983 }
1984 else if (interp_type == 1)
1985 {
1986 hypre_BoomerAMGNormalizeVecs(
1987 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level])),
1988 hypre_ParAMGDataNumSamples(amg_data), SmoothVecs);
1989
1990 hypre_BoomerAMGBuildInterpLS(NULL, CF_marker, S,
1991 coarse_pnts_global, num_functions, dof_func_data,
1992 debug_flag, trunc_factor,
1993 hypre_ParAMGDataNumSamples(amg_data), SmoothVecs, &P);
1994 }
1995 else if (interp_type == 2)
1996 {
1997 hypre_BoomerAMGBuildInterpHE(A_array[level], CF_marker,
1998 S, coarse_pnts_global, num_functions, dof_func_data,
1999 debug_flag, trunc_factor, P_max_elmts, &P);
2000 }
2001 else if (interp_type == 3 || interp_type == 15)
2002 {
2003 hypre_BoomerAMGBuildDirInterp(A_array[level], CF_marker,
2004 S, coarse_pnts_global, num_functions, dof_func_data,
2005 debug_flag, trunc_factor, P_max_elmts,
2006 interp_type, &P);
2007 }
2008 else if (interp_type == 6) /*Extended+i classical interpolation */
2009 {
2010 hypre_BoomerAMGBuildExtPIInterp(A_array[level], CF_marker,
2011 S, coarse_pnts_global, num_functions, dof_func_data,
2012 debug_flag, trunc_factor, P_max_elmts, &P);
2013 }
2014 else if (interp_type == 14) /*Extended classical interpolation */
2015 {
2016 hypre_BoomerAMGBuildExtInterp(A_array[level], CF_marker,
2017 S, coarse_pnts_global, num_functions, dof_func_data,
2018 debug_flag, trunc_factor, P_max_elmts, &P);
2019 }
2020 else if (interp_type == 16) /*Extended classical MM interpolation */
2021 {
2022 hypre_BoomerAMGBuildModExtInterp(A_array[level], CF_marker,
2023 S, coarse_pnts_global,
2024 num_functions, dof_func_data,
2025 debug_flag,
2026 trunc_factor, P_max_elmts, &P);
2027 }
2028 else if (interp_type == 17) /*Extended+i MM interpolation */
2029 {
2030 hypre_BoomerAMGBuildModExtPIInterp(A_array[level], CF_marker,
2031 S, coarse_pnts_global,
2032 num_functions, dof_func_data,
2033 debug_flag, trunc_factor, P_max_elmts, &P);
2034 }
2035 else if (interp_type == 18) /*Extended+e MM interpolation */
2036 {
2037 hypre_BoomerAMGBuildModExtPEInterp(A_array[level], CF_marker,
2038 S, coarse_pnts_global,
2039 num_functions, dof_func_data,
2040 debug_flag, trunc_factor, P_max_elmts, &P);
2041 }
2042
2043 else if (interp_type == 7) /*Extended+i (if no common C) interpolation */
2044 {
2045 hypre_BoomerAMGBuildExtPICCInterp(A_array[level], CF_marker,
2046 S, coarse_pnts_global, num_functions, dof_func_data,
2047 debug_flag, trunc_factor, P_max_elmts, &P);
2048 }
2049 else if (interp_type == 12) /*FF interpolation */
2050 {
2051 hypre_BoomerAMGBuildFFInterp(A_array[level], CF_marker,
2052 S, coarse_pnts_global, num_functions, dof_func_data,
2053 debug_flag, trunc_factor, P_max_elmts, &P);
2054 }
2055 else if (interp_type == 13) /*FF1 interpolation */
2056 {
2057 hypre_BoomerAMGBuildFF1Interp(A_array[level], CF_marker,
2058 S, coarse_pnts_global, num_functions, dof_func_data,
2059 debug_flag, trunc_factor, P_max_elmts, &P);
2060 }
2061 else if (interp_type == 8) /*Standard interpolation */
2062 {
2063 hypre_BoomerAMGBuildStdInterp(A_array[level], CF_marker,
2064 S, coarse_pnts_global, num_functions, dof_func_data,
2065 debug_flag, trunc_factor, P_max_elmts, sep_weight, &P);
2066 }
2067 else if (interp_type == 100) /* 1pt interpolation */
2068 {
2069 hypre_BoomerAMGBuildInterpOnePnt(A_array[level], CF_marker, S,
2070 coarse_pnts_global, 1, NULL,
2071 debug_flag, &P);
2072
2073 #if DEBUG_SAVE_ALL_OPS
2074 char file[256];
2075 hypre_sprintf(file, "P_%d.mtx", level);
2076 hypre_ParCSRMatrixPrintIJ(P, 1, 1, file);
2077 #endif
2078 }
2079 else if (hypre_ParAMGDataGSMG(amg_data) == 0) /* none of above choosen and not GMSMG */
2080 {
2081 if (block_mode) /* nodal interpolation */
2082 {
2083
2084 /* convert A to a block matrix if there isn't already a block
2085 matrix - there should be one already*/
2086 if (!(A_block_array[level]))
2087 {
2088 A_block_array[level] = hypre_ParCSRBlockMatrixConvertFromParCSRMatrix(
2089 A_array[level], num_functions);
2090 }
2091
2092 /* note that the current CF_marker is nodal */
2093 if (interp_type == 11)
2094 {
2095 hypre_BoomerAMGBuildBlockInterpDiag( A_block_array[level], CF_marker,
2096 SN,
2097 coarse_pnts_global, 1,
2098 NULL,
2099 debug_flag,
2100 trunc_factor, P_max_elmts,1,
2101 &P_block_array[level]);
2102
2103
2104 }
2105 else if (interp_type == 22)
2106 {
2107 hypre_BoomerAMGBuildBlockInterpRV( A_block_array[level], CF_marker,
2108 SN,
2109 coarse_pnts_global, 1,
2110 NULL,
2111 debug_flag,
2112 trunc_factor, P_max_elmts,
2113 &P_block_array[level]);
2114 }
2115 else if (interp_type == 23)
2116 {
2117 hypre_BoomerAMGBuildBlockInterpRV( A_block_array[level], CF_marker,
2118 SN,
2119 coarse_pnts_global, 1,
2120 NULL,
2121 debug_flag,
2122 trunc_factor, P_max_elmts,
2123 &P_block_array[level]);
2124 }
2125 else if (interp_type == 20)
2126 {
2127 hypre_BoomerAMGBuildBlockInterp( A_block_array[level], CF_marker,
2128 SN,
2129 coarse_pnts_global, 1,
2130 NULL,
2131 debug_flag,
2132 trunc_factor, P_max_elmts, 0,
2133 &P_block_array[level]);
2134
2135 }
2136 else if (interp_type == 21)
2137 {
2138 hypre_BoomerAMGBuildBlockInterpDiag( A_block_array[level], CF_marker,
2139 SN,
2140 coarse_pnts_global, 1,
2141 NULL,
2142 debug_flag,
2143 trunc_factor, P_max_elmts, 0,
2144 &P_block_array[level]);
2145 }
2146 else if (interp_type == 24)
2147 {
2148 hypre_BoomerAMGBuildBlockDirInterp( A_block_array[level], CF_marker,
2149 SN,
2150 coarse_pnts_global, 1,
2151 NULL,
2152 debug_flag,
2153 trunc_factor, P_max_elmts,
2154 &P_block_array[level]);
2155 }
2156
2157 else /* interp_type ==10 */
2158 {
2159
2160 hypre_BoomerAMGBuildBlockInterp( A_block_array[level], CF_marker,
2161 SN,
2162 coarse_pnts_global, 1,
2163 NULL,
2164 debug_flag,
2165 trunc_factor, P_max_elmts, 1,
2166 &P_block_array[level]);
2167
2168 }
2169
2170 /* we need to set the global number of cols in P, as this was
2171 not done in the interp
2172 (which calls the matrix create) since we didn't
2173 have the global partition */
2174 /* this has to be done before converting from block to non-block*/
2175 hypre_ParCSRBlockMatrixGlobalNumCols(P_block_array[level]) = coarse_size;
2176
2177 /* if we don't do nodal relaxation, we need a CF_array that is
2178 not nodal - right now we don't allow this to happen though*/
2179 /*
2180 if (grid_relax_type[0] < 20 )
2181 {
2182 hypre_BoomerAMGCreateScalarCF(CFN_marker, num_functions,
2183 hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(AN)),
2184 &dof_func1, &CF_marker);
2185
2186 dof_func_array[level+1] = dof_func1;
2187 hypre_TFree(CFN_marker, HYPRE_MEMORY_HOST);
2188 CF_marker_array[level] = CF_marker;
2189 }
2190 */
2191
2192 /* clean up other things */
2193 hypre_ParCSRMatrixDestroy(AN);
2194 hypre_ParCSRMatrixDestroy(SN);
2195
2196 }
2197 else /* not block mode - use default interp (interp_type = 0) */
2198 {
2199 if (nodal > -1) /* non-systems, or systems with unknown approach interpolation*/
2200 {
2201 /* if systems, do we want to use an interp. that uses the full strength matrix?*/
2202
2203 if ( (num_functions > 1) && (interp_type == 19 || interp_type == 18 || interp_type == 17 || interp_type == 16))
2204 {
2205 /* so create a second strength matrix and build interp with with num_functions = 1 */
2206 hypre_BoomerAMGCreateS(A_array[level],
2207 strong_threshold, max_row_sum,
2208 1, dof_func_data,&S2);
2209 switch (interp_type)
2210 {
2211
2212 case 19:
2213 dbg_flg = debug_flag;
2214 if (amg_print_level) dbg_flg = -debug_flag;
2215 hypre_BoomerAMGBuildInterp(A_array[level], CF_marker,
2216 S2, coarse_pnts_global, 1,
2217 dof_func_data,
2218 dbg_flg, trunc_factor, P_max_elmts, &P);
2219 break;
2220
2221 case 18:
2222 hypre_BoomerAMGBuildStdInterp(A_array[level], CF_marker,
2223 S2, coarse_pnts_global, 1, dof_func_data,
2224 debug_flag, trunc_factor, P_max_elmts, 0, &P);
2225
2226 break;
2227
2228 case 17:
2229 hypre_BoomerAMGBuildExtPIInterp(A_array[level], CF_marker,
2230 S2, coarse_pnts_global, 1, dof_func_data,
2231 debug_flag, trunc_factor, P_max_elmts, &P);
2232 break;
2233 case 16:
2234 dbg_flg = debug_flag;
2235 if (amg_print_level) dbg_flg = -debug_flag;
2236 hypre_BoomerAMGBuildInterpModUnk(A_array[level], CF_marker,
2237 S2, coarse_pnts_global, num_functions, dof_func_data,
2238 dbg_flg, trunc_factor, P_max_elmts, &P);
2239 break;
2240
2241 }
2242
2243
2244 hypre_ParCSRMatrixDestroy(S2);
2245
2246 }
2247 else /* one function only or unknown-based interpolation- */
2248 {
2249 dbg_flg = debug_flag;
2250 if (amg_print_level) dbg_flg = -debug_flag;
2251
2252 hypre_BoomerAMGBuildInterp(A_array[level], CF_marker,
2253 S, coarse_pnts_global, num_functions,
2254 dof_func_data,
2255 dbg_flg, trunc_factor, P_max_elmts, &P);
2256
2257
2258 }
2259 }
2260 }
2261 }
2262 else
2263 {
2264 hypre_BoomerAMGBuildInterpGSMG(NULL, CF_marker, S,
2265 coarse_pnts_global, num_functions, dof_func_data,
2266 debug_flag, trunc_factor, &P);
2267 }
2268 } /* end of no aggressive coarsening */
2269
2270 dof_func_array[level+1] = NULL;
2271 if (num_functions > 1 && nodal > -1 && (!block_mode) )
2272 dof_func_array[level+1] = coarse_dof_func;
2273
2274 HYPRE_ANNOTATE_REGION_END("%s", "Interpolation");
2275 } /* end of if max_levels > 1 */
2276
2277 /* if no coarse-grid, stop coarsening, and set the
2278 * coarsest solve to be a single sweep of Jacobi */
2279 if ( (coarse_size == 0) || (coarse_size == fine_size) )
2280 {
2281 HYPRE_Int *num_grid_sweeps =
2282 hypre_ParAMGDataNumGridSweeps(amg_data);
2283 HYPRE_Int **grid_relax_points =
2284 hypre_ParAMGDataGridRelaxPoints(amg_data);
2285 if (grid_relax_type[3] == 9 || grid_relax_type[3] == 99
2286 || grid_relax_type[3] == 19 || grid_relax_type[3] == 98)
2287 {
2288 grid_relax_type[3] = grid_relax_type[0];
2289 num_grid_sweeps[3] = 1;
2290 if (grid_relax_points) grid_relax_points[3][0] = 0;
2291 }
2292 if (S)
2293 hypre_ParCSRMatrixDestroy(S);
2294 if (P)
2295 hypre_ParCSRMatrixDestroy(P);
2296 if (level > 0)
2297 {
2298 /* note special case treatment of CF_marker is necessary
2299 * to do CF relaxation correctly when num_levels = 1 */
2300 hypre_IntArrayDestroy(CF_marker_array[level]);
2301 CF_marker_array[level] = NULL;
2302 hypre_ParVectorDestroy(F_array[level]);
2303 hypre_ParVectorDestroy(U_array[level]);
2304 }
2305 hypre_IntArrayDestroy(dof_func_array[level+1]);
2306 dof_func_array[level+1] = NULL;
2307
2308 break;
2309 }
2310 if (level < agg_num_levels && coarse_size < min_coarse_size)
2311 {
2312 if (S)
2313 hypre_ParCSRMatrixDestroy(S);
2314 if (P)
2315 hypre_ParCSRMatrixDestroy(P);
2316 if (level > 0)
2317 {
2318 hypre_IntArrayDestroy(CF_marker_array[level]);
2319 CF_marker_array[level] = NULL;
2320 hypre_ParVectorDestroy(F_array[level]);
2321 hypre_ParVectorDestroy(U_array[level]);
2322 }
2323 hypre_IntArrayDestroy(dof_func_array[level+1]);
2324 dof_func_array[level+1] = NULL;
2325 coarse_size = fine_size;
2326
2327 break;
2328 }
2329
2330 /*-------------------------------------------------------------
2331 * Build prolongation matrix, P, and place in P_array[level]
2332 *--------------------------------------------------------------*/
2333
2334 if (interp_refine > 0)
2335 {
2336 for (k = 0; k < interp_refine; k++)
2337 hypre_BoomerAMGRefineInterp(A_array[level],
2338 P,
2339 coarse_pnts_global,
2340 &num_functions,
2341 dof_func_data,
2342 hypre_IntArrayData(CF_marker_array[level]),
2343 level);
2344 }
2345
2346 /* Post processing of interpolation operators to incorporate
2347 smooth vectors NOTE: must pick nodal coarsening !!!
2348 (nodal is changed above to 1 if it is 0) */
2349 if (interp_vec_variant && nodal && num_interp_vectors)
2350 {
2351 /* TO DO: add option of smoothing the vectors at
2352 * coarser levels?*/
2353
2354 if (level < interp_vec_first_level)
2355 {
2356 /* coarsen the smooth vecs */
2357 hypre_BoomerAMGCoarsenInterpVectors( P,
2358 num_interp_vectors,
2359 interp_vectors_array[level],
2360 hypre_IntArrayData(CF_marker_array[level]),
2361 &interp_vectors_array[level+1],
2362 0, num_functions);
2363
2364 }
2365 /* do GM 2 and LN (3) at all levels and GM 1 only on first level */
2366 if (( interp_vec_variant > 1 && level >= interp_vec_first_level) ||
2367 (interp_vec_variant == 1 && interp_vec_first_level == level))
2368
2369 {
2370 /*if (level == 0)
2371 {
2372 hypre_ParCSRMatrixPrintIJ(A_array[0], 0, 0, "A");
2373 hypre_ParVectorPrintIJ(interp_vectors_array[0][0], 0, "rbm");
2374 }*/
2375 if (interp_vec_variant < 3) /* GM */
2376 {
2377 hypre_BoomerAMG_GMExpandInterp( A_array[level],
2378 &P,
2379 num_interp_vectors,
2380 interp_vectors_array[level],
2381 &num_functions,
2382 dof_func_data,
2383 &dof_func_array[level+1],
2384 interp_vec_variant, level,
2385 abs_q_trunc,
2386 expandp_weights,
2387 q_max,
2388 hypre_IntArrayData(CF_marker_array[level]),
2389 interp_vec_first_level);
2390 }
2391 else /* LN */
2392 {
2393 hypre_BoomerAMG_LNExpandInterp( A_array[level],
2394 &P,
2395 coarse_pnts_global,
2396 &num_functions,
2397 dof_func_data,
2398 &dof_func_array[level+1],
2399 hypre_IntArrayData(CF_marker_array[level]),
2400 level,
2401 expandp_weights,
2402 num_interp_vectors,
2403 interp_vectors_array[level],
2404 abs_q_trunc,
2405 q_max,
2406 interp_vec_first_level);
2407 }
2408
2409 if (level == interp_vec_first_level)
2410 {
2411 /* check to see if we made A bigger - this can happen
2412 * in 3D with certain coarsenings - if so, need to fix vtemp*/
2413
2414 HYPRE_Int local_sz = hypre_ParVectorActualLocalSize(Vtemp);
2415 HYPRE_Int local_P_sz = hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(P));
2416 if (local_sz < local_P_sz)
2417 {
2418 hypre_Vector *Vtemp_local = hypre_ParVectorLocalVector(Vtemp);
2419 hypre_TFree(hypre_VectorData(Vtemp_local), memory_location);
2420 hypre_VectorSize(Vtemp_local) = local_P_sz;
2421 hypre_VectorData(Vtemp_local) = hypre_CTAlloc(HYPRE_Complex, local_P_sz, memory_location);
2422 if (Ztemp)
2423 {
2424 hypre_Vector *Ztemp_local = hypre_ParVectorLocalVector(Ztemp);
2425 hypre_TFree(hypre_VectorData(Ztemp_local), memory_location);
2426 hypre_VectorSize(Ztemp_local) = local_P_sz;
2427 hypre_VectorData(Ztemp_local) = hypre_CTAlloc(HYPRE_Complex, local_P_sz, memory_location);
2428 }
2429 if (Ptemp)
2430 {
2431 hypre_Vector *Ptemp_local = hypre_ParVectorLocalVector(Ptemp);
2432 hypre_TFree(hypre_VectorData(Ptemp_local), memory_location);
2433 hypre_VectorSize(Ptemp_local) = local_P_sz;
2434 hypre_VectorData(Ptemp_local) = hypre_CTAlloc(HYPRE_Complex, local_P_sz, memory_location);
2435 }
2436 if (Rtemp)
2437 {
2438 hypre_Vector *Rtemp_local = hypre_ParVectorLocalVector(Rtemp);
2439 hypre_TFree(hypre_VectorData(Rtemp_local), memory_location);
2440 hypre_VectorSize(Rtemp_local) = local_P_sz;
2441 hypre_VectorData(Rtemp_local) = hypre_CTAlloc(HYPRE_Complex, local_P_sz, memory_location);
2442 }
2443 }
2444 /*if (hypre_ParCSRMatrixGlobalNumRows(A_array[0]) < hypre_ParCSRMatrixGlobalNumCols(P))
2445 {
2446
2447 hypre_ParVectorDestroy(Vtemp);
2448 Vtemp = NULL;
2449
2450 Vtemp = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(P),
2451 hypre_ParCSRMatrixGlobalNumCols(P),
2452 hypre_ParCSRMatrixColStarts(P));
2453 hypre_ParVectorInitialize(Vtemp);
2454 hypre_ParAMGDataVtemp(amg_data) = Vtemp;
2455 }*/
2456 }
2457 /* at the first level we have to add space for the new
2458 * unknowns in the smooth vectors */
2459 if (interp_vec_variant > 1 && level < max_levels)
2460 {
2461 HYPRE_Int expand_level = 0;
2462
2463 if (level == interp_vec_first_level)
2464 expand_level = 1;
2465
2466 hypre_BoomerAMGCoarsenInterpVectors( P,
2467 num_interp_vectors,
2468 interp_vectors_array[level],
2469 hypre_IntArrayData(CF_marker_array[level]),
2470 &interp_vectors_array[level+1],
2471 expand_level, num_functions);
2472 }
2473 } /* end apply variant */
2474 }/* end interp_vec_variant > 0 */
2475
2476 /* Improve on P with Jacobi interpolation */
2477 for (i = 0; i < post_interp_type; i++)
2478 {
2479 hypre_BoomerAMGJacobiInterp( A_array[level], &P, S,
2480 num_functions, dof_func_data,
2481 hypre_IntArrayData(CF_marker_array[level]),
2482 level, jacobi_trunc_threshold, 0.5*jacobi_trunc_threshold );
2483 }
2484
2485 dof_func_data = NULL;
2486 if (dof_func_array[level+1])
2487 {
2488 dof_func_data = hypre_IntArrayData(dof_func_array[level+1]);
2489
2490 }
2491
2492 if (!block_mode)
2493 {
2494 if (mult_addlvl > -1 && level >= mult_addlvl && level <= add_end)
2495 {
2496 hypre_Vector *d_diag = NULL;
2497
2498 if (ns == 1)
2499 {
2500 d_diag = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[level]));
2501
2502 if (add_rlx == 0)
2503 {
2504 hypre_CSRMatrix *lvl_Adiag = hypre_ParCSRMatrixDiag(A_array[level]);
2505 HYPRE_Int lvl_nrows = hypre_CSRMatrixNumRows(lvl_Adiag);
2506 HYPRE_Int *lvl_i = hypre_CSRMatrixI(lvl_Adiag);
2507 HYPRE_Real *lvl_data = hypre_CSRMatrixData(lvl_Adiag);
2508 HYPRE_Real w_inv = 1.0 / add_rlx_wt;
2509 /*HYPRE_Real w_inv = 1.0/hypre_ParAMGDataRelaxWeight(amg_data)[level];*/
2510 hypre_SeqVectorInitialize_v2(d_diag, HYPRE_MEMORY_HOST);
2511 for (i=0; i < lvl_nrows; i++)
2512 {
2513 hypre_VectorData(d_diag)[i] = lvl_data[lvl_i[i]] * w_inv;
2514 }
2515 }
2516 else
2517 {
2518 HYPRE_Real *d_diag_data = NULL;
2519
2520 hypre_ParCSRComputeL1Norms(A_array[level], 1, NULL, &d_diag_data);
2521
2522 hypre_VectorData(d_diag) = d_diag_data;
2523 hypre_SeqVectorInitialize_v2(d_diag, hypre_ParCSRMatrixMemoryLocation(A_array[level]));
2524 }
2525 }
2526
2527 HYPRE_ANNOTATE_REGION_BEGIN("%s", "RAP");
2528 if (ns == 1)
2529 {
2530 hypre_ParCSRMatrix *Q = NULL;
2531 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2532 {
2533 Q = hypre_ParCSRMatMat(A_array[level], P);
2534 hypre_ParCSRMatrixAminvDB(P, Q, hypre_VectorData(d_diag), &P_array[level]);
2535 A_H = hypre_ParCSRTMatMat(P, Q);
2536 }
2537 else
2538 {
2539 Q = hypre_ParMatmul(A_array[level], P);
2540 hypre_ParCSRMatrixAminvDB(P, Q, hypre_VectorData(d_diag), &P_array[level]);
2541 A_H = hypre_ParTMatmul(P, Q);
2542 }
2543 if (num_procs > 1)
2544 {
2545 hypre_MatvecCommPkgCreate(A_H);
2546 }
2547 /*hypre_ParCSRMatrixDestroy(P); */
2548 hypre_SeqVectorDestroy(d_diag);
2549 /* Set NonGalerkin drop tol on each level */
2550 if (level < nongalerk_num_tol) nongalerk_tol_l = nongalerk_tol[level];
2551 if (nongal_tol_array) nongalerk_tol_l = nongal_tol_array[level];
2552 if (nongalerk_tol_l > 0.0)
2553 {
2554 /* Build Non-Galerkin Coarse Grid */
2555 hypre_ParCSRMatrix *Q = NULL;
2556 hypre_BoomerAMGBuildNonGalerkinCoarseOperator(&A_H, Q,
2557 0.333*strong_threshold, max_row_sum, num_functions,
2558 dof_func_data, hypre_IntArrayData(CF_marker_array[level]),
2559 /* nongalerk_tol, sym_collapse, lump_percent, beta );*/
2560 nongalerk_tol_l, 1, 0.5, 1.0 );
2561
2562 hypre_ParCSRMatrixColStarts(P_array[level])[0] = hypre_ParCSRMatrixRowStarts(A_H)[0];
2563 hypre_ParCSRMatrixColStarts(P_array[level])[1] = hypre_ParCSRMatrixRowStarts(A_H)[1];
2564 if (!hypre_ParCSRMatrixCommPkg(A_H))
2565 hypre_MatvecCommPkgCreate(A_H);
2566 }
2567 hypre_ParCSRMatrixDestroy(Q);
2568 }
2569 else
2570 {
2571 HYPRE_Int ns_tmp = ns;
2572 hypre_ParCSRMatrix *C = NULL;
2573 hypre_ParCSRMatrix *Ptmp = NULL;
2574 /* Set NonGalerkin drop tol on each level */
2575 if (level < nongalerk_num_tol)
2576 nongalerk_tol_l = nongalerk_tol[level];
2577 if (nongal_tol_array) nongalerk_tol_l = nongal_tol_array[level];
2578
2579 if (nongalerk_tol_l > 0.0)
2580 {
2581 /* Construct AP, and then RAP */
2582 hypre_ParCSRMatrix *Q = NULL;
2583 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2584 {
2585 Q = hypre_ParCSRMatMat(A_array[level], P);
2586 A_H = hypre_ParCSRTMatMatKT(P, Q, keepTranspose);
2587 }
2588 else
2589 {
2590 Q = hypre_ParMatmul(A_array[level], P);
2591 A_H = hypre_ParTMatmul(P, Q);
2592 }
2593 if (num_procs > 1)
2594 {
2595 hypre_MatvecCommPkgCreate(A_H);
2596 }
2597
2598 /* Build Non-Galerkin Coarse Grid */
2599 hypre_BoomerAMGBuildNonGalerkinCoarseOperator(&A_H, Q,
2600 0.333*strong_threshold, max_row_sum, num_functions,
2601 dof_func_data, hypre_IntArrayData(CF_marker_array[level]),
2602 /* nongalerk_tol, sym_collapse, lump_percent, beta );*/
2603 nongalerk_tol_l, 1, 0.5, 1.0 );
2604
2605 if (!hypre_ParCSRMatrixCommPkg(A_H))
2606 {
2607 hypre_MatvecCommPkgCreate(A_H);
2608 }
2609
2610 /* Delete AP */
2611 hypre_ParCSRMatrixDestroy(Q);
2612 }
2613 else if (rap2)
2614 {
2615 /* Use two matrix products to generate A_H */
2616 hypre_ParCSRMatrix *Q = NULL;
2617 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2618 {
2619 Q = hypre_ParCSRMatMat(A_array[level], P);
2620 A_H = hypre_ParCSRTMatMatKT(P, Q, keepTranspose);
2621 }
2622 else
2623 {
2624 Q = hypre_ParMatmul(A_array[level], P);
2625 A_H = hypre_ParTMatmul(P, Q);
2626 }
2627
2628 if (num_procs > 1)
2629 {
2630 hypre_MatvecCommPkgCreate(A_H);
2631 }
2632 /* Delete AP */
2633 hypre_ParCSRMatrixDestroy(Q);
2634 }
2635 else
2636 {
2637 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2638 {
2639 A_H = hypre_ParCSRMatrixRAPKT(P, A_array[level],
2640 P, keepTranspose);
2641 }
2642 else
2643 {
2644 hypre_BoomerAMGBuildCoarseOperatorKT(P, A_array[level] , P,
2645 keepTranspose, &A_H);
2646 }
2647 }
2648
2649 if (add_rlx == 18)
2650 {
2651 C = hypre_CreateC(A_array[level], 0.0);
2652 }
2653 else
2654 {
2655 C = hypre_CreateC(A_array[level], add_rlx_wt);
2656 }
2657
2658 Ptmp = P;
2659 while (ns_tmp > 0)
2660 {
2661 Pnew = Ptmp;
2662 Ptmp = NULL;
2663 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2664 {
2665 Ptmp = hypre_ParCSRMatMat(C,Pnew);
2666 }
2667 else
2668 {
2669 Ptmp = hypre_ParMatmul(C,Pnew);
2670 }
2671 if (ns_tmp < ns)
2672 {
2673 hypre_ParCSRMatrixDestroy(Pnew);
2674 }
2675 ns_tmp--;
2676 }
2677 Pnew = Ptmp;
2678 P_array[level] = Pnew;
2679 hypre_ParCSRMatrixDestroy(C);
2680 } /* if (ns == 1) */
2681 HYPRE_ANNOTATE_REGION_END("%s", "RAP");
2682
2683 if (add_P_max_elmts || add_trunc_factor)
2684 {
2685 hypre_BoomerAMGTruncandBuild(P_array[level], add_trunc_factor,add_P_max_elmts);
2686 }
2687 /*else
2688 hypre_MatvecCommPkgCreate(P_array[level]); */
2689 hypre_ParCSRMatrixDestroy(P);
2690 }
2691 else
2692 {
2693 P_array[level] = P;
2694 /* RL: save R matrix */
2695 if (restri_type)
2696 {
2697 R_array[level] = R;
2698 }
2699 }
2700 }
2701
2702 if (S)
2703 {
2704 hypre_ParCSRMatrixDestroy(S);
2705 }
2706 S = NULL;
2707
2708 hypre_TFree(SmoothVecs, HYPRE_MEMORY_HOST);
2709
2710 if (debug_flag==1)
2711 {
2712 wall_time = time_getWallclockSeconds() - wall_time;
2713 hypre_printf("Proc = %d Level = %d Build Interp Time = %f\n",
2714 my_id,level, wall_time);
2715 fflush(NULL);
2716 }
2717
2718 /*-------------------------------------------------------------
2719 * Build coarse-grid operator, A_array[level+1] by R*A*P
2720 *--------------------------------------------------------------*/
2721
2722 HYPRE_ANNOTATE_REGION_BEGIN("%s", "RAP");
2723 if (debug_flag==1) wall_time = time_getWallclockSeconds();
2724
2725 if (block_mode)
2726 {
2727 hypre_ParCSRBlockMatrixRAP(P_block_array[level],
2728 A_block_array[level],
2729 P_block_array[level], &A_H_block);
2730
2731 hypre_ParCSRBlockMatrixSetNumNonzeros(A_H_block);
2732 hypre_ParCSRBlockMatrixSetDNumNonzeros(A_H_block);
2733 A_block_array[level+1] = A_H_block;
2734 }
2735 else if (mult_addlvl == -1 || level < mult_addlvl || level > add_end)
2736 {
2737 /* Set NonGalerkin drop tol on each level */
2738 if (level < nongalerk_num_tol)
2739 {
2740 nongalerk_tol_l = nongalerk_tol[level];
2741 }
2742 if (nongal_tol_array)
2743 {
2744 nongalerk_tol_l = nongal_tol_array[level];
2745 }
2746
2747 if (nongalerk_tol_l > 0.0)
2748 {
2749 /* Construct AP, and then RAP */
2750 hypre_ParCSRMatrix *Q = NULL;
2751 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2752 {
2753 Q = hypre_ParCSRMatMat(A_array[level],P_array[level]);
2754 A_H = hypre_ParCSRTMatMatKT(P_array[level],Q,keepTranspose);
2755 }
2756 else
2757 {
2758 Q = hypre_ParMatmul(A_array[level],P_array[level]);
2759 A_H = hypre_ParTMatmul(P_array[level],Q);
2760 }
2761 if (num_procs > 1) hypre_MatvecCommPkgCreate(A_H);
2762
2763 /* Build Non-Galerkin Coarse Grid */
2764 hypre_BoomerAMGBuildNonGalerkinCoarseOperator(&A_H, Q,
2765 0.333*strong_threshold, max_row_sum, num_functions,
2766 dof_func_data, hypre_IntArrayData(CF_marker_array[level]),
2767 /* nongalerk_tol, sym_collapse, lump_percent, beta );*/
2768 nongalerk_tol_l, 1, 0.5, 1.0 );
2769
2770 if (!hypre_ParCSRMatrixCommPkg(A_H))
2771 {
2772 hypre_MatvecCommPkgCreate(A_H);
2773 }
2774 /* Delete AP */
2775 hypre_ParCSRMatrixDestroy(Q);
2776 }
2777 else if (restri_type) /* RL: */
2778 {
2779 /* Use two matrix products to generate A_H */
2780 hypre_ParCSRMatrix *AP = NULL;
2781 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2782 {
2783 AP = hypre_ParCSRMatMat(A_array[level], P_array[level]);
2784 A_H = hypre_ParCSRMatMat(R_array[level], AP);
2785 hypre_CSRMatrixReorder(hypre_ParCSRMatrixDiag(A_H));
2786 }
2787 else
2788 {
2789 AP = hypre_ParMatmul(A_array[level], P_array[level]);
2790 A_H = hypre_ParMatmul(R_array[level], AP);
2791 }
2792 if (num_procs > 1)
2793 {
2794 hypre_MatvecCommPkgCreate(A_H);
2795 }
2796 /* Delete AP */
2797 hypre_ParCSRMatrixDestroy(AP);
2798
2799 #if DEBUG_SAVE_ALL_OPS
2800 if (level == 0)
2801 {
2802 hypre_ParCSRMatrixPrintIJ(A_array[0], 1, 1, "A_0.mtx");
2803 }
2804 char file[256];
2805 hypre_sprintf(file, "A_%d.mtx", level+1);
2806 hypre_ParCSRMatrixPrintIJ(A_H, 1, 1, file);
2807 #endif
2808 }
2809 else if (rap2)
2810 {
2811 /* Use two matrix products to generate A_H */
2812 hypre_ParCSRMatrix *Q = NULL;
2813 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2814 {
2815 Q = hypre_ParCSRMatMat(A_array[level], P_array[level]);
2816 A_H = hypre_ParCSRTMatMatKT(P_array[level], Q, keepTranspose);
2817 }
2818 else
2819 {
2820 Q = hypre_ParMatmul(A_array[level],P_array[level]);
2821 A_H = hypre_ParTMatmul(P_array[level],Q);
2822 }
2823 if (num_procs > 1)
2824 {
2825 hypre_MatvecCommPkgCreate(A_H);
2826 }
2827 /* Delete AP */
2828 hypre_ParCSRMatrixDestroy(Q);
2829 }
2830 else
2831 {
2832 /* Compute standard Galerkin coarse-grid product */
2833 if (hypre_ParAMGDataModularizedMatMat(amg_data))
2834 {
2835 A_H = hypre_ParCSRMatrixRAPKT(P_array[level], A_array[level],
2836 P_array[level], keepTranspose);
2837 }
2838 else
2839 {
2840 hypre_BoomerAMGBuildCoarseOperatorKT(P_array[level], A_array[level] ,
2841 P_array[level], keepTranspose, &A_H);
2842 }
2843
2844 if (Pnew && ns==1)
2845 {
2846 hypre_ParCSRMatrixDestroy(P);
2847 P_array[level] = Pnew;
2848 }
2849 }
2850 }
2851
2852 HYPRE_ANNOTATE_REGION_END("%s", "RAP");
2853 if (debug_flag==1)
2854 {
2855 wall_time = time_getWallclockSeconds() - wall_time;
2856 hypre_printf("Proc = %d Level = %d Build Coarse Operator Time = %f\n",
2857 my_id,level, wall_time);
2858 fflush(NULL);
2859 }
2860
2861 HYPRE_ANNOTATE_MGLEVEL_END(level);
2862 ++level;
2863 HYPRE_ANNOTATE_MGLEVEL_BEGIN(level);
2864
2865 if (!block_mode)
2866 {
2867 /* dropping in A_H */
2868 hypre_ParCSRMatrixDropSmallEntries(A_H, hypre_ParAMGDataADropTol(amg_data),
2869 hypre_ParAMGDataADropType(amg_data));
2870 /* if CommPkg for A_H was not built */
2871 if (num_procs > 1 && hypre_ParCSRMatrixCommPkg(A_H) == NULL)
2872 {
2873 hypre_MatvecCommPkgCreate(A_H);
2874 }
2875 /* NumNonzeros was set in hypre_ParCSRMatrixDropSmallEntries */
2876 if (hypre_ParAMGDataADropTol(amg_data) <= 0.0)
2877 {
2878 hypre_ParCSRMatrixSetNumNonzeros(A_H);
2879 hypre_ParCSRMatrixSetDNumNonzeros(A_H);
2880 }
2881 A_array[level] = A_H;
2882 }
2883
2884 size = ((HYPRE_Real) fine_size )*.75;
2885 if (coarsen_type > 0 && coarse_size >= (HYPRE_BigInt) size)
2886 {
2887 coarsen_type = 0;
2888 }
2889
2890 {
2891 HYPRE_Int max_thresh = hypre_max(coarse_threshold, seq_threshold);
2892 #ifdef HYPRE_USING_DSUPERLU
2893 max_thresh = hypre_max(max_thresh, dslu_threshold);
2894 #endif
2895 if ( (level == max_levels-1) || (coarse_size <= (HYPRE_BigInt) max_thresh) )
2896 {
2897 not_finished_coarsening = 0;
2898 }
2899 }
2900 } /* end of coarsening loop: while (not_finished_coarsening) */
2901 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Coarse solve");
2902
2903 /* free memory */
2904 hypre_TFree(coarse_pnts_global, HYPRE_MEMORY_HOST);
2905 hypre_TFree(coarse_pnts_global1, HYPRE_MEMORY_HOST);
2906
2907 /* redundant coarse grid solve */
2908 if ((seq_threshold >= coarse_threshold) &&
2909 (coarse_size > (HYPRE_BigInt) coarse_threshold) &&
2910 (level != max_levels-1))
2911 {
2912 hypre_seqAMGSetup(amg_data, level, coarse_threshold);
2913 }
2914 #ifdef HYPRE_USING_DSUPERLU
2915 else if ((dslu_threshold >= coarse_threshold) &&
2916 (coarse_size > (HYPRE_BigInt)coarse_threshold) &&
2917 (level != max_levels-1))
2918 {
2919 HYPRE_Solver dslu_solver;
2920 hypre_SLUDistSetup(&dslu_solver, A_array[level], amg_print_level);
2921 hypre_ParAMGDataDSLUSolver(amg_data) = dslu_solver;
2922 }
2923 #endif
2924 else if (grid_relax_type[3] == 9 ||
2925 grid_relax_type[3] == 99 ||
2926 grid_relax_type[3] == 199 ) /*use of Gaussian elimination on coarsest level */
2927 {
2928 if (coarse_size <= coarse_threshold)
2929 {
2930 hypre_GaussElimSetup(amg_data, level, grid_relax_type[3]);
2931 }
2932 else
2933 {
2934 grid_relax_type[3] = grid_relax_type[1];
2935 }
2936 }
2937 else if (grid_relax_type[3] == 19 || grid_relax_type[3] == 98) /*use of Gaussian elimination on coarsest level */
2938 {
2939 if (coarse_size > coarse_threshold)
2940 {
2941 grid_relax_type[3] = grid_relax_type[1];
2942 }
2943 }
2944 HYPRE_ANNOTATE_REGION_END("%s", "Coarse solve");
2945 HYPRE_ANNOTATE_MGLEVEL_END(level);
2946
2947 if (level > 0)
2948 {
2949 if (block_mode)
2950 {
2951 F_array[level] =
2952 hypre_ParVectorCreateFromBlock(hypre_ParCSRBlockMatrixComm(A_block_array[level]),
2953 hypre_ParCSRMatrixGlobalNumRows(A_block_array[level]),
2954 hypre_ParCSRBlockMatrixRowStarts(A_block_array[level]),
2955 hypre_ParCSRBlockMatrixBlockSize(A_block_array[level]));
2956 hypre_ParVectorInitialize(F_array[level]);
2957
2958 U_array[level] =
2959 hypre_ParVectorCreateFromBlock(hypre_ParCSRBlockMatrixComm(A_block_array[level]),
2960 hypre_ParCSRMatrixGlobalNumRows(A_block_array[level]),
2961 hypre_ParCSRBlockMatrixRowStarts(A_block_array[level]),
2962 hypre_ParCSRBlockMatrixBlockSize(A_block_array[level]));
2963
2964 hypre_ParVectorInitialize(U_array[level]);
2965 }
2966 else
2967 {
2968 F_array[level] =
2969 hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
2970 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
2971 hypre_ParCSRMatrixRowStarts(A_array[level]));
2972 hypre_ParVectorInitialize_v2(F_array[level], memory_location);
2973
2974 U_array[level] =
2975 hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
2976 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
2977 hypre_ParCSRMatrixRowStarts(A_array[level]));
2978 hypre_ParVectorInitialize_v2(U_array[level], memory_location);
2979 }
2980 }
2981
2982 /*-----------------------------------------------------------------------
2983 * enter all the stuff created, A[level], P[level], CF_marker[level],
2984 * for levels 1 through coarsest, into amg_data data structure
2985 *-----------------------------------------------------------------------*/
2986
2987 num_levels = level+1;
2988 hypre_ParAMGDataNumLevels(amg_data) = num_levels;
2989 if (hypre_ParAMGDataSmoothNumLevels(amg_data) > num_levels-1)
2990 {
2991 hypre_ParAMGDataSmoothNumLevels(amg_data) = num_levels;
2992 }
2993 smooth_num_levels = hypre_ParAMGDataSmoothNumLevels(amg_data);
2994
2995 /*-----------------------------------------------------------------------
2996 * Setup of special smoothers when needed
2997 *-----------------------------------------------------------------------*/
2998
2999 if (addlvl > -1 ||
3000 grid_relax_type[1] == 7 || grid_relax_type[2] == 7 || grid_relax_type[3] == 7 ||
3001 grid_relax_type[1] == 8 || grid_relax_type[2] == 8 || grid_relax_type[3] == 8 ||
3002 grid_relax_type[1] == 13 || grid_relax_type[2] == 13 || grid_relax_type[3] == 13 ||
3003 grid_relax_type[1] == 14 || grid_relax_type[2] == 14 || grid_relax_type[3] == 14 ||
3004 grid_relax_type[1] == 18 || grid_relax_type[2] == 18 || grid_relax_type[3] == 18)
3005 {
3006 l1_norms = hypre_CTAlloc(hypre_Vector*, num_levels, HYPRE_MEMORY_HOST);
3007 hypre_ParAMGDataL1Norms(amg_data) = l1_norms;
3008 }
3009
3010 /* Chebyshev */
3011 if (grid_relax_type[0] == 16 || grid_relax_type[1] == 16 ||
3012 grid_relax_type[2] == 16 || grid_relax_type[3] == 16)
3013 {
3014 max_eig_est = hypre_CTAlloc(HYPRE_Real, num_levels, HYPRE_MEMORY_HOST);
3015 min_eig_est = hypre_CTAlloc(HYPRE_Real, num_levels, HYPRE_MEMORY_HOST);
3016 hypre_ParAMGDataMaxEigEst(amg_data) = max_eig_est;
3017 hypre_ParAMGDataMinEigEst(amg_data) = min_eig_est;
3018 cheby_ds = hypre_CTAlloc(hypre_Vector *, num_levels, HYPRE_MEMORY_HOST);
3019 cheby_coefs = hypre_CTAlloc(HYPRE_Real *, num_levels, HYPRE_MEMORY_HOST);
3020 hypre_ParAMGDataChebyDS(amg_data) = cheby_ds;
3021 hypre_ParAMGDataChebyCoefs(amg_data) = cheby_coefs;
3022 }
3023
3024 /* CG */
3025 if (grid_relax_type[0] == 15 || grid_relax_type[1] == 15 ||
3026 grid_relax_type[2] == 15 || grid_relax_type[3] == 15)
3027 {
3028 smoother = hypre_CTAlloc(HYPRE_Solver, num_levels, HYPRE_MEMORY_HOST);
3029 hypre_ParAMGDataSmoother(amg_data) = smoother;
3030 }
3031
3032 if (addlvl == -1)
3033 {
3034 addlvl = num_levels;
3035 }
3036
3037 for (j = 0; j < addlvl; j++)
3038 {
3039 HYPRE_Real *l1_norm_data = NULL;
3040
3041 HYPRE_ANNOTATE_MGLEVEL_BEGIN(j);
3042 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Relaxation");
3043
3044 if (j < num_levels-1 &&
3045 (grid_relax_type[1] == 8 || grid_relax_type[1] == 13 || grid_relax_type[1] == 14 ||
3046 grid_relax_type[2] == 8 || grid_relax_type[2] == 13 || grid_relax_type[2] == 14))
3047 {
3048 if (relax_order)
3049 {
3050 hypre_ParCSRComputeL1Norms(A_array[j], 4, hypre_IntArrayData(CF_marker_array[j]), &l1_norm_data);
3051 }
3052 else
3053 {
3054 hypre_ParCSRComputeL1Norms(A_array[j], 4, NULL, &l1_norm_data);
3055 }
3056 }
3057 else if (j == num_levels-1 &&
3058 (grid_relax_type[3] == 8 || grid_relax_type[3] == 13 || grid_relax_type[3] == 14))
3059 {
3060 hypre_ParCSRComputeL1Norms(A_array[j], 4, NULL, &l1_norm_data);
3061 }
3062
3063 if (j < num_levels-1 && (grid_relax_type[1] == 18 || grid_relax_type[2] == 18))
3064 {
3065 if (relax_order)
3066 {
3067 hypre_ParCSRComputeL1Norms(A_array[j], 1, hypre_IntArrayData(CF_marker_array[j]), &l1_norm_data);
3068 }
3069 else
3070 {
3071 hypre_ParCSRComputeL1Norms(A_array[j], 1, NULL, &l1_norm_data);
3072 }
3073 }
3074 else if (j == num_levels-1 && grid_relax_type[3] == 18)
3075 {
3076 hypre_ParCSRComputeL1Norms(A_array[j], 1, NULL, &l1_norm_data);
3077 }
3078
3079 if (l1_norm_data)
3080 {
3081 l1_norms[j] = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[j]));
3082 hypre_VectorData(l1_norms[j]) = l1_norm_data;
3083 hypre_SeqVectorInitialize_v2(l1_norms[j], hypre_ParCSRMatrixMemoryLocation(A_array[j]));
3084 }
3085
3086 HYPRE_ANNOTATE_REGION_END("%s", "Relaxation");
3087 HYPRE_ANNOTATE_MGLEVEL_END(j);
3088 }
3089
3090 for (j = addlvl; j < hypre_min(add_end+1, num_levels) ; j++)
3091 {
3092 if (add_rlx == 18 )
3093 {
3094 HYPRE_Real *l1_norm_data = NULL;
3095
3096 HYPRE_ANNOTATE_MGLEVEL_BEGIN(j);
3097 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Relaxation");
3098
3099 hypre_ParCSRComputeL1Norms(A_array[j], 1, NULL, &l1_norm_data);
3100
3101 l1_norms[j] = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[j]));
3102 hypre_VectorData(l1_norms[j]) = l1_norm_data;
3103 hypre_SeqVectorInitialize_v2(l1_norms[j], hypre_ParCSRMatrixMemoryLocation(A_array[j]));
3104
3105 HYPRE_ANNOTATE_REGION_END("%s", "Relaxation");
3106 HYPRE_ANNOTATE_MGLEVEL_END(j);
3107 }
3108 }
3109
3110 for (j = add_end+1; j < num_levels; j++)
3111 {
3112 HYPRE_Real *l1_norm_data = NULL;
3113
3114 HYPRE_ANNOTATE_MGLEVEL_BEGIN(j);
3115 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Relaxation");
3116
3117 if (j < num_levels-1 && (grid_relax_type[1] == 8 || grid_relax_type[1] == 13 || grid_relax_type[1] == 14 ||
3118 grid_relax_type[2] == 8 || grid_relax_type[2] == 13 || grid_relax_type[2] == 14))
3119 {
3120 if (relax_order)
3121 {
3122 hypre_ParCSRComputeL1Norms(A_array[j], 4, hypre_IntArrayData(CF_marker_array[j]), &l1_norm_data);
3123 }
3124 else
3125 {
3126 hypre_ParCSRComputeL1Norms(A_array[j], 4, NULL, &l1_norm_data);
3127 }
3128 }
3129 else if ((grid_relax_type[3] == 8 || grid_relax_type[3] == 13 || grid_relax_type[3] == 14) && j == num_levels-1)
3130 {
3131 hypre_ParCSRComputeL1Norms(A_array[j], 4, NULL, &l1_norm_data);
3132 }
3133 if ((grid_relax_type[1] == 18 || grid_relax_type[2] == 18) && j < num_levels-1)
3134 {
3135 if (relax_order)
3136 {
3137 hypre_ParCSRComputeL1Norms(A_array[j], 1, hypre_IntArrayData(CF_marker_array[j]), &l1_norm_data);
3138 }
3139 else
3140 {
3141 hypre_ParCSRComputeL1Norms(A_array[j], 1, NULL, &l1_norm_data);
3142 }
3143 }
3144 else if (grid_relax_type[3] == 18 && j == num_levels-1)
3145 {
3146 hypre_ParCSRComputeL1Norms(A_array[j], 1, NULL, &l1_norm_data);
3147 }
3148
3149 if (l1_norm_data)
3150 {
3151 l1_norms[j] = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[j]));
3152 hypre_VectorData(l1_norms[j]) = l1_norm_data;
3153 hypre_SeqVectorInitialize_v2(l1_norms[j], hypre_ParCSRMatrixMemoryLocation(A_array[j]));
3154 }
3155
3156 HYPRE_ANNOTATE_REGION_END("%s", "Relaxation");
3157 HYPRE_ANNOTATE_MGLEVEL_END(j);
3158 }
3159
3160 for (j = 0; j < num_levels; j++)
3161 {
3162 HYPRE_ANNOTATE_MGLEVEL_BEGIN(j);
3163 HYPRE_ANNOTATE_REGION_BEGIN("%s", "Relaxation");
3164
3165 if (grid_relax_type[1] == 7 || grid_relax_type[2] == 7 || (grid_relax_type[3] == 7 && j == (num_levels-1)))
3166 {
3167 HYPRE_Real *l1_norm_data = NULL;
3168
3169 hypre_ParCSRComputeL1Norms(A_array[j], 5, NULL, &l1_norm_data);
3170
3171 l1_norms[j] = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[j]));
3172 hypre_VectorData(l1_norms[j]) = l1_norm_data;
3173 hypre_SeqVectorInitialize_v2(l1_norms[j], hypre_ParCSRMatrixMemoryLocation(A_array[j]));
3174 }
3175 else if (grid_relax_type[1] == 16 || grid_relax_type[2] == 16 || (grid_relax_type[3] == 16 && j== (num_levels-1)))
3176 {
3177 HYPRE_Int scale = hypre_ParAMGDataChebyScale(amg_data);
3178 /* If the full array is being considered, create the relevant temp vectors */
3179
3180 HYPRE_Int variant = hypre_ParAMGDataChebyVariant(amg_data);
3181 HYPRE_Real max_eig, min_eig = 0;
3182 HYPRE_Real *coefs = NULL;
3183 HYPRE_Int cheby_order = hypre_ParAMGDataChebyOrder(amg_data);
3184 HYPRE_Int cheby_eig_est = hypre_ParAMGDataChebyEigEst(amg_data);
3185 HYPRE_Real cheby_fraction = hypre_ParAMGDataChebyFraction(amg_data);
3186 if (cheby_eig_est)
3187 {
3188 hypre_ParCSRMaxEigEstimateCG(A_array[j], scale, cheby_eig_est,
3189 &max_eig, &min_eig);
3190 }
3191 else
3192 {
3193 hypre_ParCSRMaxEigEstimate(A_array[j], scale, &max_eig, &min_eig);
3194 }
3195 max_eig_est[j] = max_eig;
3196 min_eig_est[j] = min_eig;
3197
3198 cheby_ds[j] = hypre_SeqVectorCreate(hypre_ParCSRMatrixNumRows(A_array[j]));
3199 hypre_VectorVectorStride(cheby_ds[j]) = hypre_ParCSRMatrixNumRows(A_array[j]);
3200 hypre_VectorIndexStride(cheby_ds[j]) = 1;
3201 hypre_VectorMemoryLocation(cheby_ds[j]) = hypre_ParCSRMatrixMemoryLocation(A_array[j]);
3202
3203 hypre_ParCSRRelax_Cheby_Setup(A_array[j],
3204 max_eig,
3205 min_eig,
3206 cheby_fraction,
3207 cheby_order,
3208 scale,
3209 variant,
3210 &coefs,
3211 &hypre_VectorData(cheby_ds[j]));
3212 cheby_coefs[j] = coefs;
3213 }
3214 else if (grid_relax_type[1] == 15 || (grid_relax_type[3] == 15 && j == (num_levels-1)) )
3215 {
3216
3217 HYPRE_ParCSRPCGCreate(comm, &smoother[j]);
3218 /*HYPRE_ParCSRPCGSetup(smoother[j],
3219 (HYPRE_ParCSRMatrix) A_array[j],
3220 (HYPRE_ParVector) F_array[j],
3221 (HYPRE_ParVector) U_array[j]);*/
3222
3223 HYPRE_PCGSetTol(smoother[j], 1e-12); /* make small */
3224 HYPRE_PCGSetTwoNorm(smoother[j], 1); /* use 2-norm*/
3225
3226 HYPRE_ParCSRPCGSetup(smoother[j],
3227 (HYPRE_ParCSRMatrix) A_array[j],
3228 (HYPRE_ParVector) F_array[j],
3229 (HYPRE_ParVector) U_array[j]);
3230
3231
3232 }
3233 if (relax_weight[j] == 0.0)
3234 {
3235 hypre_ParCSRMatrixScaledNorm(A_array[j], &relax_weight[j]);
3236 if (relax_weight[j] != 0.0)
3237 relax_weight[j] = 4.0/3.0/relax_weight[j];
3238 else
3239 hypre_error_w_msg(HYPRE_ERROR_GENERIC," Warning ! Matrix norm is zero !!!");
3240 }
3241 if ((smooth_type == 6 || smooth_type == 16) && smooth_num_levels > j)
3242 {
3243 schwarz_relax_wt = hypre_ParAMGDataSchwarzRlxWeight(amg_data);
3244
3245 HYPRE_SchwarzCreate(&smoother[j]);
3246 HYPRE_SchwarzSetNumFunctions(smoother[j],num_functions);
3247 HYPRE_SchwarzSetVariant(smoother[j],
3248 hypre_ParAMGDataVariant(amg_data));
3249 HYPRE_SchwarzSetOverlap(smoother[j],
3250 hypre_ParAMGDataOverlap(amg_data));
3251 HYPRE_SchwarzSetDomainType(smoother[j],
3252 hypre_ParAMGDataDomainType(amg_data));
3253 HYPRE_SchwarzSetNonSymm(smoother[j],
3254 hypre_ParAMGDataSchwarzUseNonSymm(amg_data));
3255 if (schwarz_relax_wt > 0)
3256 HYPRE_SchwarzSetRelaxWeight(smoother[j],schwarz_relax_wt);
3257 HYPRE_SchwarzSetup(smoother[j],
3258 (HYPRE_ParCSRMatrix) A_array[j],
3259 (HYPRE_ParVector) f,
3260 (HYPRE_ParVector) u);
3261 if (schwarz_relax_wt < 0 )
3262 {
3263 num_cg_sweeps = (HYPRE_Int) (-schwarz_relax_wt);
3264 hypre_BoomerAMGCGRelaxWt(amg_data, j, num_cg_sweeps,
3265 &schwarz_relax_wt);
3266 /*hypre_printf (" schwarz weight %f \n", schwarz_relax_wt);*/
3267 HYPRE_SchwarzSetRelaxWeight(smoother[j], schwarz_relax_wt);
3268 if (hypre_ParAMGDataVariant(amg_data) > 0)
3269 {
3270 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[j]));
3271 hypre_SchwarzReScale(smoother[j], local_size, schwarz_relax_wt);
3272 }
3273 schwarz_relax_wt = 1;
3274 }
3275 }
3276 else if ((smooth_type == 9 || smooth_type == 19) && smooth_num_levels > j)
3277 {
3278 #ifdef HYPRE_MIXEDINT
3279 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"Euclid smoothing is not available in mixedint mode!");
3280 return hypre_error_flag;
3281 #endif
3282 HYPRE_EuclidCreate(comm, &smoother[j]);
3283 if (euclidfile)
3284 HYPRE_EuclidSetParamsFromFile(smoother[j],euclidfile);
3285 HYPRE_EuclidSetLevel(smoother[j],eu_level);
3286 if (eu_bj)
3287 HYPRE_EuclidSetBJ(smoother[j],eu_bj);
3288 if (eu_sparse_A)
3289 HYPRE_EuclidSetSparseA(smoother[j],eu_sparse_A);
3290 HYPRE_EuclidSetup(smoother[j],
3291 (HYPRE_ParCSRMatrix) A_array[j],
3292 (HYPRE_ParVector) F_array[j],
3293 (HYPRE_ParVector) U_array[j]);
3294 }
3295 else if ((smooth_type == 5 || smooth_type == 15) && smooth_num_levels > j)
3296 {
3297 HYPRE_ILUCreate( &smoother[j]);
3298 HYPRE_ILUSetType( smoother[j], ilu_type);
3299 HYPRE_ILUSetLocalReordering( smoother[j], ilu_reordering_type);
3300 HYPRE_ILUSetMaxIter(smoother[j], ilu_max_iter);
3301 HYPRE_ILUSetTol(smoother[j], 0.);
3302 HYPRE_ILUSetDropThreshold(smoother[j], ilu_droptol);
3303 HYPRE_ILUSetLogging(smoother[j], 0);
3304 HYPRE_ILUSetPrintLevel(smoother[j], 0);
3305 HYPRE_ILUSetLevelOfFill(smoother[j], ilu_lfil);
3306 HYPRE_ILUSetMaxNnzPerRow(smoother[j],ilu_max_row_nnz);
3307 HYPRE_ILUSetup(smoother[j],
3308 (HYPRE_ParCSRMatrix) A_array[j],
3309 (HYPRE_ParVector) F_array[j],
3310 (HYPRE_ParVector) U_array[j]);
3311 }
3312 else if ((smooth_type == 8 || smooth_type == 18) && smooth_num_levels > j)
3313 {
3314 #ifdef HYPRE_MIXEDINT
3315 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"ParaSails smoothing is not available in mixedint mode!");
3316 return hypre_error_flag;
3317 #endif
3318 HYPRE_ParCSRParaSailsCreate(comm, &smoother[j]);
3319 HYPRE_ParCSRParaSailsSetParams(smoother[j],thresh,nlevel);
3320 HYPRE_ParCSRParaSailsSetFilter(smoother[j],filter);
3321 HYPRE_ParCSRParaSailsSetSym(smoother[j],sym);
3322 HYPRE_ParCSRParaSailsSetup(smoother[j],
3323 (HYPRE_ParCSRMatrix) A_array[j],
3324 (HYPRE_ParVector) F_array[j],
3325 (HYPRE_ParVector) U_array[j]);
3326 }
3327 else if ((smooth_type == 7 || smooth_type == 17) && smooth_num_levels > j)
3328 {
3329 #ifdef HYPRE_MIXEDINT
3330 hypre_error_w_msg(HYPRE_ERROR_GENERIC,"pilut smoothing is not available in mixedint mode!");
3331 return hypre_error_flag;
3332 #endif
3333 HYPRE_ParCSRPilutCreate(comm, &smoother[j]);
3334 HYPRE_ParCSRPilutSetup(smoother[j],
3335 (HYPRE_ParCSRMatrix) A_array[j],
3336 (HYPRE_ParVector) F_array[j],
3337 (HYPRE_ParVector) U_array[j]);
3338 HYPRE_ParCSRPilutSetDropTolerance(smoother[j],drop_tol);
3339 HYPRE_ParCSRPilutSetFactorRowSize(smoother[j],max_nz_per_row);
3340 }
3341 else if ( (j < num_levels-1) || ((j == num_levels-1) &&
3342 (grid_relax_type[3] != 9 && grid_relax_type[3] != 99 &&
3343 grid_relax_type[3] != 19 && grid_relax_type[3] != 98) && coarse_size > 9) )
3344 {
3345 if (relax_weight[j] < 0)
3346 {
3347 num_cg_sweeps = (HYPRE_Int) (-relax_weight[j]);
3348 hypre_BoomerAMGCGRelaxWt(amg_data, j, num_cg_sweeps, &relax_weight[j]);
3349 }
3350 if (omega[j] < 0)
3351 {
3352 num_cg_sweeps = (HYPRE_Int) (-omega[j]);
3353 hypre_BoomerAMGCGRelaxWt(amg_data, j, num_cg_sweeps, &omega[j]);
3354 }
3355 }
3356
3357 HYPRE_ANNOTATE_REGION_END("%s", "Relaxation");
3358 HYPRE_ANNOTATE_MGLEVEL_END(j);
3359 } /* end of levels loop */
3360
3361 if (amg_logging > 1)
3362 {
3363 Residual_array = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[0]),
3364 hypre_ParCSRMatrixGlobalNumRows(A_array[0]),
3365 hypre_ParCSRMatrixRowStarts(A_array[0]) );
3366 hypre_ParVectorInitialize_v2(Residual_array, memory_location);
3367 hypre_ParAMGDataResidual(amg_data) = Residual_array;
3368 }
3369 else
3370 {
3371 hypre_ParAMGDataResidual(amg_data) = NULL;
3372 }
3373
3374 if (simple > -1 && simple < num_levels)
3375 {
3376 hypre_CreateDinv(amg_data);
3377 }
3378 else if ( (mult_additive > -1 && mult_additive < num_levels) ||
3379 (additive > -1 && additive < num_levels) )
3380 {
3381 hypre_CreateLambda(amg_data);
3382 }
3383
3384 /*-----------------------------------------------------------------------
3385 * Print some stuff
3386 *-----------------------------------------------------------------------*/
3387
3388 if (amg_print_level == 1 || amg_print_level == 3)
3389 {
3390 hypre_BoomerAMGSetupStats(amg_data,A);
3391 }
3392
3393 /* print out CF info to plot grids in matlab (see 'tools/AMGgrids.m') */
3394
3395 if (hypre_ParAMGDataPlotGrids(amg_data))
3396 {
3397 HYPRE_Int *CF, *CFc, *itemp;
3398 FILE* fp;
3399 char filename[256];
3400 HYPRE_Int coorddim = hypre_ParAMGDataCoordDim (amg_data);
3401 float *coordinates = hypre_ParAMGDataCoordinates (amg_data);
3402
3403 if (!coordinates) coorddim=0;
3404
3405 if (block_mode)
3406 {
3407 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRBlockMatrixDiag(A_block_array[0]));
3408 }
3409 else
3410 {
3411 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
3412 }
3413
3414 CF = hypre_CTAlloc(HYPRE_Int, local_size, HYPRE_MEMORY_HOST);
3415 CFc = hypre_CTAlloc(HYPRE_Int, local_size, HYPRE_MEMORY_HOST);
3416
3417 for (level = (num_levels - 2); level >= 0; level--)
3418 {
3419 /* swap pointers */
3420 itemp = CFc;
3421 CFc = CF;
3422 CF = itemp;
3423 if (block_mode)
3424 {
3425 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRBlockMatrixDiag(A_block_array[level]));
3426 }
3427 else
3428 {
3429 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A_array[level]));
3430 }
3431
3432 /* copy CF_marker to the host if needed */
3433 hypre_IntArray *CF_marker_host;
3434 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
3435 {
3436 CF_marker_host = hypre_IntArrayCloneDeep_v2(CF_marker_array[level], HYPRE_MEMORY_HOST);
3437 }
3438 else
3439 {
3440 CF_marker_host = CF_marker_array[level];
3441 }
3442 CF_marker = hypre_IntArrayData(CF_marker_host);
3443
3444 for (i = 0, j = 0; i < local_size; i++)
3445 {
3446 /* if a C-point */
3447 CF[i] = 0;
3448 if (CF_marker[i] > -1)
3449 {
3450 CF[i] = CFc[j] + 1;
3451 j++;
3452 }
3453 }
3454
3455 /* copy back to device and destroy host copy */
3456 if (hypre_GetActualMemLocation(hypre_IntArrayMemoryLocation(CF_marker_array[level])) == hypre_MEMORY_DEVICE)
3457 {
3458 hypre_IntArrayCopy(CF_marker_host, CF_marker_array[level]);
3459 hypre_IntArrayDestroy(CF_marker_host);
3460 }
3461 }
3462 if (block_mode)
3463 {
3464 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRBlockMatrixDiag(A_block_array[0]));
3465 }
3466 else
3467 {
3468 local_size = hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A));
3469 }
3470 hypre_sprintf (filename,"%s.%05d",hypre_ParAMGDataPlotFileName (amg_data),my_id);
3471 fp = fopen(filename, "w");
3472
3473 for (i = 0; i < local_size; i++)
3474 {
3475 for (j = 0; j < coorddim; j++)
3476 {
3477 hypre_fprintf (fp,"% f ",coordinates[coorddim*i+j]);
3478 }
3479 hypre_fprintf(fp, "%d\n", CF[i]);
3480 }
3481 fclose(fp);
3482
3483 hypre_TFree(CF, HYPRE_MEMORY_HOST);
3484 hypre_TFree(CFc, HYPRE_MEMORY_HOST);
3485 }
3486
3487 /* print out matrices on all levels */
3488 #if DEBUG
3489 {
3490 char filename[256];
3491
3492 if (block_mode)
3493 {
3494 hypre_ParCSRMatrix *temp_A;
3495
3496 for (level = 0; level < num_levels; level++)
3497 {
3498 hypre_sprintf(filename, "BoomerAMG.out.A_blk.%02d.ij", level);
3499 temp_A = hypre_ParCSRBlockMatrixConvertToParCSRMatrix(
3500 A_block_array[level]);
3501 hypre_ParCSRMatrixPrintIJ(temp_A, 0, 0, filename);
3502 hypre_ParCSRMatrixDestroy(temp_A);
3503 }
3504
3505 }
3506 else
3507 {
3508 for (level = 0; level < num_levels; level++)
3509 {
3510 hypre_sprintf(filename, "BoomerAMG.out.A.%02d.ij", level);
3511 hypre_ParCSRMatrixPrintIJ(A_array[level], 0, 0, filename);
3512 }
3513 for (level = 0; level < (num_levels-1); level++)
3514 {
3515 hypre_sprintf(filename, "BoomerAMG.out.P.%02d.ij", level);
3516 hypre_ParCSRMatrixPrintIJ(P_array[level], 0, 0, filename);
3517 }
3518 }
3519 }
3520 #endif
3521
3522 /* run compatible relaxation on all levels and print results */
3523 #if 0
3524 {
3525 hypre_ParVector *u_vec, *f_vec;
3526 HYPRE_Real *u, rho0, rho1, rho;
3527 HYPRE_Int n;
3528
3529 for (level = 0; level < (num_levels-1); level++)
3530 {
3531 u_vec = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
3532 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
3533 hypre_ParCSRMatrixRowStarts(A_array[level]));
3534 hypre_ParVectorInitialize(u_vec);
3535 f_vec = hypre_ParVectorCreate(hypre_ParCSRMatrixComm(A_array[level]),
3536 hypre_ParCSRMatrixGlobalNumRows(A_array[level]),
3537 hypre_ParCSRMatrixRowStarts(A_array[level]));
3538 hypre_ParVectorInitialize(f_vec);
3539
3540 hypre_ParVectorSetRandomValues(u_vec, 99);
3541 hypre_ParVectorSetConstantValues(f_vec, 0.0);
3542
3543 /* set C-pt values to zero */
3544 n = hypre_VectorSize(hypre_ParVectorLocalVector(u_vec));
3545 u = hypre_VectorData(hypre_ParVectorLocalVector(u_vec));
3546 for (i = 0; i < n; i++)
3547 {
3548 if (CF_marker_array[level][i] == 1)
3549 {
3550 u[i] = 0.0;
3551 }
3552 }
3553
3554 rho1 = hypre_ParVectorInnerProd(u_vec, u_vec);
3555 for (i = 0; i < 5; i++)
3556 {
3557 rho0 = rho1;
3558 hypre_BoomerAMGRelax(A_array[level], f_vec, CF_marker_array[level],
3559 grid_relax_type[0], -1,
3560 relax_weight[level], omega[level], l1_norms[level],
3561 u_vec, Vtemp, Ztemp);
3562 rho1 = hypre_ParVectorInnerProd(u_vec,u_vec);
3563 rho = sqrt(rho1/rho0);
3564 if (rho < 0.01)
3565 {
3566 break;
3567 }
3568 }
3569
3570 hypre_ParVectorDestroy(u_vec);
3571 hypre_ParVectorDestroy(f_vec);
3572
3573 if (my_id == 0)
3574 {
3575 hypre_printf("level = %d, rhocr = %f\n", level, rho);
3576 }
3577 }
3578 }
3579 #endif
3580
3581 HYPRE_ANNOTATE_FUNC_END;
3582
3583 return(hypre_error_flag);
3584 }
3585