Actual source code: vhyp.c


  2: /*
  3:     Creates hypre ijvector from PETSc vector
  4: */

  6: #include <petsc/private/vecimpl.h>
  7: #include <../src/vec/vec/impls/hypre/vhyp.h>
  8: #include <HYPRE.h>

 10: PetscErrorCode VecHYPRE_IJVectorCreate(PetscLayout map, VecHYPRE_IJVector *ij)
 11: {
 12:   VecHYPRE_IJVector nij;

 14:   PetscFunctionBegin;
 15:   PetscCall(PetscNew(&nij));
 16:   PetscCall(PetscLayoutSetUp(map));
 17:   PetscCallExternal(HYPRE_IJVectorCreate, map->comm, map->rstart, map->rend - 1, &nij->ij);
 18:   PetscCallExternal(HYPRE_IJVectorSetObjectType, nij->ij, HYPRE_PARCSR);
 19: #if defined(PETSC_HAVE_HYPRE_DEVICE)
 20:   PetscCallExternal(HYPRE_IJVectorInitialize_v2, nij->ij, HYPRE_MEMORY_DEVICE);
 21: #else
 22:   PetscCallExternal(HYPRE_IJVectorInitialize, nij->ij);
 23: #endif
 24:   PetscCallExternal(HYPRE_IJVectorAssemble, nij->ij);
 25:   *ij = nij;
 26:   PetscFunctionReturn(PETSC_SUCCESS);
 27: }

 29: PetscErrorCode VecHYPRE_IJVectorDestroy(VecHYPRE_IJVector *ij)
 30: {
 31:   PetscFunctionBegin;
 32:   if (!*ij) PetscFunctionReturn(PETSC_SUCCESS);
 33:   PetscCheck(!(*ij)->pvec, PetscObjectComm((PetscObject)((*ij)->pvec)), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
 34:   PetscCallExternal(HYPRE_IJVectorDestroy, (*ij)->ij);
 35:   PetscCall(PetscFree(*ij));
 36:   PetscFunctionReturn(PETSC_SUCCESS);
 37: }

 39: PetscErrorCode VecHYPRE_IJVectorCopy(Vec v, VecHYPRE_IJVector ij)
 40: {
 41:   const PetscScalar *array;

 43:   PetscFunctionBegin;
 44: #if defined(PETSC_HAVE_HYPRE_DEVICE)
 45:   PetscCallExternal(HYPRE_IJVectorInitialize_v2, ij->ij, HYPRE_MEMORY_DEVICE);
 46: #else
 47:   PetscCallExternal(HYPRE_IJVectorInitialize, ij->ij);
 48: #endif
 49:   PetscCall(VecGetArrayRead(v, &array));
 50:   PetscCallExternal(HYPRE_IJVectorSetValues, ij->ij, v->map->n, NULL, (HYPRE_Complex *)array);
 51:   PetscCall(VecRestoreArrayRead(v, &array));
 52:   PetscCallExternal(HYPRE_IJVectorAssemble, ij->ij);
 53:   PetscFunctionReturn(PETSC_SUCCESS);
 54: }

 56: /*
 57:     Replaces the address where the HYPRE vector points to its data with the address of
 58:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
 59:   Allows use to get the data into a HYPRE vector without the cost of memcopies
 60: */
 61: #define VecHYPRE_ParVectorReplacePointer(b, newvalue, savedvalue) \
 62:   { \
 63:     hypre_ParVector *par_vector   = (hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector *)b)); \
 64:     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
 65:     savedvalue                    = local_vector->data; \
 66:     local_vector->data            = newvalue; \
 67:   }

 69: /*
 70:   This routine access the pointer to the raw data of the "v" to be passed to HYPRE
 71:    - rw values indicate the type of access : 0 -> read, 1 -> write, 2 -> read-write
 72:    - hmem is the location HYPRE is expecting
 73:    - the function returns a pointer to the data (ptr) and the corresponding restore
 74:   Could be extended to VECKOKKOS if we had a way to access the raw pointer to device data.
 75: */
 76: static inline PetscErrorCode VecGetArrayForHYPRE(Vec v, int rw, HYPRE_MemoryLocation hmem, PetscScalar **ptr, PetscErrorCode (**res)(Vec, PetscScalar **))
 77: {
 78:   PetscMemType mtype;
 79:   MPI_Comm     comm;

 81:   PetscFunctionBegin;
 82: #if !defined(PETSC_HAVE_HYPRE_DEVICE)
 83:   hmem = HYPRE_MEMORY_HOST; /* this is just a convenience because HYPRE_MEMORY_HOST and HYPRE_MEMORY_DEVICE are the same in this case */
 84: #endif
 85:   *ptr = NULL;
 86:   *res = NULL;
 87:   PetscCall(PetscObjectGetComm((PetscObject)v, &comm));
 88:   switch (rw) {
 89:   case 0: /* read */
 90:     if (hmem == HYPRE_MEMORY_HOST) {
 91:       PetscCall(VecGetArrayRead(v, (const PetscScalar **)ptr));
 92:       *res = (PetscErrorCode(*)(Vec, PetscScalar **))VecRestoreArrayRead;
 93:     } else {
 94:       PetscCall(VecGetArrayReadAndMemType(v, (const PetscScalar **)ptr, &mtype));
 95:       PetscCheck(PetscMemTypeDevice(mtype), comm, PETSC_ERR_ARG_WRONG, "HYPRE_MEMORY_DEVICE expects a device vector. You need to enable PETSc device support, for example, in some cases, -vec_type cuda");
 96:       *res = (PetscErrorCode(*)(Vec, PetscScalar **))VecRestoreArrayReadAndMemType;
 97:     }
 98:     break;
 99:   case 1: /* write */
100:     if (hmem == HYPRE_MEMORY_HOST) {
101:       PetscCall(VecGetArrayWrite(v, ptr));
102:       *res = VecRestoreArrayWrite;
103:     } else {
104:       PetscCall(VecGetArrayWriteAndMemType(v, (PetscScalar **)ptr, &mtype));
105:       PetscCheck(PetscMemTypeDevice(mtype), comm, PETSC_ERR_ARG_WRONG, "HYPRE_MEMORY_DEVICE expects a device vector. You need to enable PETSc device support, for example, in some cases, -vec_type cuda");
106:       *res = VecRestoreArrayWriteAndMemType;
107:     }
108:     break;
109:   case 2: /* read/write */
110:     if (hmem == HYPRE_MEMORY_HOST) {
111:       PetscCall(VecGetArray(v, ptr));
112:       *res = VecRestoreArray;
113:     } else {
114:       PetscCall(VecGetArrayAndMemType(v, (PetscScalar **)ptr, &mtype));
115:       PetscCheck(PetscMemTypeDevice(mtype), comm, PETSC_ERR_ARG_WRONG, "HYPRE_MEMORY_DEVICE expects a device vector. You need to enable PETSc device support, for example, in some cases, -vec_type cuda");
116:       *res = VecRestoreArrayAndMemType;
117:     }
118:     break;
119:   default:
120:     SETERRQ(comm, PETSC_ERR_SUP, "Unhandled case %d", rw);
121:   }
122:   PetscFunctionReturn(PETSC_SUCCESS);
123: }

