source: CIVL/examples/mpi-omp/AMG2013/sstruct_mv/HYPRE_sstruct_vector.c

main
Last change on this file 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: 22.4 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_SStructVector interface
18 *
19 *****************************************************************************/
20
21#include "headers.h"
22#include "sstruct_mv.h"
23
24/*--------------------------------------------------------------------------
25 *--------------------------------------------------------------------------*/
26
27int
28HYPRE_SStructVectorCreate( MPI_Comm comm,
29 HYPRE_SStructGrid grid,
30 HYPRE_SStructVector *vector_ptr )
31{
32 int ierr = 0;
33
34 hypre_SStructVector *vector;
35 int nparts;
36 hypre_SStructPVector **pvectors;
37 MPI_Comm pcomm;
38 hypre_SStructPGrid *pgrid;
39 int part;
40
41 vector = hypre_TAlloc(hypre_SStructVector, 1);
42
43 hypre_SStructVectorComm(vector) = comm;
44 hypre_SStructVectorNDim(vector) = hypre_SStructGridNDim(grid);
45 hypre_SStructGridRef(grid, &hypre_SStructVectorGrid(vector));
46 hypre_SStructVectorObjectType(vector) = HYPRE_SSTRUCT;
47 nparts = hypre_SStructGridNParts(grid);
48 hypre_SStructVectorNParts(vector) = nparts;
49 pvectors = hypre_TAlloc(hypre_SStructPVector *, nparts);
50 for (part = 0; part < nparts; part++)
51 {
52 pgrid = hypre_SStructGridPGrid(grid, part);
53 pcomm = hypre_SStructPGridComm(pgrid);
54 ierr = hypre_SStructPVectorCreate(pcomm, pgrid, &pvectors[part]);
55 }
56 hypre_SStructVectorPVectors(vector) = pvectors;
57 hypre_SStructVectorIJVector(vector) = NULL;
58
59 /* GEC1002 initializing to NULL */
60
61 hypre_SStructVectorDataIndices(vector) = NULL;
62 hypre_SStructVectorData(vector) = NULL;
63
64 /* GEC1002 moving the creation of the ijvector the the initialize part
65 * ilower = hypre_SStructGridStartRank(grid);
66 * iupper = ilower + hypre_SStructGridLocalSize(grid) - 1;
67 * HYPRE_IJVectorCreate(comm, ilowergh, iuppergh,
68 * &hypre_SStructVectorIJVector(vector)); */
69
70 hypre_SStructVectorIJVector(vector) = NULL;
71 hypre_SStructVectorParVector(vector) = NULL;
72 hypre_SStructVectorGlobalSize(vector) = 0;
73 hypre_SStructVectorRefCount(vector) = 1;
74 hypre_SStructVectorDataSize(vector) = 0;
75 hypre_SStructVectorObjectType(vector) = HYPRE_SSTRUCT;
76
77 *vector_ptr = vector;
78
79 return ierr;
80}
81
82/*--------------------------------------------------------------------------
83 *--------------------------------------------------------------------------*/
84
85int
86HYPRE_SStructVectorDestroy( HYPRE_SStructVector vector )
87{
88 int ierr = 0;
89
90 int nparts;
91 hypre_SStructPVector **pvectors;
92 int part;
93 int vector_type = hypre_SStructVectorObjectType(vector);
94
95 /* GEC1002 destroying dataindices and data in vector */
96
97 if (vector)
98 {
99 hypre_SStructVectorRefCount(vector) --;
100 if (hypre_SStructVectorRefCount(vector) == 0)
101 {
102 HYPRE_SStructGridDestroy(hypre_SStructVectorGrid(vector));
103 nparts = hypre_SStructVectorNParts(vector);
104 pvectors = hypre_SStructVectorPVectors(vector);
105 for (part = 0; part < nparts; part++)
106 {
107 hypre_SStructPVectorDestroy(pvectors[part]);
108 }
109 hypre_TFree(pvectors);
110 HYPRE_IJVectorDestroy(hypre_SStructVectorIJVector(vector));
111 /* GEC1002 the ijdestroy takes care of the data when the
112 * vector is type HYPRE_SSTRUCT. This is a result that the
113 * ijvector does not use the owndata flag in the data structure
114 * unlike the structvector */
115
116 /* GEC if data has been allocated then free the pointer */
117 hypre_TFree(hypre_SStructVectorDataIndices(vector));
118
119 if (hypre_SStructVectorData(vector) && (vector_type == HYPRE_PARCSR))
120 {
121 hypre_TFree(hypre_SStructVectorData(vector));
122 }
123
124 hypre_TFree(vector);
125 }
126 }
127
128 return ierr;
129}
130
131/*---------------------------------------------------------
132 * GEC1002 changes to initialize the vector with a data chunk
133 * that includes all the part,var pieces instead of just svector-var
134 * pieces. In case of pure unstruct-variables (ucvar), which are at the
135 * end of each part, we might need to modify initialize shell vector
136 * ----------------------------------------------------------*/
137
138int
139HYPRE_SStructVectorInitialize( HYPRE_SStructVector vector )
140{
141 int ierr = 0;
142 int datasize;
143 int nvars ;
144 int nparts = hypre_SStructVectorNParts(vector) ;
145 int var,part ;
146 double *data ;
147 double *pdata ;
148 double *sdata ;
149 hypre_SStructPVector *pvector;
150 hypre_StructVector *svector;
151 int *dataindices;
152 int *pdataindices;
153 int vector_type = hypre_SStructVectorObjectType(vector);
154 hypre_SStructGrid *grid = hypre_SStructVectorGrid(vector);
155 MPI_Comm comm = hypre_SStructVectorComm(vector);
156 HYPRE_IJVector ijvector;
157
158 /* GEC0902 addition of variables for ilower and iupper */
159 HYPRE_BigInt ilower, iupper;
160 hypre_ParVector *par_vector;
161 hypre_Vector *parlocal_vector;
162
163
164 /* GEC0902 getting the datasizes and indices we need */
165
166 hypre_SStructVectorInitializeShell(vector);
167
168 datasize = hypre_SStructVectorDataSize(vector);
169
170 data = hypre_CTAlloc(double, datasize);
171
172 dataindices = hypre_SStructVectorDataIndices(vector);
173
174 hypre_SStructVectorData(vector) = data;
175
176 for (part = 0; part < nparts; part++)
177 {
178 pvector = hypre_SStructVectorPVector(vector,part);
179 pdataindices = hypre_SStructPVectorDataIndices(pvector);
180 /* shift-num = dataindices[part]; */
181 pdata = data + dataindices[part];
182 nvars = hypre_SStructPVectorNVars(pvector);
183
184 for (var = 0; var < nvars; var++)
185 {
186 svector = hypre_SStructPVectorSVector(pvector, var);
187 /* shift-pnum = pdataindices[var]; */
188 sdata = pdata + pdataindices[var];
189
190 /* GEC1002 initialization of inside data pointer of a svector
191 * because no data is alloced, we make sure the flag is zero. This
192 * affects the destroy */
193 hypre_StructVectorInitializeData(svector, sdata);
194 hypre_StructVectorDataAlloced(svector) = 0;
195 }
196
197 }
198
199 /* GEC1002 this is now the creation of the ijmatrix and the initialization
200 * by checking the type of the vector */
201
202 if(vector_type == HYPRE_PARCSR )
203 {
204 ilower = hypre_SStructGridStartRank(grid);
205 iupper = ilower + (HYPRE_BigInt)(hypre_SStructGridLocalSize(grid) - 1);
206 }
207
208 if(vector_type == HYPRE_SSTRUCT || vector_type == HYPRE_STRUCT)
209 {
210 ilower = hypre_SStructGridGhstartRank(grid);
211 iupper = ilower + (HYPRE_BigInt)(hypre_SStructGridGhlocalSize(grid) - 1);
212 }
213
214 HYPRE_IJVectorCreate(comm, ilower, iupper,
215 &hypre_SStructVectorIJVector(vector));
216
217 /* GEC1002, once the partitioning is done, it is time for the actual
218 * initialization */
219
220
221 /* u-vector: the type is for the parvector inside the ijvector */
222
223 ijvector = hypre_SStructVectorIJVector(vector);
224
225 ierr = HYPRE_IJVectorSetObjectType(ijvector, HYPRE_PARCSR);
226
227 ierr += HYPRE_IJVectorInitialize(ijvector);
228
229
230 /* GEC1002 for HYPRE_SSTRUCT type of vector, we do not need data allocated inside
231 * the parvector piece of the structure. We make that pointer within the localvector
232 * to point to the outside "data". Before redirecting the
233 * local pointer to point to the true data chunk for HYPRE_SSTRUCT: we destroy and assign.
234 * We now have two entries of the data structure pointing to the same chunk
235 * if we have a HYPRE_SSTRUCT vector
236 * We do not need the IJVectorInitializePar, we have to undoit for the SStruct case
237 * in a sense it is a desinitializepar */
238
239 if (vector_type == HYPRE_SSTRUCT || vector_type == HYPRE_STRUCT)
240 {
241 par_vector = hypre_IJVectorObject(ijvector);
242 parlocal_vector = hypre_ParVectorLocalVector(par_vector);
243 hypre_TFree(hypre_VectorData(parlocal_vector));
244 hypre_VectorData(parlocal_vector) = data ;
245 }
246
247 return ierr;
248}
249
250/*--------------------------------------------------------------------------
251 *--------------------------------------------------------------------------*/
252
253int
254HYPRE_SStructVectorSetValues( HYPRE_SStructVector vector,
255 int part,
256 int *index,
257 int var,
258 double *value )
259{
260 int ierr = 0;
261 int ndim = hypre_SStructVectorNDim(vector);
262 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
263 hypre_Index cindex;
264
265 hypre_CopyToCleanIndex(index, ndim, cindex);
266
267 if (var < hypre_SStructPVectorNVars(pvector))
268 {
269 ierr = hypre_SStructPVectorSetValues(pvector, cindex, var, value, 0);
270 }
271 else
272 {
273 /* TODO */
274 }
275
276 return ierr;
277}
278
279/*--------------------------------------------------------------------------
280 *--------------------------------------------------------------------------*/
281
282int
283HYPRE_SStructVectorSetBoxValues( HYPRE_SStructVector vector,
284 int part,
285 int *ilower,
286 int *iupper,
287 int var,
288 double *values )
289{
290 int ierr = 0;
291 int ndim = hypre_SStructVectorNDim(vector);
292 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
293 hypre_Index cilower;
294 hypre_Index ciupper;
295
296 hypre_CopyToCleanIndex(ilower, ndim, cilower);
297 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
298
299 ierr = hypre_SStructPVectorSetBoxValues(pvector, cilower, ciupper,
300 var, values, 0);
301
302 return ierr;
303}
304
305/*--------------------------------------------------------------------------
306 *--------------------------------------------------------------------------*/
307
308int
309HYPRE_SStructVectorAddToValues( HYPRE_SStructVector vector,
310 int part,
311 int *index,
312 int var,
313 double *value )
314{
315 int ierr = 0;
316 int ndim = hypre_SStructVectorNDim(vector);
317 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
318 hypre_Index cindex;
319
320 hypre_CopyToCleanIndex(index, ndim, cindex);
321
322 if (var < hypre_SStructPVectorNVars(pvector))
323 {
324 ierr = hypre_SStructPVectorSetValues(pvector, cindex, var, value, 1);
325 }
326 else
327 {
328 /* TODO */
329 }
330
331 return ierr;
332}
333
334/*--------------------------------------------------------------------------
335 *--------------------------------------------------------------------------*/
336
337int
338HYPRE_SStructVectorAddToBoxValues( HYPRE_SStructVector vector,
339 int part,
340 int *ilower,
341 int *iupper,
342 int var,
343 double *values )
344{
345 int ierr = 0;
346 int ndim = hypre_SStructVectorNDim(vector);
347 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
348 hypre_Index cilower;
349 hypre_Index ciupper;
350
351 hypre_CopyToCleanIndex(ilower, ndim, cilower);
352 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
353
354 ierr = hypre_SStructPVectorSetBoxValues(pvector, cilower, ciupper,
355 var, values, 1);
356
357 return ierr;
358}
359
360/*--------------------------------------------------------------------------
361 *--------------------------------------------------------------------------*/
362
363int
364HYPRE_SStructVectorAssemble( HYPRE_SStructVector vector )
365{
366 int ierr = 0;
367 int nparts = hypre_SStructVectorNParts(vector);
368 HYPRE_IJVector ijvector = hypre_SStructVectorIJVector(vector);
369 int part;
370
371 for (part = 0; part < nparts; part++)
372 {
373 hypre_SStructPVectorAssemble(hypre_SStructVectorPVector(vector, part));
374 }
375
376#if 0
377 /* ZTODO: construct comm_pkg for communications between parts */
378 hypre_CommPkgDestroy(comm_pkgs[var]);
379 comm_pkgs[var] =
380 hypre_CommPkgCreate(send_boxes, recv_boxes,
381 unit_stride, unit_stride,
382 hypre_StructVectorDataSpace(svectors[var]),
383 hypre_StructVectorDataSpace(svectors[var]),
384 send_processes, recv_processes, 1,
385 hypre_StructVectorComm(svectors[var]),
386 hypre_StructGridPeriodic(sgrid));
387#endif
388
389 /* u-vector */
390 ierr = HYPRE_IJVectorAssemble(ijvector);
391
392 HYPRE_IJVectorGetObject(ijvector,
393 (void **) &hypre_SStructVectorParVector(vector));
394
395 /* if the object type is parcsr, then convert the sstruct vector which has ghost
396 layers to a parcsr vector without ghostlayers. */
397 if (hypre_SStructVectorObjectType(vector) == HYPRE_PARCSR)
398 {
399 hypre_SStructVectorParConvert(vector,
400 &hypre_SStructVectorParVector(vector));
401 }
402
403 return ierr;
404}
405
406/*--------------------------------------------------------------------------
407 *--------------------------------------------------------------------------*/
408
409int
410HYPRE_SStructVectorGather( HYPRE_SStructVector vector )
411{
412 int ierr = 0;
413 int nparts = hypre_SStructVectorNParts(vector);
414 int part;
415
416 /* GEC1102 we change the name of the restore-->parrestore */
417
418 if (hypre_SStructVectorObjectType(vector) == HYPRE_PARCSR)
419 {
420 hypre_SStructVectorParRestore(vector, hypre_SStructVectorParVector(vector));
421 }
422
423 for (part = 0; part < nparts; part++)
424 {
425 hypre_SStructPVectorGather(hypre_SStructVectorPVector(vector, part));
426 }
427
428#if 0
429 /* ZTODO: gather data from other parts */
430 {
431 int nvars = hypre_SStructPVectorNVars(pvector);
432 hypre_StructVector **svectors = hypre_SStructPVectorSVectors(pvector);
433 hypre_CommPkg **comm_pkgs = hypre_SStructPVectorCommPkgs(pvector);
434 hypre_CommHandle *comm_handle;
435 int var;
436
437 for (var = 0; var < nvars; var++)
438 {
439 if (comm_pkgs[var] != NULL)
440 {
441 hypre_InitializeCommunication(comm_pkgs[var],
442 hypre_StructVectorData(svectors[var]),
443 hypre_StructVectorData(svectors[var]),
444 &comm_handle);
445 hypre_FinalizeCommunication(comm_handle);
446 }
447 }
448 }
449#endif
450
451 return ierr;
452}
453
454/*--------------------------------------------------------------------------
455 *--------------------------------------------------------------------------*/
456
457int
458HYPRE_SStructVectorGetValues( HYPRE_SStructVector vector,
459 int part,
460 int *index,
461 int var,
462 double *value )
463{
464 int ierr = 0;
465 int ndim = hypre_SStructVectorNDim(vector);
466 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
467 hypre_Index cindex;
468
469 hypre_CopyToCleanIndex(index, ndim, cindex);
470
471 if (var < hypre_SStructPVectorNVars(pvector))
472 {
473 ierr = hypre_SStructPVectorGetValues(pvector, cindex, var, value);
474 }
475 else
476 {
477 /* TODO */
478 }
479
480 return ierr;
481}
482
483/*--------------------------------------------------------------------------
484 *--------------------------------------------------------------------------*/
485
486int
487HYPRE_SStructVectorGetBoxValues(HYPRE_SStructVector vector,
488 int part,
489 int *ilower,
490 int *iupper,
491 int var,
492 double *values )
493{
494 int ierr = 0;
495 int ndim = hypre_SStructVectorNDim(vector);
496 hypre_SStructPVector *pvector = hypre_SStructVectorPVector(vector, part);
497 hypre_Index cilower;
498 hypre_Index ciupper;
499
500 hypre_CopyToCleanIndex(ilower, ndim, cilower);
501 hypre_CopyToCleanIndex(iupper, ndim, ciupper);
502
503 ierr = hypre_SStructPVectorGetBoxValues(pvector, cilower, ciupper,
504 var, values);
505
506 return ierr;
507}
508
509/*--------------------------------------------------------------------------
510 *--------------------------------------------------------------------------*/
511
512int
513HYPRE_SStructVectorSetConstantValues( HYPRE_SStructVector vector,
514 double value )
515{
516 int ierr = 0;
517 hypre_SStructPVector *pvector;
518 int part;
519 int nparts = hypre_SStructVectorNParts(vector);
520
521 for ( part = 0; part < nparts; part++ )
522 {
523 pvector = hypre_SStructVectorPVector( vector, part );
524 ierr += hypre_SStructPVectorSetConstantValues( pvector, value );
525 };
526
527 return ierr;
528}
529
530/*--------------------------------------------------------------------------
531 *--------------------------------------------------------------------------*/
532
533int
534HYPRE_SStructVectorSetObjectType( HYPRE_SStructVector vector,
535 int type )
536{
537 int ierr = 0;
538
539 /* this implements only HYPRE_PARCSR, which is always available */
540 hypre_SStructVectorObjectType(vector) = type;
541
542 return ierr;
543}
544
545/*--------------------------------------------------------------------------
546 *--------------------------------------------------------------------------*/
547
548int
549HYPRE_SStructVectorGetObject( HYPRE_SStructVector vector,
550 void **object )
551{
552 int ierr = 0;
553 int vector_type = hypre_SStructVectorObjectType(vector);
554 hypre_SStructPVector *pvector;
555 hypre_StructVector *svector;
556
557 int part, var;
558
559 /* GEC1102 in case the vector is HYPRE_SSTRUCT */
560
561 if (vector_type == HYPRE_SSTRUCT)
562 {
563 hypre_SStructVectorConvert(vector, (hypre_ParVector **) object);
564 }
565 else if (vector_type == HYPRE_PARCSR)
566 {
567 *object= hypre_SStructVectorParVector(vector);
568 }
569 else if (vector_type == HYPRE_STRUCT)
570 {
571 /* only one part & one variable */
572 part= 0;
573 var = 0;
574 pvector= hypre_SStructVectorPVector(vector, part);
575 svector= hypre_SStructPVectorSVector(pvector, var);
576 *object = svector;
577 }
578
579 return ierr;
580}
581
582/*--------------------------------------------------------------------------
583 *--------------------------------------------------------------------------*/
584
585int
586HYPRE_SStructVectorPrint( const char *filename,
587 HYPRE_SStructVector vector,
588 int all )
589{
590 int ierr = 0;
591 int nparts = hypre_SStructVectorNParts(vector);
592 int part;
593 char new_filename[255];
594
595 for (part = 0; part < nparts; part++)
596 {
597 sprintf(new_filename, "%s.%02d", filename, part);
598 hypre_SStructPVectorPrint(new_filename,
599 hypre_SStructVectorPVector(vector, part),
600 all);
601 }
602
603 return ierr;
604}
605
606/******************************************************************************
607 * Zero out the ghostlayer values of a vector. Needed when a vector is
608 * re-assembled with new values and addtovalues from off_procs are triggered.
609 *****************************************************************************/
610int
611HYPRE_SStructVectorClearGhostValues(HYPRE_SStructVector x)
612{
613 return hypre_SStructVectorClearGhostValues((hypre_SStructVector *)x);
614}
615
616
617/******************************************************************************
618 * copy x to y, y should already exist and be the same size
619 *****************************************************************************/
620int
621HYPRE_SStructVectorCopy( HYPRE_SStructVector x,
622 HYPRE_SStructVector y )
623{
624 return hypre_SStructCopy( (hypre_SStructVector *)x,
625 (hypre_SStructVector *)y );
626}
627
628/******************************************************************************
629 * y = a*y, for vector y and scalar a
630 *****************************************************************************/
631int
632HYPRE_SStructVectorScale( double alpha, HYPRE_SStructVector y )
633{
634 return hypre_SStructScale( alpha, (hypre_SStructVector *)y );
635}
636
637/******************************************************************************
638 * inner or dot product, result = < x, y >
639 *****************************************************************************/
640int
641HYPRE_SStructInnerProd( HYPRE_SStructVector x,
642 HYPRE_SStructVector y,
643 double *result )
644{
645 return hypre_SStructInnerProd( (hypre_SStructVector *)x,
646 (hypre_SStructVector *)y,
647 result );
648}
649
650/******************************************************************************
651 * y = y + alpha*x for vectors y, x and scalar alpha
652 *****************************************************************************/
653int
654HYPRE_SStructAxpy( double alpha,
655 HYPRE_SStructVector x,
656 HYPRE_SStructVector y )
657{
658 return hypre_SStructAxpy( alpha,
659 (hypre_SStructVector *)x,
660 (hypre_SStructVector *)y );
661}
Note: See TracBrowser for help on using the repository browser.