Actual source code: mpiviennacl.cxx
1: /*
2: This file contains routines for Parallel vector operations.
3: */
4: #include <petscconf.h>
5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>
6: #include <../src/vec/vec/impls/seq/seqviennacl/viennaclvecimpl.h>
8: /*MC
9: VECVIENNACL - VECVIENNACL = "viennacl" - A `VECSEQVIENNACL` on a single-process communicator, and `VECMPIVIENNACL` otherwise.
11: Options Database Keys:
12: . -vec_type viennacl - sets the vector type to V`ECVIENNACL` during a call to `VecSetFromOptions()`
14: Level: beginner
16: .seealso: `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECSEQVIENNACL`, `VECMPIVIENNACL`, `VECSTANDARD`, `VecType`, `VecCreateMPI()`
17: M*/
19: static PetscErrorCode VecDestroy_MPIViennaCL(Vec v)
20: {
21: PetscFunctionBegin;
22: try {
23: if (v->spptr) {
24: delete ((Vec_ViennaCL *)v->spptr)->GPUarray_allocated;
25: delete (Vec_ViennaCL *)v->spptr;
26: }
27: } catch (std::exception const &ex) {
28: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_LIB, "ViennaCL error: %s", ex.what());
29: }
30: PetscCall(VecDestroy_MPI(v));
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode VecNorm_MPIViennaCL(Vec xin, NormType type, PetscReal *z)
35: {
36: PetscFunctionBegin;
37: if (type == NORM_2 || type == NORM_FROBENIUS) {
38: PetscCall(VecNorm_SeqViennaCL(xin, NORM_2, z));
39: *z *= *z;
40: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
41: *z = PetscSqrtReal(*z);
42: } else if (type == NORM_1) {
43: /* Find the local part */
44: PetscCall(VecNorm_SeqViennaCL(xin, NORM_1, z));
45: /* Find the global max */
46: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 1, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
47: } else if (type == NORM_INFINITY) {
48: /* Find the local max */
49: PetscCall(VecNorm_SeqViennaCL(xin, NORM_INFINITY, z));
50: /* Find the global max */
51: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 1, MPIU_REAL, MPIU_MAX, PetscObjectComm((PetscObject)xin)));
52: } else if (type == NORM_1_AND_2) {
53: PetscCall(VecNorm_SeqViennaCL(xin, NORM_1, z));
54: PetscCall(VecNorm_SeqViennaCL(xin, NORM_2, z + 1));
55: z[1] = z[1] * z[1];
56: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 2, MPIU_REAL, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
57: z[1] = PetscSqrtReal(z[1]);
58: }
59: PetscFunctionReturn(PETSC_SUCCESS);
60: }
62: static PetscErrorCode VecDot_MPIViennaCL(Vec xin, Vec yin, PetscScalar *z)
63: {
64: PetscFunctionBegin;
65: PetscCall(VecDot_SeqViennaCL(xin, yin, z));
66: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
67: PetscFunctionReturn(PETSC_SUCCESS);
68: }
70: static PetscErrorCode VecTDot_MPIViennaCL(Vec xin, Vec yin, PetscScalar *z)
71: {
72: PetscFunctionBegin;
73: PetscCall(VecTDot_SeqViennaCL(xin, yin, z));
74: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, 1, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
75: PetscFunctionReturn(PETSC_SUCCESS);
76: }
78: static PetscErrorCode VecMDot_MPIViennaCL(Vec xin, PetscInt nv, const Vec y[], PetscScalar *z)
79: {
80: PetscFunctionBegin;
81: PetscCall(VecMDot_SeqViennaCL(xin, nv, y, z));
82: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, z, nv, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)xin)));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: /*MC
87: VECMPIVIENNACL - VECMPIVIENNACL = "mpiviennacl" - The basic parallel vector, modified to use ViennaCL
89: Options Database Keys:
90: . -vec_type mpiviennacl - sets the vector type to VECMPIVIENNACL during a call to VecSetFromOptions()
92: Level: beginner
94: .seealso: `VecCreate()`, `VecSetType()`, `VecSetFromOptions()`, `VecCreateMPIWithArray()`, `VECMPI`, `VecType`, `VecCreateMPI()`
95: M*/
97: static PetscErrorCode VecDuplicate_MPIViennaCL(Vec win, Vec *v)
98: {
99: Vec_MPI *vw, *w = (Vec_MPI *)win->data;
100: PetscScalar *array;
102: PetscFunctionBegin;
103: PetscCall(VecCreate(PetscObjectComm((PetscObject)win), v));
104: PetscCall(PetscLayoutReference(win->map, &(*v)->map));
106: PetscCall(VecCreate_MPI_Private(*v, PETSC_FALSE, w->nghost, 0));
107: vw = (Vec_MPI *)(*v)->data;
108: (*v)->ops[0] = win->ops[0];
110: /* save local representation of the parallel vector (and scatter) if it exists */
111: if (w->localrep) {
112: PetscCall(VecGetArray(*v, &array));
113: PetscCall(VecCreateSeqWithArray(PETSC_COMM_SELF, 1, win->map->n + w->nghost, array, &vw->localrep));
114: vw->localrep->ops[0] = w->localrep->ops[0];
115: PetscCall(VecRestoreArray(*v, &array));
116: vw->localupdate = w->localupdate;
117: PetscCall(PetscObjectReference((PetscObject)vw->localupdate));
118: }
120: /* New vector should inherit stashing property of parent */
121: (*v)->stash.donotstash = win->stash.donotstash;
122: (*v)->stash.ignorenegidx = win->stash.ignorenegidx;
124: /* change type_name appropriately */
125: PetscCall(PetscObjectChangeTypeName((PetscObject)*v, VECMPIVIENNACL));
127: PetscCall(PetscObjectListDuplicate(((PetscObject)win)->olist, &((PetscObject)*v)->olist));
128: PetscCall(PetscFunctionListDuplicate(((PetscObject)win)->qlist, &((PetscObject)*v)->qlist));
129: (*v)->map->bs = win->map->bs;
130: (*v)->bstash.bs = win->bstash.bs;
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode VecDotNorm2_MPIViennaCL(Vec s, Vec t, PetscScalar *dp, PetscScalar *nm)
135: {
136: PetscScalar sum[2];
138: PetscFunctionBegin;
139: PetscCall(VecDotNorm2_SeqViennaCL(s, t, sum, sum + 1));
140: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, sum, 2, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)s)));
141: *dp = sum[0];
142: *nm = sum[1];
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: static PetscErrorCode VecBindToCPU_MPIViennaCL(Vec vv, PetscBool bind)
147: {
148: PetscFunctionBegin;
149: vv->boundtocpu = bind;
151: if (bind) {
152: PetscCall(VecViennaCLCopyFromGPU(vv));
153: vv->offloadmask = PETSC_OFFLOAD_CPU; /* since the CPU code will likely change values in the vector */
154: vv->ops->dotnorm2 = NULL;
155: vv->ops->waxpy = VecWAXPY_Seq;
156: vv->ops->dot = VecDot_MPI;
157: vv->ops->mdot = VecMDot_MPI;
158: vv->ops->tdot = VecTDot_MPI;
159: vv->ops->norm = VecNorm_MPI;
160: vv->ops->scale = VecScale_Seq;
161: vv->ops->copy = VecCopy_Seq;
162: vv->ops->set = VecSet_Seq;
163: vv->ops->swap = VecSwap_Seq;
164: vv->ops->axpy = VecAXPY_Seq;
165: vv->ops->axpby = VecAXPBY_Seq;
166: vv->ops->maxpy = VecMAXPY_Seq;
167: vv->ops->aypx = VecAYPX_Seq;
168: vv->ops->axpbypcz = VecAXPBYPCZ_Seq;
169: vv->ops->pointwisemult = VecPointwiseMult_Seq;
170: vv->ops->setrandom = VecSetRandom_Seq;
171: vv->ops->placearray = VecPlaceArray_Seq;
172: vv->ops->replacearray = VecReplaceArray_Seq;
173: vv->ops->resetarray = VecResetArray_Seq;
174: vv->ops->dot_local = VecDot_Seq;
175: vv->ops->tdot_local = VecTDot_Seq;
176: vv->ops->norm_local = VecNorm_Seq;
177: vv->ops->mdot_local = VecMDot_Seq;
178: vv->ops->pointwisedivide = VecPointwiseDivide_Seq;
179: vv->ops->getlocalvector = NULL;
180: vv->ops->restorelocalvector = NULL;
181: vv->ops->getlocalvectorread = NULL;
182: vv->ops->restorelocalvectorread = NULL;
183: vv->ops->getarraywrite = NULL;
184: } else {
185: vv->ops->dotnorm2 = VecDotNorm2_MPIViennaCL;
186: vv->ops->waxpy = VecWAXPY_SeqViennaCL;
187: vv->ops->duplicate = VecDuplicate_MPIViennaCL;
188: vv->ops->dot = VecDot_MPIViennaCL;
189: vv->ops->mdot = VecMDot_MPIViennaCL;
190: vv->ops->tdot = VecTDot_MPIViennaCL;
191: vv->ops->norm = VecNorm_MPIViennaCL;
192: vv->ops->scale = VecScale_SeqViennaCL;
193: vv->ops->copy = VecCopy_SeqViennaCL;
194: vv->ops->set = VecSet_SeqViennaCL;
195: vv->ops->swap = VecSwap_SeqViennaCL;
196: vv->ops->axpy = VecAXPY_SeqViennaCL;
197: vv->ops->axpby = VecAXPBY_SeqViennaCL;
198: vv->ops->maxpy = VecMAXPY_SeqViennaCL;
199: vv->ops->aypx = VecAYPX_SeqViennaCL;
200: vv->ops->axpbypcz = VecAXPBYPCZ_SeqViennaCL;
201: vv->ops->pointwisemult = VecPointwiseMult_SeqViennaCL;
202: vv->ops->setrandom = VecSetRandom_SeqViennaCL;
203: vv->ops->dot_local = VecDot_SeqViennaCL;
204: vv->ops->tdot_local = VecTDot_SeqViennaCL;
205: vv->ops->norm_local = VecNorm_SeqViennaCL;
206: vv->ops->mdot_local = VecMDot_SeqViennaCL;
207: vv->ops->destroy = VecDestroy_MPIViennaCL;
208: vv->ops->pointwisedivide = VecPointwiseDivide_SeqViennaCL;
209: vv->ops->placearray = VecPlaceArray_SeqViennaCL;
210: vv->ops->replacearray = VecReplaceArray_SeqViennaCL;
211: vv->ops->resetarray = VecResetArray_SeqViennaCL;
212: vv->ops->getarraywrite = VecGetArrayWrite_SeqViennaCL;
213: vv->ops->getarray = VecGetArray_SeqViennaCL;
214: vv->ops->restorearray = VecRestoreArray_SeqViennaCL;
215: }
216: vv->ops->duplicatevecs = VecDuplicateVecs_Default;
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: PETSC_EXTERN PetscErrorCode VecCreate_MPIViennaCL(Vec vv)
221: {
222: PetscFunctionBegin;
223: PetscCall(PetscLayoutSetUp(vv->map));
224: PetscCall(VecViennaCLAllocateCheck(vv));
225: PetscCall(VecCreate_MPIViennaCL_Private(vv, PETSC_FALSE, 0, ((Vec_ViennaCL *)vv->spptr)->GPUarray));
226: PetscCall(VecViennaCLAllocateCheckHost(vv));
227: PetscCall(VecSet(vv, 0.0));
228: PetscCall(VecSet_Seq(vv, 0.0));
229: vv->offloadmask = PETSC_OFFLOAD_BOTH;
230: PetscFunctionReturn(PETSC_SUCCESS);
231: }
233: PETSC_EXTERN PetscErrorCode VecCreate_ViennaCL(Vec v)
234: {
235: PetscMPIInt size;
237: PetscFunctionBegin;
238: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)v), &size));
239: if (size == 1) PetscCall(VecSetType(v, VECSEQVIENNACL));
240: else PetscCall(VecSetType(v, VECMPIVIENNACL));
241: PetscFunctionReturn(PETSC_SUCCESS);
242: }
244: /*@C
245: VecCreateMPIViennaCLWithArray - Creates a parallel, array-style vector,
246: where the user provides the viennacl vector to store the vector values.
248: Collective
250: Input Parameters:
251: + comm - the MPI communicator to use
252: . bs - block size, same meaning as in `VecSetBlockSize()`
253: . n - local vector length, cannot be `PETSC_DECIDE`
254: . N - global vector length (or `PETSC_DECIDE` to have calculated)
255: - array - the user provided GPU array to store the vector values
257: Output Parameter:
258: . vv - the vector
260: Level: intermediate
262: Notes:
263: Use `VecDuplicate()` or `VecDuplicateVecs(`) to form additional vectors of the
264: same type as an existing vector.
266: If the user-provided array is `NULL`, then `VecViennaCLPlaceArray()` can be used
267: at a later stage to SET the array for storing the vector values.
269: PETSc does NOT free the array when the vector is destroyed via `VecDestroy()`.
270: The user should not free the array until the vector is destroyed.
272: .seealso: `VecCreateSeqViennaCLWithArray()`, `VecCreateMPIWithArray()`, `VecCreateSeqWithArray()`,
273: `VecCreate()`, `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecViennaCLPlaceArray()`
274: @*/
275: PetscErrorCode VecCreateMPIViennaCLWithArray(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, const ViennaCLVector *array, Vec *vv) PeNS
276: {
277: PetscFunctionBegin;
278: PetscCheck(n != PETSC_DECIDE, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Must set local size of vector");
279: PetscCall(PetscSplitOwnership(comm, &n, &N));
280: PetscCall(VecCreate(comm, vv));
281: PetscCall(VecSetSizes(*vv, n, N));
282: PetscCall(VecSetBlockSize(*vv, bs));
283: PetscCall(VecCreate_MPIViennaCL_Private(*vv, PETSC_FALSE, 0, array));
284: PetscFunctionReturn(PETSC_SUCCESS);
285: }
287: /*@C
288: VecCreateMPIViennaCLWithArrays - Creates a parallel, array-style vector,
289: where the user provides the ViennaCL vector to store the vector values.
291: Collective
293: Input Parameters:
294: + comm - the MPI communicator to use
295: . bs - block size, same meaning as VecSetBlockSize()
296: . n - local vector length, cannot be PETSC_DECIDE
297: . N - global vector length (or PETSC_DECIDE to have calculated)
298: . cpuarray - the user provided CPU array to store the vector values
299: - viennaclvec - ViennaCL vector where the Vec entries are to be stored on the device.
301: Output Parameter:
302: . vv - the vector
304: Notes:
305: If both cpuarray and viennaclvec are provided, the caller must ensure that
306: the provided arrays have identical values.
308: Use VecDuplicate() or VecDuplicateVecs() to form additional vectors of the
309: same type as an existing vector.
311: PETSc does NOT free the provided arrays when the vector is destroyed via
312: VecDestroy(). The user should not free the array until the vector is
313: destroyed.
315: Level: intermediate
317: .seealso: `VecCreateSeqViennaCLWithArrays()`, `VecCreateMPIWithArray()`,
318: `VecCreate()`, `VecDuplicate()`, `VecDuplicateVecs()`, `VecCreateGhost()`,
319: `VecCreateMPI()`, `VecCreateGhostWithArray()`, `VecViennaCLPlaceArray()`,
320: `VecPlaceArray()`, `VecCreateMPICUDAWithArrays()`,
321: `VecViennaCLAllocateCheckHost()`
322: @*/
323: PetscErrorCode VecCreateMPIViennaCLWithArrays(MPI_Comm comm, PetscInt bs, PetscInt n, PetscInt N, const PetscScalar cpuarray[], const ViennaCLVector *viennaclvec, Vec *vv) PeNS
324: {
325: PetscFunctionBegin;
326: PetscCall(VecCreateMPIViennaCLWithArray(comm, bs, n, N, viennaclvec, vv));
327: if (cpuarray && viennaclvec) {
328: Vec_MPI *s = (Vec_MPI *)((*vv)->data);
329: s->array = (PetscScalar *)cpuarray;
330: (*vv)->offloadmask = PETSC_OFFLOAD_BOTH;
331: } else if (cpuarray) {
332: Vec_MPI *s = (Vec_MPI *)((*vv)->data);
333: s->array = (PetscScalar *)cpuarray;
334: (*vv)->offloadmask = PETSC_OFFLOAD_CPU;
335: } else if (viennaclvec) {
336: (*vv)->offloadmask = PETSC_OFFLOAD_GPU;
337: } else {
338: (*vv)->offloadmask = PETSC_OFFLOAD_UNALLOCATED;
339: }
340: PetscFunctionReturn(PETSC_SUCCESS);
341: }
343: PetscErrorCode VecCreate_MPIViennaCL_Private(Vec vv, PetscBool alloc, PetscInt nghost, const ViennaCLVector *array)
344: {
345: Vec_ViennaCL *vecviennacl;
347: PetscFunctionBegin;
348: PetscCall(VecCreate_MPI_Private(vv, PETSC_FALSE, 0, 0));
349: PetscCall(PetscObjectChangeTypeName((PetscObject)vv, VECMPIVIENNACL));
351: PetscCall(VecBindToCPU_MPIViennaCL(vv, PETSC_FALSE));
352: vv->ops->bindtocpu = VecBindToCPU_MPIViennaCL;
354: if (alloc && !array) {
355: PetscCall(VecViennaCLAllocateCheck(vv));
356: PetscCall(VecViennaCLAllocateCheckHost(vv));
357: PetscCall(VecSet(vv, 0.0));
358: PetscCall(VecSet_Seq(vv, 0.0));
359: vv->offloadmask = PETSC_OFFLOAD_BOTH;
360: }
361: if (array) {
362: if (!vv->spptr) vv->spptr = new Vec_ViennaCL;
363: vecviennacl = (Vec_ViennaCL *)vv->spptr;
364: vecviennacl->GPUarray_allocated = 0;
365: vecviennacl->GPUarray = (ViennaCLVector *)array;
366: vv->offloadmask = PETSC_OFFLOAD_UNALLOCATED;
367: }
368: PetscFunctionReturn(PETSC_SUCCESS);
369: }