source: CIVL/examples/mpi-omp/AMG2013/sstruct_mv/HYPRE_sstruct_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: 31.8 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 * HYPRE_SStructMatrix interface
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22
23/*--------------------------------------------------------------------------
24 *--------------------------------------------------------------------------*/
25
26int
27HYPRE_SStructMatrixCreate( MPI_Comm comm,
28 HYPRE_SStructGraph graph,
29 HYPRE_SStructMatrix *matrix_ptr )
30{
31 int ierr = 0;
32
33 /* GEC1202 grid not needed */
34 /* hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph); */
35 hypre_SStructStencil ***stencils = hypre_SStructGraphStencils(graph);
36
37 hypre_SStructMatrix *matrix;
38 int ***splits;
39 int nparts;
40 hypre_SStructPMatrix **pmatrices;
41 int ***symmetric;
42
43 hypre_SStructPGrid *pgrid;
44 int nvars;
45
46 int stencil_size;
47 int *stencil_vars;
48 int pstencil_size;
49
50 HYPRE_SStructVariable vitype, vjtype;
51 int part, vi, vj, i;
52 int size;
53
54 matrix = hypre_TAlloc(hypre_SStructMatrix, 1);
55
56 hypre_SStructMatrixComm(matrix) = comm;
57 hypre_SStructMatrixNDim(matrix) = hypre_SStructGraphNDim(graph);
58 hypre_SStructGraphRef(graph, &hypre_SStructMatrixGraph(matrix));
59
60 /* compute S/U-matrix split */
61 nparts = hypre_SStructGraphNParts(graph);
62 hypre_SStructMatrixNParts(matrix) = nparts;
63 splits = hypre_TAlloc(int **, nparts);
64 hypre_SStructMatrixSplits(matrix) = splits;
65 pmatrices = hypre_TAlloc(hypre_SStructPMatrix *, nparts);
66 hypre_SStructMatrixPMatrices(matrix) = pmatrices;
67 symmetric = hypre_TAlloc(int **, nparts);
68 hypre_SStructMatrixSymmetric(matrix) = symmetric;
69 for (part = 0; part < nparts; part++)
70 {
71 pgrid = hypre_SStructGraphPGrid(graph, part);
72 nvars = hypre_SStructPGridNVars(pgrid);
73 splits[part] = hypre_TAlloc(int *, nvars);
74 symmetric[part] = hypre_TAlloc(int *, nvars);
75 for (vi = 0; vi < nvars; vi++)
76 {
77 stencil_size = hypre_SStructStencilSize(stencils[part][vi]);
78 stencil_vars = hypre_SStructStencilVars(stencils[part][vi]);
79 pstencil_size = 0;
80 splits[part][vi] = hypre_TAlloc(int, stencil_size);
81 symmetric[part][vi] = hypre_TAlloc(int, nvars);
82 for (i = 0; i < stencil_size; i++)
83 {
84 vj = stencil_vars[i];
85 vitype = hypre_SStructPGridVarType(pgrid, vi);
86 vjtype = hypre_SStructPGridVarType(pgrid, vj);
87 if (vjtype == vitype)
88 {
89 splits[part][vi][i] = pstencil_size;
90 pstencil_size++;
91 }
92 else
93 {
94 splits[part][vi][i] = -1;
95 }
96 }
97 for (vj = 0; vj < nvars; vj++)
98 {
99 symmetric[part][vi][vj] = 0;
100 }
101 }
102 }
103
104 /* GEC0902 move the IJ creation to the initialization phase
105 * ilower = hypre_SStructGridGhstartRank(grid);
106 * iupper = ilower + hypre_SStructGridGhlocalSize(grid) - 1;
107 * HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper,
108 * &hypre_SStructMatrixIJMatrix(matrix)); */
109
110 hypre_SStructMatrixIJMatrix(matrix) = NULL;
111 hypre_SStructMatrixParCSRMatrix(matrix) = NULL;
112
113 size = 0;
114 for (part = 0; part < nparts; part++)
115 {
116 pgrid = hypre_SStructGraphPGrid(graph, part);
117 nvars = hypre_SStructPGridNVars(pgrid);
118 for (vi = 0; vi < nvars; vi++)
119 {
120 size = hypre_max(size, hypre_SStructStencilSize(stencils[part][vi]));
121 }
122 }
123 hypre_SStructMatrixEntriesSize(matrix) = size;
124 hypre_SStructMatrixSEntries(matrix) = hypre_TAlloc(int, size);
125 hypre_SStructMatrixUEntries(matrix) = hypre_TAlloc(int, size);
126 hypre_SStructMatrixTmpColCoords(matrix) = NULL;
127 hypre_SStructMatrixTmpCoeffs(matrix) = NULL;
128
129 hypre_SStructMatrixNSSymmetric(matrix) = 0;
130 hypre_SStructMatrixGlobalSize(matrix) = 0;
131 hypre_SStructMatrixRefCount(matrix) = 1;
132
133 /* GEC0902 setting the default of the object_type to HYPRE_SSTRUCT */
134
135 hypre_SStructMatrixObjectType(matrix) = HYPRE_SSTRUCT;
136
137 *matrix_ptr = matrix;
138
139 return ierr;
140}
141
142/*--------------------------------------------------------------------------
143 *--------------------------------------------------------------------------*/
144
145int
146HYPRE_SStructMatrixDestroy( HYPRE_SStructMatrix matrix )
147{
148 int ierr = 0;
149
150 hypre_SStructGraph *graph;
151 int ***splits;
152 int nparts;
153 hypre_SStructPMatrix **pmatrices;
154 int ***symmetric;
155 hypre_SStructPGrid *pgrid;
156 int nvars;
157 int part, var;
158
159 if (matrix)
160 {
161 hypre_SStructMatrixRefCount(matrix) --;
162 if (hypre_SStructMatrixRefCount(matrix) == 0)
163 {
164 graph = hypre_SStructMatrixGraph(matrix);
165 splits = hypre_SStructMatrixSplits(matrix);
166 nparts = hypre_SStructMatrixNParts(matrix);
167 pmatrices = hypre_SStructMatrixPMatrices(matrix);
168 symmetric = hypre_SStructMatrixSymmetric(matrix);
169 for (part = 0; part < nparts; part++)
170 {
171 pgrid = hypre_SStructGraphPGrid(graph, part);
172 nvars = hypre_SStructPGridNVars(pgrid);
173 for (var = 0; var < nvars; var++)
174 {
175 hypre_TFree(splits[part][var]);
176 hypre_TFree(symmetric[part][var]);
177 }
178 hypre_TFree(splits[part]);
179 hypre_TFree(symmetric[part]);
180 hypre_SStructPMatrixDestroy(pmatrices[part]);
181 }
182 HYPRE_SStructGraphDestroy(graph);
183 hypre_TFree(splits);
184 hypre_TFree(pmatrices);
185 hypre_TFree(symmetric);
186 HYPRE_IJMatrixDestroy(hypre_SStructMatrixIJMatrix(matrix));
187 hypre_TFree(hypre_SStructMatrixSEntries(matrix));
188 hypre_TFree(hypre_SStructMatrixUEntries(matrix));
189 hypre_TFree(hypre_SStructMatrixTmpColCoords(matrix));
190 hypre_TFree(hypre_SStructMatrixTmpCoeffs(matrix));
191 hypre_TFree(matrix);
192 }
193 }
194
195 return ierr;
196}
197
198/*--------------------------------------------------------------------------
199 *--------------------------------------------------------------------------*/
200
201int
202HYPRE_SStructMatrixInitialize( HYPRE_SStructMatrix matrix )
203{
204 int ierr = 0;
205
206 int nparts = hypre_SStructMatrixNParts(matrix);
207 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
208 hypre_SStructPMatrix **pmatrices = hypre_SStructMatrixPMatrices(matrix);
209 int ***symmetric = hypre_SStructMatrixSymmetric(matrix);
210 hypre_SStructStencil ***stencils = hypre_SStructGraphStencils(graph);
211 int *split;
212
213 MPI_Comm pcomm;
214 hypre_SStructPGrid *pgrid;
215 hypre_SStructStencil **pstencils;
216 int nvars;
217
218 int stencil_size;
219 hypre_Index *stencil_shape;
220 int *stencil_vars;
221 int pstencil_ndim;
222 int pstencil_size;
223
224 int part, var, i;
225
226 /* GEC0902 addition of variables for ilower and iupper */
227 MPI_Comm comm;
228 hypre_SStructGrid *grid;
229 HYPRE_BigInt ilower, iupper;
230 int matrix_type = hypre_SStructMatrixObjectType(matrix);
231
232 /* S-matrix */
233 for (part = 0; part < nparts; part++)
234 {
235 pgrid = hypre_SStructGraphPGrid(graph, part);
236 nvars = hypre_SStructPGridNVars(pgrid);
237 pstencils = hypre_TAlloc(hypre_SStructStencil *, nvars);
238 for (var = 0; var < nvars; var++)
239 {
240 split = hypre_SStructMatrixSplit(matrix, part, var);
241 stencil_size = hypre_SStructStencilSize(stencils[part][var]);
242 stencil_shape = hypre_SStructStencilShape(stencils[part][var]);
243 stencil_vars = hypre_SStructStencilVars(stencils[part][var]);
244 pstencil_ndim = hypre_SStructStencilNDim(stencils[part][var]);
245 pstencil_size = 0;
246 for (i = 0; i < stencil_size; i++)
247 {
248 if (split[i] > -1)
249 {
250 pstencil_size++;
251 }
252 }
253 HYPRE_SStructStencilCreate(pstencil_ndim, pstencil_size,
254 &pstencils[var]);
255 for (i = 0; i < stencil_size; i++)
256 {
257 if (split[i] > -1)
258 {
259 HYPRE_SStructStencilSetEntry(pstencils[var], split[i],
260 stencil_shape[i],
261 stencil_vars[i]);
262 }
263 }
264 }
265 pcomm = hypre_SStructPGridComm(pgrid);
266 hypre_SStructPMatrixCreate(pcomm, pgrid, pstencils, &pmatrices[part]);
267 for (var = 0; var < nvars; var++)
268 {
269 for (i = 0; i < nvars; i++)
270 {
271 hypre_SStructPMatrixSetSymmetric(pmatrices[part], var, i,
272 symmetric[part][var][i]);
273 }
274 }
275 hypre_SStructPMatrixInitialize(pmatrices[part]);
276 }
277
278 /* U-matrix */
279
280 /* GEC0902 knowing the kind of matrix we can create the IJMATRIX with the
281 * the right dimension (HYPRE_PARCSR without ghosts) */
282
283 grid = hypre_SStructGraphGrid(graph);
284 comm = hypre_SStructMatrixComm(matrix);
285
286 if(matrix_type == HYPRE_PARCSR)
287 {
288 ilower = hypre_SStructGridStartRank(grid);
289 iupper = ilower + (HYPRE_BigInt)(hypre_SStructGridLocalSize(grid) - 1);
290 }
291
292 if(matrix_type == HYPRE_SSTRUCT || matrix_type == HYPRE_STRUCT)
293 {
294 ilower = hypre_SStructGridGhstartRank(grid);
295 iupper = ilower + (HYPRE_BigInt)(hypre_SStructGridGhlocalSize(grid) - 1);
296 }
297
298 HYPRE_IJMatrixCreate(comm, ilower, iupper, ilower, iupper,
299 &hypre_SStructMatrixIJMatrix(matrix));
300
301
302 /* U-matrix */
303 hypre_SStructUMatrixInitialize(matrix);
304
305 return ierr;
306}
307
308/*--------------------------------------------------------------------------
309 *--------------------------------------------------------------------------*/
310
311int
312HYPRE_SStructMatrixSetValues( HYPRE_SStructMatrix matrix,
313 int part,
314 int *index,
315 int var,
316 int nentries,
317 int *entries,
318 double *values )
319{
320 int ierr = 0;
321 int ndim = hypre_SStructMatrixNDim(matrix);
322 int *Sentries;
323 int *Uentries;
324 int nSentries;
325 int nUentries;
326 hypre_SStructPMatrix *pmatrix;
327 hypre_Index cindex;
328
329 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
330 &nSentries, &Sentries,
331 &nUentries, &Uentries);
332
333 hypre_CopyToCleanIndex(index, ndim, cindex);
334
335 /* S-matrix */
336 if (nSentries > 0)
337 {
338 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
339 hypre_SStructPMatrixSetValues(pmatrix, cindex, var,
340 nSentries, Sentries, values, 0);
341 }
342 /* U-matrix */
343 if (nUentries > 0)
344 {
345 hypre_SStructUMatrixSetValues(matrix, part, cindex, var,
346 nUentries, Uentries, values, 0);
347 }
348
349 return ierr;
350}
351
352/*--------------------------------------------------------------------------
353 *--------------------------------------------------------------------------*/
354
355int
356HYPRE_SStructMatrixSetBoxValues( HYPRE_SStructMatrix matrix,
357 int part,
358 int *ilower,
359 int *iupper,
360 int var,
361 int nentries,
362 int *entries,
363 double *values )
364{
365 int ierr = 0;
366 int ndim = hypre_SStructMatrixNDim(matrix);
367 int *Sentries;
368 int *Uentries;
369 int nSentries;
370 int nUentries;
371 hypre_SStructPMatrix *pmatrix;
372 hypre_Index cilower;
373 hypre_Index ciupper;
374
375 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
376 &nSentries, &Sentries,
377 &nUentries, &Uentries);
378
379 hypre_CopyToCleanIndex(ilower, ndim, cilower);
380 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
381
382 /* S-matrix */
383 if (nSentries > 0)
384 {
385 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
386 hypre_SStructPMatrixSetBoxValues(pmatrix, cilower, ciupper, var,
387 nSentries, Sentries, values, 0);
388 }
389 /* U-matrix */
390 if (nUentries > 0)
391 {
392 hypre_SStructUMatrixSetBoxValues(matrix, part, cilower, ciupper, var,
393 nUentries, Uentries, values, 0);
394 }
395
396 return ierr;
397}
398
399/*--------------------------------------------------------------------------
400 *--------------------------------------------------------------------------*/
401int
402HYPRE_SStructMatrixGetValues( HYPRE_SStructMatrix matrix,
403 int part,
404 int *index,
405 int var,
406 int nentries,
407 int *entries,
408 double *values )
409{
410 int ierr = 0;
411 int ndim = hypre_SStructMatrixNDim(matrix);
412 int *Sentries;
413 int *Uentries;
414 int nSentries;
415 int nUentries;
416 hypre_SStructPMatrix *pmatrix;
417 hypre_Index cindex;
418
419 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
420 &nSentries, &Sentries,
421 &nUentries, &Uentries);
422
423 hypre_CopyToCleanIndex(index, ndim, cindex);
424
425 /* S-matrix */
426 if (nSentries > 0)
427 {
428 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
429 hypre_SStructPMatrixSetValues(pmatrix, cindex, var,
430 nSentries, Sentries, values, -1);
431 }
432 /* U-matrix */
433 if (nUentries > 0)
434 {
435 hypre_SStructUMatrixSetValues(matrix, part, cindex, var,
436 nUentries, Uentries, values, -1);
437 }
438
439 return ierr;
440}
441
442/*--------------------------------------------------------------------------
443 *--------------------------------------------------------------------------*/
444
445int
446HYPRE_SStructMatrixGetBoxValues( HYPRE_SStructMatrix matrix,
447 int part,
448 int *ilower,
449 int *iupper,
450 int var,
451 int nentries,
452 int *entries,
453 double *values )
454{
455 int ierr = 0;
456 int ndim = hypre_SStructMatrixNDim(matrix);
457 int *Sentries;
458 int *Uentries;
459 int nSentries;
460 int nUentries;
461 hypre_SStructPMatrix *pmatrix;
462 hypre_Index cilower;
463 hypre_Index ciupper;
464 int action;
465
466 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
467 &nSentries, &Sentries,
468 &nUentries, &Uentries);
469
470 hypre_CopyToCleanIndex(ilower, ndim, cilower);
471 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
472
473 action= -2; /* action < -1: get values */
474 /* S-matrix */
475 if (nSentries > 0)
476 {
477 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
478 hypre_SStructPMatrixSetBoxValues(pmatrix, cilower, ciupper, var,
479 nSentries, Sentries, values, action);
480 }
481 /* U-matrix */
482 if (nUentries > 0)
483 {
484 hypre_SStructUMatrixSetBoxValues(matrix, part, cilower, ciupper, var,
485 nUentries, Uentries, values, action);
486 }
487
488 return ierr;
489}
490/*--------------------------------------------------------------------------
491 *--------------------------------------------------------------------------*/
492
493int
494HYPRE_SStructMatrixAddToValues( HYPRE_SStructMatrix matrix,
495 int part,
496 int *index,
497 int var,
498 int nentries,
499 int *entries,
500 double *values )
501{
502 int ierr = 0;
503 int ndim = hypre_SStructMatrixNDim(matrix);
504 int *Sentries;
505 int *Uentries;
506 int nSentries;
507 int nUentries;
508 hypre_SStructPMatrix *pmatrix;
509 hypre_Index cindex;
510
511 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
512 &nSentries, &Sentries,
513 &nUentries, &Uentries);
514
515 hypre_CopyToCleanIndex(index, ndim, cindex);
516
517 /* S-matrix */
518 if (nSentries > 0)
519 {
520 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
521 hypre_SStructPMatrixSetValues(pmatrix, cindex, var,
522 nSentries, Sentries, values, 1);
523 }
524 /* U-matrix */
525 if (nUentries > 0)
526 {
527 hypre_SStructUMatrixSetValues(matrix, part, cindex, var,
528 nUentries, Uentries, values, 1);
529 }
530
531 return ierr;
532}
533
534/*--------------------------------------------------------------------------
535 *--------------------------------------------------------------------------*/
536
537int
538HYPRE_SStructMatrixAddToBoxValues( HYPRE_SStructMatrix matrix,
539 int part,
540 int *ilower,
541 int *iupper,
542 int var,
543 int nentries,
544 int *entries,
545 double *values )
546{
547 int ierr = 0;
548 int ndim = hypre_SStructMatrixNDim(matrix);
549 int *Sentries;
550 int *Uentries;
551 int nSentries;
552 int nUentries;
553 hypre_SStructPMatrix *pmatrix;
554 hypre_Index cilower;
555 hypre_Index ciupper;
556
557 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
558 &nSentries, &Sentries,
559 &nUentries, &Uentries);
560
561 hypre_CopyToCleanIndex(ilower, ndim, cilower);
562 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
563
564 /* S-matrix */
565 if (nSentries > 0)
566 {
567 pmatrix = hypre_SStructMatrixPMatrix(matrix, part);
568 hypre_SStructPMatrixSetBoxValues(pmatrix, cilower, ciupper, var,
569 nSentries, Sentries, values, 1);
570 }
571 /* U-matrix */
572 if (nUentries > 0)
573 {
574 hypre_SStructUMatrixSetBoxValues(matrix, part, cilower, ciupper, var,
575 nUentries, Uentries, values, 1);
576 }
577
578 return ierr;
579}
580
581/*--------------------------------------------------------------------------
582 *--------------------------------------------------------------------------*/
583
584int
585HYPRE_SStructMatrixAssemble( HYPRE_SStructMatrix matrix )
586{
587 int ierr = 0;
588
589 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
590 int nparts = hypre_SStructMatrixNParts(matrix);
591 hypre_SStructPMatrix **pmatrices = hypre_SStructMatrixPMatrices(matrix);
592 hypre_SStructGrid *grid = hypre_SStructGraphGrid(graph);
593 int **nvneighbors = hypre_SStructGridNVNeighbors(grid);
594 hypre_SStructNeighbor ***vneighbors = hypre_SStructGridVNeighbors(grid);
595
596 hypre_SStructPMatrix *pmatrix;
597 hypre_SStructStencil *stencil;
598 hypre_Index *shape;
599 int *smap;
600 int *vars;
601 hypre_StructMatrix *smatrix;
602 hypre_StructGrid *sgrid;
603 hypre_SStructNeighbor *vneighbor;
604
605 hypre_Box *box, *sbox, *ibox;
606 hypre_IndexRef offset;
607
608 int *entries;
609 int *Sentries;
610 int *Uentries;
611 int nSentries;
612 int nUentries;
613
614 double *values = NULL;
615
616 int nvars, nentries;
617 int part, var, entry, sentry, b, sb;
618
619 /*------------------------------------------------------
620 * Move off-part couplings (described by neighbor info)
621 * from S-matrix structure into U-matrix structure.
622 *------------------------------------------------------*/
623
624 box = hypre_BoxCreate();
625 ibox = hypre_BoxCreate();
626
627 nentries = hypre_SStructMatrixEntriesSize(matrix);
628 entries = hypre_TAlloc(int, nentries);
629 for (entry = 0; entry < nentries; entry++)
630 {
631 entries[entry] = entry;
632 }
633
634 for (part = 0; part < nparts; part++)
635 {
636 pmatrix = pmatrices[part];
637
638 nvars = hypre_SStructPMatrixNVars(pmatrix);
639 for (var = 0; var < nvars; var++)
640 {
641 stencil = hypre_SStructPMatrixStencil(pmatrix, var);
642 smap = hypre_SStructPMatrixSMap(pmatrix, var);
643 shape = hypre_SStructStencilShape(stencil);
644 vars = hypre_SStructStencilVars(stencil);
645 nentries = hypre_SStructStencilSize(stencil);
646
647 hypre_SStructMatrixSplitEntries(matrix, part, var, nentries, entries,
648 &nSentries, &Sentries,
649 &nUentries, &Uentries);
650
651 for (entry = 0; entry < nSentries; entry++)
652 {
653 smatrix = hypre_SStructPMatrixSMatrix(pmatrix, var,
654 vars[entries[entry]]);
655 sentry = smap[entries[entry]];
656
657 /* Shift/intersect neighbor box and move values */
658 for (b = 0; b < nvneighbors[part][var]; b++)
659 {
660 vneighbor = &vneighbors[part][var][b];
661 hypre_CopyBox(hypre_SStructNeighborBox(vneighbor), box);
662
663 /* shift box by stencil offset */
664 offset = shape[entry];
665 hypre_BoxIMinX(box) -= hypre_IndexX(offset);
666 hypre_BoxIMinY(box) -= hypre_IndexY(offset);
667 hypre_BoxIMinZ(box) -= hypre_IndexZ(offset);
668 hypre_BoxIMaxX(box) -= hypre_IndexX(offset);
669 hypre_BoxIMaxY(box) -= hypre_IndexY(offset);
670 hypre_BoxIMaxZ(box) -= hypre_IndexZ(offset);
671
672 sgrid = hypre_StructMatrixGrid(smatrix);
673 hypre_ForStructGridBoxI(sb, sgrid)
674 {
675 sbox = hypre_StructGridBox(sgrid, sb);
676 hypre_IntersectBoxes(box, sbox, ibox);
677
678 if (hypre_BoxVolume(ibox))
679 {
680 values = hypre_TReAlloc(values, double,
681 hypre_BoxVolume(ibox));
682
683 /* move matrix values from S-matrix to U-matrix */
684 hypre_StructMatrixSetBoxValues(smatrix, ibox,
685 1, &sentry, values, -1);
686
687 hypre_SStructUMatrixSetBoxValues(matrix, part,
688 hypre_BoxIMin(ibox),
689 hypre_BoxIMax(ibox),
690 var, 1, &entry,
691 values, 1);
692 }
693 }
694 }
695 }
696 }
697 }
698
699 hypre_TFree(entries);
700 hypre_TFree(values);
701 hypre_BoxDestroy(box);
702 hypre_BoxDestroy(ibox);
703
704 for (part = 0; part < nparts; part++)
705 {
706 hypre_SStructPMatrixAssemble(pmatrices[part]);
707 }
708
709 /* U-matrix */
710 hypre_SStructUMatrixAssemble(matrix);
711
712 return ierr;
713}
714
715/*--------------------------------------------------------------------------
716 * NOTE: Should set things up so that this information can be passed
717 * immediately to the PMatrix. Unfortunately, the PMatrix is
718 * currently not created until the SStructMatrix is initialized.
719 *--------------------------------------------------------------------------*/
720
721int
722HYPRE_SStructMatrixSetSymmetric( HYPRE_SStructMatrix matrix,
723 int part,
724 int var,
725 int to_var,
726 int symmetric )
727{
728 int ierr = 0;
729
730 int ***msymmetric = hypre_SStructMatrixSymmetric(matrix);
731 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
732 hypre_SStructPGrid *pgrid;
733
734 int pstart = part;
735 int psize = 1;
736 int vstart = var;
737 int vsize = 1;
738 int tstart = to_var;
739 int tsize = 1;
740 int p, v, t;
741
742 if (part == -1)
743 {
744 pstart = 0;
745 psize = hypre_SStructMatrixNParts(matrix);
746 }
747
748 for (p = pstart; p < psize; p++)
749 {
750 pgrid = hypre_SStructGraphPGrid(graph, p);
751 if (var == -1)
752 {
753 vstart = 0;
754 vsize = hypre_SStructPGridNVars(pgrid);
755 }
756 if (to_var == -1)
757 {
758 tstart = 0;
759 tsize = hypre_SStructPGridNVars(pgrid);
760 }
761
762 for (v = vstart; v < vsize; v++)
763 {
764 for (t = tstart; t < tsize; t++)
765 {
766 msymmetric[p][v][t] = symmetric;
767 }
768 }
769 }
770
771 return ierr;
772}
773
774/*--------------------------------------------------------------------------
775 *--------------------------------------------------------------------------*/
776
777int
778HYPRE_SStructMatrixSetNSSymmetric( HYPRE_SStructMatrix matrix,
779 int symmetric )
780{
781 int ierr = 0;
782
783 hypre_SStructMatrixNSSymmetric(matrix) = symmetric;
784
785 return ierr;
786}
787
788/*--------------------------------------------------------------------------
789 *--------------------------------------------------------------------------*/
790
791int
792HYPRE_SStructMatrixSetObjectType( HYPRE_SStructMatrix matrix,
793 int type )
794{
795 int ierr = 0;
796
797 hypre_SStructGraph *graph = hypre_SStructMatrixGraph(matrix);
798 int ***splits = hypre_SStructMatrixSplits(matrix);
799 int nparts = hypre_SStructMatrixNParts(matrix);
800 hypre_SStructStencil ***stencils = hypre_SStructGraphStencils(graph);
801
802 hypre_SStructPGrid *pgrid;
803 int nvars;
804 int stencil_size;
805 int part, var, i;
806
807 hypre_SStructMatrixObjectType(matrix) = type ;
808
809 /* RDF: This and all other modifications to 'split' really belong
810 * in the Initialize routine */
811 if (type != HYPRE_SSTRUCT && type != HYPRE_STRUCT)
812 {
813 for (part = 0; part < nparts; part++)
814 {
815 pgrid = hypre_SStructGraphPGrid(graph, part);
816 nvars = hypre_SStructPGridNVars(pgrid);
817 for (var = 0; var < nvars; var++)
818 {
819 stencil_size = hypre_SStructStencilSize(stencils[part][var]);
820 for (i = 0; i < stencil_size; i++)
821 {
822 splits[part][var][i] = -1;
823 }
824 }
825 }
826 }
827
828 return ierr;
829}
830
831/*--------------------------------------------------------------------------
832 *--------------------------------------------------------------------------*/
833
834int
835HYPRE_SStructMatrixGetObject( HYPRE_SStructMatrix matrix,
836 void **object )
837{
838 int ierr = 0;
839
840 int type = hypre_SStructMatrixObjectType(matrix);
841 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
842 hypre_SStructPMatrix *pA;
843 hypre_StructMatrix *sA;
844 int part, var;
845
846
847 if (type == HYPRE_PARCSR)
848 {
849 HYPRE_IJMatrixGetObject(ijmatrix, object);
850 }
851
852 else if (type == HYPRE_SSTRUCT)
853 {
854 *object= matrix;
855 }
856
857 else if (type == HYPRE_STRUCT)
858 {
859 /* only one part & one variable */
860 part= 0;
861 pA = hypre_SStructMatrixPMatrix(matrix, part);
862 var = 0;
863 sA = hypre_SStructPMatrixSMatrix(pA, var, var);
864 *object= sA;
865 }
866
867 return ierr;
868}
869
870/*--------------------------------------------------------------------------
871 *--------------------------------------------------------------------------*/
872
873int
874HYPRE_SStructMatrixGetObject2( HYPRE_SStructMatrix matrix,
875 void **object )
876{
877 int ierr = 0;
878
879 int type = hypre_SStructMatrixObjectType(matrix);
880 HYPRE_IJMatrix ijmatrix = hypre_SStructMatrixIJMatrix(matrix);
881 hypre_SStructPMatrix *pA;
882 hypre_StructMatrix *sA;
883 int part, var;
884
885
886 if (type == HYPRE_PARCSR)
887 {
888 /* only difference from ..GetObject: here returns an IJMatrix, not a ParCSRMatrix */
889 *object = ijmatrix;
890 }
891
892 else if (type == HYPRE_SSTRUCT)
893 {
894 *object= matrix;
895 }
896
897 else if (type == HYPRE_STRUCT)
898 {
899 /* only one part & one variable */
900 part= 0;
901 pA = hypre_SStructMatrixPMatrix(matrix, part);
902 var = 0;
903 sA = hypre_SStructPMatrixSMatrix(pA, var, var);
904 *object= sA;
905 }
906
907 return ierr;
908}
909
910/*--------------------------------------------------------------------------
911 *--------------------------------------------------------------------------*/
912
913int
914HYPRE_SStructMatrixPrint( const char *filename,
915 HYPRE_SStructMatrix matrix,
916 int all )
917{
918 int ierr = 0;
919 int nparts = hypre_SStructMatrixNParts(matrix);
920 int part;
921 char new_filename[255];
922
923 for (part = 0; part < nparts; part++)
924 {
925 sprintf(new_filename, "%s.%02d", filename, part);
926 hypre_SStructPMatrixPrint(new_filename,
927 hypre_SStructMatrixPMatrix(matrix, part),
928 all);
929 }
930
931 /* U-matrix */
932 sprintf(new_filename, "%s.UMatrix", filename);
933 HYPRE_IJMatrixPrint(hypre_SStructMatrixIJMatrix(matrix), new_filename);
934
935 return ierr;
936}
937
938/*--------------------------------------------------------------------------
939 * HYPRE_StructMatrixMatvec
940 *--------------------------------------------------------------------------*/
941
942int
943HYPRE_SStructMatrixMatvec( double alpha,
944 HYPRE_SStructMatrix A,
945 HYPRE_SStructVector x,
946 double beta,
947 HYPRE_SStructVector y )
948{
949 return ( hypre_SStructMatvec( alpha, (hypre_SStructMatrix *) A,
950 (hypre_SStructVector *) x, beta,
951 (hypre_SStructVector *) y) );
952}
953
Note: See TracBrowser for help on using the repository browser.