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