source: CIVL/examples/mpi-omp/AMG2013/struct_mv/struct_matrix.c@ beab7f2

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

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

git-svn-id: ​svn://vsl.cis.udel.edu/civl/trunk@5704 fb995dde-84ed-4084-dfe6-e5aef3e2452c
  • Property mode set to 100644
File size: 74.0 KB
Line 
1/*BHEADER**********************************************************************
2 * Copyright (c) 2008, Lawrence Livermore National Security, LLC.
3 * Produced at the Lawrence Livermore National Laboratory.
4 * This file is part of HYPRE. See file COPYRIGHT for details.
5 *
6 * HYPRE is free software; you can redistribute it and/or modify it under the
7 * terms of the GNU Lesser General Public License (as published by the Free
8 * Software Foundation) version 2.1 dated February 1999.
9 *
10 * $Revision: 2.4 $
11 ***********************************************************************EHEADER*/
12
13
14/******************************************************************************
15 *
16 * Member functions for hypre_StructMatrix class.
17 *
18 *****************************************************************************/
19
20#include "headers.h"
21
22/*--------------------------------------------------------------------------
23 * hypre_StructMatrixExtractPointerByIndex
24 * Returns pointer to data for stencil entry coresponding to
25 * `index' in `matrix'. If the index does not exist in the matrix's
26 * stencil, the NULL pointer is returned.
27 *--------------------------------------------------------------------------*/
28
29double *
30hypre_StructMatrixExtractPointerByIndex( hypre_StructMatrix *matrix,
31 int b,
32 hypre_Index index )
33{
34 hypre_StructStencil *stencil;
35 int rank;
36
37 stencil = hypre_StructMatrixStencil(matrix);
38 rank = hypre_StructStencilElementRank( stencil, index );
39
40 if ( rank >= 0 )
41 return hypre_StructMatrixBoxData(matrix, b, rank);
42 else
43 return NULL; /* error - invalid index */
44}
45
46/*--------------------------------------------------------------------------
47 * hypre_StructMatrixCreate
48 *--------------------------------------------------------------------------*/
49
50hypre_StructMatrix *
51hypre_StructMatrixCreate( MPI_Comm comm,
52 hypre_StructGrid *grid,
53 hypre_StructStencil *user_stencil )
54{
55 hypre_StructMatrix *matrix;
56
57 int ndim = hypre_StructGridDim(grid);
58 int i;
59
60 matrix = hypre_CTAlloc(hypre_StructMatrix, 1);
61
62 hypre_StructMatrixComm(matrix) = comm;
63 hypre_StructGridRef(grid, &hypre_StructMatrixGrid(matrix));
64 hypre_StructMatrixUserStencil(matrix) =
65 hypre_StructStencilRef(user_stencil);
66 hypre_StructMatrixDataAlloced(matrix) = 1;
67 hypre_StructMatrixRefCount(matrix) = 1;
68
69 /* set defaults */
70 hypre_StructMatrixSymmetric(matrix) = 0;
71 hypre_StructMatrixConstantCoefficient(matrix) = 0;
72 for (i = 0; i < 6; i++)
73 {
74 hypre_StructMatrixNumGhost(matrix)[i] = 0;
75 hypre_StructMatrixAddNumGhost(matrix)[i]= 0;
76 }
77
78 for (i = 0; i < ndim; i++)
79 {
80 hypre_StructMatrixNumGhost(matrix)[2*i] = 1;
81 hypre_StructMatrixNumGhost(matrix)[2*i+1] = 1;
82
83 hypre_StructMatrixAddNumGhost(matrix)[2*i]= 1;
84 hypre_StructMatrixAddNumGhost(matrix)[2*i+1]= 1;
85 }
86
87 hypre_StructMatrixOffProcAdd(matrix) = 0;
88
89 return matrix;
90}
91
92/*--------------------------------------------------------------------------
93 * hypre_StructMatrixRef
94 *--------------------------------------------------------------------------*/
95
96hypre_StructMatrix *
97hypre_StructMatrixRef( hypre_StructMatrix *matrix )
98{
99 hypre_StructMatrixRefCount(matrix) ++;
100
101 return matrix;
102}
103
104/*--------------------------------------------------------------------------
105 * hypre_StructMatrixDestroy
106 *--------------------------------------------------------------------------*/
107
108int
109hypre_StructMatrixDestroy( hypre_StructMatrix *matrix )
110{
111 int i;
112 int ierr = 0;
113
114 if (matrix)
115 {
116 hypre_StructMatrixRefCount(matrix) --;
117 if (hypre_StructMatrixRefCount(matrix) == 0)
118 {
119 if (hypre_StructMatrixDataAlloced(matrix))
120 {
121 hypre_SharedTFree(hypre_StructMatrixData(matrix));
122 }
123 hypre_CommPkgDestroy(hypre_StructMatrixCommPkg(matrix));
124
125 hypre_ForBoxI(i, hypre_StructMatrixDataSpace(matrix))
126 hypre_TFree(hypre_StructMatrixDataIndices(matrix)[i]);
127 hypre_TFree(hypre_StructMatrixDataIndices(matrix));
128
129 hypre_BoxArrayDestroy(hypre_StructMatrixDataSpace(matrix));
130
131 hypre_TFree(hypre_StructMatrixSymmElements(matrix));
132 hypre_StructStencilDestroy(hypre_StructMatrixUserStencil(matrix));
133 hypre_StructStencilDestroy(hypre_StructMatrixStencil(matrix));
134 hypre_StructGridDestroy(hypre_StructMatrixGrid(matrix));
135
136 hypre_TFree(matrix);
137 }
138 }
139
140 return ierr;
141}
142
143/*--------------------------------------------------------------------------
144 * hypre_StructMatrixInitializeShell
145 *--------------------------------------------------------------------------*/
146
147int
148hypre_StructMatrixInitializeShell( hypre_StructMatrix *matrix )
149{
150 int ierr = 0;
151
152 hypre_StructGrid *grid;
153
154 hypre_StructStencil *user_stencil;
155 hypre_StructStencil *stencil;
156 hypre_Index *stencil_shape;
157 int stencil_size;
158 int num_values;
159 int *symm_elements;
160 int constant_coefficient;
161
162 int *num_ghost;
163 int extra_ghost[] = {0, 0, 0, 0, 0, 0};
164
165 hypre_BoxArray *data_space;
166 hypre_BoxArray *boxes;
167 hypre_Box *box;
168 hypre_Box *data_box;
169
170 int **data_indices;
171 int data_size;
172 int data_box_volume;
173
174 int i, j, d;
175
176 grid = hypre_StructMatrixGrid(matrix);
177
178 /*-----------------------------------------------------------------------
179 * Set up stencil and num_values:
180 *
181 * If the matrix is symmetric, then the stencil is a "symmetrized"
182 * version of the user's stencil. If the matrix is not symmetric,
183 * then the stencil is the same as the user's stencil.
184 *
185 * The `symm_elements' array is used to determine what data is
186 * explicitely stored (symm_elements[i] < 0) and what data does is
187 * not explicitely stored (symm_elements[i] >= 0), but is instead
188 * stored as the transpose coefficient at a neighboring grid point.
189 *-----------------------------------------------------------------------*/
190
191 if (hypre_StructMatrixStencil(matrix) == NULL)
192 {
193 user_stencil = hypre_StructMatrixUserStencil(matrix);
194
195 if (hypre_StructMatrixSymmetric(matrix))
196 {
197 /* store only symmetric stencil entry data */
198 hypre_StructStencilSymmetrize(user_stencil, &stencil, &symm_elements);
199 num_values = ( hypre_StructStencilSize(stencil) + 1 ) / 2;
200 }
201 else
202 {
203 /* store all stencil entry data */
204 stencil = hypre_StructStencilRef(user_stencil);
205 num_values = hypre_StructStencilSize(stencil);
206 symm_elements = hypre_TAlloc(int, num_values);
207 for (i = 0; i < num_values; i++)
208 {
209 symm_elements[i] = -1;
210 }
211 }
212
213 hypre_StructMatrixStencil(matrix) = stencil;
214 hypre_StructMatrixSymmElements(matrix) = symm_elements;
215 hypre_StructMatrixNumValues(matrix) = num_values;
216 }
217
218 /*-----------------------------------------------------------------------
219 * Set ghost-layer size for symmetric storage
220 * - All stencil coeffs are to be available at each point in the
221 * grid, as well as in the user-specified ghost layer.
222 *-----------------------------------------------------------------------*/
223
224 num_ghost = hypre_StructMatrixNumGhost(matrix);
225 stencil = hypre_StructMatrixStencil(matrix);
226 stencil_shape = hypre_StructStencilShape(stencil);
227 stencil_size = hypre_StructStencilSize(stencil);
228 symm_elements = hypre_StructMatrixSymmElements(matrix);
229
230 for (i = 0; i < stencil_size; i++)
231 {
232 if (symm_elements[i] >= 0)
233 {
234 for (d = 0; d < 3; d++)
235 {
236 extra_ghost[2*d] =
237 hypre_max(extra_ghost[2*d], -hypre_IndexD(stencil_shape[i], d));
238 extra_ghost[2*d + 1] =
239 hypre_max(extra_ghost[2*d + 1], hypre_IndexD(stencil_shape[i], d));
240 }
241 }
242 }
243
244 for (d = 0; d < 3; d++)
245 {
246 num_ghost[2*d] += extra_ghost[2*d];
247 num_ghost[2*d + 1] += extra_ghost[2*d + 1];
248 }
249
250 /*-----------------------------------------------------------------------
251 * Set up data_space
252 *-----------------------------------------------------------------------*/
253
254 if (hypre_StructMatrixDataSpace(matrix) == NULL)
255 {
256 boxes = hypre_StructGridBoxes(grid);
257 data_space = hypre_BoxArrayCreate(hypre_BoxArraySize(boxes));
258
259 hypre_ForBoxI(i, boxes)
260 {
261 box = hypre_BoxArrayBox(boxes, i);
262 data_box = hypre_BoxArrayBox(data_space, i);
263
264 hypre_CopyBox(box, data_box);
265 for (d = 0; d < 3; d++)
266 {
267 hypre_BoxIMinD(data_box, d) -= num_ghost[2*d];
268 hypre_BoxIMaxD(data_box, d) += num_ghost[2*d + 1];
269 }
270 }
271
272 hypre_StructMatrixDataSpace(matrix) = data_space;
273 }
274
275 /*-----------------------------------------------------------------------
276 * Set up data_indices array and data-size
277 *-----------------------------------------------------------------------*/
278
279 if (hypre_StructMatrixDataIndices(matrix) == NULL)
280 {
281 data_space = hypre_StructMatrixDataSpace(matrix);
282 data_indices = hypre_CTAlloc(int *, hypre_BoxArraySize(data_space));
283 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
284
285 data_size = 0;
286 if ( constant_coefficient==0 )
287 {
288 hypre_ForBoxI(i, data_space)
289 {
290 data_box = hypre_BoxArrayBox(data_space, i);
291 data_box_volume = hypre_BoxVolume(data_box);
292
293 data_indices[i] = hypre_CTAlloc(int, stencil_size);
294
295 /* set pointers for "stored" coefficients */
296 for (j = 0; j < stencil_size; j++)
297 {
298 if (symm_elements[j] < 0)
299 {
300 data_indices[i][j] = data_size;
301 data_size += data_box_volume;
302 }
303 }
304
305 /* set pointers for "symmetric" coefficients */
306 for (j = 0; j < stencil_size; j++)
307 {
308 if (symm_elements[j] >= 0)
309 {
310 data_indices[i][j] = data_indices[i][symm_elements[j]] +
311 hypre_BoxOffsetDistance(data_box, stencil_shape[j]);
312 }
313 }
314 }
315 }
316 else if ( constant_coefficient==1 )
317 {
318
319 hypre_ForBoxI(i, data_space)
320 {
321 data_box = hypre_BoxArrayBox(data_space, i);
322 data_box_volume = hypre_BoxVolume(data_box);
323
324 data_indices[i] = hypre_CTAlloc(int, stencil_size);
325
326 /* set pointers for "stored" coefficients */
327 for (j = 0; j < stencil_size; j++)
328 {
329 if (symm_elements[j] < 0)
330 {
331 data_indices[i][j] = data_size;
332 ++data_size;
333 }
334 }
335
336 /* set pointers for "symmetric" coefficients */
337 for (j = 0; j < stencil_size; j++)
338 {
339 if (symm_elements[j] >= 0)
340 {
341 data_indices[i][j] = data_indices[i][symm_elements[j]];
342 }
343 }
344 }
345 }
346 else
347 {
348 hypre_assert( constant_coefficient == 2 );
349 data_size += stencil_size; /* all constant coefficients at the beginning */
350 /* ... this allocates a little more space than is absolutely necessary */
351 hypre_ForBoxI(i, data_space)
352 {
353 data_box = hypre_BoxArrayBox(data_space, i);
354 data_box_volume = hypre_BoxVolume(data_box);
355
356 data_indices[i] = hypre_CTAlloc(int, stencil_size);
357
358 /* set pointers for "stored" coefficients */
359 for (j = 0; j < stencil_size; j++)
360 {
361 if (symm_elements[j] < 0)
362 {
363 if (
364 hypre_IndexX(stencil_shape[j])==0 &&
365 hypre_IndexY(stencil_shape[j])==0 &&
366 hypre_IndexZ(stencil_shape[j])==0 ) /* diagonal, variable coefficient */
367 {
368 data_indices[i][j] = data_size;
369 data_size += data_box_volume;
370 }
371 else /* off-diagonal, constant coefficient */
372 {
373 data_indices[i][j] = j;
374 }
375 }
376 }
377
378 /* set pointers for "symmetric" coefficients */
379 for (j = 0; j < stencil_size; j++)
380 {
381 if (symm_elements[j] >= 0)
382 {
383 if (
384 hypre_IndexX(stencil_shape[j])==0 &&
385 hypre_IndexY(stencil_shape[j])==0 &&
386 hypre_IndexZ(stencil_shape[j])==0 ) /* diagonal, variable coefficient */
387 {
388 data_indices[i][j] = data_indices[i][symm_elements[j]] +
389 hypre_BoxOffsetDistance(data_box, stencil_shape[j]);
390 }
391 else /* off-diagonal, constant coefficient */
392 {
393 data_indices[i][j] = data_indices[i][symm_elements[j]];
394 }
395 }
396 }
397 }
398 }
399
400 hypre_StructMatrixDataIndices(matrix) = data_indices;
401 hypre_StructMatrixDataSize(matrix) = data_size;
402 }
403
404 /*-----------------------------------------------------------------------
405 * Set total number of nonzero coefficients
406 * For constant coefficients, this is unrelated to the amount of data
407 * actually stored.
408 *-----------------------------------------------------------------------*/
409
410 hypre_StructMatrixGlobalSize(matrix) =
411 hypre_StructGridGlobalSize(grid) * stencil_size;
412
413 /*-----------------------------------------------------------------------
414 * Return
415 *-----------------------------------------------------------------------*/
416
417 return ierr;
418}
419
420/*--------------------------------------------------------------------------
421 * hypre_StructMatrixInitializeData
422 *--------------------------------------------------------------------------*/
423
424int
425hypre_StructMatrixInitializeData( hypre_StructMatrix *matrix,
426 double *data )
427{
428 int ierr = 0;
429
430 hypre_BoxArray *data_boxes;
431 hypre_Box *data_box;
432 hypre_Index loop_size;
433 hypre_Index index;
434 hypre_IndexRef start;
435 hypre_Index stride;
436 double *datap;
437 int datai;
438 int constant_coefficient;
439 int i;
440 int loopi, loopj, loopk;
441
442 hypre_StructMatrixData(matrix) = data;
443 hypre_StructMatrixDataAlloced(matrix) = 0;
444 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
445
446 /*-------------------------------------------------
447 * If the matrix has a diagonal, set these entries
448 * to 1 inside the grid_boxes. Ghostvalues will
449 * be set to 1 in the assembly. This reduces the
450 * complexity of many computations by eliminating
451 * divide-by-zero in the ghost region.
452 *-------------------------------------------------*/
453
454 hypre_SetIndex(index, 0, 0, 0);
455 hypre_SetIndex(stride, 1, 1, 1);
456
457 data_boxes = hypre_StructMatrixDataSpace(matrix);
458 hypre_ForBoxI(i, data_boxes)
459 {
460 datap = hypre_StructMatrixExtractPointerByIndex(matrix, i, index);
461
462 if (datap)
463 {
464 data_box = hypre_BoxArrayBox(data_boxes, i);
465 start = hypre_BoxIMin(data_box);
466
467 if ( constant_coefficient==1 )
468 {
469 datai = hypre_CCBoxIndexRank(data_box,start);
470 datap[datai] = 1.0;
471 }
472 else
473 /* either fully variable coefficient matrix, constant_coefficient=0,
474 or variable diagonal (otherwise constant) coefficient matrix,
475 constant_coefficient=2 */
476 {
477 hypre_BoxGetSize(data_box, loop_size);
478
479 hypre_BoxLoop1Begin(loop_size,
480 data_box, start, stride, datai);
481
482 hypre_BoxLoop1For(loopi, loopj, loopk, datai)
483 {
484 datap[datai] = 1.0;
485 }
486 hypre_BoxLoop1End(datai);
487 }
488 }
489 }
490
491 return ierr;
492}
493
494/*--------------------------------------------------------------------------
495 * hypre_StructMatrixInitialize
496 *--------------------------------------------------------------------------*/
497int
498hypre_StructMatrixInitialize( hypre_StructMatrix *matrix )
499{
500 int ierr = 0;
501
502 double *data;
503
504 ierr = hypre_StructMatrixInitializeShell(matrix);
505
506 data = hypre_StructMatrixData(matrix);
507 data = hypre_SharedCTAlloc(double, hypre_StructMatrixDataSize(matrix));
508 hypre_StructMatrixInitializeData(matrix, data);
509 hypre_StructMatrixDataAlloced(matrix) = 1;
510
511 return ierr;
512}
513
514/*--------------------------------------------------------------------------
515 * (action > 0): add-to values
516 * (action = 0): set values
517 * (action < 0): get values
518 * should not be called to set a constant-coefficient part of the matrix,
519 * call hypre_StructMatrixSetConstantValues instead
520 *--------------------------------------------------------------------------*/
521
522int
523hypre_StructMatrixSetValues( hypre_StructMatrix *matrix,
524 hypre_Index grid_index,
525 int num_stencil_indices,
526 int *stencil_indices,
527 double *values,
528 int action )
529{
530 int ierr = 0;
531 MPI_Comm comm= hypre_StructMatrixComm(matrix);
532
533 hypre_BoxArray *boxes;
534 hypre_Box *box;
535 hypre_Index center_index;
536 hypre_StructStencil *stencil;
537 int center_rank;
538 int constant_coefficient;
539
540 double *matp;
541
542 int i, s, found;
543 int true = 1;
544 int false= 0;
545
546 int nprocs;
547
548 MPI_Comm_size(comm, &nprocs );
549
550 boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
551 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
552
553 found= true; /* index found will be set to false later on. This
554 eliminates the constant_coefficient= 1 case correctly. */
555 hypre_ForBoxI(i, boxes)
556 {
557 box = hypre_BoxArrayBox(boxes, i);
558
559 if ( constant_coefficient==1 )
560 {
561 ++ierr; /* call SetConstantValues instead */
562 if (action > 0)
563 {
564 for (s = 0; s < num_stencil_indices; s++)
565 {
566 matp = hypre_StructMatrixBoxData(matrix, i,
567 stencil_indices[s]);
568 *matp += values[s];
569 }
570 }
571 else if (action > -1)
572 {
573 for (s = 0; s < num_stencil_indices; s++)
574 {
575 matp = hypre_StructMatrixBoxData(matrix, i,
576 stencil_indices[s]);
577 *matp = values[s];
578 }
579 }
580 else /* action < 0 */
581 {
582 for (s = 0; s < num_stencil_indices; s++)
583 {
584 matp = hypre_StructMatrixBoxData(matrix, i,
585 stencil_indices[s]);
586 values[s] = *matp;
587 }
588 }
589 }
590 else if ( constant_coefficient==2 )
591 {
592 hypre_SetIndex(center_index, 0, 0, 0);
593 stencil = hypre_StructMatrixStencil(matrix);
594 center_rank = hypre_StructStencilElementRank( stencil, center_index );
595
596 found= false;
597 if ( action > 0 )
598 {
599 for (s = 0; s < num_stencil_indices; s++)
600 {
601 if ( stencil_indices[s] == center_rank )
602 { /* center (diagonal), like constant_coefficient==0 */
603 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
604 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
605 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
606 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
607 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
608 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
609 {
610 matp = hypre_StructMatrixBoxDataValue(matrix, i,
611 stencil_indices[s],
612 grid_index);
613 *matp += values[s];
614 found= true;
615 }
616 }
617 else
618 { /* non-center, like constant_coefficient==1 */
619 ++ierr; /* should have called SetConstantValues */
620 matp = hypre_StructMatrixBoxData(matrix, i,
621 stencil_indices[s]);
622 *matp += values[s];
623 found= true;
624 }
625 }
626 }
627 else if ( action > -1 )
628 {
629 for (s = 0; s < num_stencil_indices; s++)
630 {
631 if ( stencil_indices[s] == center_rank )
632 { /* center (diagonal), like constant_coefficient==0 */
633 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
634 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
635 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
636 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
637 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
638 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
639 {
640 matp = hypre_StructMatrixBoxDataValue(matrix, i,
641 stencil_indices[s],
642 grid_index);
643 *matp = values[s];
644 found= true;
645 }
646 }
647 else
648 { /* non-center, like constant_coefficient==1 */
649 ++ierr; /* should have called SetConstantValues */
650 matp = hypre_StructMatrixBoxData(matrix, i,
651 stencil_indices[s]);
652 *matp += values[s];
653 found= true;
654 }
655 }
656 }
657 else /* action<0 */
658 {
659 found= true; /* no need to set-off proc for get values */
660 for (s = 0; s < num_stencil_indices; s++)
661 {
662 if ( stencil_indices[s] == center_rank )
663 { /* center (diagonal), like constant_coefficient==0 */
664 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
665 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
666 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
667 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
668 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
669 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
670 {
671 matp = hypre_StructMatrixBoxDataValue(matrix, i,
672 stencil_indices[s],
673 grid_index);
674 *matp += values[s];
675 }
676 }
677 else
678 { /* non-center, like constant_coefficient==1 */
679 ++ierr; /* should have called SetConstantValues */
680 matp = hypre_StructMatrixBoxData(matrix, i,
681 stencil_indices[s]);
682 values[s] = *matp;
683 }
684 }
685 }
686 }
687 else
688 /* variable coefficient, constant_coefficient=0 */
689 {
690 found= false;
691 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
692 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
693 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
694 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
695 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
696 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
697 {
698 if (action > 0)
699 {
700 for (s = 0; s < num_stencil_indices; s++)
701 {
702 matp = hypre_StructMatrixBoxDataValue(matrix, i,
703 stencil_indices[s],
704 grid_index);
705 *matp += values[s];
706 }
707 }
708 else if (action > -1)
709 {
710 for (s = 0; s < num_stencil_indices; s++)
711 {
712 matp = hypre_StructMatrixBoxDataValue(matrix, i,
713 stencil_indices[s],
714 grid_index);
715 *matp = values[s];
716 }
717 }
718 else
719 {
720 for (s = 0; s < num_stencil_indices; s++)
721 {
722 matp = hypre_StructMatrixBoxDataValue(matrix, i,
723 stencil_indices[s],
724 grid_index);
725 values[s] = *matp;
726 }
727 }
728
729 found= true;
730 }
731 }
732 }
733
734 /* if the index was not found on this proc and the user wants to ADD
735 an off_proc value, see if it is along the ghostlayer. */
736 if ((!found) && (action > 0) && (nprocs > 1))
737 {
738 hypre_Box *orig_box;
739
740 int *add_num_ghost= hypre_StructMatrixAddNumGhost(matrix);
741 int j;
742
743 hypre_ForBoxI(i, boxes)
744 {
745 orig_box = hypre_BoxArrayBox(boxes, i);
746 box = hypre_BoxDuplicate(orig_box);
747 for (j= 0; j< 3; j++)
748 {
749 hypre_BoxIMin(box)[j]-= add_num_ghost[2*j];
750 hypre_BoxIMax(box)[j]+= add_num_ghost[2*j+1];
751 }
752
753 if ( constant_coefficient==2 ) /* only center and action > 0 */
754 {
755 hypre_SetIndex(center_index, 0, 0, 0);
756 stencil = hypre_StructMatrixStencil(matrix);
757 center_rank = hypre_StructStencilElementRank( stencil, center_index );
758
759 /* center (diagonal), like constant_coefficient==0 */
760 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
761 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
762 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
763 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
764 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
765 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
766 {
767 matp = hypre_StructMatrixBoxDataValue(matrix, i,
768 stencil_indices[center_rank],
769 grid_index);
770 *matp += values[s];
771 found= true;
772 }
773 }
774
775 else
776 /* variable coefficient, constant_coefficient=0 */
777 {
778 if ((hypre_IndexX(grid_index) >= hypre_BoxIMinX(box)) &&
779 (hypre_IndexX(grid_index) <= hypre_BoxIMaxX(box)) &&
780 (hypre_IndexY(grid_index) >= hypre_BoxIMinY(box)) &&
781 (hypre_IndexY(grid_index) <= hypre_BoxIMaxY(box)) &&
782 (hypre_IndexZ(grid_index) >= hypre_BoxIMinZ(box)) &&
783 (hypre_IndexZ(grid_index) <= hypre_BoxIMaxZ(box)) )
784 {
785 for (s = 0; s < num_stencil_indices; s++)
786 {
787 matp = hypre_StructMatrixBoxDataValue(matrix, i,
788 stencil_indices[s],
789 grid_index);
790 *matp += values[s];
791 }
792 found= true;
793 }
794 }
795 if (found) break;
796 }
797
798 /* set OffProcAdd for communication */
799 if (found)
800 {
801 hypre_StructMatrixOffProcAdd(matrix)= 1;
802 }
803 else
804 {
805 printf("not found- grid_index off the extended matrix grid\n");
806 }
807 } /* if ((!found) && (action > 0)) */
808
809 return(ierr);
810}
811
812/*--------------------------------------------------------------------------
813 * (action > 0): add-to values
814 * (action = 0): set values
815 * (action =-1): get values and zero out
816 * (action <-1): get values
817 * should not be called to set a constant-coefficient part of the matrix,
818 * call hypre_StructMatrixSetConstantValues instead
819 *--------------------------------------------------------------------------*/
820
821int
822hypre_StructMatrixSetBoxValues( hypre_StructMatrix *matrix,
823 hypre_Box *value_box,
824 int num_stencil_indices,
825 int *stencil_indices,
826 double *values,
827 int action )
828{
829 int ierr = 0;
830 MPI_Comm comm = hypre_StructMatrixComm(matrix);
831 int *add_num_ghost= hypre_StructVectorAddNumGhost(matrix);
832
833 hypre_BoxArray *grid_boxes;
834 hypre_Box *grid_box;
835 hypre_BoxArray *box_array1, *box_array2, *tmp_box_array;
836 hypre_BoxArray *value_boxarray;
837 hypre_BoxArrayArray *box_aarray;
838 hypre_Box *box, *tmp_box, *orig_box, *vbox;
839 hypre_Index center_index;
840 hypre_StructStencil *stencil;
841 int center_rank;
842
843 hypre_BoxArray *data_space;
844 hypre_Box *data_box;
845 hypre_IndexRef data_start;
846 hypre_Index data_stride;
847 int datai;
848 double *datap;
849 int constant_coefficient;
850
851 hypre_Box *dval_box;
852 hypre_Index dval_start;
853 hypre_Index dval_stride;
854 int dvali;
855
856 hypre_Index loop_size;
857
858 int i, j, k, s, vol_vbox, vol_iboxes, vol_offproc;
859 int loopi, loopj, loopk;
860
861 int nprocs;
862
863 MPI_Comm_size(comm, &nprocs );
864
865 /*-----------------------------------------------------------------------
866 * Set up `box_array' by intersecting `box' with the grid boxes
867 *-----------------------------------------------------------------------*/
868
869 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
870
871 /* Find the intersecting boxes of the grid with value_box. Record the
872 volumes of the intersections for possible off_proc settings. */
873 vol_vbox = hypre_BoxVolume(value_box);
874 vol_iboxes= 0;
875
876 grid_boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
877 box_array1 = hypre_BoxArrayCreate(hypre_BoxArraySize(grid_boxes));
878 box = hypre_BoxCreate();
879 hypre_ForBoxI(i, grid_boxes)
880 {
881 grid_box = hypre_BoxArrayBox(grid_boxes, i);
882 hypre_IntersectBoxes(value_box, grid_box, box);
883 hypre_CopyBox(box, hypre_BoxArrayBox(box_array1, i));
884 vol_iboxes+= hypre_BoxVolume(box);
885 }
886
887 vol_offproc= 0;
888 if ((vol_vbox > vol_iboxes) && (action > 0) && (nprocs > 1)) /* only for addto values */
889 {
890 box_aarray= hypre_BoxArrayArrayCreate(hypre_BoxArraySize(grid_boxes));
891
892 /* to prevent overlapping intersected boxes, we subtract the intersected
893 boxes from value_box. This requires a box_array structure. */
894 value_boxarray= hypre_BoxArrayCreate(0);
895 hypre_AppendBox(value_box, value_boxarray);
896
897 hypre_ForBoxI(i, grid_boxes)
898 {
899 tmp_box_array= hypre_BoxArrayCreate(0);
900
901 /* get ghostlayer boxes */
902 orig_box= hypre_BoxArrayBox(grid_boxes, i);
903 tmp_box = hypre_BoxDuplicate(orig_box );
904 for (j= 0; j< 3; j++)
905 {
906 hypre_BoxIMin(tmp_box)[j]-= add_num_ghost[2*j];
907 hypre_BoxIMax(tmp_box)[j]+= add_num_ghost[2*j+1];
908 }
909
910 /* get ghostlayer boxes */
911 hypre_SubtractBoxes(tmp_box, orig_box, tmp_box_array);
912 hypre_BoxDestroy(tmp_box);
913
914 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
915 /* intersect the value_box and the ghostlayer boxes */
916 hypre_ForBoxI(j, tmp_box_array)
917 {
918 tmp_box= hypre_BoxArrayBox(tmp_box_array, j);
919 hypre_ForBoxI(k, value_boxarray)
920 {
921 vbox= hypre_BoxArrayBox(value_boxarray, k);
922 hypre_IntersectBoxes(vbox, tmp_box, box);
923 hypre_AppendBox(box, box_array2);
924
925 vol_offproc+= hypre_BoxVolume(box);
926 }
927 }
928
929 /* eliminate intersected boxes so that we do not get overlapping */
930 hypre_SubtractBoxArrays(value_boxarray, box_array2, tmp_box_array);
931 hypre_BoxArrayDestroy(tmp_box_array);
932 } /* hypre_ForBoxI(i, grid_boxes) */
933
934 /* if vol_offproc= 0, trying to set values away from ghostlayer */
935 if (!vol_offproc)
936 {
937 hypre_BoxArrayArrayDestroy(box_aarray);
938 }
939 else
940 {
941 /* set OffProcAdd for communication, i.e., off-proc values
942 must be communicated */
943 hypre_StructMatrixOffProcAdd(matrix)= 1;
944 }
945 hypre_BoxArrayDestroy(value_boxarray);
946
947 }
948 hypre_BoxDestroy(box);
949
950 /*-----------------------------------------------------------------------
951 * Set the matrix coefficients
952 *-----------------------------------------------------------------------*/
953
954 if (box_array1)
955 {
956 data_space = hypre_StructMatrixDataSpace(matrix);
957 hypre_SetIndex(data_stride, 1, 1, 1);
958
959 dval_box = hypre_BoxDuplicate(value_box);
960 hypre_BoxIMinD(dval_box, 0) *= num_stencil_indices;
961 hypre_BoxIMaxD(dval_box, 0) *= num_stencil_indices;
962 hypre_BoxIMaxD(dval_box, 0) += num_stencil_indices - 1;
963 hypre_SetIndex(dval_stride, num_stencil_indices, 1, 1);
964
965 hypre_ForBoxI(i, box_array1)
966 {
967 box = hypre_BoxArrayBox(box_array1, i);
968 data_box = hypre_BoxArrayBox(data_space, i);
969
970 /* if there was an intersection */
971 if (box)
972 {
973 data_start = hypre_BoxIMin(box);
974 hypre_CopyIndex(data_start, dval_start);
975 hypre_IndexD(dval_start, 0) *= num_stencil_indices;
976
977 if ( constant_coefficient==2 )
978 {
979 hypre_SetIndex(center_index, 0, 0, 0);
980 stencil = hypre_StructMatrixStencil(matrix);
981 center_rank = hypre_StructStencilElementRank( stencil, center_index );
982 }
983
984 for (s = 0; s < num_stencil_indices; s++)
985 {
986 datap = hypre_StructMatrixBoxData(matrix, i,
987 stencil_indices[s]);
988
989 if ( constant_coefficient==1 ||
990 ( constant_coefficient==2 && stencil_indices[s]!=center_rank ))
991 /* datap has only one data point for a given i and s */
992 {
993 ++ierr; /* should have called SetConstantValues */
994 hypre_BoxGetSize(box, loop_size);
995
996 if (action > 0)
997 {
998 datai = hypre_CCBoxIndexRank(data_box,data_start);
999 dvali = hypre_BoxIndexRank(dval_box,dval_start);
1000 datap[datai] += values[dvali];
1001 }
1002 else if (action > -1)
1003 {
1004 datai = hypre_CCBoxIndexRank(data_box,data_start);
1005 dvali = hypre_BoxIndexRank(dval_box,dval_start);
1006 datap[datai] = values[dvali];
1007 }
1008 else
1009 {
1010 datai = hypre_CCBoxIndexRank(data_box,data_start);
1011 dvali = hypre_BoxIndexRank(dval_box,dval_start);
1012 values[dvali] = datap[datai];
1013 if (action == -1)
1014 {
1015 datap[datai] = 0;
1016 }
1017 }
1018
1019 }
1020 else /* variable coefficient: constant_coefficient==0
1021 or diagonal with constant_coefficient==2 */
1022 {
1023 hypre_BoxGetSize(box, loop_size);
1024
1025 if (action > 0)
1026 {
1027 hypre_BoxLoop2Begin(loop_size,
1028 data_box,data_start,data_stride,datai,
1029 dval_box,dval_start,dval_stride,dvali);
1030
1031 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
1032 {
1033 datap[datai] += values[dvali];
1034 }
1035 hypre_BoxLoop2End(datai, dvali);
1036 }
1037 else if (action > -1)
1038 {
1039 hypre_BoxLoop2Begin(loop_size,
1040 data_box,data_start,data_stride,datai,
1041 dval_box,dval_start,dval_stride,dvali);
1042
1043 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
1044 {
1045 datap[datai] = values[dvali];
1046 }
1047 hypre_BoxLoop2End(datai, dvali);
1048 }
1049 else
1050 {
1051 hypre_BoxLoop2Begin(loop_size,
1052 data_box,data_start,data_stride,datai,
1053 dval_box,dval_start,dval_stride,dvali);
1054
1055 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
1056 {
1057 values[dvali] = datap[datai];
1058 if (action == -1)
1059 {
1060 datap[datai] = 0;
1061 }
1062 }
1063 hypre_BoxLoop2End(datai, dvali);
1064 }
1065 }
1066
1067 hypre_IndexD(dval_start, 0) ++;
1068 }
1069 }
1070 }
1071
1072 hypre_BoxDestroy(dval_box);
1073 }
1074
1075 hypre_BoxArrayDestroy(box_array1);
1076
1077 if (vol_offproc) /* only adding values, action > 0 */
1078 {
1079 data_space= hypre_StructMatrixDataSpace(matrix);
1080 hypre_SetIndex(data_stride, 1, 1, 1);
1081
1082 dval_box = hypre_BoxDuplicate(value_box);
1083 hypre_BoxIMinD(dval_box, 0) *= num_stencil_indices;
1084 hypre_BoxIMaxD(dval_box, 0) *= num_stencil_indices;
1085 hypre_BoxIMaxD(dval_box, 0) += num_stencil_indices - 1;
1086 hypre_SetIndex(dval_stride, num_stencil_indices, 1, 1);
1087
1088 hypre_ForBoxI(i, data_space)
1089 {
1090 data_box = hypre_BoxArrayBox(data_space, i);
1091 box_array2= hypre_BoxArrayArrayBoxArray(box_aarray, i);
1092
1093 hypre_ForBoxI(j, box_array2)
1094 {
1095 box= hypre_BoxArrayBox(box_array2, j);
1096
1097 /* if there was an intersection */
1098 if (box)
1099 {
1100 data_start = hypre_BoxIMin(box);
1101 hypre_CopyIndex(data_start, dval_start);
1102 hypre_IndexD(dval_start, 0) *= num_stencil_indices;
1103
1104 if ( constant_coefficient==2 )
1105 {
1106 hypre_SetIndex(center_index, 0, 0, 0);
1107 stencil= hypre_StructMatrixStencil(matrix);
1108 center_rank=
1109 hypre_StructStencilElementRank(stencil, center_index);
1110 }
1111
1112 for (s = 0; s < num_stencil_indices; s++)
1113 {
1114 datap = hypre_StructMatrixBoxData(matrix, i,
1115 stencil_indices[s]);
1116
1117 /* variable coefficient: constant_coefficient==0
1118 or diagonal with constant_coefficient==2 */
1119 if ( constant_coefficient==0 ||
1120 ( constant_coefficient==2 && stencil_indices[s]==center_rank ))
1121 {
1122 hypre_BoxGetSize(box, loop_size);
1123
1124 hypre_BoxLoop2Begin(loop_size,
1125 data_box,data_start,data_stride,datai,
1126 dval_box,dval_start,dval_stride,dvali);
1127
1128 hypre_BoxLoop2For(loopi, loopj, loopk, datai, dvali)
1129 {
1130 datap[datai] += values[dvali];
1131 }
1132 hypre_BoxLoop2End(datai, dvali);
1133
1134 } /* else variable coefficient */
1135 hypre_IndexD(dval_start, 0)++;
1136 } /* for (s = 0; s < num_stencil_indices; s++) */
1137 } /* if (box) */
1138 } /* hypre_ForBoxI(j, box_array2) */
1139 } /* hypre_ForBoxI(i, data_space) */
1140
1141 hypre_BoxDestroy(dval_box);
1142 hypre_BoxArrayArrayDestroy(box_aarray);
1143 }
1144
1145 return(ierr);
1146}
1147
1148
1149/*--------------------------------------------------------------------------
1150 * (action > 0): add-to values
1151 * (action = 0): set values
1152 * (action =-1): get values and zero out (not implemented, just gets values)
1153 * (action <-1): get values
1154 * should be called to set a constant-coefficient part of the matrix
1155 *--------------------------------------------------------------------------*/
1156
1157int
1158hypre_StructMatrixSetConstantValues( hypre_StructMatrix *matrix,
1159 int num_stencil_indices,
1160 int *stencil_indices,
1161 double *values,
1162 int action )
1163{
1164 int ierr = 0;
1165 hypre_BoxArray *boxes;
1166 hypre_Box *box;
1167 hypre_Index center_index;
1168 hypre_StructStencil *stencil;
1169 int center_rank;
1170 int constant_coefficient;
1171
1172 double *matp;
1173
1174 int i, s;
1175
1176 boxes = hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
1177 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
1178
1179 if ( constant_coefficient==1 )
1180 {
1181 hypre_ForBoxI(i, boxes)
1182 {
1183 box = hypre_BoxArrayBox(boxes, i);
1184 if (action > 0)
1185 {
1186 for (s = 0; s < num_stencil_indices; s++)
1187 {
1188 matp = hypre_StructMatrixBoxData(matrix, i,
1189 stencil_indices[s]);
1190 *matp += values[s];
1191 }
1192 }
1193 else if (action > -1)
1194 {
1195 for (s = 0; s < num_stencil_indices; s++)
1196 {
1197 matp = hypre_StructMatrixBoxData(matrix, i,
1198 stencil_indices[s]);
1199 *matp = values[s];
1200 }
1201 }
1202 else /* action < 0 */
1203 {
1204 for (s = 0; s < num_stencil_indices; s++)
1205 {
1206 matp = hypre_StructMatrixBoxData(matrix, i,
1207 stencil_indices[s]);
1208 values[s] = *matp;
1209 }
1210 }
1211 }
1212 }
1213 else if ( constant_coefficient==2 )
1214 {
1215 hypre_SetIndex(center_index, 0, 0, 0);
1216 stencil = hypre_StructMatrixStencil(matrix);
1217 center_rank = hypre_StructStencilElementRank( stencil, center_index );
1218 if ( action > 0 )
1219 {
1220 for (s = 0; s < num_stencil_indices; s++)
1221 {
1222 if ( stencil_indices[s] == center_rank )
1223 { /* center (diagonal), like constant_coefficient==0
1224 We consider it an error, but do the best we can. */
1225 ++ierr;
1226 hypre_ForBoxI(i, boxes)
1227 {
1228 box = hypre_BoxArrayBox(boxes, i);
1229 ierr += hypre_StructMatrixSetBoxValues(
1230 matrix, box,
1231 num_stencil_indices, stencil_indices,
1232 values, action );
1233 }
1234 }
1235 else
1236 { /* non-center, like constant_coefficient==1 */
1237 matp = hypre_StructMatrixBoxData(matrix, 0,
1238 stencil_indices[s]);
1239 *matp += values[s];
1240 }
1241 }
1242 }
1243 else if ( action > -1 )
1244 {
1245 for (s = 0; s < num_stencil_indices; s++)
1246 {
1247 if ( stencil_indices[s] == center_rank )
1248 { /* center (diagonal), like constant_coefficient==0
1249 We consider it an error, but do the best we can. */
1250 ++ierr;
1251 hypre_ForBoxI(i, boxes)
1252 {
1253 box = hypre_BoxArrayBox(boxes, i);
1254 ierr += hypre_StructMatrixSetBoxValues(
1255 matrix, box,
1256 num_stencil_indices, stencil_indices,
1257 values, action );
1258 }
1259 }
1260 else
1261 { /* non-center, like constant_coefficient==1 */
1262 matp = hypre_StructMatrixBoxData(matrix, 0,
1263 stencil_indices[s]);
1264 *matp += values[s];
1265 }
1266 }
1267 }
1268 else /* action<0 */
1269 {
1270 for (s = 0; s < num_stencil_indices; s++)
1271 {
1272 if ( stencil_indices[s] == center_rank )
1273 { /* center (diagonal), like constant_coefficient==0
1274 We consider it an error, but do the best we can. */
1275 ++ierr;
1276 hypre_ForBoxI(i, boxes)
1277 {
1278 box = hypre_BoxArrayBox(boxes, i);
1279 ierr += hypre_StructMatrixSetBoxValues(
1280 matrix, box,
1281 num_stencil_indices, stencil_indices,
1282 values, -1);
1283 }
1284 }
1285 else
1286 { /* non-center, like constant_coefficient==1 */
1287 matp = hypre_StructMatrixBoxData(matrix, 0,
1288 stencil_indices[s]);
1289 values[s] = *matp;
1290 }
1291 }
1292 }
1293 }
1294 else /* constant_coefficient==0 */
1295 {
1296 /* We consider this an error, but do the best we can. */
1297 ++ierr;
1298 hypre_ForBoxI(i, boxes)
1299 {
1300 box = hypre_BoxArrayBox(boxes, i);
1301 ierr += hypre_StructMatrixSetBoxValues(
1302 matrix, box,
1303 num_stencil_indices, stencil_indices,
1304 values, action );
1305 }
1306 }
1307 return ierr;
1308}
1309
1310/*--------------------------------------------------------------------------
1311 * hypre_StructMatrixAssemble:
1312 * Before assembling the matrix, all matrix values added from off_procs
1313 * are communicated to the correct proc. However, note that since
1314 * the comm_pkg is created from an "inverted" comm_info derived from the
1315 * matrix, not all the communicated data is valid (i.e., we did not mark
1316 * which values were actually set off_proc). Some of them are invalid
1317 * communicated data. These data values are zero or one (centre). The
1318 * centre will be checked against 1 (not the best solution).
1319 *--------------------------------------------------------------------------*/
1320
1321int
1322hypre_StructMatrixAssemble( hypre_StructMatrix *matrix )
1323{
1324 int ierr = 0;
1325
1326 int *num_ghost = hypre_StructMatrixNumGhost(matrix);
1327
1328 int comm_num_values, mat_num_values, constant_coefficient;
1329 int stencil_size;
1330 hypre_StructStencil *stencil;
1331
1332 hypre_CommInfo *comm_info;
1333 hypre_CommPkg *comm_pkg;
1334
1335 hypre_CommHandle *comm_handle;
1336 int data_initial_offset = 0;
1337 double *matrix_data = hypre_StructMatrixData(matrix);
1338 double *matrix_data_comm = matrix_data;
1339
1340 int sum_OffProcAdd;
1341 int OffProcAdd= hypre_StructMatrixOffProcAdd(matrix);
1342
1343 hypre_Box *data_box;
1344 hypre_Box *box;
1345 hypre_Index loop_size;
1346 hypre_IndexRef start;
1347 hypre_Index unit_stride;
1348
1349 int i, j;
1350 int loopi, loopj, loopk;
1351
1352 /* add_values may be off-proc. Communication needed, which is triggered
1353 if one of the OffProcAdd is 1 */
1354 sum_OffProcAdd= 0;
1355 MPI_Allreduce(&OffProcAdd, &sum_OffProcAdd, 1, MPI_INT, MPI_SUM,
1356 hypre_StructMatrixComm(matrix));
1357
1358 /* If matrix_data has an initial segment which is not mesh-based,
1359 it will not need to be communicated between processors, so
1360 matrix_data_comm will be set to point to the mesh-based part
1361 of the data */
1362
1363 /*-----------------------------------------------------------------------
1364 * If the CommPkg has not been set up, set it up
1365 *
1366 * The matrix data array is assumed to have two segments - an initial
1367 * segment of data constant over all space, followed by a segment with
1368 * comm_num_values matrix entries for each mesh element. The mesh-dependent
1369 * data is, of course, the only part relevent to communications.
1370 * For constant_coefficient==0, all the data is mesh-dependent.
1371 * For constant_coefficient==1, all data is constant.
1372 * For constant_coefficient==2, both segments are non-null.
1373 *-----------------------------------------------------------------------*/
1374
1375 constant_coefficient = hypre_StructMatrixConstantCoefficient( matrix );
1376
1377 mat_num_values = hypre_StructMatrixNumValues(matrix);
1378
1379 if ( constant_coefficient==0 )
1380 {
1381 comm_num_values = mat_num_values;
1382 }
1383 else if ( constant_coefficient==1 )
1384 {
1385 comm_num_values = 0;
1386 }
1387 else /* constant_coefficient==2 */
1388 {
1389 comm_num_values = 1;
1390 stencil = hypre_StructMatrixStencil(matrix);
1391 stencil_size = hypre_StructStencilSize(stencil);
1392 data_initial_offset = stencil_size;
1393 matrix_data_comm = &( matrix_data[data_initial_offset] );
1394 }
1395
1396 comm_pkg = hypre_StructMatrixCommPkg(matrix);
1397
1398 if (!comm_pkg)
1399 {
1400 hypre_CreateCommInfoFromNumGhost(hypre_StructMatrixGrid(matrix),
1401 num_ghost, &comm_info);
1402 hypre_CommPkgCreate(comm_info,
1403 hypre_StructMatrixDataSpace(matrix),
1404 hypre_StructMatrixDataSpace(matrix),
1405 comm_num_values,
1406 hypre_StructMatrixComm(matrix), &comm_pkg);
1407
1408 hypre_StructMatrixCommPkg(matrix) = comm_pkg;
1409
1410 /* inverse communication pattern if OffProcAdd values.
1411 Note that we cannot use the comm_info above because the
1412 add_num_ghost is generally different from num_ghost (add_num_ghost
1413 depends on the variable type. */
1414 if (sum_OffProcAdd)
1415 {
1416 /* since the off_proc add_values are located on the ghostlayer, we
1417 need "inverse" communication. */
1418 hypre_CommInfo *add_comm_info;
1419 hypre_CommInfo *inv_comm_info;
1420 hypre_CommPkg *inv_comm_pkg;
1421 int *add_num_ghost=
1422 hypre_StructMatrixAddNumGhost(matrix);
1423
1424 hypre_BoxArrayArray *send_boxes;
1425 hypre_BoxArrayArray *recv_boxes;
1426 int **send_procs;
1427 int **recv_procs;
1428 int **send_rboxnums;
1429 int **recv_rboxnums;
1430 hypre_BoxArrayArray *send_rboxes;
1431 hypre_BoxArray *box_array, *recv_array;
1432
1433 double *data;
1434 hypre_CommHandle *inv_comm_handle;
1435
1436 double *data_matrix;
1437 double *comm_data;
1438 int center_rank;
1439
1440 int s, xi;
1441
1442 hypre_CreateCommInfoFromNumGhost(hypre_StructMatrixGrid(matrix),
1443 add_num_ghost, &add_comm_info);
1444 /* inverse CommInfo achieved by switching send & recv structures of
1445 add_comm_info */
1446 send_boxes=
1447 hypre_BoxArrayArrayDuplicate(hypre_CommInfoRecvBoxes(add_comm_info));
1448 recv_boxes=
1449 hypre_BoxArrayArrayDuplicate(hypre_CommInfoSendBoxes(add_comm_info));
1450 send_rboxes= hypre_BoxArrayArrayDuplicate(send_boxes);
1451
1452 send_procs= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(send_boxes));
1453 recv_procs= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(recv_boxes));
1454 send_rboxnums= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(send_boxes));
1455 recv_rboxnums= hypre_CTAlloc(int *, hypre_BoxArrayArraySize(recv_boxes));
1456
1457 hypre_ForBoxArrayI(i, send_boxes)
1458 {
1459 box_array= hypre_BoxArrayArrayBoxArray(send_boxes, i);
1460 send_procs[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
1461 memcpy(send_procs[i], hypre_CommInfoRecvProcesses(add_comm_info)[i],
1462 hypre_BoxArraySize(box_array)*sizeof(int));
1463
1464 send_rboxnums[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
1465 memcpy(send_rboxnums[i], hypre_CommInfoRecvRBoxnums(add_comm_info)[i],
1466 hypre_BoxArraySize(box_array)*sizeof(int));
1467 }
1468
1469 hypre_ForBoxArrayI(i, recv_boxes)
1470 {
1471 box_array= hypre_BoxArrayArrayBoxArray(recv_boxes, i);
1472 recv_procs[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
1473 memcpy(recv_procs[i], hypre_CommInfoSendProcesses(add_comm_info)[i],
1474 hypre_BoxArraySize(box_array)*sizeof(int));
1475
1476 recv_rboxnums[i]= hypre_CTAlloc(int, hypre_BoxArraySize(box_array));
1477 memcpy(recv_rboxnums[i], hypre_CommInfoSendRBoxnums(add_comm_info)[i],
1478 hypre_BoxArraySize(box_array)*sizeof(int));
1479 }
1480
1481 hypre_CommInfoCreate(send_boxes, recv_boxes, send_procs, recv_procs,
1482 send_rboxnums, recv_rboxnums, send_rboxes,
1483 &inv_comm_info);
1484
1485 hypre_CommPkgCreate(inv_comm_info,
1486 hypre_StructMatrixDataSpace(matrix),
1487 hypre_StructMatrixDataSpace(matrix),
1488 comm_num_values,
1489 hypre_StructMatrixComm(matrix),
1490 &inv_comm_pkg);
1491
1492 /* communicate the add value entries */
1493 data= hypre_CTAlloc(double, hypre_StructMatrixDataSize(matrix));
1494 hypre_InitializeCommunication(inv_comm_pkg,
1495 hypre_StructMatrixData(matrix),
1496 data,
1497 &inv_comm_handle);
1498 hypre_FinalizeCommunication(inv_comm_handle);
1499
1500 /* this proc will recved data in it's send_boxes of add_comm_info, or
1501 equivalently, in the recv_boxes of inv_comm_info. Since inv_comm_info
1502 has already been destroyed, we use the send_boxes of add_comm_info */
1503 stencil = hypre_StructMatrixStencil(matrix);
1504 stencil_size = hypre_StructStencilSize(stencil);
1505 hypre_SetIndex(unit_stride, 0, 0, 0);
1506 center_rank = hypre_StructStencilElementRank(stencil, unit_stride);
1507 hypre_SetIndex(unit_stride, 1, 1, 1);
1508
1509 box_array= hypre_StructGridBoxes(hypre_StructMatrixGrid(matrix));
1510 hypre_ForBoxI(i, box_array)
1511 {
1512 recv_array=
1513 hypre_BoxArrayArrayBoxArray(hypre_CommInfoSendBoxes(add_comm_info), i);
1514
1515 data_box= hypre_BoxArrayBox(hypre_StructMatrixDataSpace(matrix), i);
1516 for (s= 0; s< stencil_size; s++)
1517 {
1518 data_matrix= hypre_StructMatrixBoxData(matrix, i, s);
1519 comm_data=(data + hypre_StructMatrixDataIndices(matrix)[i][s]);
1520
1521 hypre_ForBoxI(j, recv_array)
1522 {
1523 box = hypre_BoxArrayBox(recv_array, j);
1524 start= hypre_BoxIMin(box);
1525 hypre_BoxGetSize(box, loop_size);
1526
1527 /* only adding offproc values. */
1528 if (s != center_rank)
1529 {
1530 hypre_BoxLoop1Begin(loop_size,
1531 data_box, start, unit_stride, xi)
1532
1533 hypre_BoxLoop1For(loopi, loopj, loopk, xi)
1534 {
1535 data_matrix[xi] += comm_data[xi];
1536 }
1537 hypre_BoxLoop1End(xi);
1538 }
1539
1540 else
1541 {
1542 hypre_BoxLoop1Begin(loop_size,
1543 data_box, start, unit_stride, xi)
1544
1545 hypre_BoxLoop1For(loopi, loopj, loopk, xi)
1546 {
1547 if (comm_data[xi] != 1.0)
1548 {
1549 data_matrix[xi] += comm_data[xi];
1550 }
1551 }
1552 hypre_BoxLoop1End(xi);
1553 }
1554
1555 } /* hypre_ForBoxI(j, recv_array) */
1556 } /* for (s= 0; s< stencil_size; s++) */
1557 } /* hypre_ForBoxI(i, box_array) */
1558
1559 hypre_TFree(data);
1560 hypre_CommPkgDestroy(inv_comm_pkg);
1561 hypre_CommInfoDestroy(add_comm_info);
1562 }
1563 }
1564
1565 /*-----------------------------------------------------------------------
1566 * Update the ghost data
1567 * This takes care of the communication needs of all known functions
1568 * referencing the matrix.
1569 *
1570 * At present this is the only place where matrix data gets communicated.
1571 * However, comm_pkg is kept as long as the matrix is, in case some
1572 * future version hypre has a use for it - e.g. if the user replaces
1573 * a matrix with a very similar one, we may not want to recompute comm_pkg.
1574 *-----------------------------------------------------------------------*/
1575
1576 if ( constant_coefficient!=1 )
1577 {
1578 hypre_InitializeCommunication( comm_pkg,
1579 matrix_data_comm,
1580 matrix_data_comm,
1581 &comm_handle );
1582 hypre_FinalizeCommunication( comm_handle );
1583 }
1584
1585 return(ierr);
1586}
1587
1588/*--------------------------------------------------------------------------
1589 * hypre_StructMatrixSetNumGhost
1590 *--------------------------------------------------------------------------*/
1591
1592int
1593hypre_StructMatrixSetNumGhost( hypre_StructMatrix *matrix,
1594 int *num_ghost )
1595{
1596 int ierr = 0;
1597 int i;
1598
1599 for (i = 0; i < 6; i++)
1600 hypre_StructMatrixNumGhost(matrix)[i] = num_ghost[i];
1601
1602 return ierr;
1603}
1604
1605/*--------------------------------------------------------------------------
1606 * hypre_StructMatrixSetConstantCoefficient
1607 * deprecated in user interface, in favor of SetConstantEntries.
1608 * left here for internal use
1609 *--------------------------------------------------------------------------*/
1610
1611int
1612hypre_StructMatrixSetConstantCoefficient( hypre_StructMatrix *matrix,
1613 int constant_coefficient )
1614{
1615 int ierr = 0;
1616
1617 hypre_StructMatrixConstantCoefficient(matrix) = constant_coefficient;
1618
1619 return ierr;
1620}
1621
1622/*--------------------------------------------------------------------------
1623 * hypre_StructMatrixSetConstantEntries
1624 * - nentries is the number of array entries
1625 * - Each int entries[i] is an index into the shape array of the stencil of the
1626 * matrix
1627 * In the present version, only three possibilites are recognized:
1628 * - no entries constant (constant_coefficient==0)
1629 * - all entries constant (constant_coefficient==1)
1630 * - all but the diagonal entry constant (constant_coefficient==2)
1631 * If something else is attempted, this function will return a nonzero error.
1632 * In the present version, if this function is called more than once, only
1633 * the last call will take effect.
1634 *--------------------------------------------------------------------------*/
1635
1636
1637int hypre_StructMatrixSetConstantEntries( hypre_StructMatrix *matrix,
1638 int nentries,
1639 int *entries )
1640{
1641 /* We make an array offdconst corresponding to the stencil's shape array,
1642 and use "entries" to fill it with flags - 1 for constant, 0 otherwise.
1643 By counting the nonzeros in offdconst, and by checking whether its
1644 diagonal entry is nonzero, we can distinguish among the three
1645 presently legal values of constant_coefficient, and detect input errors.
1646 We do not need to treat duplicates in "entries" as an error condition.
1647 */
1648 int ierr = 0;
1649 hypre_StructStencil *stencil = hypre_StructMatrixUserStencil(matrix);
1650 /* ... Stencil doesn't exist yet */
1651 int stencil_size = hypre_StructStencilSize(stencil);
1652 int *offdconst = hypre_CTAlloc(int, stencil_size);
1653 /* ... note: CTAlloc initializes to 0 (normally it works by calling calloc) */
1654 int nconst = 0;
1655 int constant_coefficient, diag_rank;
1656 hypre_Index diag_index;
1657 int i, j;
1658
1659 for ( i=0; i<nentries; ++i )
1660 {
1661 offdconst[ entries[i] ] = 1;
1662 }
1663 for ( j=0; j<stencil_size; ++j )
1664 {
1665 nconst += offdconst[j];
1666 }
1667 if ( nconst<=0 ) constant_coefficient=0;
1668 else if ( nconst>=stencil_size ) constant_coefficient=1;
1669 else
1670 {
1671 hypre_SetIndex(diag_index, 0, 0, 0);
1672 diag_rank = hypre_StructStencilElementRank( stencil, diag_index );
1673 if ( offdconst[diag_rank]==0 )
1674 {
1675 constant_coefficient=2;
1676 if ( nconst!=(stencil_size-1) ) ++ierr;
1677 }
1678 else
1679 {
1680 constant_coefficient=0;
1681 ++ierr;
1682 }
1683 }
1684
1685 ierr += hypre_StructMatrixSetConstantCoefficient( matrix, constant_coefficient );
1686
1687 hypre_TFree(offdconst);
1688 return ierr;
1689}
1690
1691/*--------------------------------------------------------------------------
1692 * hypre_StructMatrixPrint
1693 *--------------------------------------------------------------------------*/
1694
1695int
1696hypre_StructMatrixPrint( const char *filename,
1697 hypre_StructMatrix *matrix,
1698 int all )
1699{
1700 int ierr = 0;
1701
1702 FILE *file;
1703 char new_filename[255];
1704
1705 hypre_StructGrid *grid;
1706 hypre_BoxArray *boxes;
1707
1708 hypre_StructStencil *stencil;
1709 hypre_Index *stencil_shape;
1710 int stencil_size;
1711 hypre_Index center_index;
1712
1713 int num_values;
1714
1715 hypre_BoxArray *data_space;
1716
1717 int *symm_elements;
1718
1719 int i, j;
1720 int constant_coefficient;
1721 int center_rank;
1722 int myid;
1723
1724 constant_coefficient = hypre_StructMatrixConstantCoefficient(matrix);
1725
1726 /*----------------------------------------
1727 * Open file
1728 *----------------------------------------*/
1729
1730#ifdef HYPRE_USE_PTHREADS
1731#if MPI_Comm_rank == hypre_thread_MPI_Comm_rank
1732#undef MPI_Comm_rank
1733#endif
1734#endif
1735
1736 MPI_Comm_rank(hypre_StructMatrixComm(matrix), &myid);
1737
1738 sprintf(new_filename, "%s.%05d", filename, myid);
1739
1740 if ((file = fopen(new_filename, "w")) == NULL)
1741 {
1742 printf("Error: can't open output file %s\n", new_filename);
1743 exit(1);
1744 }
1745
1746 /*----------------------------------------
1747 * Print header info
1748 *----------------------------------------*/
1749
1750 fprintf(file, "StructMatrix\n");
1751
1752 fprintf(file, "\nSymmetric: %d\n", hypre_StructMatrixSymmetric(matrix));
1753 fprintf(file, "\nConstantCoefficient: %d\n", hypre_StructMatrixConstantCoefficient(matrix));
1754
1755 /* print grid info */
1756 fprintf(file, "\nGrid:\n");
1757 grid = hypre_StructMatrixGrid(matrix);
1758 hypre_StructGridPrint(file, grid);
1759
1760 /* print stencil info */
1761 fprintf(file, "\nStencil:\n");
1762 stencil = hypre_StructMatrixStencil(matrix);
1763 stencil_shape = hypre_StructStencilShape(stencil);
1764
1765 num_values = hypre_StructMatrixNumValues(matrix);
1766 symm_elements = hypre_StructMatrixSymmElements(matrix);
1767 fprintf(file, "%d\n", num_values);
1768 stencil_size = hypre_StructStencilSize(stencil);
1769 j = 0;
1770 for (i=0; i<stencil_size; i++)
1771 {
1772 if (symm_elements[i] < 0)
1773 {
1774 fprintf(file, "%d: %d %d %d\n", j++,
1775 hypre_IndexX(stencil_shape[i]),
1776 hypre_IndexY(stencil_shape[i]),
1777 hypre_IndexZ(stencil_shape[i]));
1778 }
1779 }
1780
1781 /*----------------------------------------
1782 * Print data
1783 *----------------------------------------*/
1784
1785 data_space = hypre_StructMatrixDataSpace(matrix);
1786
1787 if (all)
1788 boxes = data_space;
1789 else
1790 boxes = hypre_StructGridBoxes(grid);
1791
1792 fprintf(file, "\nData:\n");
1793 if ( constant_coefficient==1 )
1794 {
1795 hypre_PrintCCBoxArrayData(file, boxes, data_space, num_values,
1796 hypre_StructMatrixData(matrix));
1797 }
1798 else if ( constant_coefficient==2 )
1799 {
1800 hypre_SetIndex(center_index, 0, 0, 0);
1801 center_rank = hypre_StructStencilElementRank( stencil, center_index );
1802
1803 hypre_PrintCCVDBoxArrayData(file, boxes, data_space, num_values,
1804 center_rank, stencil_size, symm_elements,
1805 hypre_StructMatrixData(matrix));
1806 }
1807 else
1808 {
1809 hypre_PrintBoxArrayData(file, boxes, data_space, num_values,
1810 hypre_StructMatrixData(matrix));
1811 }
1812
1813 /*----------------------------------------
1814 * Close file
1815 *----------------------------------------*/
1816
1817 fflush(file);
1818 fclose(file);
1819
1820 return ierr;
1821}
1822
1823/*--------------------------------------------------------------------------
1824 * hypre_StructMatrixMigrate
1825 *--------------------------------------------------------------------------*/
1826
1827int
1828hypre_StructMatrixMigrate( hypre_StructMatrix *from_matrix,
1829 hypre_StructMatrix *to_matrix )
1830{
1831 hypre_CommInfo *comm_info;
1832 hypre_CommPkg *comm_pkg;
1833 hypre_CommHandle *comm_handle;
1834
1835 int ierr = 0;
1836 int constant_coefficient, comm_num_values;
1837 int stencil_size, mat_num_values;
1838 hypre_StructStencil *stencil;
1839 int data_initial_offset = 0;
1840 double *matrix_data_from = hypre_StructMatrixData(from_matrix);
1841 double *matrix_data_to = hypre_StructMatrixData(to_matrix);
1842 double *matrix_data_comm_from = matrix_data_from;
1843 double *matrix_data_comm_to = matrix_data_to;
1844
1845 /*------------------------------------------------------
1846 * Set up hypre_CommPkg
1847 *------------------------------------------------------*/
1848
1849 constant_coefficient = hypre_StructMatrixConstantCoefficient( from_matrix );
1850 hypre_assert( constant_coefficient == hypre_StructMatrixConstantCoefficient( to_matrix ) );
1851
1852 mat_num_values = hypre_StructMatrixNumValues(from_matrix);
1853 hypre_assert( mat_num_values = hypre_StructMatrixNumValues(to_matrix) );
1854
1855 if ( constant_coefficient==0 )
1856 {
1857 comm_num_values = mat_num_values;
1858 }
1859 else if ( constant_coefficient==1 )
1860 {
1861 comm_num_values = 0;
1862 }
1863 else /* constant_coefficient==2 */
1864 {
1865 comm_num_values = 1;
1866 stencil = hypre_StructMatrixStencil(from_matrix);
1867 stencil_size = hypre_StructStencilSize(stencil);
1868 hypre_assert(stencil_size ==
1869 hypre_StructStencilSize( hypre_StructMatrixStencil(to_matrix) ) );
1870 data_initial_offset = stencil_size;
1871 matrix_data_comm_from = &( matrix_data_from[data_initial_offset] );
1872 matrix_data_comm_to = &( matrix_data_to[data_initial_offset] );
1873 }
1874
1875 hypre_CreateCommInfoFromGrids(hypre_StructMatrixGrid(from_matrix),
1876 hypre_StructMatrixGrid(to_matrix),
1877 &comm_info);
1878 hypre_CommPkgCreate(comm_info,
1879 hypre_StructMatrixDataSpace(from_matrix),
1880 hypre_StructMatrixDataSpace(to_matrix),
1881 comm_num_values,
1882 hypre_StructMatrixComm(from_matrix), &comm_pkg);
1883 /* is this correct for periodic? */
1884
1885 /*-----------------------------------------------------------------------
1886 * Migrate the matrix data
1887 *-----------------------------------------------------------------------*/
1888
1889 if ( constant_coefficient!=1 )
1890 {
1891 hypre_InitializeCommunication( comm_pkg,
1892 matrix_data_comm_from,
1893 matrix_data_comm_to,
1894 &comm_handle );
1895 hypre_FinalizeCommunication( comm_handle );
1896 }
1897
1898 return ierr;
1899}
1900
1901/*--------------------------------------------------------------------------
1902 * hypre_StructMatrixRead
1903 *--------------------------------------------------------------------------*/
1904
1905hypre_StructMatrix *
1906hypre_StructMatrixRead( MPI_Comm comm,
1907 const char *filename,
1908 int *num_ghost )
1909{
1910 FILE *file;
1911 char new_filename[255];
1912
1913 hypre_StructMatrix *matrix;
1914
1915 hypre_StructGrid *grid;
1916 hypre_BoxArray *boxes;
1917 int dim;
1918
1919 hypre_StructStencil *stencil;
1920 hypre_Index *stencil_shape;
1921 int stencil_size, real_stencil_size;
1922
1923 int num_values;
1924
1925 hypre_BoxArray *data_space;
1926
1927 int symmetric;
1928 int constant_coefficient;
1929
1930 int i, idummy;
1931
1932 int myid;
1933
1934 /*----------------------------------------
1935 * Open file
1936 *----------------------------------------*/
1937
1938#ifdef HYPRE_USE_PTHREADS
1939#if MPI_Comm_rank == hypre_thread_MPI_Comm_rank
1940#undef MPI_Comm_rank
1941#endif
1942#endif
1943
1944 MPI_Comm_rank(comm, &myid );
1945
1946 sprintf(new_filename, "%s.%05d", filename, myid);
1947
1948 if ((file = fopen(new_filename, "r")) == NULL)
1949 {
1950 printf("Error: can't open output file %s\n", new_filename);
1951 exit(1);
1952 }
1953
1954 /*----------------------------------------
1955 * Read header info
1956 *----------------------------------------*/
1957
1958 fscanf(file, "StructMatrix\n");
1959
1960 fscanf(file, "\nSymmetric: %d\n", &symmetric);
1961 fscanf(file, "\nConstantCoefficient: %d\n", &constant_coefficient);
1962
1963 /* read grid info */
1964 fscanf(file, "\nGrid:\n");
1965 hypre_StructGridRead(comm,file,&grid);
1966
1967 /* read stencil info */
1968 fscanf(file, "\nStencil:\n");
1969 dim = hypre_StructGridDim(grid);
1970 fscanf(file, "%d\n", &stencil_size);
1971 if (symmetric) { real_stencil_size = 2*stencil_size-1; }
1972 else { real_stencil_size = stencil_size; }
1973 /* ... real_stencil_size is the stencil size of the matrix after it's fixed up
1974 by the call (if any) of hypre_StructStencilSymmetrize from
1975 hypre_StructMatrixInitializeShell.*/
1976 stencil_shape = hypre_CTAlloc(hypre_Index, stencil_size);
1977 for (i = 0; i < stencil_size; i++)
1978 {
1979 fscanf(file, "%d: %d %d %d\n", &idummy,
1980 &hypre_IndexX(stencil_shape[i]),
1981 &hypre_IndexY(stencil_shape[i]),
1982 &hypre_IndexZ(stencil_shape[i]));
1983 }
1984 stencil = hypre_StructStencilCreate(dim, stencil_size, stencil_shape);
1985
1986 /*----------------------------------------
1987 * Initialize the matrix
1988 *----------------------------------------*/
1989
1990 matrix = hypre_StructMatrixCreate(comm, grid, stencil);
1991 hypre_StructMatrixSymmetric(matrix) = symmetric;
1992 hypre_StructMatrixConstantCoefficient(matrix) = constant_coefficient;
1993 hypre_StructMatrixSetNumGhost(matrix, num_ghost);
1994 hypre_StructMatrixInitialize(matrix);
1995
1996 /*----------------------------------------
1997 * Read data
1998 *----------------------------------------*/
1999
2000 boxes = hypre_StructGridBoxes(grid);
2001 data_space = hypre_StructMatrixDataSpace(matrix);
2002 num_values = hypre_StructMatrixNumValues(matrix);
2003
2004 fscanf(file, "\nData:\n");
2005 if ( constant_coefficient==0 )
2006 {
2007 hypre_ReadBoxArrayData(file, boxes, data_space, num_values,
2008 hypre_StructMatrixData(matrix));
2009 }
2010 else
2011 {
2012 hypre_assert( constant_coefficient<=2 );
2013 hypre_ReadBoxArrayData_CC( file, boxes, data_space,
2014 stencil_size, real_stencil_size,
2015 constant_coefficient,
2016 hypre_StructMatrixData(matrix));
2017 }
2018
2019 /*----------------------------------------
2020 * Assemble the matrix
2021 *----------------------------------------*/
2022
2023 hypre_StructMatrixAssemble(matrix);
2024
2025 /*----------------------------------------
2026 * Close file
2027 *----------------------------------------*/
2028
2029 fclose(file);
2030
2031 return matrix;
2032}
2033
Note: See TracBrowser for help on using the repository browser.