1 /******************************************************************************
2 * Copyright 1998-2019 Lawrence Livermore National Security, LLC and other
3 * HYPRE Project Developers. See the top-level COPYRIGHT file for details.
4 *
5 * SPDX-License-Identifier: (Apache-2.0 OR MIT)
6 ******************************************************************************/
7
8 /******************************************************************************
9 * OpenMP Problems
10 *
11 * Need to fix the way these variables are set and incremented in loops:
12 * vals
13 *
14 ******************************************************************************/
15
16 #include "_hypre_sstruct_ls.h"
17 #include "_hypre_struct_mv.hpp"
18 #include "fac.h"
19
20 #define MapStencilRank(stencil, rank) \
21 { \
22 HYPRE_Int ii,jj,kk; \
23 ii = hypre_IndexX(stencil); \
24 jj = hypre_IndexY(stencil); \
25 kk = hypre_IndexZ(stencil); \
26 if (ii==-1) \
27 ii=2; \
28 if (jj==-1) \
29 jj=2; \
30 if (kk==-1) \
31 kk=2; \
32 rank = ii + 3*jj + 9*kk; \
33 }
34
35 #define InverseMapStencilRank(rank, stencil) \
36 { \
37 HYPRE_Int ij,ii,jj,kk; \
38 ij = (rank%9); \
39 ii = (ij%3); \
40 jj = (ij-ii)/3; \
41 kk = (rank-3*jj-ii)/9; \
42 if (ii==2) \
43 ii= -1; \
44 if (jj==2) \
45 jj= -1; \
46 if (kk==2) \
47 kk= -1; \
48 hypre_SetIndex3(stencil, ii, jj, kk); \
49 }
50
51
52 #define AbsStencilShape(stencil, abs_shape) \
53 { \
54 HYPRE_Int ii,jj,kk; \
55 ii = hypre_IndexX(stencil); \
56 jj = hypre_IndexY(stencil); \
57 kk = hypre_IndexZ(stencil); \
58 abs_shape= hypre_abs(ii) + hypre_abs(jj) + hypre_abs(kk); \
59 }
60
61 /*--------------------------------------------------------------------------
62 * hypre_AMR_FCoarsen: Coarsen the fbox and f/c connections. Forms the
63 * coarse operator by averaging neighboring connections in the refinement
64 * patch.
65 *--------------------------------------------------------------------------*/
66
67 HYPRE_Int
hypre_AMR_FCoarsen(hypre_SStructMatrix * A,hypre_SStructMatrix * fac_A,hypre_SStructPMatrix * A_crse,hypre_Index refine_factors,HYPRE_Int level)68 hypre_AMR_FCoarsen( hypre_SStructMatrix * A,
69 hypre_SStructMatrix * fac_A,
70 hypre_SStructPMatrix * A_crse,
71 hypre_Index refine_factors,
72 HYPRE_Int level )
73
74 {
75 hypre_Box fine_box;
76 hypre_Box intersect_box;
77
78 MPI_Comm comm = hypre_SStructMatrixComm(A);
79
80 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(A);
81 HYPRE_Int graph_type = hypre_SStructGraphObjectType(graph);
82 hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph);
83 HYPRE_IJMatrix ij_A = hypre_SStructMatrixIJMatrix(A);
84 HYPRE_Int matrix_type= hypre_SStructMatrixObjectType(A);
85 HYPRE_Int ndim = hypre_SStructMatrixNDim(A);
86
87 hypre_SStructPMatrix *A_pmatrix = hypre_SStructMatrixPMatrix(fac_A, level);
88
89 hypre_StructMatrix *smatrix_var;
90 hypre_StructStencil *stencils, *stencils_last;
91 HYPRE_Int stencil_size, stencil_last_size;
92 hypre_Index stencil_shape_i, stencil_last_shape_i;
93 hypre_Index loop_size;
94 hypre_Box loop_box;
95 HYPRE_Real **a_ptrs;
96 hypre_Box *A_dbox;
97
98 HYPRE_Int part_crse= level-1;
99 HYPRE_Int part_fine= level;
100
101 hypre_StructMatrix *crse_smatrix;
102 HYPRE_Real *crse_ptr;
103 HYPRE_Real **crse_ptrs;
104 hypre_Box *crse_dbox;
105
106 hypre_StructGrid *cgrid;
107 hypre_BoxArray *cgrid_boxes;
108 hypre_Box *cgrid_box;
109 hypre_Index cstart;
110 hypre_Index fstart, fend;
111 hypre_Index stridec, stridef;
112
113 hypre_StructGrid *fgrid;
114 hypre_BoxArray *fgrid_boxes;
115 hypre_Box *fgrid_box;
116 hypre_BoxArray ***fgrid_crse_extents;
117 hypre_BoxArray ***fbox_interior;
118 hypre_BoxArrayArray ***fbox_bdy;
119 HYPRE_Int ***interior_fboxi;
120 HYPRE_Int ***bdy_fboxi;
121 HYPRE_Int ***cboxi_fboxes;
122 HYPRE_Int **cboxi_fcnt;
123
124 hypre_BoxArray *fbox_interior_ci, *fbox_bdy_ci_fi;
125 hypre_BoxArrayArray *fbox_bdy_ci;
126 HYPRE_Int *interior_fboxi_ci;
127 HYPRE_Int *bdy_fboxi_ci;
128
129 HYPRE_Int centre;
130
131 hypre_BoxArray *data_space;
132
133 HYPRE_Int ci, fi, arrayi;
134 HYPRE_Int max_stencil_size= 27;
135 HYPRE_Int trueV = 1;
136 HYPRE_Int falseV= 0;
137 HYPRE_Int found, sort;
138 HYPRE_Int stencil_marker;
139 HYPRE_Int *stencil_ranks, *rank_stencils;
140 HYPRE_Int *stencil_contrib_cnt;
141 HYPRE_Int **stencil_contrib_i;
142 HYPRE_Real **weight_contrib_i;
143 HYPRE_Real weights[4]= {1.0, 0.25, 0.125, 0.0625};
144 HYPRE_Real sum;
145 HYPRE_Int abs_stencil_shape;
146 hypre_Box **shift_box;
147 hypre_Box coarse_cell_box;
148 HYPRE_Int volume_coarse_cell_box;
149 HYPRE_Int *volume_shift_box;
150 HYPRE_Int max_contribut_size, stencil_i;
151 HYPRE_BigInt startrank, rank;
152 HYPRE_Real *vals, *vals2;
153
154 HYPRE_Int i, j, k, l, m, n, ll, kk, jj;
155 HYPRE_Int nvars, var1, var2, var2_start;
156 HYPRE_Int iA_shift_z, iA_shift_zy, iA_shift_zyx;
157
158 hypre_Index lindex;
159 hypre_Index index1, index2;
160 hypre_Index index_temp;
161
162 HYPRE_Int **box_graph_indices;
163 HYPRE_Int *box_graph_cnts;
164 HYPRE_Int *box_ranks, *box_ranks_cnt, *box_to_ranks_cnt;
165 HYPRE_Int *cdata_space_ranks, *box_starts, *box_ends;
166 HYPRE_Int *box_connections;
167 HYPRE_Int **coarse_contrib_Uv;
168 HYPRE_Int *fine_interface_ranks;
169 HYPRE_Int nUventries= hypre_SStructGraphNUVEntries(graph);
170 HYPRE_Int *iUventries = hypre_SStructGraphIUVEntries(graph);
171 hypre_SStructUVEntry **Uventries = hypre_SStructGraphUVEntries(graph);
172 hypre_SStructUVEntry *Uventry;
173 HYPRE_Int nUentries, cnt1;
174 hypre_Index index, *cindex, *Uv_cindex;
175 HYPRE_Int box_array_size, cbox_array_size;
176
177 HYPRE_Int nrows;
178 HYPRE_BigInt to_rank;
179 HYPRE_Int *ncols;
180 HYPRE_BigInt *rows, *cols;
181 HYPRE_Int **interface_max_stencil_ranks;
182 HYPRE_Int **interface_max_stencil_cnt;
183 HYPRE_Int **interface_rank_stencils;
184 HYPRE_Int **interface_stencil_ranks;
185 HYPRE_Int *coarse_stencil_cnt;
186 HYPRE_Real *stencil_vals;
187 HYPRE_Int *common_rank_stencils, *common_stencil_ranks;
188 HYPRE_Int *common_stencil_i;
189 hypre_BoxManEntry *boxman_entry;
190
191 HYPRE_Int *temp1, *temp2;
192 HYPRE_Real *temp3;
193 HYPRE_Real sum_contrib, scaling;
194
195 HYPRE_Int **OffsetA;
196
197 HYPRE_Int *parents;
198 HYPRE_Int *parents_cnodes;
199
200 HYPRE_Int myid;
201
202 hypre_MPI_Comm_rank(comm, &myid);
203
204 hypre_BoxInit(&fine_box, ndim);
205 hypre_BoxInit(&intersect_box, ndim);
206 hypre_BoxInit(&loop_box, ndim);
207 hypre_BoxInit(&coarse_cell_box, ndim);
208
209 /*--------------------------------------------------------------------------
210 * Task: Coarsen the fbox and f/c connections to form the coarse grid
211 * operator inside the fgrid.
212 *--------------------------------------------------------------------------*/
213
214 if (graph_type == HYPRE_SSTRUCT)
215 {
216 startrank = hypre_SStructGridGhstartRank(grid);
217 }
218 if (graph_type == HYPRE_PARCSR)
219 {
220 startrank = hypre_SStructGridStartRank(grid);
221 }
222
223 /*--------------------------------------------------------------------------
224 * Fine grid strides by the refinement factors.
225 *--------------------------------------------------------------------------*/
226 hypre_SetIndex3(stridec, 1, 1, 1);
227 for (i= 0; i< ndim; i++)
228 {
229 stridef[i]= refine_factors[i];
230 }
231 for (i= ndim; i< 3; i++)
232 {
233 stridef[i]= 1;
234 }
235
236 /*--------------------------------------------------------------------------
237 * Scaling for averaging row sum.
238 *--------------------------------------------------------------------------*/
239 scaling= 1.0;
240 for (i= 0; i< ndim-2; i++)
241 {
242 scaling*= refine_factors[0];
243 }
244
245 /*--------------------------------------------------------------------------
246 * Determine the coarsened fine grid- fgrid_crse_extents.
247 * These are between fpart= level and cpart= (level-1). The
248 * fgrid_crse_extents will be indexed by cboxes- the boxarray of coarsened
249 * fboxes FULLY in a given cbox.
250 *
251 * Also, determine the interior and boundary boxes of each fbox. Having
252 * these will allow us to determine the f/c interface nodes without
253 * extensive checking. These are also indexed by the cboxes.
254 * fgrid_interior- for each cbox, we have a collection of child fboxes,
255 * each leading to an interior=> boxarray
256 * fgrid_bdy - for each cbox, we have a collection of child fboxes,
257 * each leading to a boxarray of bdies=> boxarrayarray.
258 * Because we need to know the fbox id for these boxarray/boxarrayarray,
259 * we will need one for each fbox.
260 *
261 * And, determine which cboxes contain a given fbox. That is, given a
262 * fbox, find all cboxes that contain a chunk of it.
263 *--------------------------------------------------------------------------*/
264 nvars = hypre_SStructPMatrixNVars(A_pmatrix);
265
266 fgrid_crse_extents = hypre_TAlloc(hypre_BoxArray **, nvars, HYPRE_MEMORY_HOST);
267 fbox_interior = hypre_TAlloc(hypre_BoxArray **, nvars, HYPRE_MEMORY_HOST);
268 fbox_bdy = hypre_TAlloc(hypre_BoxArrayArray **, nvars, HYPRE_MEMORY_HOST);
269 interior_fboxi = hypre_TAlloc(HYPRE_Int **, nvars, HYPRE_MEMORY_HOST);
270 bdy_fboxi = hypre_TAlloc(HYPRE_Int **, nvars, HYPRE_MEMORY_HOST);
271 cboxi_fboxes = hypre_TAlloc(HYPRE_Int **, nvars, HYPRE_MEMORY_HOST);
272 cboxi_fcnt = hypre_TAlloc(HYPRE_Int *, nvars, HYPRE_MEMORY_HOST);
273
274 for (var1= 0; var1< nvars; var1++)
275 {
276 cgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_crse), var1);
277 cgrid_boxes= hypre_StructGridBoxes(cgrid);
278 fgrid_crse_extents[var1]= hypre_TAlloc(hypre_BoxArray *,
279 hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
280 fbox_interior[var1]= hypre_TAlloc(hypre_BoxArray *,
281 hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
282 fbox_bdy[var1] = hypre_TAlloc(hypre_BoxArrayArray *,
283 hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
284 interior_fboxi[var1]= hypre_TAlloc(HYPRE_Int *, hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
285 bdy_fboxi[var1] = hypre_TAlloc(HYPRE_Int *, hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
286
287 fgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_pmatrix), var1);
288 fgrid_boxes= hypre_StructGridBoxes(fgrid);
289
290 cboxi_fboxes[var1]= hypre_CTAlloc(HYPRE_Int *, hypre_BoxArraySize(fgrid_boxes), HYPRE_MEMORY_HOST);
291 cboxi_fcnt[var1] = hypre_CTAlloc(HYPRE_Int , hypre_BoxArraySize(fgrid_boxes), HYPRE_MEMORY_HOST);
292
293 /*-----------------------------------------------------------------------
294 * Determine the fine grid boxes that are underlying a coarse grid box.
295 * Coarsen the indices to determine the looping extents of these
296 * boxes. Also, find the looping extents for the extended coarsened
297 * boxes, and the interior and boundary extents of a fine_grid box.
298 * The fine_grid boxes must be adjusted so that only the coarse nodes
299 * inside these boxes are included. Only the lower bound needs to be
300 * adjusted.
301 *-----------------------------------------------------------------------*/
302 hypre_ForBoxI(ci, cgrid_boxes)
303 {
304 cgrid_box= hypre_BoxArrayBox(cgrid_boxes, ci);
305 hypre_CopyIndex(hypre_BoxIMin(cgrid_box), cstart);
306
307 cnt1= 0;
308 temp1= hypre_CTAlloc(HYPRE_Int, hypre_BoxArraySize(fgrid_boxes), HYPRE_MEMORY_HOST);
309
310 hypre_ClearIndex(index_temp);
311 hypre_ForBoxI(fi, fgrid_boxes)
312 {
313 fgrid_box= hypre_BoxArrayBox(fgrid_boxes, fi);
314 hypre_CopyIndex(hypre_BoxIMin(fgrid_box), fstart);
315 for (i= 0; i< ndim; i++)
316 {
317 j= fstart[i]%refine_factors[i];
318 if (j)
319 {
320 fstart[i]+= refine_factors[i] - j;
321 }
322 }
323
324 hypre_StructMapFineToCoarse(fstart, index_temp,
325 refine_factors, hypre_BoxIMin(&fine_box));
326 hypre_StructMapFineToCoarse(hypre_BoxIMax(fgrid_box), index_temp,
327 refine_factors, hypre_BoxIMax(&fine_box));
328
329 hypre_IntersectBoxes(&fine_box, cgrid_box, &intersect_box);
330 if (hypre_BoxVolume(&intersect_box) > 0)
331 {
332 temp1[cnt1++]= fi;
333 }
334 }
335
336 fgrid_crse_extents[var1][ci]= hypre_BoxArrayCreate(cnt1, ndim);
337 fbox_interior[var1][ci] = hypre_BoxArrayCreate(cnt1, ndim);
338 fbox_bdy[var1][ci] = hypre_BoxArrayArrayCreate(cnt1, ndim);
339 interior_fboxi[var1][ci] = hypre_CTAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
340 bdy_fboxi[var1][ci] = hypre_CTAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
341
342 for (fi= 0; fi< cnt1; fi++)
343 {
344 fgrid_box= hypre_BoxArrayBox(fgrid_boxes, temp1[fi]);
345 hypre_CopyIndex(hypre_BoxIMin(fgrid_box), fstart);
346 hypre_CopyIndex(hypre_BoxIMax(fgrid_box), fend);
347
348 /*--------------------------------------------------------------------
349 * record which sides will be adjusted- fstart adjustments will
350 * decrease the box size, whereas fend adjustments will increase the
351 * box size. Since we fstart decreases the box size, we cannot
352 * have an f/c interface at an adjusted fstart end. fend may
353 * correspond to an f/c interface whether it has been adjusted or not.
354 *--------------------------------------------------------------------*/
355 hypre_SetIndex3(index1, 1, 1, 1);
356 for (i= 0; i< ndim; i++)
357 {
358 j= fstart[i]%refine_factors[i];
359 if (j)
360 {
361 fstart[i]+= refine_factors[i] - j;
362 index1[i] = 0;
363 }
364
365 j= fend[i]%refine_factors[i];
366 if (refine_factors[i]-1 - j)
367 {
368 fend[i] +=(refine_factors[i]-1) - j;
369 }
370 }
371
372 hypre_StructMapFineToCoarse(fstart, index_temp,
373 refine_factors, hypre_BoxIMin(&fine_box));
374 hypre_StructMapFineToCoarse(hypre_BoxIMax(fgrid_box), index_temp,
375 refine_factors, hypre_BoxIMax(&fine_box));
376 hypre_IntersectBoxes(&fine_box, cgrid_box, &intersect_box);
377
378 hypre_CopyBox(&intersect_box,
379 hypre_BoxArrayBox(fgrid_crse_extents[var1][ci], fi));
380
381 /*--------------------------------------------------------------------
382 * adjust the fine intersect_box so that we get the interior and
383 * boundaries separately.
384 *--------------------------------------------------------------------*/
385 hypre_StructMapCoarseToFine(hypre_BoxIMin(&intersect_box), index_temp,
386 refine_factors, hypre_BoxIMin(&fine_box));
387
388 /* the following index2 shift for ndim<3 is no problem since
389 refine_factors[j]= 1 for j>=ndim. */
390 hypre_SetIndex3(index2, refine_factors[0]-1, refine_factors[1]-1,
391 refine_factors[2]-1);
392 hypre_StructMapCoarseToFine(hypre_BoxIMax(&intersect_box), index2,
393 refine_factors, hypre_BoxIMax(&fine_box));
394
395 hypre_SetIndex3(index2, 1, 1, 1);
396 hypre_CopyBox(&fine_box, &loop_box);
397 for (i= 0; i< ndim; i++)
398 {
399 hypre_BoxIMin(&loop_box)[i]+= refine_factors[i]*index1[i];
400 hypre_BoxIMax(&loop_box)[i]-= refine_factors[i]*index2[i];
401 }
402 hypre_CopyBox(&loop_box,
403 hypre_BoxArrayBox(fbox_interior[var1][ci], fi));
404 interior_fboxi[var1][ci][fi]= temp1[fi];
405
406 hypre_SubtractBoxes(&fine_box, &loop_box,
407 hypre_BoxArrayArrayBoxArray(fbox_bdy[var1][ci], fi));
408 bdy_fboxi[var1][ci][fi]= temp1[fi];
409 }
410 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
411
412 } /* hypre_ForBoxI(ci, cgrid_boxes) */
413
414 /*--------------------------------------------------------------------
415 * Determine the cboxes that contain a chunk of a given fbox.
416 *--------------------------------------------------------------------*/
417 hypre_ForBoxI(fi, fgrid_boxes)
418 {
419 fgrid_box= hypre_BoxArrayBox(fgrid_boxes, fi);
420 hypre_CopyIndex(hypre_BoxIMin(fgrid_box), fstart);
421 for (i= 0; i< ndim; i++)
422 {
423 j= fstart[i]%refine_factors[i];
424 if (j)
425 {
426 fstart[i]+= refine_factors[i] - j;
427 }
428 }
429
430 hypre_StructMapFineToCoarse(fstart, index_temp,
431 refine_factors, hypre_BoxIMin(&fine_box));
432 hypre_StructMapFineToCoarse(hypre_BoxIMax(fgrid_box), index_temp,
433 refine_factors, hypre_BoxIMax(&fine_box));
434
435 temp1= hypre_CTAlloc(HYPRE_Int, hypre_BoxArraySize(cgrid_boxes), HYPRE_MEMORY_HOST);
436 hypre_ForBoxI(i, cgrid_boxes)
437 {
438 cgrid_box= hypre_BoxArrayBox(cgrid_boxes, i);
439 hypre_IntersectBoxes(&fine_box, cgrid_box, &intersect_box);
440 if (hypre_BoxVolume(&intersect_box) > 0)
441 {
442 temp1[cboxi_fcnt[var1][fi]]= i;
443 cboxi_fcnt[var1][fi]++;
444 }
445 }
446
447 cboxi_fboxes[var1][fi]= hypre_TAlloc(HYPRE_Int, cboxi_fcnt[var1][fi], HYPRE_MEMORY_HOST);
448 for (i= 0; i< cboxi_fcnt[var1][fi]; i++)
449 {
450 cboxi_fboxes[var1][fi][i]= temp1[i];
451 }
452 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
453 }
454 } /* for (var1= 0; var1< nvars; var1++) */
455
456 /*--------------------------------------------------------------------------
457 * STEP 1:
458 * COMPUTE THE COARSE LEVEL OPERATOR INSIDE OF A REFINED BOX.
459 *
460 * We assume that the coarse and fine grid variables are of the same type.
461 *
462 * Coarse stencils in the refinement patches are obtained by averaging the
463 * fine grid coefficients. Since we are assuming cell-centred discretization,
464 * we apply a weighted averaging of ONLY the fine grid coefficients along
465 * interfaces of adjacent agglomerated coarse cells.
466 *
467 * Since the stencil pattern is assumed arbitrary, we must determine the
468 * stencil pattern of each var1-var2 struct_matrix to get the correct
469 * contributing stencil coefficients, averaging weights, etc.
470 *--------------------------------------------------------------------------*/
471
472 /*--------------------------------------------------------------------------
473 * Agglomerated coarse cell info. These are needed in defining the looping
474 * extents for averaging- i.e., we loop over extents determined by the
475 * size of the agglomerated coarse cell.
476 * Note that the agglomerated coarse cell is constructed correctly for
477 * any dimensions (1, 2, or 3).
478 *--------------------------------------------------------------------------*/
479 hypre_ClearIndex(index_temp);
480 hypre_CopyIndex(index_temp, hypre_BoxIMin(&coarse_cell_box));
481 hypre_SetIndex3(index_temp, refine_factors[0]-1, refine_factors[1]-1,
482 refine_factors[2]-1 );
483 hypre_CopyIndex(index_temp, hypre_BoxIMax(&coarse_cell_box));
484
485 volume_coarse_cell_box= hypre_BoxVolume(&coarse_cell_box);
486
487
488 /*--------------------------------------------------------------------------
489 * Offsets in y & z directions for refinement patches. These will be used
490 * for pointing to correct coarse stencil location.
491 *--------------------------------------------------------------------------*/
492 OffsetA = hypre_CTAlloc(HYPRE_Int *, 2, HYPRE_MEMORY_HOST);
493 for (i= 0; i< 2; i++)
494 {
495 OffsetA[i]= hypre_CTAlloc(HYPRE_Int, refine_factors[i+1], HYPRE_MEMORY_HOST);
496 }
497
498 /*--------------------------------------------------------------------------
499 * Stencil contribution cnts, weights, etc are computed only if we have
500 * a new stencil pattern. If the pattern is the same, the previously
501 * computed stencil contribution cnts, weights, etc can be used.
502 *
503 * Mark the stencil_marker so that the first time the stencil is non-null,
504 * the stencil contribution cnts, weights, etc are computed.
505 *--------------------------------------------------------------------------*/
506 stencil_marker= trueV;
507 for (var1= 0; var1< nvars; var1++)
508 {
509 cgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_crse), var1);
510 cgrid_boxes= hypre_StructGridBoxes(cgrid);
511
512 fgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_pmatrix), var1);
513 fgrid_boxes= hypre_StructGridBoxes(fgrid);
514
515
516 for (var2= 0; var2< nvars; var2++)
517 {
518 stencils= hypre_SStructPMatrixSStencil(A_crse, var1, var2);
519 if (stencils != NULL)
520 {
521 stencil_size= hypre_StructStencilSize(stencils);
522
523 /*-----------------------------------------------------------------
524 * When stencil_marker== true, form the stencil contributions cnts,
525 * weights, etc. This occurs for the first non-null stencil or
526 * when the stencil shape of the current non-null stencil has a
527 * different stencil shape from that of the latest non-null stencil.
528 *
529 * But when stencil_marker== false, we must check to see if we
530 * need new stencil contributions cnts, weights, etc. Thus, find
531 * the latest non-null stencil for comparison.
532 *-----------------------------------------------------------------*/
533 if (stencil_marker == falseV)
534 {
535 /* search for the first previous non-null stencil */
536 found = falseV;
537 var2_start= var2-1;
538 for (j= var1; j>= 0; j--)
539 {
540 for (i= var2_start; i>= 0; i--)
541 {
542 stencils_last= hypre_SStructPMatrixSStencil(A_crse, j, i);
543 if (stencils_last != NULL)
544 {
545 found= trueV;
546 break;
547 }
548 }
549 if (found)
550 {
551 break;
552 }
553 else
554 {
555 var2_start= nvars-1;
556 }
557 }
558
559 /*--------------------------------------------------------------
560 * Compare the stencil shape.
561 *--------------------------------------------------------------*/
562 stencil_last_size= hypre_StructStencilSize(stencils_last);
563 if (stencil_last_size != stencil_size)
564 {
565 stencil_marker= trueV;
566 break;
567 }
568 else
569 {
570 found= falseV;
571 for (i= 0; i< stencil_size; i++)
572 {
573 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
574 stencil_shape_i);
575 hypre_CopyIndex(hypre_StructStencilElement(stencils_last,i),
576 stencil_last_shape_i);
577
578 hypre_SetIndex3(index_temp,
579 stencil_shape_i[0]-stencil_last_shape_i[0],
580 stencil_shape_i[1]-stencil_last_shape_i[1],
581 stencil_shape_i[2]-stencil_last_shape_i[2]);
582
583 AbsStencilShape(index_temp, abs_stencil_shape);
584 if (abs_stencil_shape)
585 {
586 found= trueV;
587 stencil_marker= trueV;
588 hypre_TFree(stencil_contrib_cnt, HYPRE_MEMORY_HOST);
589 hypre_TFree(stencil_ranks, HYPRE_MEMORY_HOST);
590 for (i= 0; i< stencil_size; i++)
591 {
592 hypre_BoxDestroy(shift_box[i]);
593 }
594 hypre_TFree(shift_box, HYPRE_MEMORY_HOST);
595 hypre_TFree(volume_shift_box, HYPRE_MEMORY_HOST);
596 hypre_TFree(vals, HYPRE_MEMORY_HOST);
597
598 for (j= 1; j< max_stencil_size; j++)
599 {
600 stencil_i= rank_stencils[j];
601 if (stencil_i != -1)
602 {
603 hypre_TFree(stencil_contrib_i[stencil_i], HYPRE_MEMORY_HOST);
604 hypre_TFree(weight_contrib_i[stencil_i], HYPRE_MEMORY_HOST);
605 }
606 }
607 hypre_TFree(stencil_contrib_i, HYPRE_MEMORY_HOST);
608 hypre_TFree(weight_contrib_i, HYPRE_MEMORY_HOST);
609 hypre_TFree(rank_stencils, HYPRE_MEMORY_HOST);
610 }
611
612 if (found)
613 {
614 break;
615 }
616 } /* for (i= 0; i< stencil_size; i++) */
617 } /* else */
618 } /* if (stencil_marker == false) */
619
620 /*-----------------------------------------------------------------
621 * If stencil_marker==true, form the contribution structures.
622 * Since the type of averaging is determined by the stencil shapes,
623 * we need a ranking of the stencil shape to allow for easy
624 * determination.
625 *
626 * top: 14 12 13 centre: 5 3 4 bottom 23 21 22
627 * 11 9 10 2 0 1 20 18 19
628 * 17 15 16 8 6 7 26 24 25
629 *
630 * for stencil of max. size 27.
631 *
632 * stencil_contrib_cnt[i]= no. of fine stencils averaged to
633 * form stencil entry i.
634 * stencil_contrib_i[i] = rank of fine stencils contributing
635 * to form stencil entry i.
636 * weight_contrib_i[i] = array of weights for weighting
637 * the contributions to stencil entry i.
638 * stencil_ranks[i] = rank of stencil entry i.
639 * rank_stencils[i] = stencil entry of rank i.
640 *-----------------------------------------------------------------*/
641
642 if (stencil_marker == trueV)
643 {
644
645 /* mark stencil_marker for the next stencil */
646 stencil_marker= falseV;
647
648 stencil_contrib_cnt= hypre_CTAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
649 stencil_contrib_i = hypre_TAlloc(HYPRE_Int *, stencil_size, HYPRE_MEMORY_HOST);
650 weight_contrib_i = hypre_TAlloc(HYPRE_Real *, stencil_size, HYPRE_MEMORY_HOST);
651 stencil_ranks = hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
652 rank_stencils = hypre_TAlloc(HYPRE_Int, max_stencil_size, HYPRE_MEMORY_HOST);
653 shift_box = hypre_TAlloc(hypre_Box *, stencil_size, HYPRE_MEMORY_HOST);
654 volume_shift_box = hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
655
656 for (i= 0; i< max_stencil_size; i++)
657 {
658 rank_stencils[i]= -1;
659 if (i < stencil_size)
660 {
661 stencil_ranks[i]= -1;
662 }
663 }
664
665 /*-----------------------------------------------------------------
666 * Get mappings between stencil entries and ranks and vice versa;
667 * fine grid looping extents for averaging of the fine coefficients;
668 * and the number of fine grid values to be averaged.
669 * Note that the shift_boxes are constructed correctly for any
670 * dimensions. For j>=ndim,
671 * hypre_BoxIMin(shift_box[i])[j]=hypre_BoxIMax(shift_box[i])[j]= 0.
672 *-----------------------------------------------------------------*/
673 for (i= 0; i< stencil_size; i++)
674 {
675 shift_box[i] = hypre_BoxCreate(ndim);
676 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
677 stencil_shape_i);
678 MapStencilRank(stencil_shape_i, j);
679 stencil_ranks[i]= j;
680 rank_stencils[stencil_ranks[i]] = i;
681
682 hypre_SetIndex3(hypre_BoxIMin(shift_box[i]),
683 (refine_factors[0]-1)*stencil_shape_i[0],
684 (refine_factors[1]-1)*stencil_shape_i[1],
685 (refine_factors[2]-1)*stencil_shape_i[2]);
686
687 hypre_AddIndexes(hypre_BoxIMin(shift_box[i]),
688 hypre_BoxIMax(&coarse_cell_box), 3,
689 hypre_BoxIMax(shift_box[i]));
690
691 hypre_IntersectBoxes(&coarse_cell_box, shift_box[i], shift_box[i]);
692
693 volume_shift_box[i]= hypre_BoxVolume(shift_box[i]);
694 }
695
696 /*-----------------------------------------------------------------
697 * Derive the contribution info.
698 * The above rank table is used to determine the direction indices.
699 * Weight construction procedure valid for any dimensions.
700 *-----------------------------------------------------------------*/
701
702 /* east */
703 stencil_i= rank_stencils[1];
704 if (stencil_i != -1)
705 {
706 stencil_contrib_cnt[stencil_i]++;
707 for (i= 4; i<= 7; i+=3)
708 {
709 if (rank_stencils[i] != -1) /* ne or se */
710 stencil_contrib_cnt[stencil_i]++;
711 }
712
713 if (ndim > 2)
714 {
715 for (j= 1; j<= 2; j++)
716 {
717 for (i= 1; i<= 7; i+=3)
718 {
719 if (rank_stencils[j*9+i] != -1) /* bottom or top planes */
720 stencil_contrib_cnt[stencil_i]++;
721 }
722 }
723 }
724 max_contribut_size= stencil_contrib_cnt[stencil_i];
725 }
726
727 /* fill up the east contribution stencil indices */
728 if (stencil_i != -1)
729 {
730 stencil_contrib_i[stencil_i]=
731 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
732 weight_contrib_i[stencil_i] =
733 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
734 sum= 0.0;
735 k= 0;
736
737 stencil_contrib_i[stencil_i][k]= stencil_i;
738 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
739 abs_stencil_shape );
740 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
741 sum+= weights[abs_stencil_shape];
742
743 for (i= 4; i<= 7; i+=3)
744 {
745 if (rank_stencils[i] != -1)
746 {
747 stencil_contrib_i[stencil_i][k] = rank_stencils[i];
748
749 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
750 abs_stencil_shape );
751 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
752 sum+= weights[abs_stencil_shape];
753 }
754 }
755
756 if (ndim > 2)
757 {
758 for (j= 1; j<= 2; j++)
759 {
760 for (i= 1; i<= 7; i+=3)
761 {
762 if (rank_stencils[j*9+i] != -1)
763 {
764 stencil_contrib_i[stencil_i][k] = rank_stencils[j*9+i];
765
766 AbsStencilShape(
767 hypre_StructStencilElement(stencils,rank_stencils[j*9+i]),
768 abs_stencil_shape );
769 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
770 sum+= weights[abs_stencil_shape];
771 }
772 }
773 }
774 }
775
776 for (i= 0; i< k ; i++)
777 {
778 weight_contrib_i[stencil_i][i]/= sum;
779 }
780 }
781
782
783 /* west */
784 stencil_i= rank_stencils[2];
785 if (stencil_i != -1)
786 {
787 stencil_contrib_cnt[stencil_i]++;
788 for (i= 5; i<= 8; i+=3)
789 {
790 if (rank_stencils[i] != -1) /* nw or sw */
791 stencil_contrib_cnt[stencil_i]++;
792 }
793
794 if (ndim > 2)
795 {
796 for (j= 1; j<= 2; j++)
797 {
798 for (i= 2; i<= 8; i+=3)
799 {
800 if (rank_stencils[j*9+i] != -1) /* bottom or top planes */
801 stencil_contrib_cnt[stencil_i]++;
802 }
803 }
804 }
805 max_contribut_size= hypre_max( max_contribut_size,
806 stencil_contrib_cnt[stencil_i] );
807 }
808
809 if (stencil_i != -1)
810 {
811 stencil_contrib_i[stencil_i]=
812 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
813 weight_contrib_i[stencil_i] =
814 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
815 sum= 0.0;
816 k= 0;
817
818 stencil_contrib_i[stencil_i][k]= stencil_i;
819 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
820 abs_stencil_shape );
821 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
822 sum+= weights[abs_stencil_shape];
823
824 for (i= 5; i<= 8; i+=3)
825 {
826 if (rank_stencils[i] != -1)
827 {
828 stencil_contrib_i[stencil_i][k] = rank_stencils[i];
829
830 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
831 abs_stencil_shape );
832 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
833 sum+= weights[abs_stencil_shape];
834 }
835 }
836
837 if (ndim > 2)
838 {
839 for (j= 1; j<= 2; j++)
840 {
841 for (i= 2; i<= 8; i+=3)
842 {
843 if (rank_stencils[j*9+i] != -1)
844 {
845 stencil_contrib_i[stencil_i][k] = rank_stencils[j*9+i];
846
847 AbsStencilShape(
848 hypre_StructStencilElement(stencils,rank_stencils[j*9+i]),
849 abs_stencil_shape );
850 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
851 sum+= weights[abs_stencil_shape];
852 }
853 }
854 }
855 }
856
857 for (i= 0; i< k ; i++)
858 {
859 weight_contrib_i[stencil_i][i]/= sum;
860 }
861 }
862
863
864 /* north */
865 stencil_i= rank_stencils[3];
866 if (stencil_i != -1)
867 {
868 stencil_contrib_cnt[stencil_i]++;
869 for (i= 4; i<= 5; i++)
870 {
871 if (rank_stencils[i] != -1) /* ne or nw */
872 stencil_contrib_cnt[stencil_i]++;
873 }
874
875 if (ndim > 2)
876 {
877 for (j= 1; j<= 2; j++)
878 {
879 for (i= 3; i<= 5; i++)
880 {
881 if (rank_stencils[j*9+i] != -1) /* bottom or top planes */
882 stencil_contrib_cnt[stencil_i]++;
883 }
884 }
885 }
886 max_contribut_size= hypre_max( max_contribut_size,
887 stencil_contrib_cnt[stencil_i] );
888 }
889
890 if (stencil_i != -1)
891 {
892 stencil_contrib_i[stencil_i]=
893 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
894 weight_contrib_i[stencil_i] =
895 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
896 sum= 0.0;
897 k= 0;
898
899 stencil_contrib_i[stencil_i][k]= stencil_i;
900 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
901 abs_stencil_shape );
902 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
903 sum+= weights[abs_stencil_shape];
904
905 for (i= 4; i<= 5; i++)
906 {
907 if (rank_stencils[i] != -1)
908 {
909 stencil_contrib_i[stencil_i][k] = rank_stencils[i];
910
911 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
912 abs_stencil_shape );
913 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
914 sum+= weights[abs_stencil_shape];
915 }
916 }
917
918 if (ndim > 2)
919 {
920 for (j= 1; j<= 2; j++)
921 {
922 for (i= 3; i<= 5; i++)
923 {
924 if (rank_stencils[j*9+i] != -1)
925 {
926 stencil_contrib_i[stencil_i][k] = rank_stencils[j*9+i];
927
928 AbsStencilShape(
929 hypre_StructStencilElement(stencils,rank_stencils[j*9+i]),
930 abs_stencil_shape );
931 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
932 sum+= weights[abs_stencil_shape];
933 }
934 }
935 }
936 }
937
938 for (i= 0; i< k ; i++)
939 {
940 weight_contrib_i[stencil_i][i]/= sum;
941 }
942 }
943
944 /* south */
945 stencil_i= rank_stencils[6];
946 if (stencil_i != -1)
947 {
948 stencil_contrib_cnt[stencil_i]++;
949 for (i= 7; i<= 8; i++)
950 {
951 if (rank_stencils[i] != -1) /* ne or nw */
952 stencil_contrib_cnt[stencil_i]++;
953 }
954
955 if (ndim > 2)
956 {
957 for (j= 1; j<= 2; j++)
958 {
959 for (i= 6; i<= 8; i++)
960 {
961 if (rank_stencils[j*9+i] != -1) /* bottom or top planes */
962 stencil_contrib_cnt[stencil_i]++;
963 }
964 }
965 }
966 max_contribut_size= hypre_max( max_contribut_size,
967 stencil_contrib_cnt[stencil_i] );
968 }
969
970
971 if (stencil_i != -1)
972 {
973 stencil_contrib_i[stencil_i]=
974 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
975 weight_contrib_i[stencil_i] =
976 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
977 sum= 0.0;
978 k= 0;
979
980 stencil_contrib_i[stencil_i][k]= stencil_i;
981 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
982 abs_stencil_shape );
983 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
984 sum+= weights[abs_stencil_shape];
985
986 for (i= 7; i<= 8; i++)
987 {
988 if (rank_stencils[i] != -1)
989 {
990 stencil_contrib_i[stencil_i][k] = rank_stencils[i];
991
992 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
993 abs_stencil_shape );
994 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
995 sum+= weights[abs_stencil_shape];
996 }
997 }
998
999 if (ndim > 2)
1000 {
1001 for (j= 1; j<= 2; j++)
1002 {
1003 for (i= 6; i<= 8; i++)
1004 {
1005 if (rank_stencils[j*9+i] != -1)
1006 {
1007 stencil_contrib_i[stencil_i][k] = rank_stencils[j*9+i];
1008
1009 AbsStencilShape(
1010 hypre_StructStencilElement(stencils,rank_stencils[j*9+i]),
1011 abs_stencil_shape );
1012 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1013 sum+= weights[abs_stencil_shape];
1014 }
1015 }
1016 }
1017 }
1018
1019 for (i= 0; i< k ; i++)
1020 {
1021 weight_contrib_i[stencil_i][i]/= sum;
1022 }
1023 }
1024
1025 /*-----------------------------------------------------------------
1026 * If only 2-d, extract the corner indices.
1027 *-----------------------------------------------------------------*/
1028 if (ndim == 2)
1029 {
1030 /* corners: ne & nw */
1031 for (i= 4; i<= 5; i++)
1032 {
1033 stencil_i= rank_stencils[i];
1034 if (stencil_i != -1)
1035 {
1036 stencil_contrib_cnt[stencil_i]++;
1037 stencil_contrib_i[stencil_i]= hypre_TAlloc(HYPRE_Int, 1, HYPRE_MEMORY_HOST);
1038 weight_contrib_i[stencil_i] = hypre_TAlloc(HYPRE_Real, 1, HYPRE_MEMORY_HOST);
1039 stencil_contrib_i[stencil_i][0]= stencil_i;
1040 weight_contrib_i[stencil_i][0] = weights[0];
1041 }
1042 }
1043
1044 /* corners: se & sw */
1045 for (i= 7; i<= 8; i++)
1046 {
1047 stencil_i= rank_stencils[i];
1048 if (stencil_i != -1)
1049 {
1050 stencil_contrib_cnt[stencil_i]++;
1051 stencil_contrib_i[stencil_i]= hypre_TAlloc(HYPRE_Int, 1, HYPRE_MEMORY_HOST);
1052 weight_contrib_i[stencil_i] = hypre_TAlloc(HYPRE_Real, 1, HYPRE_MEMORY_HOST);
1053 stencil_contrib_i[stencil_i][0]= stencil_i;
1054 weight_contrib_i[stencil_i][0] = weights[0];
1055 }
1056 }
1057 }
1058
1059 /*-----------------------------------------------------------------
1060 * Additional directions for 3-dim case
1061 *-----------------------------------------------------------------*/
1062 if (ndim > 2)
1063 {
1064 /* sides: top */
1065 stencil_i= rank_stencils[9];
1066 if (stencil_i != -1)
1067 {
1068 stencil_contrib_cnt[stencil_i]++;
1069 for (i=1; i<= 8; i++)
1070 {
1071 if (rank_stencils[9+i] != -1)
1072 stencil_contrib_cnt[stencil_i]++;
1073 }
1074 max_contribut_size= hypre_max( max_contribut_size,
1075 stencil_contrib_cnt[stencil_i] );
1076 }
1077
1078 if (stencil_i != -1)
1079 {
1080 stencil_contrib_i[stencil_i]=
1081 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1082 weight_contrib_i[stencil_i] =
1083 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1084 sum= 0.0;
1085 k= 0;
1086
1087 stencil_contrib_i[stencil_i][k]= stencil_i;
1088 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1089 abs_stencil_shape );
1090 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1091 sum+= weights[abs_stencil_shape];
1092
1093 for (i=1; i<= 8; i++)
1094 {
1095 if (rank_stencils[9+i] != -1)
1096 {
1097 stencil_contrib_i[stencil_i][k]= rank_stencils[9+i];
1098
1099 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[9+i]),
1100 abs_stencil_shape );
1101 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1102 sum+= weights[abs_stencil_shape];
1103 }
1104 }
1105
1106 for (i= 0; i< k ; i++)
1107 {
1108 weight_contrib_i[stencil_i][i]/= sum;
1109 }
1110 }
1111
1112 /* sides: bottom */
1113 stencil_i= rank_stencils[18];
1114 if (stencil_i != -1)
1115 {
1116 stencil_contrib_cnt[stencil_i]++;
1117 for (i=1; i<= 8; i++)
1118 {
1119 if (rank_stencils[18+i] != -1)
1120 stencil_contrib_cnt[stencil_i]++;
1121 }
1122 max_contribut_size= hypre_max( max_contribut_size,
1123 stencil_contrib_cnt[stencil_i] );
1124 }
1125
1126 if (stencil_i != -1)
1127 {
1128 stencil_contrib_i[stencil_i]=
1129 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1130 weight_contrib_i[stencil_i] =
1131 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1132 sum= 0.0;
1133 k= 0;
1134
1135 stencil_contrib_i[stencil_i][k]= stencil_i;
1136 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1137 abs_stencil_shape );
1138 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1139 sum+= weights[abs_stencil_shape];
1140
1141 for (i=1; i<= 8; i++)
1142 {
1143 if (rank_stencils[18+i] != -1)
1144 {
1145 stencil_contrib_i[stencil_i][k]= rank_stencils[18+i];
1146
1147 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[18+i]),
1148 abs_stencil_shape );
1149 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1150 sum+= weights[abs_stencil_shape];
1151 }
1152 }
1153
1154 for (i= 0; i< k ; i++)
1155 {
1156 weight_contrib_i[stencil_i][i]/= sum;
1157 }
1158 }
1159
1160 /* edges: cne */
1161 stencil_i= rank_stencils[4];
1162 if (stencil_i != -1)
1163 {
1164 stencil_contrib_cnt[stencil_i]++;
1165 for (j=1; j<= 2; j++)
1166 {
1167 if (rank_stencils[j*9+4] != -1) /* bottom or top planes */
1168 stencil_contrib_cnt[stencil_i]++;
1169 }
1170 max_contribut_size= hypre_max( max_contribut_size,
1171 stencil_contrib_cnt[stencil_i] );
1172 }
1173
1174 if (stencil_i != -1)
1175 {
1176 stencil_contrib_i[stencil_i]=
1177 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1178 weight_contrib_i[stencil_i] =
1179 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1180 sum= 0.0;
1181 k= 0;
1182
1183 stencil_contrib_i[stencil_i][k]= stencil_i;
1184 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1185 abs_stencil_shape );
1186 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1187 sum+= weights[abs_stencil_shape];
1188
1189 for (j=1; j<= 2; j++)
1190 {
1191 if (rank_stencils[j*9+4] != -1)
1192 {
1193 stencil_contrib_i[stencil_i][k]= rank_stencils[j*9+4];
1194
1195 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[j*9+4]),
1196 abs_stencil_shape );
1197 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1198 sum+= weights[abs_stencil_shape];
1199 }
1200 }
1201
1202 for (i= 0; i< k ; i++)
1203 {
1204 weight_contrib_i[stencil_i][i]/= sum;
1205 }
1206 }
1207
1208 /* edges: cse */
1209 stencil_i= rank_stencils[7];
1210 if (stencil_i != -1)
1211 {
1212 stencil_contrib_cnt[stencil_i]++;
1213 for (j=1; j<= 2; j++)
1214 {
1215 if (rank_stencils[j*9+7] != -1) /* bottom or top planes */
1216 stencil_contrib_cnt[stencil_i]++;
1217 }
1218 max_contribut_size= hypre_max( max_contribut_size,
1219 stencil_contrib_cnt[stencil_i] );
1220 }
1221
1222 if (stencil_i != -1)
1223 {
1224 stencil_contrib_i[stencil_i]=
1225 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1226 weight_contrib_i[stencil_i] =
1227 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1228 sum= 0.0;
1229 k= 0;
1230
1231 stencil_contrib_i[stencil_i][k]= stencil_i;
1232 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1233 abs_stencil_shape );
1234 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1235 sum+= weights[abs_stencil_shape];
1236
1237 for (j=1; j<= 2; j++)
1238 {
1239 if (rank_stencils[j*9+7] != -1)
1240 {
1241 stencil_contrib_i[stencil_i][k]= rank_stencils[j*9+7];
1242
1243 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[j*9+7]),
1244 abs_stencil_shape );
1245 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1246 sum+= weights[abs_stencil_shape];
1247 }
1248 }
1249
1250 for (i= 0; i< k ; i++)
1251 {
1252 weight_contrib_i[stencil_i][i]/= sum;
1253 }
1254 }
1255
1256 /* edges: cnw */
1257 stencil_i= rank_stencils[5];
1258 if (stencil_i != -1)
1259 {
1260 stencil_contrib_cnt[stencil_i]++;
1261 for (j=1; j<= 2; j++)
1262 {
1263 if (rank_stencils[j*9+5] != -1) /* bottom or top planes */
1264 stencil_contrib_cnt[stencil_i]++;
1265 }
1266 max_contribut_size= hypre_max( max_contribut_size,
1267 stencil_contrib_cnt[stencil_i] );
1268 }
1269
1270 if (stencil_i != -1)
1271 {
1272 stencil_contrib_i[stencil_i]=
1273 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1274 weight_contrib_i[stencil_i] =
1275 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1276 sum= 0.0;
1277 k= 0;
1278
1279 stencil_contrib_i[stencil_i][k]= stencil_i;
1280 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1281 abs_stencil_shape );
1282 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1283 sum+= weights[abs_stencil_shape];
1284
1285 for (j=1; j<= 2; j++)
1286 {
1287 if (rank_stencils[j*9+5] != -1)
1288 {
1289 stencil_contrib_i[stencil_i][k]= rank_stencils[j*9+5];
1290
1291 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[j*9+5]),
1292 abs_stencil_shape );
1293 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1294 sum+= weights[abs_stencil_shape];
1295 }
1296 }
1297
1298 for (i= 0; i< k ; i++)
1299 {
1300 weight_contrib_i[stencil_i][i]/= sum;
1301 }
1302 }
1303
1304 /* edges: csw */
1305 stencil_i= rank_stencils[8];
1306 if (stencil_i != -1)
1307 {
1308 stencil_contrib_cnt[stencil_i]++;
1309 for (j=1; j<= 2; j++)
1310 {
1311 if (rank_stencils[j*9+8] != -1) /* bottom or top planes */
1312 stencil_contrib_cnt[stencil_i]++;
1313 }
1314 max_contribut_size= hypre_max( max_contribut_size,
1315 stencil_contrib_cnt[stencil_i] );
1316 }
1317
1318 if (stencil_i != -1)
1319 {
1320 stencil_contrib_i[stencil_i]=
1321 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1322 weight_contrib_i[stencil_i] =
1323 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1324 sum= 0.0;
1325 k= 0;
1326
1327 stencil_contrib_i[stencil_i][k]= stencil_i;
1328 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1329 abs_stencil_shape );
1330 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1331 sum+= weights[abs_stencil_shape];
1332
1333 for (j=1; j<= 2; j++)
1334 {
1335 if (rank_stencils[j*9+8] != -1)
1336 {
1337 stencil_contrib_i[stencil_i][k]= rank_stencils[j*9+8];
1338
1339 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[j*9+8]),
1340 abs_stencil_shape );
1341 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1342 sum+= weights[abs_stencil_shape];
1343 }
1344 }
1345
1346 for (i= 0; i< k ; i++)
1347 {
1348 weight_contrib_i[stencil_i][i]/= sum;
1349 }
1350 }
1351
1352 /* edges: top east */
1353 stencil_i= rank_stencils[10];
1354 if (stencil_i != -1)
1355 {
1356 stencil_contrib_cnt[stencil_i]++;
1357 for (i=3; i<= 6; i+=3)
1358 {
1359 if (rank_stencils[10+i] != -1)
1360 stencil_contrib_cnt[stencil_i]++;
1361 }
1362 max_contribut_size= hypre_max( max_contribut_size,
1363 stencil_contrib_cnt[stencil_i] );
1364 }
1365
1366 if (stencil_i != -1)
1367 {
1368 stencil_contrib_i[stencil_i]=
1369 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1370 weight_contrib_i[stencil_i] =
1371 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1372 sum= 0.0;
1373 k= 0;
1374
1375 stencil_contrib_i[stencil_i][k]= stencil_i;
1376 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1377 abs_stencil_shape );
1378 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1379 sum+= weights[abs_stencil_shape];
1380
1381 for (i=3; i<= 6; i+=3)
1382 {
1383 if (rank_stencils[10+i] != -1)
1384 {
1385 stencil_contrib_i[stencil_i][k]= rank_stencils[10+i];
1386
1387 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[10+i]),
1388 abs_stencil_shape );
1389 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1390 sum+= weights[abs_stencil_shape];
1391 }
1392 }
1393
1394 for (i= 0; i< k ; i++)
1395 {
1396 weight_contrib_i[stencil_i][i]/= sum;
1397 }
1398 }
1399
1400 /* edges: top west */
1401 stencil_i= rank_stencils[11];
1402 if (stencil_i != -1)
1403 {
1404 stencil_contrib_cnt[stencil_i]++;
1405 for (i=3; i<= 6; i+=3)
1406 {
1407 if (rank_stencils[11+i] != -1)
1408 stencil_contrib_cnt[stencil_i]++;
1409 }
1410 max_contribut_size= hypre_max( max_contribut_size,
1411 stencil_contrib_cnt[stencil_i] );
1412 }
1413
1414 if (stencil_i != -1)
1415 {
1416 stencil_contrib_i[stencil_i]=
1417 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1418 weight_contrib_i[stencil_i] =
1419 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1420 sum= 0.0;
1421 k= 0;
1422
1423 stencil_contrib_i[stencil_i][k]= stencil_i;
1424 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1425 abs_stencil_shape );
1426 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1427 sum+= weights[abs_stencil_shape];
1428
1429 for (i=3; i<= 6; i+=3)
1430 {
1431 if (rank_stencils[11+i] != -1)
1432 {
1433 stencil_contrib_i[stencil_i][k]= rank_stencils[11+i];
1434
1435 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[11+i]),
1436 abs_stencil_shape );
1437 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1438 sum+= weights[abs_stencil_shape];
1439 }
1440 }
1441
1442 for (i= 0; i< k ; i++)
1443 {
1444 weight_contrib_i[stencil_i][i]/= sum;
1445 }
1446 }
1447
1448 /* edges: top north */
1449 stencil_i= rank_stencils[12];
1450 if (stencil_i != -1)
1451 {
1452 stencil_contrib_cnt[stencil_i]++;
1453 for (i=13; i<= 14; i++)
1454 {
1455 if (rank_stencils[i] != -1)
1456 stencil_contrib_cnt[stencil_i]++;
1457 }
1458 max_contribut_size= hypre_max( max_contribut_size,
1459 stencil_contrib_cnt[stencil_i] );
1460 }
1461
1462 if (stencil_i != -1)
1463 {
1464 stencil_contrib_i[stencil_i]=
1465 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1466 weight_contrib_i[stencil_i] =
1467 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1468 sum= 0.0;
1469 k= 0;
1470
1471 stencil_contrib_i[stencil_i][k]= stencil_i;
1472 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1473 abs_stencil_shape );
1474 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1475 sum+= weights[abs_stencil_shape];
1476
1477 for (i=13; i<= 14; i++)
1478 {
1479 if (rank_stencils[i] != -1)
1480 {
1481 stencil_contrib_i[stencil_i][k]= rank_stencils[i];
1482
1483 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
1484 abs_stencil_shape );
1485 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1486 sum+= weights[abs_stencil_shape];
1487 }
1488 }
1489
1490 for (i= 0; i< k ; i++)
1491 {
1492 weight_contrib_i[stencil_i][i]/= sum;
1493 }
1494 }
1495
1496 /* edges: top south*/
1497 stencil_i= rank_stencils[15];
1498 if (stencil_i != -1)
1499 {
1500 stencil_contrib_cnt[stencil_i]++;
1501 for (i=16; i<= 17; i++)
1502 {
1503 if (rank_stencils[i] != -1)
1504 stencil_contrib_cnt[stencil_i]++;
1505 }
1506 max_contribut_size= hypre_max( max_contribut_size,
1507 stencil_contrib_cnt[stencil_i] );
1508 }
1509
1510 if (stencil_i != -1)
1511 {
1512 stencil_contrib_i[stencil_i]=
1513 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1514 weight_contrib_i[stencil_i] =
1515 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1516 sum= 0.0;
1517 k= 0;
1518
1519 stencil_contrib_i[stencil_i][k]= stencil_i;
1520 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1521 abs_stencil_shape );
1522 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1523 sum+= weights[abs_stencil_shape];
1524
1525 for (i=16; i<= 17; i++)
1526 {
1527 if (rank_stencils[i] != -1)
1528 {
1529 stencil_contrib_i[stencil_i][k]= rank_stencils[i];
1530
1531 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
1532 abs_stencil_shape );
1533 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1534 sum+= weights[abs_stencil_shape];
1535 }
1536 }
1537
1538 for (i= 0; i< k ; i++)
1539 {
1540 weight_contrib_i[stencil_i][i]/= sum;
1541 }
1542 }
1543
1544 /* edges: bottom east */
1545 stencil_i= rank_stencils[19];
1546 if (stencil_i != -1)
1547 {
1548 stencil_contrib_cnt[stencil_i]++;
1549 for (i=3; i<= 6; i+=3)
1550 {
1551 if (rank_stencils[19+i] != -1)
1552 stencil_contrib_cnt[stencil_i]++;
1553 }
1554 max_contribut_size= hypre_max( max_contribut_size,
1555 stencil_contrib_cnt[stencil_i] );
1556 }
1557
1558 if (stencil_i != -1)
1559 {
1560 stencil_contrib_i[stencil_i]=
1561 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1562 weight_contrib_i[stencil_i] =
1563 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1564 sum= 0.0;
1565 k= 0;
1566
1567 stencil_contrib_i[stencil_i][k]= stencil_i;
1568 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1569 abs_stencil_shape );
1570 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1571 sum+= weights[abs_stencil_shape];
1572
1573 for (i=3; i<= 6; i+=3)
1574 {
1575 if (rank_stencils[19+i] != -1)
1576 {
1577 stencil_contrib_i[stencil_i][k]= rank_stencils[19+i];
1578
1579 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[19+i]),
1580 abs_stencil_shape );
1581 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1582 sum+= weights[abs_stencil_shape];
1583 }
1584 }
1585
1586 for (i= 0; i< k ; i++)
1587 {
1588 weight_contrib_i[stencil_i][i]/= sum;
1589 }
1590 }
1591
1592 /* edges: bottom west */
1593 stencil_i= rank_stencils[20];
1594 if (stencil_i != -1)
1595 {
1596 stencil_contrib_cnt[stencil_i]++;
1597 for (i=3; i<= 6; i+=3)
1598 {
1599 if (rank_stencils[20+i] != -1)
1600 stencil_contrib_cnt[stencil_i]++;
1601 }
1602 max_contribut_size= hypre_max( max_contribut_size,
1603 stencil_contrib_cnt[stencil_i] );
1604 }
1605
1606 if (stencil_i != -1)
1607 {
1608 stencil_contrib_i[stencil_i]=
1609 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1610 weight_contrib_i[stencil_i] =
1611 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1612 sum= 0.0;
1613 k= 0;
1614
1615 stencil_contrib_i[stencil_i][k]= stencil_i;
1616 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1617 abs_stencil_shape );
1618 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1619 sum+= weights[abs_stencil_shape];
1620
1621 for (i=3; i<= 6; i+=3)
1622 {
1623 if (rank_stencils[20+i] != -1)
1624 {
1625 stencil_contrib_i[stencil_i][k]= rank_stencils[20+i];
1626
1627 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[20+i]),
1628 abs_stencil_shape );
1629 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1630 sum+= weights[abs_stencil_shape];
1631 }
1632 }
1633
1634 for (i= 0; i< k ; i++)
1635 {
1636 weight_contrib_i[stencil_i][i]/= sum;
1637 }
1638 }
1639
1640 /* edges: bottom north */
1641 stencil_i= rank_stencils[21];
1642 if (stencil_i != -1)
1643 {
1644 stencil_contrib_cnt[stencil_i]++;
1645 for (i=22; i<= 23; i++)
1646 {
1647 if (rank_stencils[i] != -1)
1648 stencil_contrib_cnt[stencil_i]++;
1649 }
1650 max_contribut_size= hypre_max( max_contribut_size,
1651 stencil_contrib_cnt[stencil_i] );
1652 }
1653
1654 if (stencil_i != -1)
1655 {
1656 stencil_contrib_i[stencil_i]=
1657 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1658 weight_contrib_i[stencil_i] =
1659 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1660 sum= 0.0;
1661 k= 0;
1662
1663 stencil_contrib_i[stencil_i][k]= stencil_i;
1664 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1665 abs_stencil_shape );
1666 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1667 sum+= weights[abs_stencil_shape];
1668
1669 for (i=22; i<= 23; i++)
1670 {
1671 if (rank_stencils[i] != -1)
1672 {
1673 stencil_contrib_i[stencil_i][k]= rank_stencils[i];
1674
1675 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
1676 abs_stencil_shape );
1677 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1678 sum+= weights[abs_stencil_shape];
1679 }
1680 }
1681
1682 for (i= 0; i< k ; i++)
1683 {
1684 weight_contrib_i[stencil_i][i]/= sum;
1685 }
1686 }
1687
1688 /* edges: bottom south*/
1689 stencil_i= rank_stencils[24];
1690 if (stencil_i != -1)
1691 {
1692 stencil_contrib_cnt[stencil_i]++;
1693 for (i=25; i<= 26; i++)
1694 {
1695 if (rank_stencils[i] != -1)
1696 stencil_contrib_cnt[stencil_i]++;
1697 }
1698 max_contribut_size= hypre_max( max_contribut_size,
1699 stencil_contrib_cnt[stencil_i] );
1700 }
1701
1702 if (stencil_i != -1)
1703 {
1704 stencil_contrib_i[stencil_i]=
1705 hypre_TAlloc(HYPRE_Int, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1706 weight_contrib_i[stencil_i] =
1707 hypre_TAlloc(HYPRE_Real, stencil_contrib_cnt[stencil_i], HYPRE_MEMORY_HOST);
1708 sum= 0.0;
1709 k= 0;
1710
1711 stencil_contrib_i[stencil_i][k]= stencil_i;
1712 AbsStencilShape( hypre_StructStencilElement(stencils,stencil_i),
1713 abs_stencil_shape );
1714 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1715 sum+= weights[abs_stencil_shape];
1716
1717 for (i=25; i<= 26; i++)
1718 {
1719 if (rank_stencils[i] != -1)
1720 {
1721 stencil_contrib_i[stencil_i][k]= rank_stencils[i];
1722
1723 AbsStencilShape(hypre_StructStencilElement(stencils,rank_stencils[i]),
1724 abs_stencil_shape );
1725 weight_contrib_i[stencil_i][k++] = weights[abs_stencil_shape];
1726 sum+= weights[abs_stencil_shape];
1727 }
1728 }
1729
1730 for (i= 0; i< k ; i++)
1731 {
1732 weight_contrib_i[stencil_i][i]/= sum;
1733 }
1734 }
1735
1736 /* corners*/
1737 for (j= 1; j<= 2; j++)
1738 {
1739 for (i= 4; i<= 5; i++)
1740 {
1741 stencil_i= rank_stencils[9*j+i];
1742 if (stencil_i != -1)
1743 {
1744 stencil_contrib_cnt[stencil_i]++;
1745 stencil_contrib_i[stencil_i]= hypre_TAlloc(HYPRE_Int, 1, HYPRE_MEMORY_HOST);
1746 weight_contrib_i[stencil_i] = hypre_TAlloc(HYPRE_Real, 1, HYPRE_MEMORY_HOST);
1747 stencil_contrib_i[stencil_i][0]= stencil_i;
1748 weight_contrib_i[stencil_i][0] = weights[0];
1749 }
1750 }
1751 for (i= 7; i<= 8; i++)
1752 {
1753 stencil_i= rank_stencils[9*j+i];
1754 if (stencil_i != -1)
1755 {
1756 stencil_contrib_cnt[stencil_i]++;
1757 stencil_contrib_i[stencil_i]= hypre_TAlloc(HYPRE_Int, 1, HYPRE_MEMORY_HOST);
1758 weight_contrib_i[stencil_i] = hypre_TAlloc(HYPRE_Real, 1, HYPRE_MEMORY_HOST);
1759 stencil_contrib_i[stencil_i][0]= stencil_i;
1760 weight_contrib_i[stencil_i][0] = weights[0];
1761 }
1762 }
1763 }
1764
1765 } /* if ndim > 2 */
1766 /*-----------------------------------------------------------------
1767 * Allocate for the temporary vector used in computing the
1768 * averages.
1769 *-----------------------------------------------------------------*/
1770 vals= hypre_CTAlloc(HYPRE_Real, max_contribut_size, HYPRE_MEMORY_HOST);
1771
1772 /*-----------------------------------------------------------------
1773 * coarse grid stencil contributor structures have been formed.
1774 *-----------------------------------------------------------------*/
1775 } /* if (stencil_marker == true) */
1776
1777 /*---------------------------------------------------------------------
1778 * Loop over gridboxes to average stencils
1779 *---------------------------------------------------------------------*/
1780 smatrix_var = hypre_SStructPMatrixSMatrix(A_pmatrix, var1, var2);
1781 crse_smatrix= hypre_SStructPMatrixSMatrix(A_crse, var1, var2);
1782
1783 /*---------------------------------------------------------------------
1784 * data ptrs to extract and fill in data.
1785 *---------------------------------------------------------------------*/
1786 a_ptrs = hypre_TAlloc(HYPRE_Real *, stencil_size, HYPRE_MEMORY_HOST);
1787 crse_ptrs= hypre_TAlloc(HYPRE_Real *, stencil_size, HYPRE_MEMORY_HOST);
1788
1789 hypre_ForBoxI(ci, cgrid_boxes)
1790 {
1791 cgrid_box= hypre_BoxArrayBox(cgrid_boxes, ci);
1792 fbox_interior_ci = fbox_interior[var1][ci];
1793 fbox_bdy_ci = fbox_bdy[var1][ci];
1794 interior_fboxi_ci= interior_fboxi[var1][ci];
1795 bdy_fboxi_ci = bdy_fboxi[var1][ci];
1796
1797 crse_dbox= hypre_BoxArrayBox(hypre_StructMatrixDataSpace(crse_smatrix),
1798 ci);
1799 /*------------------------------------------------------------------
1800 * grab the correct coarse grid pointers. These are the parent base
1801 * grids.
1802 *------------------------------------------------------------------*/
1803 for (i= 0; i< stencil_size; i++)
1804 {
1805 hypre_CopyIndex(hypre_StructStencilElement(stencils, i), stencil_shape_i);
1806 crse_ptrs[i]= hypre_StructMatrixExtractPointerByIndex(crse_smatrix,
1807 ci,
1808 stencil_shape_i);
1809 }
1810 /*------------------------------------------------------------------
1811 * Loop over the interior of each patch inside cgrid_box.
1812 *------------------------------------------------------------------*/
1813 hypre_ForBoxI(fi, fbox_interior_ci)
1814 {
1815 fgrid_box= hypre_BoxArrayBox(fbox_interior_ci, fi);
1816 /*--------------------------------------------------------------
1817 * grab the fine grid ptrs & create the offsets for the fine
1818 * grid ptrs.
1819 *--------------------------------------------------------------*/
1820 A_dbox= hypre_BoxArrayBox(hypre_StructMatrixDataSpace(smatrix_var),
1821 interior_fboxi_ci[fi]);
1822 for (i= 0; i< stencil_size; i++)
1823 {
1824 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
1825 stencil_shape_i);
1826 a_ptrs[i]=
1827 hypre_StructMatrixExtractPointerByIndex(smatrix_var,
1828 interior_fboxi_ci[fi],
1829 stencil_shape_i);
1830 }
1831
1832 /*---------------------------------------------------------------
1833 * Compute the offsets for pointing to the correct data.
1834 * Note that for 1-d, OffsetA[j][i]= 0. Therefore, this ptr
1835 * will be correct for 1-d.
1836 *---------------------------------------------------------------*/
1837 for (j= 0; j< 2; j++)
1838 {
1839 OffsetA[j][0]= 0;
1840 for (i= 1; i< refine_factors[j+1]; i++)
1841 {
1842 if (j == 0)
1843 {
1844 hypre_SetIndex3(index_temp, 0, i, 0);
1845 }
1846 else
1847 {
1848 hypre_SetIndex3(index_temp, 0, 0, i);
1849 }
1850 OffsetA[j][i] = hypre_BoxOffsetDistance(A_dbox, index_temp);
1851 }
1852 }
1853
1854 hypre_CopyIndex(hypre_BoxIMin(fgrid_box), fstart);
1855 hypre_CopyIndex(hypre_BoxIMax(fgrid_box), fend);
1856
1857 /* coarsen the interior patch box*/
1858 hypre_ClearIndex(index_temp);
1859 hypre_StructMapFineToCoarse(fstart, index_temp, stridef,
1860 hypre_BoxIMin(&fine_box));
1861 hypre_StructMapFineToCoarse(fend, index_temp, stridef,
1862 hypre_BoxIMax(&fine_box));
1863
1864 hypre_CopyIndex(hypre_BoxIMin(&fine_box), cstart);
1865
1866 /*----------------------------------------------------------------
1867 * Loop over interior grid box.
1868 *----------------------------------------------------------------*/
1869
1870 hypre_BoxGetSize(&fine_box, loop_size);
1871
1872 hypre_SerialBoxLoop2Begin(ndim, loop_size,
1873 A_dbox, fstart, stridef, iA,
1874 crse_dbox, cstart, stridec, iAc);
1875 {
1876 for (i= 0; i< stencil_size; i++)
1877 {
1878 rank= stencil_ranks[i];
1879
1880 /*------------------------------------------------------------
1881 * Loop over refinement agglomeration making up a coarse cell
1882 * when a non-centre stencil.
1883 *------------------------------------------------------------*/
1884 if (rank)
1885 {
1886 /*--------------------------------------------------------
1887 * Loop over refinement agglomeration extents making up a
1888 * a coarse cell.
1889 *--------------------------------------------------------*/
1890 hypre_CopyIndex(hypre_BoxIMin(shift_box[i]), index1);
1891 hypre_CopyIndex(hypre_BoxIMax(shift_box[i]), index2);
1892
1893 for (m= 0; m< stencil_contrib_cnt[i]; m++)
1894 {
1895 vals[m]= 0.0;
1896 }
1897
1898 /*--------------------------------------------------------
1899 * For 1-d, index1[l]= index2[l]= 0, l>=1. So
1900 * iA_shift_zyx= j,
1901 * which is correct. Similarly, 2-d is correct.
1902 *--------------------------------------------------------*/
1903 for (l= index1[2]; l<= index2[2]; l++)
1904 {
1905 iA_shift_z= iA + OffsetA[1][l];
1906 for (k= index1[1]; k<= index2[1]; k++)
1907 {
1908 iA_shift_zy= iA_shift_z + OffsetA[0][k];
1909 for (j= index1[0]; j<= index2[0]; j++)
1910 {
1911 iA_shift_zyx= iA_shift_zy + j;
1912
1913 for (m= 0; m< stencil_contrib_cnt[i]; m++)
1914 {
1915 stencil_i= stencil_contrib_i[i][m];
1916 vals[m]+= a_ptrs[stencil_i][iA_shift_zyx];
1917 }
1918 }
1919 }
1920 }
1921 /*----------------------------------------------------------
1922 * average & weight the contributions and place into coarse
1923 * stencil entry.
1924 *----------------------------------------------------------*/
1925 crse_ptrs[i][iAc]= 0.0;
1926 for (m= 0; m< stencil_contrib_cnt[i]; m++)
1927 {
1928 crse_ptrs[i][iAc]+= vals[m]*weight_contrib_i[i][m];
1929 }
1930 crse_ptrs[i][iAc]/= volume_shift_box[i];
1931
1932 } /* if (rank) */
1933 } /* for i */
1934
1935 /*------------------------------------------------------------------
1936 * centre stencil:
1937 * The centre stencil is computed so that the row sum is equal to
1938 * the sum of the row sums of the fine matrix. Uses the computed
1939 * coarse off-diagonal stencils.
1940 *
1941 * No fine-coarse interface for the interior boxes.
1942 *------------------------------------------------------------------*/
1943 hypre_CopyIndex(hypre_BoxIMin(&coarse_cell_box), index1);
1944 hypre_CopyIndex(hypre_BoxIMax(&coarse_cell_box), index2);
1945
1946 sum= 0.0;
1947 for (l= index1[2]; l<= index2[2]; l++)
1948 {
1949 iA_shift_z= iA + OffsetA[1][l];
1950 for (k= index1[1]; k<= index2[1]; k++)
1951 {
1952 iA_shift_zy= iA_shift_z + OffsetA[0][k];
1953 for (j= index1[0]; j<= index2[0]; j++)
1954 {
1955 iA_shift_zyx= iA_shift_zy + j;
1956 for (m= 0; m< stencil_size; m++)
1957 {
1958 sum+= a_ptrs[m][iA_shift_zyx];
1959 }
1960 }
1961 }
1962 }
1963
1964 /*---------------------------------------------------------------
1965 * coarse centre coefficient- when away from the fine-coarse
1966 * interface, the centre coefficient is the sum of the
1967 * off-diagonal components.
1968 *---------------------------------------------------------------*/
1969 sum /= scaling;
1970 for (m= 0; m< stencil_size; m++)
1971 {
1972 rank= stencil_ranks[m];
1973 if (rank)
1974 {
1975 sum-= crse_ptrs[m][iAc];
1976 }
1977 }
1978 crse_ptrs[ rank_stencils[0] ][iAc]= sum;
1979 }
1980 hypre_SerialBoxLoop2End(iA, iAc);
1981 } /* end hypre_ForBoxI(fi, fbox_interior_ci) */
1982
1983 /*------------------------------------------------------------------
1984 * Loop over the boundaries of each patch inside cgrid_box.
1985 *------------------------------------------------------------------*/
1986 hypre_ForBoxArrayI(arrayi, fbox_bdy_ci)
1987 {
1988 fbox_bdy_ci_fi= hypre_BoxArrayArrayBoxArray(fbox_bdy_ci, arrayi);
1989 hypre_ForBoxI(fi, fbox_bdy_ci_fi)
1990 {
1991 fgrid_box= hypre_BoxArrayBox(fbox_bdy_ci_fi, fi);
1992
1993 /*-----------------------------------------------------------
1994 * grab the fine grid ptrs & create the offsets for the fine
1995 * grid ptrs.
1996 *-----------------------------------------------------------*/
1997 A_dbox= hypre_BoxArrayBox(hypre_StructMatrixDataSpace(smatrix_var),
1998 bdy_fboxi_ci[arrayi]);
1999 for (i= 0; i< stencil_size; i++)
2000 {
2001 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
2002 stencil_shape_i);
2003 a_ptrs[i]=
2004 hypre_StructMatrixExtractPointerByIndex(smatrix_var,
2005 bdy_fboxi_ci[arrayi],
2006 stencil_shape_i);
2007 }
2008
2009 /*--------------------------------------------------------------
2010 * Compute the offsets for pointing to the correct data.
2011 *--------------------------------------------------------------*/
2012 for (j= 0; j< 2; j++)
2013 {
2014 OffsetA[j][0]= 0;
2015 for (i= 1; i< refine_factors[j+1]; i++)
2016 {
2017 if (j == 0)
2018 {
2019 hypre_SetIndex3(index_temp, 0, i, 0);
2020 }
2021 else
2022 {
2023 hypre_SetIndex3(index_temp, 0, 0, i);
2024 }
2025 OffsetA[j][i] = hypre_BoxOffsetDistance(A_dbox, index_temp);
2026 }
2027 }
2028
2029 hypre_CopyIndex(hypre_BoxIMin(fgrid_box), fstart);
2030 hypre_CopyIndex(hypre_BoxIMax(fgrid_box), fend);
2031
2032 /* coarsen the patch box*/
2033 hypre_ClearIndex(index_temp);
2034 hypre_StructMapFineToCoarse(fstart, index_temp, stridef,
2035 hypre_BoxIMin(&fine_box));
2036 hypre_StructMapFineToCoarse(fend, index_temp, stridef,
2037 hypre_BoxIMax(&fine_box));
2038
2039 hypre_CopyIndex(hypre_BoxIMin(&fine_box), cstart);
2040
2041 /*--------------------------------------------------------------
2042 * Loop over boundary grid box.
2043 *--------------------------------------------------------------*/
2044
2045 hypre_BoxGetSize(&fine_box, loop_size);
2046
2047 hypre_SerialBoxLoop2Begin(ndim, loop_size,
2048 A_dbox, fstart, stridef, iA,
2049 crse_dbox, cstart, stridec, iAc);
2050 {
2051 zypre_BoxLoopGetIndex(lindex);
2052 for (i= 0; i< stencil_size; i++)
2053 {
2054 rank= stencil_ranks[i];
2055
2056 /*--------------------------------------------------------
2057 * Loop over refinement agglomeration making up a coarse
2058 * cell when a non-centre stencil.
2059 *--------------------------------------------------------*/
2060 if (rank)
2061 {
2062 /*-----------------------------------------------------
2063 * Loop over refinement agglomeration extents making up
2064 * a coarse cell.
2065 *-----------------------------------------------------*/
2066 hypre_CopyIndex(hypre_BoxIMin(shift_box[i]), index1);
2067 hypre_CopyIndex(hypre_BoxIMax(shift_box[i]), index2);
2068
2069 for (m= 0; m< stencil_contrib_cnt[i]; m++)
2070 {
2071 vals[m]= 0.0;
2072 }
2073
2074 for (l= index1[2]; l<= index2[2]; l++)
2075 {
2076 iA_shift_z= iA + OffsetA[1][l];
2077 for (k= index1[1]; k<= index2[1]; k++)
2078 {
2079 iA_shift_zy= iA_shift_z + OffsetA[0][k];
2080 for (j= index1[0]; j<= index2[0]; j++)
2081 {
2082 iA_shift_zyx= iA_shift_zy + j;
2083
2084 for (m= 0; m< stencil_contrib_cnt[i]; m++)
2085 {
2086 stencil_i= stencil_contrib_i[i][m];
2087 vals[m]+= a_ptrs[stencil_i][iA_shift_zyx];
2088 }
2089 }
2090 }
2091 }
2092 /*---------------------------------------------------------
2093 * average & weight the contributions and place into coarse
2094 * stencil entry.
2095 *---------------------------------------------------------*/
2096 crse_ptrs[i][iAc]= 0.0;
2097 for (m= 0; m< stencil_contrib_cnt[i]; m++)
2098 {
2099 crse_ptrs[i][iAc]+= vals[m]*weight_contrib_i[i][m];
2100 }
2101 crse_ptrs[i][iAc]/= volume_shift_box[i];
2102
2103 } /* if (rank) */
2104 } /* for i */
2105
2106 /*---------------------------------------------------------------
2107 * centre stencil:
2108 * The centre stencil is computed so that the row sum is equal to
2109 * th sum of the row sums of the fine matrix. Uses the computed
2110 * coarse off-diagonal stencils.
2111 *
2112 * Along the fine-coarse interface, we need to add the unstructured
2113 * connections.
2114 *---------------------------------------------------------------*/
2115 hypre_CopyIndex(hypre_BoxIMin(&coarse_cell_box), index1);
2116 hypre_CopyIndex(hypre_BoxIMax(&coarse_cell_box), index2);
2117
2118 temp3= hypre_CTAlloc(HYPRE_Real, volume_coarse_cell_box, HYPRE_MEMORY_HOST);
2119
2120 /*---------------------------------------------------------------
2121 * iA_shift_zyx is computed correctly for 1 & 2-d. Also,
2122 * ll= 0 for 2-d, and ll= kk= 0 for 1-d. Correct ptrs.
2123 *---------------------------------------------------------------*/
2124 for (l= index1[2]; l<= index2[2]; l++)
2125 {
2126 iA_shift_z= iA + OffsetA[1][l];
2127 ll = l*refine_factors[1]*refine_factors[0];
2128 for (k= index1[1]; k<= index2[1]; k++)
2129 {
2130 iA_shift_zy= iA_shift_z + OffsetA[0][k];
2131 kk = ll + k*refine_factors[0];
2132 for (j= index1[0]; j<= index2[0]; j++)
2133 {
2134 iA_shift_zyx= iA_shift_zy + j;
2135 jj = kk + j;
2136 for (m= 0; m< stencil_size; m++)
2137 {
2138 temp3[jj]+= a_ptrs[m][iA_shift_zyx];
2139 }
2140 }
2141 }
2142 }
2143
2144 /*------------------------------------------------------------
2145 * extract all unstructured connections. Note that we extract
2146 * from sstruct_matrix A, which already has been assembled.
2147 *------------------------------------------------------------*/
2148 if (nUventries > 0)
2149 {
2150 temp2= hypre_CTAlloc(HYPRE_Int, volume_coarse_cell_box, HYPRE_MEMORY_HOST);
2151 cnt1= 0;
2152 for (l= index1[2]; l<= index2[2]; l++)
2153 {
2154 ll= l*refine_factors[1]*refine_factors[0];
2155 for (k= index1[1]; k<= index2[1]; k++)
2156 {
2157 kk= ll + k*refine_factors[0];
2158 for (j= index1[0]; j<= index2[0]; j++)
2159 {
2160 jj= kk+ j;
2161
2162 hypre_SetIndex3(index_temp,
2163 j+lindex[0]*stridef[0],
2164 k+lindex[1]*stridef[1],
2165 l+lindex[2]*stridef[2]);
2166 hypre_AddIndexes(fstart, index_temp, 3, index_temp);
2167
2168 hypre_SStructGridFindBoxManEntry(grid, part_fine, index_temp,
2169 var1, &boxman_entry);
2170 hypre_SStructBoxManEntryGetGlobalRank(boxman_entry, index_temp,
2171 &rank, matrix_type);
2172
2173 found= falseV;
2174 i= hypre_SStructGraphIUVEntry(graph, 0);
2175 m= hypre_SStructGraphIUVEntry(graph, nUventries-1);
2176 if ((rank-startrank) >= i && (rank-startrank) <= m)
2177 {
2178 found= trueV;
2179 }
2180
2181 if (found)
2182 {
2183 Uventry= hypre_SStructGraphUVEntry(graph, rank-startrank);
2184
2185 if (Uventry != NULL)
2186 {
2187 nUentries= hypre_SStructUVEntryNUEntries(Uventry);
2188
2189 m= 0;
2190 for (i= 0; i< nUentries; i++)
2191 {
2192 if (hypre_SStructUVEntryToPart(Uventry, i)==part_crse)
2193 {
2194 m++;
2195 }
2196 } /* for (i= 0; i< nUentries; i++) */
2197
2198 temp2[jj]= m;
2199 cnt1 += m;
2200
2201 } /* if (Uventry != NULL) */
2202 } /* if (found) */
2203
2204 } /* for (j= index1[0]; j<= index2[0]; j++) */
2205 } /* for (k= index1[1]; k<= index2[1]; k++) */
2206 } /* for (l= index1[2]; l<= index2[2]; l++) */
2207
2208 ncols= hypre_TAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
2209 for (l= 0; l< cnt1; l++)
2210 {
2211 ncols[l]= 1;
2212 }
2213
2214 rows = hypre_TAlloc(HYPRE_BigInt, cnt1, HYPRE_MEMORY_HOST);
2215 cols = hypre_TAlloc(HYPRE_BigInt, cnt1, HYPRE_MEMORY_HOST);
2216 vals2= hypre_CTAlloc(HYPRE_Real, cnt1, HYPRE_MEMORY_HOST);
2217
2218 cnt1= 0;
2219 for (l= index1[2]; l<= index2[2]; l++)
2220 {
2221 ll= l*refine_factors[1]*refine_factors[0];
2222 for (k= index1[1]; k<= index2[1]; k++)
2223 {
2224 kk= ll + k*refine_factors[0];
2225 for (j= index1[0]; j<= index2[0]; j++)
2226 {
2227 jj= kk+ j;
2228
2229 hypre_SetIndex3(index_temp,
2230 j+lindex[0]*stridef[0],
2231 k+lindex[1]*stridef[1],
2232 l+lindex[2]*stridef[2]);
2233 hypre_AddIndexes(fstart, index_temp, 3, index_temp);
2234
2235 hypre_SStructGridFindBoxManEntry(grid, part_fine, index_temp,
2236 var1, &boxman_entry);
2237 hypre_SStructBoxManEntryGetGlobalRank(boxman_entry, index_temp,
2238 &rank, matrix_type);
2239
2240 found= falseV;
2241 if (nUventries > 0)
2242 {
2243 i= hypre_SStructGraphIUVEntry(graph, 0);
2244 m= hypre_SStructGraphIUVEntry(graph, nUventries-1);
2245 if ((HYPRE_Int)(rank-startrank) >= i && (HYPRE_Int)(rank-startrank) <= m)
2246 {
2247 found= trueV;
2248 }
2249 }
2250
2251 if (found)
2252 {
2253 Uventry= hypre_SStructGraphUVEntry(graph, (HYPRE_Int)(rank-startrank));
2254
2255 if (Uventry != NULL)
2256 {
2257 nUentries= hypre_SStructUVEntryNUEntries(Uventry);
2258 for (i= 0; i< nUentries; i++)
2259 {
2260 if (hypre_SStructUVEntryToPart(Uventry, i)==part_crse)
2261 {
2262 rows[cnt1]= rank;
2263 cols[cnt1++]= hypre_SStructUVEntryToRank(Uventry, i);
2264 }
2265
2266 } /* for (i= 0; i< nUentries; i++) */
2267 } /* if (Uventry != NULL) */
2268 } /* if (found) */
2269
2270 } /* for (j= index1[0]; j<= index2[0]; j++) */
2271 } /* for (k= index1[1]; k<= index2[1]; k++) */
2272 } /* for (l= index1[2]; l<= index2[2]; l++) */
2273
2274 HYPRE_IJMatrixGetValues(ij_A, cnt1, ncols, rows, cols, vals2);
2275
2276 cnt1= 0;
2277 for (l= index1[2]; l<= index2[2]; l++)
2278 {
2279 ll= l*refine_factors[1]*refine_factors[0];
2280 for (k= index1[1]; k<= index2[1]; k++)
2281 {
2282 kk= ll + k*refine_factors[0];
2283 for (j= index1[0]; j<= index2[0]; j++)
2284 {
2285 jj= kk+ j;
2286 for (m= 0; m< temp2[jj]; m++)
2287 {
2288 temp3[jj]+= vals2[cnt1];
2289 cnt1++;
2290 }
2291 temp2[jj]= 0; /* zero off for next time */
2292 } /* for (j= index1[0]; j<= index2[0]; j++) */
2293 } /* for (k= index1[1]; k<= index2[1]; k++) */
2294 } /* for (l= index1[2]; l<= index2[2]; l++) */
2295
2296 hypre_TFree(ncols, HYPRE_MEMORY_HOST);
2297 hypre_TFree(rows, HYPRE_MEMORY_HOST);
2298 hypre_TFree(cols, HYPRE_MEMORY_HOST);
2299 hypre_TFree(vals2, HYPRE_MEMORY_HOST);
2300 hypre_TFree(temp2, HYPRE_MEMORY_HOST);
2301
2302 } /* if Uventries > 0 */
2303
2304 sum= 0.0;
2305 for (l= index1[2]; l<= index2[2]; l++)
2306 {
2307 ll= l*refine_factors[1]*refine_factors[0];
2308 for (k= index1[1]; k<= index2[1]; k++)
2309 {
2310 kk= ll + k*refine_factors[0];
2311 for (j= index1[0]; j<= index2[0]; j++)
2312 {
2313 jj= kk + j;
2314 sum+= temp3[jj];
2315 }
2316 }
2317 }
2318
2319 sum /= scaling;
2320 crse_ptrs[ rank_stencils[0] ][iAc]= sum;
2321
2322 hypre_TFree(temp3, HYPRE_MEMORY_HOST);
2323
2324 }
2325 hypre_SerialBoxLoop2End(iA, iAc);
2326
2327 } /* hypre_ForBoxI(fi, fbox_bdy_ci_fi) */
2328 } /* hypre_ForBoxArrayI(arrayi, fbox_bdy_ci) */
2329 } /* hypre_ForBoxI(ci, cgrid_boxes) */
2330
2331 hypre_TFree(a_ptrs, HYPRE_MEMORY_HOST);
2332 hypre_TFree(crse_ptrs, HYPRE_MEMORY_HOST);
2333
2334 } /* if (stencils != NULL) */
2335 } /* end var2 */
2336 } /* end var1 */
2337
2338 if (stencil_contrib_cnt)
2339 hypre_TFree(stencil_contrib_cnt, HYPRE_MEMORY_HOST);
2340 if (stencil_ranks)
2341 hypre_TFree(stencil_ranks, HYPRE_MEMORY_HOST);
2342 if (volume_shift_box)
2343 hypre_TFree(volume_shift_box, HYPRE_MEMORY_HOST);
2344 if (vals)
2345 hypre_TFree(vals, HYPRE_MEMORY_HOST);
2346
2347 if (shift_box)
2348 {
2349 for (j= 0; j< stencil_size; j++)
2350 {
2351 if (shift_box[j])
2352 hypre_BoxDestroy(shift_box[j]);
2353 }
2354 hypre_TFree(shift_box, HYPRE_MEMORY_HOST);
2355 }
2356
2357 if (stencil_contrib_i)
2358 {
2359 for (j= 1; j< max_stencil_size; j++)
2360 {
2361 stencil_i= rank_stencils[j];
2362 if (stencil_i != -1)
2363 {
2364 if (stencil_contrib_i[stencil_i])
2365 hypre_TFree(stencil_contrib_i[stencil_i], HYPRE_MEMORY_HOST);
2366 }
2367 }
2368 hypre_TFree(stencil_contrib_i, HYPRE_MEMORY_HOST);
2369 }
2370
2371 if (weight_contrib_i)
2372 {
2373 for (j= 1; j< max_stencil_size; j++)
2374 {
2375 stencil_i= rank_stencils[j];
2376 if (stencil_i != -1)
2377 {
2378 if (weight_contrib_i[stencil_i])
2379 hypre_TFree(weight_contrib_i[stencil_i], HYPRE_MEMORY_HOST);
2380 }
2381 }
2382 hypre_TFree(weight_contrib_i, HYPRE_MEMORY_HOST);
2383 }
2384
2385 if (rank_stencils)
2386 hypre_TFree(rank_stencils, HYPRE_MEMORY_HOST);
2387
2388 if (OffsetA)
2389 {
2390 for (j= 0; j< 2; j++)
2391 {
2392 if (OffsetA[j])
2393 hypre_TFree(OffsetA[j], HYPRE_MEMORY_HOST);
2394 }
2395 hypre_TFree(OffsetA, HYPRE_MEMORY_HOST);
2396 }
2397
2398 /*--------------------------------------------------------------------------
2399 * STEP 2:
2400 *
2401 * Interface coarsening: fine-to-coarse connections. We are
2402 * assuming that only like-variables couple along interfaces.
2403 *
2404 * The task is to coarsen all the fine-to-coarse unstructured
2405 * connections and to compute coarse coefficients along the
2406 * interfaces (coarse-to-fine coefficients are obtained from these
2407 * computed values assuming symmetry). This involves
2408 * 1) scanning over the graph entries to find the locations of
2409 * the unstructure connections;
2410 * 2) determining the stencil shape of the coarsened connections;
2411 * 3) averaging the unstructured coefficients to compute
2412 * coefficient entries for the interface stencils;
2413 * 4) determining the weights of the interface stencil coefficients
2414 * to construct the structured coarse grid matrix along the
2415 * interfaces.
2416 *
2417 * We perform this task by
2418 * 1) scanning over the graph entries to group the locations
2419 * of the fine-to-coarse connections wrt the boxes of the
2420 * fine grid. Temporary vectors storing the Uventries indices
2421 * and the number of connections for each box will be created;
2422 * 2) for each fine grid box, group the fine-to-coarse connections
2423 * with respect to the connected coarse nodes. Temporary vectors
2424 * storing the Uventry indices and the Uentry indices for each
2425 * coarse node will be created (i.e., for a fixed coarse grid node,
2426 * record the fine node Uventries indices that connect to this
2427 * coarse node and Uentry index of the Uventry that contains
2428 * this coarse node.). The grouping is accomplished comparing the
2429 * ranks of the coarse nodes;
2430 * 3) using the Uventries and Uentry indices for each coarse node,
2431 * "coarsen" the fine grid connections to this coarse node to
2432 * create interface stencils (wrt to the coarse nodes- i.e.,
2433 * the centre of the stencil is at a coarse node). Also, find
2434 * the IJ rows and columns corresponding to all the fine-to-coarse
2435 * connections in a box, and extract the unstructured coefficients;
2436 * 4) looping over all coarse grid nodes connected to a fixed fine box,
2437 * compute the arithmetically averaged interface stencils;
2438 * 5) compare the underlying coarse grid structured stencil shape
2439 * to the interface stencil shape to determine how to weight the
2440 * averaged interface stencil coefficients.
2441 *
2442 * EXCEPTION: A NODE CAN CONTAIN ONLY UNSTRUCTURED CONNECTIONS
2443 * BETWEEN ONLY TWO AMR LEVELS- I.E., WE CANNOT HAVE A NODE THAT
2444 * IS ON THE INTERFACE OF MORE THAN TWO AMR LEVELS. CHANGES TO
2445 * HANDLE THIS LATTER CASE WILL INVOLVE THE SEARCH FOR f/c
2446 * CONNECTIONS.
2447 *-----------------------------------------------------------------*/
2448 if (nUventries > 0)
2449 {
2450 nvars = hypre_SStructPMatrixNVars(A_pmatrix);
2451
2452 for (var1= 0; var1< nvars; var1++)
2453 {
2454 /*-----------------------------------------------------------------
2455 * Yank out the structured stencils for this variable (only like
2456 * variables considered) and find their ranks.
2457 *-----------------------------------------------------------------*/
2458 stencils = hypre_SStructPMatrixSStencil(A_crse, var1, var1);
2459 stencil_size= hypre_StructStencilSize(stencils);
2460
2461 stencil_ranks= hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
2462 rank_stencils= hypre_TAlloc(HYPRE_Int, max_stencil_size, HYPRE_MEMORY_HOST);
2463 for (i= 0; i< stencil_size; i++)
2464 {
2465 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
2466 stencil_shape_i);
2467 MapStencilRank( stencil_shape_i, stencil_ranks[i] );
2468 rank_stencils[ stencil_ranks[i] ] = i;
2469 }
2470 /*-----------------------------------------------------------------
2471 * qsort the ranks into ascending order
2472 *-----------------------------------------------------------------*/
2473 hypre_qsort0(stencil_ranks, 0, stencil_size-1);
2474
2475 crse_smatrix= hypre_SStructPMatrixSMatrix(A_crse, var1, var1);
2476 cgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_crse), var1);
2477 cgrid_boxes= hypre_StructGridBoxes(cgrid);
2478
2479 fgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_pmatrix), var1);
2480 fgrid_boxes= hypre_StructGridBoxes(fgrid);
2481
2482 box_starts= hypre_CTAlloc(HYPRE_Int, hypre_BoxArraySize(fgrid_boxes), HYPRE_MEMORY_HOST);
2483 box_ends = hypre_CTAlloc(HYPRE_Int, hypre_BoxArraySize(fgrid_boxes), HYPRE_MEMORY_HOST);
2484 hypre_SStructGraphFindSGridEndpts(graph, part_fine, var1, myid,
2485 0, box_starts);
2486 hypre_SStructGraphFindSGridEndpts(graph, part_fine, var1, myid,
2487 1, box_ends);
2488
2489 /*-----------------------------------------------------------------
2490 * Step 1: scanning over the graph entries to group the locations
2491 * of the unstructured connections wrt to fine grid boxes.
2492 *
2493 * Count the components that couple for each box.
2494 *
2495 * box_graph_indices[fi]= array of Uventries indices in box fi.
2496 * box_graph_cnts[fi] = number of Uventries in box fi.
2497 * cdata_space_rank[ci] = begin offset rank of coarse data_space
2498 * box ci.
2499 *-----------------------------------------------------------------*/
2500 box_array_size = hypre_BoxArraySize(fgrid_boxes);
2501 cbox_array_size = hypre_BoxArraySize(cgrid_boxes);
2502 box_graph_indices= hypre_CTAlloc(HYPRE_Int *, box_array_size, HYPRE_MEMORY_HOST);
2503 box_graph_cnts = hypre_CTAlloc(HYPRE_Int , box_array_size, HYPRE_MEMORY_HOST);
2504
2505 data_space = hypre_StructMatrixDataSpace(crse_smatrix);
2506 cdata_space_ranks= hypre_CTAlloc(HYPRE_Int, cbox_array_size, HYPRE_MEMORY_HOST);
2507 cdata_space_ranks[0]= 0;
2508 for (i= 1; i< cbox_array_size; i++)
2509 {
2510 cdata_space_ranks[i]= cdata_space_ranks[i-1]+
2511 hypre_BoxVolume(hypre_BoxArrayBox(data_space, i-1));
2512 }
2513
2514 /*-----------------------------------------------------------------
2515 * Scanning obtained by searching iUventries between the start
2516 * and end of a fine box. Binary search used to find the interval
2517 * between these two endpts. Index (-1) returned if no interval
2518 * bounds found. Note that if start has positive index, then end
2519 * must have a positive index also.
2520 *-----------------------------------------------------------------*/
2521 for (fi= 0; fi< box_array_size; fi++)
2522 {
2523 i= hypre_LowerBinarySearch(iUventries, box_starts[fi], nUventries);
2524 if (i >= 0)
2525 {
2526 j= hypre_UpperBinarySearch(iUventries, box_ends[fi], nUventries);
2527 box_graph_indices[fi]= hypre_TAlloc(HYPRE_Int, j-i+1, HYPRE_MEMORY_HOST);
2528
2529 for (k= 0; k< (j-i+1); k++)
2530 {
2531 Uventry= hypre_SStructGraphUVEntry(graph,
2532 iUventries[i+k]);
2533
2534 for (m= 0; m< hypre_SStructUVEntryNUEntries(Uventry); m++)
2535 {
2536 if (hypre_SStructUVEntryToPart(Uventry, m) == part_crse)
2537 {
2538 box_graph_indices[fi][box_graph_cnts[fi]]= iUventries[i+k];
2539 box_graph_cnts[fi]++;
2540 break;
2541 }
2542 } /* for (m= 0; m< hypre_SStructUVEntryNUEntries(Uventry); m++) */
2543 } /* for (k= 0; k< (j-i+1); k++) */
2544 } /* if (i >= 0) */
2545 } /* for (fi= 0; fi< box_array_size; fi++) */
2546
2547 /*-----------------------------------------------------------------
2548 * Step 2:
2549 * Determine and group the fine-to-coarse connections in a box.
2550 * Grouped according to the coarsened fine grid interface nodes.
2551 *
2552 * box_ranks = ranks of coarsened fine grid interface
2553 * nodes.
2554 * box_connections = counter for the distinct coarsened fine
2555 * grid interface nodes. This can be
2556 * used to group all the Uventries of a
2557 * coarsened fine grid node.
2558 * cindex[l] = the hypre_Index of coarsen node l.
2559 * parents_cnodes[l] = parent box that contains the coarsened
2560 * fine grid interface node l.
2561 * fine_interface_ranks[l]= rank of coarsened fine grid interface
2562 * node l.
2563 * box_ranks_cnt[l] = counter for no. of Uventries for
2564 * coarsened node l.
2565 * coarse_contrib_Uv[l] = Uventry indices for Uventries that
2566 * contain fine-to-coarse connections of
2567 * coarse node l.
2568 *-----------------------------------------------------------------*/
2569 for (fi= 0; fi< box_array_size; fi++)
2570 {
2571 /*-------------------------------------------------------------
2572 * Determine the coarse data ptrs corresponding to fine box fi.
2573 * These are needed in assigning the averaged unstructured
2574 * coefficients.
2575 *
2576 * Determine how many distinct coarse grid nodes are in the
2577 * unstructured connection for a given box. Each node has a
2578 * structures.
2579 *
2580 * temp1 & temp2 are linked lists vectors used for grouping the
2581 * Uventries for a given coarse node.
2582 *-------------------------------------------------------------*/
2583 box_ranks = hypre_TAlloc(HYPRE_Int, box_graph_cnts[fi], HYPRE_MEMORY_HOST);
2584 box_connections = hypre_TAlloc(HYPRE_Int, box_graph_cnts[fi], HYPRE_MEMORY_HOST);
2585 parents = hypre_TAlloc(HYPRE_Int, box_graph_cnts[fi], HYPRE_MEMORY_HOST);
2586 temp1 = hypre_CTAlloc(HYPRE_Int, box_graph_cnts[fi]+1, HYPRE_MEMORY_HOST);
2587 temp2 = hypre_CTAlloc(HYPRE_Int, box_graph_cnts[fi], HYPRE_MEMORY_HOST);
2588 Uv_cindex = hypre_TAlloc(hypre_Index, box_graph_cnts[fi], HYPRE_MEMORY_HOST);
2589
2590 /*-------------------------------------------------------------
2591 * determine the parent box of this fgrid_box.
2592 *-------------------------------------------------------------*/
2593 hypre_ClearIndex(index_temp);
2594 for (i= 0; i < box_graph_cnts[fi]; i++)
2595 {
2596 Uventry = Uventries[box_graph_indices[fi][i]];
2597
2598 /*-------------------------------------------------------------
2599 * Coarsen the fine grid interface nodes and then get their
2600 * ranks. The correct coarse grid is needed to determine the
2601 * correct data_box.
2602 * Save the rank of the coarsened index & the parent box id.
2603 *-------------------------------------------------------------*/
2604 hypre_CopyIndex(hypre_SStructUVEntryIndex(Uventry), index);
2605 hypre_StructMapFineToCoarse(index, index_temp, stridef, Uv_cindex[i]);
2606 hypre_BoxSetExtents(&fine_box, Uv_cindex[i], Uv_cindex[i]);
2607
2608 for (j= 0; j< cboxi_fcnt[var1][fi]; j++)
2609 {
2610 ci= cboxi_fboxes[var1][fi][j];
2611 cgrid_box= hypre_BoxArrayBox(cgrid_boxes, ci);
2612 hypre_IntersectBoxes(&fine_box, cgrid_box, &intersect_box);
2613 if (hypre_BoxVolume(&intersect_box) > 0)
2614 {
2615 break;
2616 }
2617 }
2618
2619 parents[i] = ci;
2620 box_ranks[i]= cdata_space_ranks[ci] +
2621 hypre_BoxIndexRank(hypre_BoxArrayBox(data_space, ci),
2622 Uv_cindex[i]);
2623 }
2624
2625 /*---------------------------------------------------------------
2626 * Determine and "group" the Uventries using the box_ranks.
2627 * temp2 stores the Uventries indices for a coarsen node.
2628 *---------------------------------------------------------------*/
2629 cnt1= 0;
2630 j = 0;
2631 temp1[cnt1]= j;
2632
2633 for (i= 0; i< box_graph_cnts[fi]; i++)
2634 {
2635 if (box_ranks[i] != -1)
2636 {
2637 k = box_ranks[i];
2638 box_connections[i]= cnt1;
2639 temp2[j++] = box_graph_indices[fi][i];
2640
2641 for (l= i+1; l< box_graph_cnts[fi]; l++)
2642 {
2643 if (box_ranks[l] == k)
2644 {
2645 box_connections[l]= cnt1;
2646 temp2[j++] = box_graph_indices[fi][l];
2647 box_ranks[l] =-1;
2648 }
2649 }
2650 cnt1++;
2651 temp1[cnt1]= j;
2652 }
2653 }
2654
2655 /*-----------------------------------------------------------------
2656 * Store the graph entry info and other index info for each coarse
2657 * grid node.
2658 *-----------------------------------------------------------------*/
2659 parents_cnodes = hypre_TAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
2660 fine_interface_ranks= hypre_TAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
2661 box_ranks_cnt = hypre_CTAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
2662 coarse_contrib_Uv = hypre_TAlloc(HYPRE_Int *, cnt1, HYPRE_MEMORY_HOST);
2663 cindex = hypre_TAlloc(hypre_Index, cnt1, HYPRE_MEMORY_HOST);
2664
2665 for (i= 0; i< box_graph_cnts[fi]; i++)
2666 {
2667 if (box_ranks[i] != -1)
2668 {
2669 j = box_connections[i];
2670 parents_cnodes[j] = parents[i];
2671 fine_interface_ranks[j]=
2672 hypre_BoxIndexRank(hypre_BoxArrayBox(data_space, parents[i]),
2673 Uv_cindex[i]);
2674 hypre_CopyIndex(Uv_cindex[i], cindex[j]);
2675
2676 box_ranks_cnt[j] = temp1[j+1] - temp1[j];
2677 coarse_contrib_Uv[j] = hypre_TAlloc(HYPRE_Int, box_ranks_cnt[j], HYPRE_MEMORY_HOST);
2678
2679 l = temp1[j];
2680 for (k= 0; k< box_ranks_cnt[j]; k++)
2681 {
2682 coarse_contrib_Uv[j][k]= temp2[l+k];
2683 }
2684 }
2685 }
2686
2687 if (box_ranks)
2688 hypre_TFree(box_ranks, HYPRE_MEMORY_HOST);
2689 if (box_connections)
2690 hypre_TFree(box_connections, HYPRE_MEMORY_HOST);
2691 if (parents)
2692 hypre_TFree(parents, HYPRE_MEMORY_HOST);
2693 if (temp1)
2694 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
2695 if (temp2)
2696 hypre_TFree(temp2, HYPRE_MEMORY_HOST);
2697 if (Uv_cindex)
2698 hypre_TFree(Uv_cindex, HYPRE_MEMORY_HOST);
2699
2700 /*------------------------------------------------------------------------
2701 * Step 3:
2702 * Create the interface stencils.
2703 *
2704 * interface_max_stencil_ranks[i] = stencil_shape rank for each coarse
2705 * Uentry connection of coarsened node
2706 * i (i.e., the stencil_shape ranks of
2707 * the interface stencils at node i).
2708 * interface_max_stencil_cnt[i][m]= counter for number of Uentries
2709 * that describes a connection which
2710 * coarsens into stencil_shape rank m.
2711 * coarse_stencil_cnts[i] = counter for the no. of distinct
2712 * interface stencil_shapes (i.e., the
2713 * no. entries of the interface stencil).
2714 * interface_stencil_ranks[i][l] = stencil_shape rank for interface
2715 * stencil entry l, for coarse node i.
2716 * interface_rank_stencils[i][j] = interface stencil entry for
2717 * stencil_shape rank j, for node i.
2718 *------------------------------------------------------------------------*/
2719
2720 /*-----------------------------------------------------------------
2721 * Extract rows & cols info for extracting data from IJ matrix.
2722 * Extract for all connections for a box.
2723 *-----------------------------------------------------------------*/
2724 hypre_ClearIndex(index_temp);
2725
2726 nrows= 0;
2727 box_to_ranks_cnt= hypre_CTAlloc(HYPRE_Int, cnt1, HYPRE_MEMORY_HOST);
2728 for (i= 0; i< cnt1; i++)
2729 {
2730 for (j= 0; j< box_ranks_cnt[i]; j++)
2731 {
2732 Uventry = Uventries[ coarse_contrib_Uv[i][j] ];
2733 for (k= 0; k< hypre_SStructUVEntryNUEntries(Uventry); k++)
2734 {
2735 if (hypre_SStructUVEntryToPart(Uventry, k) == part_crse)
2736 {
2737 box_to_ranks_cnt[i]++;
2738 }
2739 }
2740 }
2741 nrows+= box_to_ranks_cnt[i];
2742 }
2743
2744 ncols= hypre_TAlloc(HYPRE_Int, nrows, HYPRE_MEMORY_HOST);
2745 for (i= 0; i< nrows; i++)
2746 {
2747 ncols[i]= 1;
2748 }
2749
2750 rows= hypre_TAlloc(HYPRE_BigInt, nrows, HYPRE_MEMORY_HOST);
2751 cols= hypre_TAlloc(HYPRE_BigInt, nrows, HYPRE_MEMORY_HOST);
2752 vals= hypre_CTAlloc(HYPRE_Real, nrows, HYPRE_MEMORY_HOST);
2753
2754 interface_max_stencil_ranks= hypre_TAlloc(HYPRE_Int *, cnt1, HYPRE_MEMORY_HOST);
2755 interface_max_stencil_cnt = hypre_TAlloc(HYPRE_Int *, cnt1, HYPRE_MEMORY_HOST);
2756 interface_rank_stencils = hypre_TAlloc(HYPRE_Int *, cnt1, HYPRE_MEMORY_HOST);
2757 interface_stencil_ranks = hypre_TAlloc(HYPRE_Int *, cnt1, HYPRE_MEMORY_HOST);
2758 coarse_stencil_cnt = hypre_CTAlloc(HYPRE_Int , cnt1, HYPRE_MEMORY_HOST);
2759
2760 k= 0;
2761 for (i= 0; i< cnt1; i++)
2762 {
2763 /*-----------------------------------------------------------------
2764 * for each coarse interface node, we get a stencil. We compute only
2765 * the ranks assuming a maximum size stencil of 27.
2766 *-----------------------------------------------------------------*/
2767 interface_max_stencil_ranks[i]= hypre_TAlloc(HYPRE_Int, box_to_ranks_cnt[i], HYPRE_MEMORY_HOST);
2768 interface_max_stencil_cnt[i] = hypre_CTAlloc(HYPRE_Int, max_stencil_size, HYPRE_MEMORY_HOST);
2769
2770 /*-----------------------------------------------------------------
2771 * conjugate the coarse node index for determining the stencil
2772 * shapes for the Uentry connections.
2773 *-----------------------------------------------------------------*/
2774 hypre_CopyIndex(cindex[i], index1);
2775 hypre_SetIndex3(index1, -index1[0], -index1[1], -index1[2]);
2776
2777 n= 0;
2778 for (j= 0; j< box_ranks_cnt[i]; j++)
2779 {
2780 /*--------------------------------------------------------------
2781 * extract the row rank for a given Uventry. Note that these
2782 * are the ranks in the grid of A. Therefore, we grab the index
2783 * from the nested_graph Uventry to determine the global rank.
2784 * With the rank, find the corresponding Uventry of the graph
2785 * of A. The to_ranks now can be extracted out.
2786 *--------------------------------------------------------------*/
2787 Uventry = Uventries[ coarse_contrib_Uv[i][j] ];
2788 hypre_CopyIndex(hypre_SStructUVEntryIndex(Uventry), index);
2789
2790 hypre_SStructGridFindBoxManEntry(grid, part_fine, index, var1, &boxman_entry);
2791 hypre_SStructBoxManEntryGetGlobalRank(boxman_entry, index, &rank, matrix_type);
2792
2793 Uventry= hypre_SStructGraphUVEntry(graph, rank-startrank);
2794 nUentries= hypre_SStructUVEntryNUEntries(Uventry);
2795
2796 for (l= 0; l< nUentries; l++)
2797 {
2798 if (hypre_SStructUVEntryToPart(Uventry, l) == part_crse)
2799 {
2800 to_rank = hypre_SStructUVEntryToRank(Uventry, l);
2801 rows[k] = rank;
2802 cols[k++]= to_rank;
2803
2804 /*---------------------------------------------------------
2805 * compute stencil shape for this Uentry.
2806 *---------------------------------------------------------*/
2807 hypre_CopyIndex( hypre_SStructUVEntryToIndex(Uventry,l),
2808 index );
2809 hypre_AddIndexes(index, index1, 3, index2);
2810
2811 MapStencilRank(index2, m);
2812 interface_max_stencil_ranks[i][n++]= m;
2813 interface_max_stencil_cnt[i][m]++;
2814 }
2815 }
2816 }
2817 hypre_TFree(coarse_contrib_Uv[i], HYPRE_MEMORY_HOST);
2818
2819 /*-----------------------------------------------------------------
2820 * Determine only the distinct stencil ranks for coarse node i.
2821 *-----------------------------------------------------------------*/
2822 l= 0;
2823 for (j= 0; j< max_stencil_size; j++)
2824 {
2825 if (interface_max_stencil_cnt[i][j])
2826 {
2827 l++;
2828 }
2829 }
2830
2831 coarse_stencil_cnt[i]= l;
2832 interface_stencil_ranks[i]= hypre_TAlloc(HYPRE_Int, l, HYPRE_MEMORY_HOST);
2833 interface_rank_stencils[i]= hypre_TAlloc(HYPRE_Int, max_stencil_size, HYPRE_MEMORY_HOST);
2834
2835 /*-----------------------------------------------------------------
2836 * For each stencil rank, assign one of the stencil_shape_i index.
2837 *-----------------------------------------------------------------*/
2838 l= 0;
2839 for (j= 0; j< max_stencil_size; j++)
2840 {
2841 if (interface_max_stencil_cnt[i][j])
2842 {
2843 interface_rank_stencils[i][j]= l;
2844 interface_stencil_ranks[i][l]= j;
2845 l++;
2846 }
2847 }
2848 } /* for (i= 0; i< cnt1; i++) */
2849
2850 hypre_TFree(coarse_contrib_Uv, HYPRE_MEMORY_HOST);
2851 hypre_TFree(box_ranks_cnt, HYPRE_MEMORY_HOST);
2852 hypre_TFree(cindex, HYPRE_MEMORY_HOST);
2853
2854 /*-----------------------------------------------------------------
2855 * Extract data from IJ matrix
2856 *-----------------------------------------------------------------*/
2857 HYPRE_IJMatrixGetValues(ij_A, nrows, ncols, rows, cols, vals);
2858
2859 hypre_TFree(ncols, HYPRE_MEMORY_HOST);
2860 hypre_TFree(rows, HYPRE_MEMORY_HOST);
2861 hypre_TFree(cols, HYPRE_MEMORY_HOST);
2862
2863 /*-----------------------------------------------------------------
2864 * Steps 4 & 5:
2865 * Compute the arithmetically averaged interface stencils,
2866 * and determine the interface stencil weights.
2867 *
2868 * stencil_vals[l] = averaged stencil coeff for interface
2869 * stencil entry l.
2870 * common_rank_stencils = final structured coarse stencil entries
2871 * for the stencil_shapes that the
2872 * interface stencils must collapse to.
2873 * common_stencil_ranks = final structured coarse stencil_shape
2874 * ranks for the stencil_shapes that the
2875 * interface stencils must collapse to.
2876 * common_stencil_i = stencil entry of the interface stencil
2877 * corresponding to the common
2878 * stencil_shape.
2879 *-----------------------------------------------------------------*/
2880 k= 0;
2881 for (i= 0; i< cnt1; i++)
2882 {
2883 stencil_vals= hypre_CTAlloc(HYPRE_Real, coarse_stencil_cnt[i], HYPRE_MEMORY_HOST);
2884
2885 /*-----------------------------------------------------------------
2886 * Compute the arithmetic stencil averages for coarse node i.
2887 *-----------------------------------------------------------------*/
2888 for (j= 0; j< box_to_ranks_cnt[i]; j++)
2889 {
2890 m= interface_max_stencil_ranks[i][j];
2891 l= interface_rank_stencils[i][m];
2892 stencil_vals[l]+= vals[k]/interface_max_stencil_cnt[i][m];
2893 k++;
2894 }
2895 hypre_TFree(interface_max_stencil_ranks[i], HYPRE_MEMORY_HOST);
2896 hypre_TFree(interface_max_stencil_cnt[i], HYPRE_MEMORY_HOST);
2897 hypre_TFree(interface_rank_stencils[i], HYPRE_MEMORY_HOST);
2898
2899 /*-----------------------------------------------------------------
2900 * Determine which stencil has to be formed. This is accomplished
2901 * by comparing the coarse grid stencil ranks with the computed
2902 * interface stencil ranks. We qsort (if there are more than one
2903 * rank) the ranks to give quick comparisons. Note that we need
2904 * to swap the elements of stencil_vals & fine_interface_ranks[i]'s
2905 * accordingly.
2906 *-----------------------------------------------------------------*/
2907
2908 sort= falseV;
2909 for (j= 0; j< (coarse_stencil_cnt[i]-1); j++)
2910 {
2911 if (interface_stencil_ranks[i][j] > interface_stencil_ranks[i][j+1])
2912 {
2913 sort= trueV;
2914 break;
2915 }
2916 }
2917
2918 if ( (coarse_stencil_cnt[i]>1) && (sort==trueV) )
2919 {
2920 temp1= hypre_TAlloc(HYPRE_Int, coarse_stencil_cnt[i], HYPRE_MEMORY_HOST);
2921 for (j= 0; j< coarse_stencil_cnt[i]; j++)
2922 {
2923 temp1[j]= j;
2924 }
2925
2926 hypre_qsort1(interface_stencil_ranks[i], (HYPRE_Real *) temp1, 0,
2927 coarse_stencil_cnt[i]-1);
2928
2929 /*---------------------------------------------------------------
2930 * swap the stencil_vals to agree with the rank swapping.
2931 *---------------------------------------------------------------*/
2932 temp3 = hypre_TAlloc(HYPRE_Real, coarse_stencil_cnt[i], HYPRE_MEMORY_HOST);
2933 for (j=0; j< coarse_stencil_cnt[i]; j++)
2934 {
2935 m = temp1[j];
2936 temp3[j] = stencil_vals[m];
2937 }
2938 for (j=0; j< coarse_stencil_cnt[i]; j++)
2939 {
2940 stencil_vals[j]= temp3[j];
2941 }
2942
2943 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
2944 hypre_TFree(temp3, HYPRE_MEMORY_HOST);
2945 }
2946
2947 /*-----------------------------------------------------------------
2948 * Compute the weights for the averaged stencil contributions.
2949 * We need to convert the ranks back to stencil_shapes and then
2950 * find the abs of the stencil shape.
2951 *-----------------------------------------------------------------*/
2952 temp3= hypre_TAlloc(HYPRE_Real, coarse_stencil_cnt[i], HYPRE_MEMORY_HOST);
2953 for (j=0; j< coarse_stencil_cnt[i]; j++)
2954 {
2955 InverseMapStencilRank(interface_stencil_ranks[i][j], index_temp);
2956 AbsStencilShape(index_temp, abs_stencil_shape);
2957 temp3[j]= weights[abs_stencil_shape];
2958 }
2959
2960 /*-----------------------------------------------------------------
2961 * Compare the coarse stencil and the interface stencil and
2962 * extract the common stencil shapes.
2963 * WE ARE ASSUMING THAT THE COARSE INTERFACE STENCIL HAS SOME
2964 * COMMON STENCIL SHAPE WITH THE COARSE STENCIL.
2965 *-----------------------------------------------------------------*/
2966 common_rank_stencils= hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
2967 common_stencil_ranks= hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
2968 common_stencil_i = hypre_TAlloc(HYPRE_Int, stencil_size, HYPRE_MEMORY_HOST);
2969
2970 l= 0;
2971 m= 0;
2972 for (j= 0; j< stencil_size; j++)
2973 {
2974 while( (l < coarse_stencil_cnt[i])
2975 && (stencil_ranks[j] > interface_stencil_ranks[i][l]) )
2976 {
2977 l++;
2978 }
2979
2980 if (l >= coarse_stencil_cnt[i])
2981 {
2982 break;
2983 }
2984 /*--------------------------------------------------------------
2985 * Check if a common stencil shape rank has been found.
2986 *--------------------------------------------------------------*/
2987 if ( (stencil_ranks[j] == interface_stencil_ranks[i][l])
2988 && (l < coarse_stencil_cnt[i]) )
2989 {
2990 common_rank_stencils[m]= rank_stencils[ stencil_ranks[j] ];
2991 common_stencil_ranks[m]= stencil_ranks[j];
2992 common_stencil_i[m++] = l;
2993 l++;
2994 }
2995 }
2996 /*-----------------------------------------------------------------
2997 * Find the contribution and weights for the averaged stencils.
2998 *-----------------------------------------------------------------*/
2999 for (j= 0; j< m; j++)
3000 {
3001 hypre_CopyIndex(hypre_StructStencilElement(
3002 stencils, common_rank_stencils[j]),
3003 stencil_shape_i);
3004 AbsStencilShape(stencil_shape_i, abs_stencil_shape);
3005
3006 crse_ptr= hypre_StructMatrixExtractPointerByIndex(crse_smatrix,
3007 parents_cnodes[i],
3008 stencil_shape_i);
3009
3010 /*-----------------------------------------------------------------
3011 * For a compact stencil (e.g., -1 <= hypre_Index[i] <= 1, i= 0-2),
3012 * the value of abs_stencil_shape can be used to determine the
3013 * stencil:
3014 * abs_stencil_shape= 3 only corners in 3-d
3015 * 2 corners in 2-d; or the centre plane
3016 * in 3-d, or e,w,n,s of the bottom
3017 * or top plane in 3-d
3018 * 1 e,w in 1-d; or e,w,n,s in 2-d;
3019 * or the centre plane in 3-d,
3020 * or c of the bottom or top plane
3021 * in 3-d
3022 * 0 c in 1-d, 2-d, or 3-d.
3023 *-----------------------------------------------------------------*/
3024
3025 switch(abs_stencil_shape)
3026 {
3027 case 3: /* corners of 3-d stencil */
3028
3029 l= common_stencil_i[j];
3030 crse_ptr[fine_interface_ranks[i]]= stencil_vals[l];
3031
3032 break;
3033
3034
3035 case 2: /* corners in 2-d or edges in 3-d */
3036
3037 if (ndim ==2)
3038 {
3039 l= common_stencil_i[j];
3040 crse_ptr[fine_interface_ranks[i]]= stencil_vals[l];
3041 }
3042
3043 else if (ndim == 3)
3044 {
3045 /*----------------------------------------------------------
3046 * The edge values are weighted sums of the averaged
3047 * coefficients. The weights and averaged coefficients must
3048 * be found. The contributions are found using the stencil
3049 * ranks and the stencil ordering
3050 * top: 14 12 13 centre: 5 3 4 bottom 23 21 22
3051 * 11 9 10 2 0 1 20 18 19
3052 * 17 15 16 8 6 7 26 24 25
3053 *----------------------------------------------------------*/
3054 l = common_stencil_ranks[j];
3055 temp1= hypre_TAlloc(HYPRE_Int, 2, HYPRE_MEMORY_HOST);
3056
3057 switch(l)
3058 {
3059 case 4: /* centre plane ne */
3060
3061 temp1[0]= 13;
3062 temp1[1]= 22;
3063 break;
3064
3065 case 5: /* centre plane nw */
3066
3067 temp1[0]= 14;
3068 temp1[1]= 23;
3069 break;
3070
3071 case 7: /* centre plane se */
3072
3073 temp1[0]= 16;
3074 temp1[1]= 25;
3075 break;
3076
3077 case 8: /* centre plane sw */
3078
3079 temp1[0]= 17;
3080 temp1[1]= 26;
3081 break;
3082
3083 case 10: /* top plane e */
3084
3085 temp1[0]= 13;
3086 temp1[1]= 16;
3087 break;
3088
3089 case 11: /* top plane w */
3090
3091 temp1[0]= 14;
3092 temp1[1]= 17;
3093 break;
3094
3095 case 12: /* top plane n */
3096
3097 temp1[0]= 13;
3098 temp1[1]= 14;
3099 break;
3100
3101 case 15: /* top plane s */
3102
3103 temp1[0]= 16;
3104 temp1[1]= 17;
3105 break;
3106
3107 case 19: /* bottom plane e */
3108
3109 temp1[0]= 22;
3110 temp1[1]= 25;
3111 break;
3112
3113 case 20: /* bottom plane w */
3114
3115 temp1[0]= 23;
3116 temp1[1]= 26;
3117 break;
3118
3119 case 21: /* bottom plane n */
3120
3121 temp1[0]= 22;
3122 temp1[1]= 23;
3123 break;
3124
3125 case 24: /* bottom plane s */
3126
3127 temp1[0]= 25;
3128 temp1[1]= 26;
3129 break;
3130 }
3131
3132
3133 /*-------------------------------------------------------
3134 * Add up the weighted contributions of the interface
3135 * stencils. This involves searching the ranks of
3136 * interface_stencil_ranks. The weights must be averaged.
3137 *-------------------------------------------------------*/
3138
3139 l= common_stencil_i[j];
3140 sum= temp3[l];
3141 sum_contrib= sum*stencil_vals[l];
3142
3143 n= 1;
3144 for (l= 0; l< 2; l++)
3145 {
3146 while ( (n < coarse_stencil_cnt[i])
3147 &&(interface_stencil_ranks[i][n] < temp1[l]) )
3148 {
3149 n++;
3150 }
3151
3152 if (n >= coarse_stencil_cnt[i])
3153 {
3154 break;
3155 }
3156
3157 if (interface_stencil_ranks[i][n] == temp1[l])
3158 {
3159 sum+= temp3[n];
3160 sum_contrib+= temp3[n]*stencil_vals[n];
3161 n++;
3162 }
3163 }
3164
3165 sum_contrib/= sum; /* average out the weights */
3166 l= common_stencil_i[j];
3167 crse_ptr[fine_interface_ranks[i]]= sum_contrib;
3168
3169 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
3170
3171 } /* else if (ndim == 3) */
3172
3173 break;
3174
3175 case 1: /* e,w in 1-d, or edges in 2-d, or faces in 3-d */
3176
3177 if (ndim == 1)
3178 {
3179 l= common_stencil_i[j];
3180 crse_ptr[fine_interface_ranks[i]]= stencil_vals[l];
3181 }
3182
3183 else if (ndim == 2)
3184 {
3185 l = common_stencil_ranks[j];
3186 temp1= hypre_TAlloc(HYPRE_Int, 2, HYPRE_MEMORY_HOST);
3187
3188 switch(l)
3189 {
3190 case 1: /* e */
3191
3192 temp1[0]= 4;
3193 temp1[1]= 7;
3194 break;
3195
3196 case 2: /* w */
3197
3198 temp1[0]= 5;
3199 temp1[1]= 8;
3200 break;
3201
3202 case 3: /* n */
3203
3204 temp1[0]= 4;
3205 temp1[1]= 5;
3206 break;
3207
3208 case 6: /* s */
3209
3210 temp1[0]= 7;
3211 temp1[1]= 8;
3212 break;
3213 }
3214
3215 /*-------------------------------------------------------
3216 * Add up the weighted contributions of the interface
3217 * stencils.
3218 *-------------------------------------------------------*/
3219
3220 l= common_stencil_i[j];
3221 sum= temp3[l];
3222 sum_contrib= sum*stencil_vals[l];
3223
3224 n= 1;
3225 for (l= 0; l< 2; l++)
3226 {
3227 while ( (n < coarse_stencil_cnt[i])
3228 &&(interface_stencil_ranks[i][n] < temp1[l]) )
3229 {
3230 n++;
3231 }
3232
3233 if (n >= coarse_stencil_cnt[i])
3234 {
3235 break;
3236 }
3237
3238 if (interface_stencil_ranks[i][n] == temp1[l])
3239 {
3240 sum+= temp3[n];
3241 sum_contrib+= temp3[n]*stencil_vals[n];
3242 n++;
3243 }
3244 }
3245
3246 sum_contrib/= sum; /* average out the weights */
3247 l= common_stencil_i[j];
3248 crse_ptr[fine_interface_ranks[i]]= sum_contrib;
3249
3250 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
3251
3252 } /* else if (ndim == 2) */
3253
3254 else /* 3-d */
3255 {
3256 l = common_stencil_ranks[j];
3257 temp1= hypre_TAlloc(HYPRE_Int, 8, HYPRE_MEMORY_HOST);
3258
3259 switch(l)
3260 {
3261 case 1: /* centre plane e */
3262
3263 temp1[0]= 4;
3264 temp1[1]= 7;
3265 temp1[2]= 10;
3266 temp1[3]= 13;
3267 temp1[4]= 16;
3268 temp1[5]= 19;
3269 temp1[6]= 22;
3270 temp1[7]= 25;
3271 break;
3272
3273 case 2: /* centre plane w */
3274
3275 temp1[0]= 5;
3276 temp1[1]= 8;
3277 temp1[2]= 11;
3278 temp1[3]= 14;
3279 temp1[4]= 17;
3280 temp1[5]= 20;
3281 temp1[6]= 23;
3282 temp1[7]= 26;
3283 break;
3284
3285 case 3: /* centre plane n */
3286
3287 temp1[0]= 4;
3288 temp1[1]= 5;
3289 temp1[2]= 12;
3290 temp1[3]= 13;
3291 temp1[4]= 14;
3292 temp1[5]= 21;
3293 temp1[6]= 22;
3294 temp1[7]= 23;
3295 break;
3296
3297 case 6: /* centre plane s */
3298
3299 temp1[0]= 7;
3300 temp1[1]= 8;
3301 temp1[2]= 15;
3302 temp1[3]= 16;
3303 temp1[4]= 17;
3304 temp1[5]= 24;
3305 temp1[6]= 25;
3306 temp1[7]= 26;
3307 break;
3308
3309 case 9: /* top plane c */
3310
3311 for (n= 0; n< 8; n++)
3312 {
3313 temp1[n]= 10+n;
3314 }
3315 break;
3316
3317 case 18: /* bottom plane c */
3318
3319 for (n= 0; n< 8; n++)
3320 {
3321 temp1[n]= 19+n;
3322 }
3323 break;
3324
3325 }
3326
3327 /*-------------------------------------------------------
3328 * Add up the weighted contributions of the interface
3329 * stencils.
3330 *-------------------------------------------------------*/
3331
3332 l= common_stencil_i[j];
3333 sum= temp3[l];
3334 sum_contrib= sum*stencil_vals[l];
3335
3336 n= 1;
3337 for (l= 0; l< 8; l++)
3338 {
3339 while ( (n < coarse_stencil_cnt[i])
3340 && (interface_stencil_ranks[i][n] < temp1[l]) )
3341 {
3342 n++;
3343 }
3344
3345 if (n >= coarse_stencil_cnt[i])
3346 {
3347 break;
3348 }
3349
3350 if (interface_stencil_ranks[i][n] == temp1[l])
3351 {
3352 sum+= temp3[n];
3353 sum_contrib+= temp3[n]*stencil_vals[n];
3354 n++;
3355 }
3356 }
3357
3358 sum_contrib/= sum; /* average out the weights */
3359 l= common_stencil_i[j];
3360 crse_ptr[fine_interface_ranks[i]]= sum_contrib;
3361
3362 hypre_TFree(temp1, HYPRE_MEMORY_HOST);
3363
3364 } /* else */
3365
3366 break;
3367
3368 } /* switch(abs_stencil_shape) */
3369 } /* for (j= 0; j< m; j++) */
3370
3371 hypre_TFree(interface_stencil_ranks[i], HYPRE_MEMORY_HOST);
3372
3373 hypre_TFree(stencil_vals, HYPRE_MEMORY_HOST);
3374 hypre_TFree(temp3, HYPRE_MEMORY_HOST);
3375 hypre_TFree(common_rank_stencils, HYPRE_MEMORY_HOST);
3376 hypre_TFree(common_stencil_ranks, HYPRE_MEMORY_HOST);
3377 hypre_TFree(common_stencil_ranks, HYPRE_MEMORY_HOST);
3378 hypre_TFree(common_stencil_i, HYPRE_MEMORY_HOST);
3379
3380 } /* for (i= 0; i< cnt1; i++) */
3381
3382 hypre_TFree(box_to_ranks_cnt, HYPRE_MEMORY_HOST);
3383 hypre_TFree(interface_max_stencil_ranks, HYPRE_MEMORY_HOST);
3384 hypre_TFree(interface_max_stencil_cnt, HYPRE_MEMORY_HOST);
3385 hypre_TFree(interface_rank_stencils, HYPRE_MEMORY_HOST);
3386 hypre_TFree(interface_stencil_ranks, HYPRE_MEMORY_HOST);
3387 hypre_TFree(coarse_stencil_cnt, HYPRE_MEMORY_HOST);
3388 hypre_TFree(fine_interface_ranks, HYPRE_MEMORY_HOST);
3389 hypre_TFree(parents_cnodes, HYPRE_MEMORY_HOST);
3390 hypre_TFree(vals, HYPRE_MEMORY_HOST);
3391
3392 /*-----------------------------------------------------------
3393 * Box fi is completed.
3394 *-----------------------------------------------------------*/
3395 } /* for (fi= 0; fi< box_array_size; fi++) */
3396
3397 hypre_TFree(stencil_ranks, HYPRE_MEMORY_HOST);
3398 hypre_TFree(rank_stencils, HYPRE_MEMORY_HOST);
3399 hypre_TFree(cdata_space_ranks, HYPRE_MEMORY_HOST);
3400 hypre_TFree(box_graph_cnts, HYPRE_MEMORY_HOST);
3401 for (i= 0; i< box_array_size; i++)
3402 {
3403 if (box_graph_indices[i])
3404 hypre_TFree(box_graph_indices[i], HYPRE_MEMORY_HOST);
3405 }
3406 hypre_TFree(box_graph_indices, HYPRE_MEMORY_HOST);
3407
3408 hypre_TFree(box_starts, HYPRE_MEMORY_HOST);
3409 hypre_TFree(box_ends, HYPRE_MEMORY_HOST);
3410 } /* for (var1= 0; var1< nvars; var1++) */
3411 } /* if (nUventries > 0) */
3412
3413
3414 /*--------------------------------------------------------------------------
3415 * STEP 3:
3416 * Coarsened f/c interface coefficients can be used to create the
3417 * centre components along the coarsened f/c nodes now. Loop over
3418 * the coarsened fbox_bdy's and set the centre stencils.
3419 *--------------------------------------------------------------------------*/
3420 hypre_ClearIndex(index_temp);
3421 for (var1= 0; var1< nvars; var1++)
3422 {
3423 /* only like variables couple. */
3424 smatrix_var = hypre_SStructPMatrixSMatrix(A_crse, var1, var1);
3425 stencils = hypre_SStructPMatrixSStencil(A_crse, var1, var1);
3426 stencil_size = hypre_StructStencilSize(stencils);
3427 a_ptrs = hypre_TAlloc(HYPRE_Real *, stencil_size, HYPRE_MEMORY_HOST);
3428
3429 rank_stencils= hypre_TAlloc(HYPRE_Int, max_stencil_size, HYPRE_MEMORY_HOST);
3430 for (i= 0; i< stencil_size; i++)
3431 {
3432 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
3433 stencil_shape_i);
3434 MapStencilRank(stencil_shape_i, rank);
3435 rank_stencils[rank]= i;
3436 }
3437 centre= rank_stencils[0];
3438
3439 cgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_crse), var1);
3440 cgrid_boxes= hypre_StructGridBoxes(cgrid);
3441
3442 hypre_ForBoxI(ci, cgrid_boxes)
3443 {
3444 A_dbox = hypre_BoxArrayBox(hypre_StructMatrixDataSpace(smatrix_var), ci);
3445 fbox_bdy_ci= fbox_bdy[var1][ci];
3446
3447 for (i= 0; i< stencil_size; i++)
3448 {
3449 hypre_CopyIndex(hypre_StructStencilElement(stencils, i),
3450 stencil_shape_i);
3451 a_ptrs[i]= hypre_StructMatrixExtractPointerByIndex(smatrix_var,
3452 ci,
3453 stencil_shape_i);
3454 }
3455
3456 /*------------------------------------------------------------------
3457 * Loop over the boundaries of each patch inside cgrid_box ci.
3458 * These patch boxes must be coarsened to get the correct extents.
3459 *------------------------------------------------------------------*/
3460 hypre_ForBoxArrayI(arrayi, fbox_bdy_ci)
3461 {
3462 fbox_bdy_ci_fi= hypre_BoxArrayArrayBoxArray(fbox_bdy_ci, arrayi);
3463 hypre_ForBoxI(fi, fbox_bdy_ci_fi)
3464 {
3465 fgrid_box= hypre_BoxArrayBox(fbox_bdy_ci_fi, fi);
3466 hypre_StructMapFineToCoarse(hypre_BoxIMin(fgrid_box), index_temp,
3467 stridef, hypre_BoxIMin(&fine_box));
3468 hypre_StructMapFineToCoarse(hypre_BoxIMax(fgrid_box), index_temp,
3469 stridef, hypre_BoxIMax(&fine_box));
3470
3471 hypre_CopyIndex(hypre_BoxIMin(&fine_box), cstart);
3472 hypre_BoxGetSize(&fine_box, loop_size);
3473
3474 #define DEVICE_VAR is_device_ptr(a_ptrs)
3475 hypre_BoxLoop1Begin(ndim, loop_size,
3476 A_dbox, cstart, stridec, iA);
3477 {
3478 HYPRE_Int i;
3479 for (i= 0; i< stencil_size; i++)
3480 {
3481 if (i != centre)
3482 {
3483 a_ptrs[centre][iA]-= a_ptrs[i][iA];
3484 }
3485 }
3486 }
3487 hypre_BoxLoop1End(iA);
3488 #undef DEVICE_VAR
3489
3490 } /* hypre_ForBoxI(fi, fbox_bdy_ci_fi) */
3491 } /* hypre_ForBoxArrayI(arrayi, fbox_bdy_ci) */
3492 } /* hypre_ForBoxI(ci, cgrid_boxes) */
3493
3494 hypre_TFree(a_ptrs, HYPRE_MEMORY_HOST);
3495 hypre_TFree(rank_stencils, HYPRE_MEMORY_HOST);
3496
3497 } /* for (var1= 0; var1< nvars; var1++) */
3498
3499 for (var1= 0; var1< nvars; var1++)
3500 {
3501 cgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_crse), var1);
3502 cgrid_boxes= hypre_StructGridBoxes(cgrid);
3503
3504 fgrid= hypre_SStructPGridSGrid(hypre_SStructPMatrixPGrid(A_pmatrix), var1);
3505 fgrid_boxes= hypre_StructGridBoxes(fgrid);
3506
3507 hypre_ForBoxI(ci, cgrid_boxes)
3508 {
3509 hypre_BoxArrayDestroy(fgrid_crse_extents[var1][ci]);
3510 hypre_BoxArrayDestroy(fbox_interior[var1][ci]);
3511 hypre_BoxArrayArrayDestroy(fbox_bdy[var1][ci]);
3512 hypre_TFree(interior_fboxi[var1][ci], HYPRE_MEMORY_HOST);
3513 hypre_TFree(bdy_fboxi[var1][ci], HYPRE_MEMORY_HOST);
3514 }
3515 hypre_TFree(fgrid_crse_extents[var1], HYPRE_MEMORY_HOST);
3516 hypre_TFree(fbox_interior[var1], HYPRE_MEMORY_HOST);
3517 hypre_TFree(fbox_bdy[var1], HYPRE_MEMORY_HOST);
3518 hypre_TFree(interior_fboxi[var1], HYPRE_MEMORY_HOST);
3519 hypre_TFree(bdy_fboxi[var1], HYPRE_MEMORY_HOST);
3520
3521 hypre_ForBoxI(fi, fgrid_boxes)
3522 {
3523 hypre_TFree(cboxi_fboxes[var1][fi], HYPRE_MEMORY_HOST);
3524 }
3525 hypre_TFree(cboxi_fboxes[var1], HYPRE_MEMORY_HOST);
3526 hypre_TFree(cboxi_fcnt[var1], HYPRE_MEMORY_HOST);
3527 }
3528 hypre_TFree(fgrid_crse_extents, HYPRE_MEMORY_HOST);
3529 hypre_TFree(fbox_interior, HYPRE_MEMORY_HOST);
3530 hypre_TFree(fbox_bdy, HYPRE_MEMORY_HOST);
3531 hypre_TFree(interior_fboxi, HYPRE_MEMORY_HOST);
3532 hypre_TFree(bdy_fboxi, HYPRE_MEMORY_HOST);
3533 hypre_TFree(cboxi_fboxes, HYPRE_MEMORY_HOST);
3534 hypre_TFree(cboxi_fcnt, HYPRE_MEMORY_HOST);
3535
3536 return 0;
3537 }
3538
3539