source: CIVL/examples/mpi-omp/AMG2013/sstruct_mv/sstruct_matrix.c@ 7d77e64

main test-branch
Last change on this file since 7d77e64 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: 36.7 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 *
17 * Member functions for hypre_SStructPMatrix class.
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23/*==========================================================================
24 * SStructPMatrix routines
25 *==========================================================================*/
26
27/*--------------------------------------------------------------------------
28 * hypre_SStructPMatrixRef
29 *--------------------------------------------------------------------------*/
30
31int
32hypre_SStructPMatrixRef( hypre_SStructPMatrix *matrix,
33 hypre_SStructPMatrix **matrix_ref )
34{
35 hypre_SStructPMatrixRefCount(matrix) ++;
36 *matrix_ref = matrix;
37
38 return hypre_error_flag;
39}
40
41/*--------------------------------------------------------------------------
42 * hypre_SStructPMatrixCreate
43 *--------------------------------------------------------------------------*/
44
45int
46hypre_SStructPMatrixCreate( MPI_Comm comm,
47 hypre_SStructPGrid *pgrid,
48 hypre_SStructStencil **stencils,
49 hypre_SStructPMatrix **pmatrix_ptr )
50{
51 hypre_SStructPMatrix *pmatrix;
52 int nvars;
53 int **smaps;
54 hypre_StructStencil ***sstencils;
55 hypre_StructMatrix ***smatrices;
56 int **symmetric;
57
58 hypre_StructStencil *sstencil;
59 int *vars;
60 hypre_Index *sstencil_shape;
61 int sstencil_size;
62 int new_dim;
63 int *new_sizes;
64 hypre_Index **new_shapes;
65 int size;
66 hypre_StructGrid *sgrid;
67
68 int vi, vj;
69 int i, j, k;
70
71 pmatrix = hypre_TAlloc(hypre_SStructPMatrix, 1);
72
73 hypre_SStructPMatrixComm(pmatrix) = comm;
74 hypre_SStructPMatrixPGrid(pmatrix) = pgrid;
75 hypre_SStructPMatrixStencils(pmatrix) = stencils;
76 nvars = hypre_SStructPGridNVars(pgrid);
77 hypre_SStructPMatrixNVars(pmatrix) = nvars;
78
79 /* create sstencils */
80 smaps = hypre_TAlloc(int *, nvars);
81 sstencils = hypre_TAlloc(hypre_StructStencil **, nvars);
82 new_sizes = hypre_TAlloc(int, nvars);
83 new_shapes = hypre_TAlloc(hypre_Index *, nvars);
84 size = 0;
85 for (vi = 0; vi < nvars; vi++)
86 {
87 sstencils[vi] = hypre_TAlloc(hypre_StructStencil *, nvars);
88 for (vj = 0; vj < nvars; vj++)
89 {
90 sstencils[vi][vj] = NULL;
91 new_sizes[vj] = 0;
92 }
93
94 sstencil = hypre_SStructStencilSStencil(stencils[vi]);
95 vars = hypre_SStructStencilVars(stencils[vi]);
96 sstencil_shape = hypre_StructStencilShape(sstencil);
97 sstencil_size = hypre_StructStencilSize(sstencil);
98
99 smaps[vi] = hypre_TAlloc(int, sstencil_size);
100 for (i = 0; i < sstencil_size; i++)
101 {
102 j = vars[i];
103 new_sizes[j]++;
104 }
105 for (vj = 0; vj < nvars; vj++)
106 {
107 if (new_sizes[vj])
108 {
109 new_shapes[vj] = hypre_TAlloc(hypre_Index, new_sizes[vj]);
110 new_sizes[vj] = 0;
111 }
112 }
113 for (i = 0; i < sstencil_size; i++)
114 {
115 j = vars[i];
116 k = new_sizes[j];
117 hypre_CopyIndex(sstencil_shape[i], new_shapes[j][k]);
118 smaps[vi][i] = k;
119 new_sizes[j]++;
120 }
121 new_dim = hypre_StructStencilDim(sstencil);
122 for (vj = 0; vj < nvars; vj++)
123 {
124 if (new_sizes[vj])
125 {
126 sstencils[vi][vj] = hypre_StructStencilCreate(new_dim,
127 new_sizes[vj],
128 new_shapes[vj]);
129 }
130 size = hypre_max(size, new_sizes[vj]);
131 }
132 }
133 hypre_SStructPMatrixSMaps(pmatrix) = smaps;
134 hypre_SStructPMatrixSStencils(pmatrix) = sstencils;
135 hypre_TFree(new_sizes);
136 hypre_TFree(new_shapes);
137
138 /* create smatrices */
139 smatrices = hypre_TAlloc(hypre_StructMatrix **, nvars);
140 for (vi = 0; vi < nvars; vi++)
141 {
142 smatrices[vi] = hypre_TAlloc(hypre_StructMatrix *, nvars);
143 for (vj = 0; vj < nvars; vj++)
144 {
145 smatrices[vi][vj] = NULL;
146 if (sstencils[vi][vj] != NULL)
147 {
148 sgrid = hypre_SStructPGridSGrid(pgrid, vi);
149 smatrices[vi][vj] =
150 hypre_StructMatrixCreate(comm, sgrid, sstencils[vi][vj]);
151 }
152 }
153 }
154 hypre_SStructPMatrixSMatrices(pmatrix) = smatrices;
155
156 /* create symmetric */
157 symmetric = hypre_TAlloc(int *, nvars);
158 for (vi = 0; vi < nvars; vi++)
159 {
160 symmetric[vi] = hypre_TAlloc(int, nvars);
161 for (vj = 0; vj < nvars; vj++)
162 {
163 symmetric[vi][vj] = 0;
164 }
165 }
166 hypre_SStructPMatrixSymmetric(pmatrix) = symmetric;
167
168 hypre_SStructPMatrixSEntriesSize(pmatrix) = size;
169 hypre_SStructPMatrixSEntries(pmatrix) = hypre_TAlloc(int, size);
170
171 hypre_SStructPMatrixRefCount(pmatrix) = 1;
172
173 *pmatrix_ptr = pmatrix;
174
175 return hypre_error_flag;
176}
177
178/*--------------------------------------------------------------------------
179 * hypre_SStructPMatrixDestroy
180 *--------------------------------------------------------------------------*/
181
182int
183hypre_SStructPMatrixDestroy( hypre_SStructPMatrix *pmatrix )
184{
185 hypre_SStructStencil **stencils;
186 int nvars;
187 int **smaps;
188 hypre_StructStencil ***sstencils;
189 hypre_StructMatrix ***smatrices;
190 int **symmetric;
191 int vi, vj;
192
193 if (pmatrix)
194 {
195 hypre_SStructPMatrixRefCount(pmatrix) --;
196 if (hypre_SStructPMatrixRefCount(pmatrix) == 0)
197 {
198 stencils = hypre_SStructPMatrixStencils(pmatrix);
199 nvars = hypre_SStructPMatrixNVars(pmatrix);
200 smaps = hypre_SStructPMatrixSMaps(pmatrix);
201 sstencils = hypre_SStructPMatrixSStencils(pmatrix);
202 smatrices = hypre_SStructPMatrixSMatrices(pmatrix);
203 symmetric = hypre_SStructPMatrixSymmetric(pmatrix);
204 for (vi = 0; vi < nvars; vi++)
205 {
206 HYPRE_SStructStencilDestroy(stencils[vi]);
207 hypre_TFree(smaps[vi]);
208 for (vj = 0; vj < nvars; vj++)
209 {
210 hypre_StructStencilDestroy(sstencils[vi][vj]);
211 hypre_StructMatrixDestroy(smatrices[vi][vj]);
212 }
213 hypre_TFree(sstencils[vi]);
214 hypre_TFree(smatrices[vi]);
215 hypre_TFree(symmetric[vi]);
216 }
217 hypre_TFree(stencils);
218 hypre_TFree(smaps);
219 hypre_TFree(sstencils);
220 hypre_TFree(smatrices);
221 hypre_TFree(symmetric);
222 hypre_TFree(hypre_SStructPMatrixSEntries(pmatrix));
223 hypre_TFree(pmatrix);
224 }
225 }
226
227 return hypre_error_flag;
228}
229
230/*--------------------------------------------------------------------------
231 * hypre_SStructPMatrixInitialize
232 *--------------------------------------------------------------------------*/
233int
234hypre_SStructPMatrixInitialize( hypre_SStructPMatrix *pmatrix )
235{
236 int nvars = hypre_SStructPMatrixNVars(pmatrix);
237 int **symmetric = hypre_SStructPMatrixSymmetric(pmatrix);
238 hypre_SStructPGrid *pgrid = hypre_SStructPMatrixPGrid(pmatrix);
239 HYPRE_SStructVariable *vartypes = hypre_SStructPGridVarTypes(pgrid);
240 int ndim = hypre_SStructPGridNDim(pgrid);
241
242 int num_ghost[6]= {1, 1, 1, 1, 1, 1};
243 hypre_StructMatrix *smatrix;
244 hypre_StructGrid *sgrid;
245
246 hypre_Index varoffset;
247 int vi, vj, d;
248
249 for (vi = 0; vi < nvars; vi++)
250 {
251 /* use variable vi add_numghost */
252 sgrid= hypre_SStructPGridSGrid(pgrid, vi);
253 hypre_SStructVariableGetOffset(vartypes[vi], ndim, varoffset);
254
255 for (vj = 0; vj < nvars; vj++)
256 {
257 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, vi, vj);
258 if (smatrix != NULL)
259 {
260 HYPRE_StructMatrixSetSymmetric(smatrix, symmetric[vi][vj]);
261 hypre_StructMatrixSetNumGhost(smatrix, num_ghost);
262
263 for (d = 0; d < 3; d++)
264 {
265 hypre_StructMatrixAddNumGhost(smatrix)[2*d]=
266 hypre_IndexD(varoffset, d);
267 hypre_StructMatrixAddNumGhost(smatrix)[2*d+1]=
268 hypre_IndexD(varoffset, d);
269 }
270
271 hypre_StructMatrixInitialize(smatrix);
272 }
273 }
274 }
275
276 return hypre_error_flag;
277}
278
279/*--------------------------------------------------------------------------
280 * hypre_SStructPMatrixSetValues
281 *--------------------------------------------------------------------------*/
282
283int
284hypre_SStructPMatrixSetValues( hypre_SStructPMatrix *pmatrix,
285 hypre_Index index,
286 int var,
287 int nentries,
288 int *entries,
289 double *values,
290 int add_to )
291{
292 hypre_SStructStencil *stencil = hypre_SStructPMatrixStencil(pmatrix, var);
293 int *smap = hypre_SStructPMatrixSMap(pmatrix, var);
294 int *vars = hypre_SStructStencilVars(stencil);
295 hypre_StructMatrix *smatrix;
296 int *sentries;
297 int i;
298
299 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, var, vars[entries[0]]);
300
301 sentries = hypre_SStructPMatrixSEntries(pmatrix);
302 for (i = 0; i < nentries; i++)
303 {
304 sentries[i] = smap[entries[i]];
305 }
306
307 hypre_StructMatrixSetValues(smatrix, index, nentries, sentries, values, add_to);
308
309 return hypre_error_flag;
310}
311
312/*--------------------------------------------------------------------------
313 * hypre_SStructPMatrixSetBoxValues
314 *--------------------------------------------------------------------------*/
315
316int
317hypre_SStructPMatrixSetBoxValues( hypre_SStructPMatrix *pmatrix,
318 hypre_Index ilower,
319 hypre_Index iupper,
320 int var,
321 int nentries,
322 int *entries,
323 double *values,
324 int add_to )
325{
326 hypre_SStructStencil *stencil = hypre_SStructPMatrixStencil(pmatrix, var);
327 int *smap = hypre_SStructPMatrixSMap(pmatrix, var);
328 int *vars = hypre_SStructStencilVars(stencil);
329 hypre_StructMatrix *smatrix;
330 hypre_Box *box;
331 int *sentries;
332 int i;
333
334 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, var, vars[entries[0]]);
335
336 box = hypre_BoxCreate();
337 hypre_CopyIndex(ilower, hypre_BoxIMin(box));
338 hypre_CopyIndex(iupper, hypre_BoxIMax(box));
339
340 sentries = hypre_SStructPMatrixSEntries(pmatrix);
341 for (i = 0; i < nentries; i++)
342 {
343 sentries[i] = smap[entries[i]];
344 }
345
346 hypre_StructMatrixSetBoxValues(smatrix, box, nentries, sentries, values, add_to);
347
348 hypre_BoxDestroy(box);
349
350 return hypre_error_flag;
351}
352
353/*--------------------------------------------------------------------------
354 * hypre_SStructPMatrixAssemble
355 *--------------------------------------------------------------------------*/
356
357int
358hypre_SStructPMatrixAssemble( hypre_SStructPMatrix *pmatrix )
359{
360 int nvars = hypre_SStructPMatrixNVars(pmatrix);
361 hypre_StructMatrix *smatrix;
362 int vi, vj;
363
364 for (vi = 0; vi < nvars; vi++)
365 {
366 for (vj = 0; vj < nvars; vj++)
367 {
368 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, vi, vj);
369 if (smatrix != NULL)
370 {
371 hypre_StructMatrixAssemble(smatrix);
372 }
373 }
374 }
375
376 return hypre_error_flag;
377}
378
379/*--------------------------------------------------------------------------
380 *--------------------------------------------------------------------------*/
381
382int
383hypre_SStructPMatrixSetSymmetric( hypre_SStructPMatrix *pmatrix,
384 int var,
385 int to_var,
386 int symmetric )
387{
388 int **pmsymmetric = hypre_SStructPMatrixSymmetric(pmatrix);
389
390 int vstart = var;
391 int vsize = 1;
392 int tstart = to_var;
393 int tsize = 1;
394 int v, t;
395
396 if (var == -1)
397 {
398 vstart = 0;
399 vsize = hypre_SStructPMatrixNVars(pmatrix);
400 }
401 if (to_var == -1)
402 {
403 tstart = 0;
404 tsize = hypre_SStructPMatrixNVars(pmatrix);
405 }
406
407 for (v = vstart; v < vsize; v++)
408 {
409 for (t = tstart; t < tsize; t++)
410 {
411 pmsymmetric[v][t] = symmetric;
412 }
413 }
414
415 return hypre_error_flag;
416}
417
418/*--------------------------------------------------------------------------
419 * hypre_SStructPMatrixPrint
420 *--------------------------------------------------------------------------*/
421
422int
423hypre_SStructPMatrixPrint( const char *filename,
424 hypre_SStructPMatrix *pmatrix,
425 int all )
426{
427 int nvars = hypre_SStructPMatrixNVars(pmatrix);
428 hypre_StructMatrix *smatrix;
429 int vi, vj;
430 char new_filename[255];
431
432 for (vi = 0; vi < nvars; vi++)
433 {
434 for (vj = 0; vj < nvars; vj++)
435 {
436 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, vi, vj);
437 if (smatrix != NULL)
438 {
439 sprintf(new_filename, "%s.%02d.%02d", filename, vi, vj);
440 hypre_StructMatrixPrint(new_filename, smatrix, all);
441 }
442 }
443 }
444
445 return hypre_error_flag;
446}
447
448/*==========================================================================
449 * SStructUMatrix routines
450 *==========================================================================*/
451
452/*--------------------------------------------------------------------------
453 * hypre_SStructUMatrixInitialize
454 *--------------------------------------------------------------------------*/
455
456int
457hypre_SStructUMatrixInitialize( hypre_SStructMatrix *matrix )
458{
459 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
460 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
461 hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph);
462 int nparts = hypre_SStructGraphNParts(graph);
463 hypre_SStructPGrid **pgrids = hypre_SStructGraphPGrids(graph);
464 hypre_SStructStencil ***stencils = hypre_SStructGraphStencils(graph);
465 int nUventries = hypre_SStructGraphNUVEntries(graph);
466 int *iUventries = hypre_SStructGraphIUVEntries(graph);
467 hypre_SStructUVEntry **Uventries = hypre_SStructGraphUVEntries(graph);
468 int **nvneighbors = hypre_SStructGridNVNeighbors(grid);
469 hypre_StructGrid *sgrid;
470 hypre_SStructStencil *stencil;
471 int *split;
472 int nvars;
473 int nrows, nnzs ;
474 int part, var, entry, i, j, k,m,b;
475 int *row_sizes;
476 int max_row_size;
477
478 int matrix_type = hypre_SStructMatrixObjectType(matrix);
479
480 hypre_Box *gridbox;
481 hypre_Box *loopbox;
482 hypre_Box *ghostbox;
483 hypre_BoxArray *boxes;
484 int *num_ghost;
485
486
487 HYPRE_IJMatrixSetObjectType(ijmatrix, HYPRE_PARCSR);
488
489 /* GEC1002 the ghlocalsize is used to set the number of rows */
490
491 if (matrix_type == HYPRE_PARCSR)
492 {
493 nrows = hypre_SStructGridLocalSize(grid);
494 }
495 if (matrix_type == HYPRE_SSTRUCT || matrix_type == HYPRE_STRUCT)
496 {
497 nrows = hypre_SStructGridGhlocalSize(grid) ;
498 }
499
500 /* set row sizes */
501 m = 0;
502 row_sizes = hypre_CTAlloc(int, nrows);
503 max_row_size = 0;
504 for (part = 0; part < nparts; part++)
505 {
506 nvars = hypre_SStructPGridNVars(pgrids[part]);
507 for (var = 0; var < nvars; var++)
508 {
509 sgrid = hypre_SStructPGridSGrid(pgrids[part], var);
510
511 stencil = stencils[part][var];
512 split = hypre_SStructMatrixSplit(matrix, part, var);
513 nnzs = 0;
514 for (entry = 0; entry < hypre_SStructStencilSize(stencil); entry++)
515 {
516 if (split[entry] == -1)
517 {
518 nnzs++;
519 }
520 }
521#if 0
522 /* TODO: For now, assume stencil is full/complete */
523 if (hypre_SStructMatrixSymmetric(matrix))
524 {
525 nnzs = 2*nnzs - 1;
526 }
527#endif
528
529 /**************/
530
531 boxes = hypre_StructGridBoxes(sgrid) ;
532 num_ghost = hypre_StructGridNumGhost(sgrid);
533 for (b = 0; b < hypre_BoxArraySize(boxes); b++)
534 {
535 gridbox = hypre_BoxArrayBox(boxes, b);
536 ghostbox = hypre_BoxCreate();
537 loopbox = hypre_BoxCreate();
538 hypre_CopyBox(gridbox,ghostbox);
539 hypre_BoxExpand(ghostbox,num_ghost);
540
541 if (matrix_type == HYPRE_SSTRUCT || matrix_type == HYPRE_STRUCT)
542 {
543 hypre_CopyBox(ghostbox,loopbox);
544 }
545 if (matrix_type == HYPRE_PARCSR)
546 {
547 hypre_CopyBox(gridbox,loopbox);
548 }
549
550 for (k = hypre_BoxIMinZ(loopbox); k <= hypre_BoxIMaxZ(loopbox); k++)
551 {
552 for (j = hypre_BoxIMinY(loopbox); j <= hypre_BoxIMaxY(loopbox); j++)
553 {
554 for (i = hypre_BoxIMinX(loopbox); i <= hypre_BoxIMaxX(loopbox); i++)
555 {
556 if ( ( ( i>=hypre_BoxIMinX(gridbox) )
557 && ( j>=hypre_BoxIMinY(gridbox) ) )
558 && ( k>=hypre_BoxIMinZ(gridbox) ) )
559 {
560 if ( ( ( i<=hypre_BoxIMaxX(gridbox) )
561 && ( j<=hypre_BoxIMaxY(gridbox) ) )
562 && ( k<=hypre_BoxIMaxZ(gridbox) ) )
563 {
564 row_sizes[m] = nnzs;
565 max_row_size = hypre_max(max_row_size, row_sizes[m]);
566 }
567 }
568 m++;
569 }
570 }
571 }
572 hypre_BoxDestroy(ghostbox);
573 hypre_BoxDestroy(loopbox);
574 }
575
576
577 if (nvneighbors[part][var])
578 {
579 max_row_size = hypre_max(max_row_size,
580 hypre_SStructStencilSize(stencil));
581 }
582
583
584 /*********************/
585 }
586 }
587
588 /* GEC0902 essentially for each UVentry we figure out how many extra columns
589 * we need to add to the rowsizes */
590
591 for (entry = 0; entry < nUventries; entry++)
592 {
593 i = iUventries[entry];
594 row_sizes[i] += hypre_SStructUVEntryNUEntries(Uventries[i]);
595 max_row_size = hypre_max(max_row_size, row_sizes[i]);
596 }
597
598 /* ZTODO: Update row_sizes based on neighbor off-part couplings */
599 HYPRE_IJMatrixSetRowSizes (ijmatrix, (const int *) row_sizes);
600
601 hypre_TFree(row_sizes);
602 hypre_SStructMatrixTmpColCoords(matrix) =
603 hypre_CTAlloc(HYPRE_BigInt, max_row_size);
604 hypre_SStructMatrixTmpCoeffs(matrix) =
605 hypre_CTAlloc(double, max_row_size);
606
607 /* GEC1002 at this point the processor has the partitioning (creation of ij) */
608
609 HYPRE_IJMatrixInitialize(ijmatrix);
610
611 return hypre_error_flag;
612}
613
614/*--------------------------------------------------------------------------
615 * hypre_SStructUMatrixSetValues
616 *
617 * (add_to > 0): add-to values
618 * (add_to = 0): set values
619 * (add_to < 0): get values
620 *--------------------------------------------------------------------------*/
621
622int
623hypre_SStructUMatrixSetValues( hypre_SStructMatrix *matrix,
624 int part,
625 hypre_Index index,
626 int var,
627 int nentries,
628 int *entries,
629 double *values,
630 int add_to )
631{
632 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
633 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
634 hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph);
635 hypre_SStructStencil *stencil = hypre_SStructGraphStencil(graph, part, var);
636 int *vars = hypre_SStructStencilVars(stencil);
637 hypre_Index *shape = hypre_SStructStencilShape(stencil);
638 int size = hypre_SStructStencilSize(stencil);
639 hypre_IndexRef offset;
640 hypre_Index to_index;
641 hypre_SStructUVEntry *Uventry;
642 hypre_BoxMapEntry *map_entry;
643 hypre_SStructMapInfo *entry_info;
644 HYPRE_BigInt row_coord;
645 HYPRE_BigInt *col_coords;
646 int ncoeffs;
647 double *coeffs;
648 int i, entry;
649 int proc, myproc;
650 /* GEC1002 the matrix type */
651 int matrix_type = hypre_SStructMatrixObjectType(matrix);
652
653 hypre_SStructGridFindMapEntry(grid, part, index, var, &map_entry);
654 if (map_entry == NULL)
655 {
656 hypre_error_in_arg(1);
657 hypre_error_in_arg(2);
658 hypre_error_in_arg(3);
659 /* RDF: This printing shouldn't be on by default */
660 printf("Warning: Attempt to set coeffs for point not in grid\n");
661 printf("hypre_SStructUMatrixSetValues call aborted for grid point\n");
662 printf(" part=%d, var=%d, index=(%d, %d, %d)\n", part, var,
663 hypre_IndexD(index,0),
664 hypre_IndexD(index,1),
665 hypre_IndexD(index,2) );
666 return hypre_error_flag;
667 }
668 else
669 {
670 hypre_BoxMapEntryGetInfo(map_entry, (void **) &entry_info);
671 }
672
673 /* Only Set values if I am the owner process; off-process AddTo and Get
674 * values are done by IJ */
675 if (!add_to)
676 {
677 hypre_SStructMapEntryGetProcess(map_entry, &proc);
678 MPI_Comm_rank(hypre_SStructGridComm(grid), &myproc);
679 if (proc != myproc)
680 {
681 return hypre_error_flag;
682 }
683 }
684
685 /* GEC1002 get the rank using the function with the type=matrixtype*/
686 hypre_SStructMapEntryGetGlobalRank(map_entry, index, &row_coord, matrix_type);
687
688
689 col_coords = hypre_SStructMatrixTmpColCoords(matrix);
690 coeffs = hypre_SStructMatrixTmpCoeffs(matrix);
691 ncoeffs = 0;
692 for (i = 0; i < nentries; i++)
693 {
694 entry = entries[i];
695
696 if (entry < size)
697 {
698 /* stencil entries */
699 offset = shape[entry];
700 hypre_IndexX(to_index) = hypre_IndexX(index) + hypre_IndexX(offset);
701 hypre_IndexY(to_index) = hypre_IndexY(index) + hypre_IndexY(offset);
702 hypre_IndexZ(to_index) = hypre_IndexZ(index) + hypre_IndexZ(offset);
703
704 hypre_SStructGridFindMapEntry(grid, part, to_index, vars[entry],
705 &map_entry);
706
707 if (map_entry != NULL)
708 {
709
710
711 hypre_SStructMapEntryGetGlobalRank(map_entry, to_index,
712 &col_coords[ncoeffs],matrix_type);
713
714
715 coeffs[ncoeffs] = values[i];
716 ncoeffs++;
717 }
718 }
719 else
720 {
721 /* non-stencil entries */
722 entry -= size;
723 hypre_SStructGraphFindUVEntry(graph, part, index, var, &Uventry);
724
725 col_coords[ncoeffs] = hypre_SStructUVEntryRank(Uventry, entry);
726 coeffs[ncoeffs] = values[i];
727 ncoeffs++;
728 }
729 }
730
731 if (add_to > 0)
732 {
733 HYPRE_IJMatrixAddToValues(ijmatrix, 1, &ncoeffs, &row_coord,
734 (const HYPRE_BigInt *) col_coords,
735 (const double *) coeffs);
736 }
737 else if (add_to > -1)
738 {
739 HYPRE_IJMatrixSetValues(ijmatrix, 1, &ncoeffs, &row_coord,
740 (const HYPRE_BigInt *) col_coords,
741 (const double *) coeffs);
742 }
743 else
744 {
745 HYPRE_IJMatrixGetValues(ijmatrix, 1, &ncoeffs, &row_coord,
746 col_coords, values);
747 }
748
749 return hypre_error_flag;
750}
751
752/*--------------------------------------------------------------------------
753 * Note: Entries must all be of type stencil or non-stencil, but not both.
754 *
755 * (add_to > 0): add-to values
756 * (add_to = 0): set values
757 * (add_to < 0): get values
758 *--------------------------------------------------------------------------*/
759
760int
761hypre_SStructUMatrixSetBoxValues( hypre_SStructMatrix *matrix,
762 int part,
763 hypre_Index ilower,
764 hypre_Index iupper,
765 int var,
766 int nentries,
767 int *entries,
768 double *values,
769 int add_to )
770{
771 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
772 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
773 hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph);
774 hypre_SStructStencil *stencil = hypre_SStructGraphStencil(graph, part, var);
775 int *vars = hypre_SStructStencilVars(stencil);
776 hypre_Index *shape = hypre_SStructStencilShape(stencil);
777 int size = hypre_SStructStencilSize(stencil);
778 hypre_IndexRef offset;
779 hypre_BoxMap *map;
780 hypre_BoxMapEntry **map_entries;
781 int nmap_entries;
782 hypre_BoxMapEntry **map_to_entries;
783 int nmap_to_entries;
784 int nrows;
785 int *ncols;
786 HYPRE_BigInt *rows;
787 HYPRE_BigInt *cols;
788 double *ijvalues;
789 hypre_Box *box;
790 hypre_Box *to_box;
791 hypre_Box *map_box;
792 hypre_Box *int_box;
793 hypre_Index index;
794 hypre_Index rs, cs;
795 int sy, sz;
796 HYPRE_BigInt row_base, col_base;
797 int val_base;
798 int e, entry, ii, jj, i, j, k;
799 int proc, myproc;
800 /* GEC1002 the matrix type */
801 int matrix_type = hypre_SStructMatrixObjectType(matrix);
802
803 box = hypre_BoxCreate();
804
805 /*------------------------------------------
806 * all stencil entries
807 *------------------------------------------*/
808
809 if (entries[0] < size)
810 {
811 to_box = hypre_BoxCreate();
812 map_box = hypre_BoxCreate();
813 int_box = hypre_BoxCreate();
814
815 hypre_CopyIndex(ilower, hypre_BoxIMin(box));
816 hypre_CopyIndex(iupper, hypre_BoxIMax(box));
817 /* ZTODO: check that this change fixes multiple-entry problem */
818 nrows = hypre_BoxVolume(box)*nentries;
819 ncols = hypre_CTAlloc(int, nrows);
820 for (i = 0; i < nrows; i++)
821 {
822 ncols[i] = 1;
823 }
824 rows = hypre_CTAlloc(HYPRE_BigInt, nrows);
825 cols = hypre_CTAlloc(HYPRE_BigInt, nrows);
826 ijvalues = hypre_CTAlloc(double, nrows);
827
828 sy = (hypre_IndexX(iupper) - hypre_IndexX(ilower) + 1);
829 sz = (hypre_IndexY(iupper) - hypre_IndexY(ilower) + 1) * sy;
830
831 map = hypre_SStructGridMap(grid, part, var);
832 hypre_BoxMapIntersect(map, ilower, iupper, &map_entries, &nmap_entries);
833
834 for (ii = 0; ii < nmap_entries; ii++)
835 {
836 /* Only Set values if I am the owner process; off-process AddTo and Get
837 * values are done by IJ */
838 if (!add_to)
839 {
840 hypre_SStructMapEntryGetProcess(map_entries[ii], &proc);
841 MPI_Comm_rank(hypre_SStructGridComm(grid), &myproc);
842 if (proc != myproc)
843 {
844 continue;
845 }
846 }
847
848 /* GEC1002 introducing the strides based on the type of the matrix */
849 hypre_SStructMapEntryGetStrides(map_entries[ii], rs, matrix_type);
850
851 hypre_CopyIndex(ilower, hypre_BoxIMin(box));
852 hypre_CopyIndex(iupper, hypre_BoxIMax(box));
853 hypre_BoxMapEntryGetExtents(map_entries[ii],
854 hypre_BoxIMin(map_box),
855 hypre_BoxIMax(map_box));
856 hypre_IntersectBoxes(box, map_box, int_box);
857 hypre_CopyBox(int_box, box);
858
859 nrows = 0;
860 for (e = 0; e < nentries; e++)
861 {
862 entry = entries[e];
863
864 hypre_CopyBox(box, to_box);
865
866 offset = shape[entry];
867 hypre_BoxIMinX(to_box) += hypre_IndexX(offset);
868 hypre_BoxIMinY(to_box) += hypre_IndexY(offset);
869 hypre_BoxIMinZ(to_box) += hypre_IndexZ(offset);
870 hypre_BoxIMaxX(to_box) += hypre_IndexX(offset);
871 hypre_BoxIMaxY(to_box) += hypre_IndexY(offset);
872 hypre_BoxIMaxZ(to_box) += hypre_IndexZ(offset);
873
874 map = hypre_SStructGridMap(grid, part, vars[entry]);
875 hypre_BoxMapIntersect(map, hypre_BoxIMin(to_box),
876 hypre_BoxIMax(to_box),
877 &map_to_entries, &nmap_to_entries );
878
879 for (jj = 0; jj < nmap_to_entries; jj++)
880 {
881
882 /* GEC1002 introducing the strides based on the type of the matrix */
883
884 hypre_SStructMapEntryGetStrides(map_to_entries[jj], cs, matrix_type);
885
886 hypre_BoxMapEntryGetExtents(map_to_entries[jj],
887 hypre_BoxIMin(map_box),
888 hypre_BoxIMax(map_box));
889 hypre_IntersectBoxes(to_box, map_box, int_box);
890
891 hypre_CopyIndex(hypre_BoxIMin(int_box), index);
892
893 /* GEC1002 introducing the rank based on the type of the matrix */
894
895 hypre_SStructMapEntryGetGlobalRank(map_to_entries[jj],
896 index, &col_base,matrix_type);
897
898 hypre_IndexX(index) -= hypre_IndexX(offset);
899 hypre_IndexY(index) -= hypre_IndexY(offset);
900 hypre_IndexZ(index) -= hypre_IndexZ(offset);
901
902 /* GEC1002 introducing the rank based on the type of the matrix */
903
904 hypre_SStructMapEntryGetGlobalRank(map_entries[ii],
905 index, &row_base,matrix_type);
906
907 hypre_IndexX(index) -= hypre_IndexX(ilower);
908 hypre_IndexY(index) -= hypre_IndexY(ilower);
909 hypre_IndexZ(index) -= hypre_IndexZ(ilower);
910 val_base = e + (hypre_IndexX(index) +
911 hypre_IndexY(index)*sy +
912 hypre_IndexZ(index)*sz) * nentries;
913
914 for (k = 0; k < hypre_BoxSizeZ(int_box); k++)
915 {
916 for (j = 0; j < hypre_BoxSizeY(int_box); j++)
917 {
918 for (i = 0; i < hypre_BoxSizeX(int_box); i++)
919 {
920 rows[nrows] = row_base + (HYPRE_BigInt)(i*rs[0] + j*rs[1] + k*rs[2]);
921 cols[nrows] = col_base + (HYPRE_BigInt)(i*cs[0] + j*cs[1] + k*cs[2]);
922 ijvalues[nrows] =
923 values[val_base + (i + j*sy + k*sz)*nentries];
924 nrows++;
925 }
926 }
927 }
928 }
929
930 hypre_TFree(map_to_entries);
931 }
932
933 /*------------------------------------------
934 * set IJ values one stencil entry at a time
935 *------------------------------------------*/
936
937 if (add_to > 0)
938 {
939 HYPRE_IJMatrixAddToValues(ijmatrix, nrows, ncols,
940 (const HYPRE_BigInt *) rows,
941 (const HYPRE_BigInt *) cols,
942 (const double *) ijvalues);
943 }
944 else if (add_to > -1)
945 {
946 HYPRE_IJMatrixSetValues(ijmatrix, nrows, ncols,
947 (const HYPRE_BigInt *) rows,
948 (const HYPRE_BigInt *) cols,
949 (const double *) ijvalues);
950 }
951 else
952 {
953 HYPRE_IJMatrixGetValues(ijmatrix, nrows, ncols, rows, cols, values);
954 }
955 }
956
957 hypre_TFree(map_entries);
958
959 hypre_TFree(ncols);
960 hypre_TFree(rows);
961 hypre_TFree(cols);
962 hypre_TFree(ijvalues);
963
964 hypre_BoxDestroy(to_box);
965 hypre_BoxDestroy(map_box);
966 hypre_BoxDestroy(int_box);
967 }
968
969 /*------------------------------------------
970 * non-stencil entries
971 *------------------------------------------*/
972
973 else
974 {
975 hypre_CopyIndex(ilower, hypre_BoxIMin(box));
976 hypre_CopyIndex(iupper, hypre_BoxIMax(box));
977
978 for (k = hypre_BoxIMinZ(box); k <= hypre_BoxIMaxZ(box); k++)
979 {
980 for (j = hypre_BoxIMinY(box); j <= hypre_BoxIMaxY(box); j++)
981 {
982 for (i = hypre_BoxIMinX(box); i <= hypre_BoxIMaxX(box); i++)
983 {
984 hypre_SetIndex(index, i, j, k);
985 hypre_SStructUMatrixSetValues(matrix, part, index, var,
986 nentries, entries, values, add_to);
987 values += nentries;
988 }
989 }
990 }
991 }
992
993 hypre_BoxDestroy(box);
994
995 return hypre_error_flag;
996}
997
998/*--------------------------------------------------------------------------
999 * hypre_SStructUMatrixAssemble
1000 *--------------------------------------------------------------------------*/
1001
1002int
1003hypre_SStructUMatrixAssemble( hypre_SStructMatrix *matrix )
1004{
1005 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
1006
1007 HYPRE_IJMatrixAssemble(ijmatrix);
1008 HYPRE_IJMatrixGetObject(ijmatrix,
1009 (void **) &hypre_SStructMatrixParCSRMatrix(matrix));
1010
1011 return hypre_error_flag;
1012}
1013
1014/*==========================================================================
1015 * SStructMatrix routines
1016 *==========================================================================*/
1017
1018/*--------------------------------------------------------------------------
1019 * hypre_SStructMatrixRef
1020 *--------------------------------------------------------------------------*/
1021
1022int
1023hypre_SStructMatrixRef( hypre_SStructMatrix *matrix,
1024 hypre_SStructMatrix **matrix_ref )
1025{
1026 hypre_SStructMatrixRefCount(matrix) ++;
1027 *matrix_ref = matrix;
1028
1029 return hypre_error_flag;
1030}
1031
1032/*--------------------------------------------------------------------------
1033 * hypre_SStructMatrixSplitEntries
1034 *--------------------------------------------------------------------------*/
1035
1036int
1037hypre_SStructMatrixSplitEntries( hypre_SStructMatrix *matrix,
1038 int part,
1039 int var,
1040 int nentries,
1041 int *entries,
1042 int *nSentries_ptr,
1043 int **Sentries_ptr,
1044 int *nUentries_ptr,
1045 int **Uentries_ptr )
1046{
1047 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
1048 int *split = hypre_SStructMatrixSplit(matrix, part, var);
1049 hypre_SStructStencil *stencil = hypre_SStructGraphStencil(graph, part, var);
1050 int entry;
1051 int i;
1052
1053 int nSentries = 0;
1054 int *Sentries = hypre_SStructMatrixSEntries(matrix);
1055 int nUentries = 0;
1056 int *Uentries = hypre_SStructMatrixUEntries(matrix);
1057
1058 for (i = 0; i < nentries; i++)
1059 {
1060 entry = entries[i];
1061 if (entry < hypre_SStructStencilSize(stencil))
1062 {
1063 /* stencil entries */
1064 if (split[entry] > -1)
1065 {
1066 Sentries[nSentries] = split[entry];
1067 nSentries++;
1068 }
1069 else
1070 {
1071 Uentries[nUentries] = entry;
1072 nUentries++;
1073 }
1074 }
1075 else
1076 {
1077 /* non-stencil entries */
1078 Uentries[nUentries] = entry;
1079 nUentries++;
1080 }
1081 }
1082
1083 *nSentries_ptr = nSentries;
1084 *Sentries_ptr = Sentries;
1085 *nUentries_ptr = nUentries;
1086 *Uentries_ptr = Uentries;
1087
1088 return hypre_error_flag;
1089}
1090
Note: See TracBrowser for help on using the repository browser.