125: #define VecHYPRE_IJVectorMemoryLocation(v) hypre_IJVectorMemoryLocation((hypre_IJVector *)(v))

127: /* Temporarily pushes the array of the data in v to ij (read access)
128:    depending on the value of the ij memory location
129:    Must be completed with a call to VecHYPRE_IJVectorPopVec */
130: PetscErrorCode VecHYPRE_IJVectorPushVecRead(VecHYPRE_IJVector ij, Vec v)
131: {
132:   HYPRE_Complex *pv;

134:   PetscFunctionBegin;
136:   PetscCheck(!ij->pvec, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
137:   PetscCheck(!ij->hv, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
138:   PetscCall(VecGetArrayForHYPRE(v, 0, VecHYPRE_IJVectorMemoryLocation(ij->ij), (PetscScalar **)&pv, &ij->restore));
139:   VecHYPRE_ParVectorReplacePointer(ij->ij, pv, ij->hv);
140:   ij->pvec = v;
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: /* Temporarily pushes the array of the data in v to ij (write access)
145:    depending on the value of the ij memory location
146:    Must be completed with a call to VecHYPRE_IJVectorPopVec */
147: PetscErrorCode VecHYPRE_IJVectorPushVecWrite(VecHYPRE_IJVector ij, Vec v)
148: {
149:   HYPRE_Complex *pv;

151:   PetscFunctionBegin;
153:   PetscCheck(!ij->pvec, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
154:   PetscCheck(!ij->hv, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
155:   PetscCall(VecGetArrayForHYPRE(v, 1, VecHYPRE_IJVectorMemoryLocation(ij->ij), (PetscScalar **)&pv, &ij->restore));
156:   VecHYPRE_ParVectorReplacePointer(ij->ij, pv, ij->hv);
157:   ij->pvec = v;
158:   PetscFunctionReturn(PETSC_SUCCESS);
159: }

161: /* Temporarily pushes the array of the data in v to ij (read/write access)
162:    depending on the value of the ij memory location
163:    Must be completed with a call to VecHYPRE_IJVectorPopVec */
164: PetscErrorCode VecHYPRE_IJVectorPushVec(VecHYPRE_IJVector ij, Vec v)
165: {
166:   HYPRE_Complex *pv;

168:   PetscFunctionBegin;
170:   PetscCheck(!ij->pvec, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
171:   PetscCheck(!ij->hv, PetscObjectComm((PetscObject)v), PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
172:   PetscCall(VecGetArrayForHYPRE(v, 2, VecHYPRE_IJVectorMemoryLocation(ij->ij), (PetscScalar **)&pv, &ij->restore));
173:   VecHYPRE_ParVectorReplacePointer(ij->ij, pv, ij->hv);
174:   ij->pvec = v;
175:   PetscFunctionReturn(PETSC_SUCCESS);
176: }

178: /* Restores the pointer data to v */
179: PetscErrorCode VecHYPRE_IJVectorPopVec(VecHYPRE_IJVector ij)
180: {
181:   HYPRE_Complex *pv;

183:   PetscFunctionBegin;
184:   PetscCheck(ij->pvec, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPushVec()");
185:   PetscCheck(ij->restore, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPushVec()");
186:   VecHYPRE_ParVectorReplacePointer(ij->ij, ij->hv, pv);
187:   PetscCall((*ij->restore)(ij->pvec, (PetscScalar **)&pv));
188:   ij->hv      = NULL;
189:   ij->pvec    = NULL;
190:   ij->restore = NULL;
191:   PetscFunctionReturn(PETSC_SUCCESS);
192: }

194: PetscErrorCode VecHYPRE_IJBindToCPU(VecHYPRE_IJVector ij, PetscBool bind)
195: {
196:   HYPRE_MemoryLocation hmem = bind ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE;
197:   hypre_ParVector     *hij;

199:   PetscFunctionBegin;
200:   PetscCheck(!ij->pvec, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
201:   PetscCheck(!ij->hv, PETSC_COMM_SELF, PETSC_ERR_ORDER, "Forgot to call VecHYPRE_IJVectorPopVec()");
202: #if !defined(PETSC_HAVE_HYPRE_DEVICE)
203:   hmem = HYPRE_MEMORY_HOST;
204: #endif
205: #if PETSC_PKG_HYPRE_VERSION_GT(2, 19, 0)
206:   if (hmem != VecHYPRE_IJVectorMemoryLocation(ij->ij)) {
207:     PetscCallExternal(HYPRE_IJVectorGetObject, ij->ij, (void **)&hij);
208:     PetscCallExternal(hypre_ParVectorMigrate, hij, hmem);
209:   }
210: #endif
211:   PetscFunctionReturn(PETSC_SUCCESS);
212: }