Actual source code: mpiviennacl.cxx

petsc-master 2019-08-18
Report Typos and Errors

  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: PetscErrorCode VecDestroy_MPIViennaCL(Vec v)
 10: {

 14:   try {
 15:     if (v->spptr) {
 16:       delete ((Vec_ViennaCL*)v->spptr)->GPUarray;
 17:       delete (Vec_ViennaCL*) v->spptr;
 18:     }
 19:   } catch(std::exception const & ex) {
 20:     SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"ViennaCL error: %s", ex.what());
 21:   }
 22:   VecDestroy_MPI(v);
 23:   return(0);
 24: }

 26: PetscErrorCode VecNorm_MPIViennaCL(Vec xin,NormType type,PetscReal *z)
 27: {
 28:   PetscReal      sum,work = 0.0;

 32:   if (type == NORM_2 || type == NORM_FROBENIUS) {
 33:     VecNorm_SeqViennaCL(xin,NORM_2,&work);
 34:     work *= work;
 35:     MPIU_Allreduce(&work,&sum,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 36:     *z    = PetscSqrtReal(sum);
 37:   } else if (type == NORM_1) {
 38:     /* Find the local part */
 39:     VecNorm_SeqViennaCL(xin,NORM_1,&work);
 40:     /* Find the global max */
 41:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 42:   } else if (type == NORM_INFINITY) {
 43:     /* Find the local max */
 44:     VecNorm_SeqViennaCL(xin,NORM_INFINITY,&work);
 45:     /* Find the global max */
 46:     MPIU_Allreduce(&work,z,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)xin));
 47:   } else if (type == NORM_1_AND_2) {
 48:     PetscReal temp[2];
 49:     VecNorm_SeqViennaCL(xin,NORM_1,temp);
 50:     VecNorm_SeqViennaCL(xin,NORM_2,temp+1);
 51:     temp[1] = temp[1]*temp[1];
 52:     MPIU_Allreduce(temp,z,2,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 53:     z[1] = PetscSqrtReal(z[1]);
 54:   }
 55:   return(0);
 56: }

 58: PetscErrorCode VecDot_MPIViennaCL(Vec xin,Vec yin,PetscScalar *z)
 59: {
 60:   PetscScalar    sum,work;

 64:   VecDot_SeqViennaCL(xin,yin,&work);
 65:   MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 66:   *z   = sum;
 67:   return(0);
 68: }

 70: PetscErrorCode VecTDot_MPIViennaCL(Vec xin,Vec yin,PetscScalar *z)
 71: {
 72:   PetscScalar    sum,work;

 76:   VecTDot_SeqViennaCL(xin,yin,&work);
 77:   MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 78:   *z   = sum;
 79:   return(0);
 80: }

 82: PetscErrorCode VecMDot_MPIViennaCL(Vec xin,PetscInt nv,const Vec y[],PetscScalar *z)
 83: {
 84:   PetscScalar    awork[128],*work = awork;

 88:   if (nv > 128) {
 89:     PetscMalloc1(nv,&work);
 90:   }
 91:   VecMDot_SeqViennaCL(xin,nv,y,work);
 92:   MPIU_Allreduce(work,z,nv,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)xin));
 93:   if (nv > 128) {
 94:     PetscFree(work);
 95:   }
 96:   return(0);
 97: }

 99: /*MC
100:    VECMPIVIENNACL - VECMPIVIENNACL = "mpiviennacl" - The basic parallel vector, modified to use ViennaCL

102:    Options Database Keys:
103: . -vec_type mpiviennacl - sets the vector type to VECMPIVIENNACL during a call to VecSetFromOptions()

105:   Level: beginner

107: .seealso: VecCreate(), VecSetType(), VecSetFromOptions(), VecCreateMPIWithArray(), VECMPI, VecType, VecCreateMPI(), VecCreateMPI()
108: M*/


111: PetscErrorCode VecDuplicate_MPIViennaCL(Vec win,Vec *v)
112: {
114:   Vec_MPI        *vw,*w = (Vec_MPI*)win->data;
115:   PetscScalar    *array;

118:   VecCreate(PetscObjectComm((PetscObject)win),v);
119:   PetscLayoutReference(win->map,&(*v)->map);

121:   VecCreate_MPI_Private(*v,PETSC_FALSE,w->nghost,0);
122:   vw   = (Vec_MPI*)(*v)->data;
123:   PetscMemcpy((*v)->ops,win->ops,sizeof(struct _VecOps));

125:   /* save local representation of the parallel vector (and scatter) if it exists */
126:   if (w->localrep) {
127:     VecGetArray(*v,&array);
128:     VecCreateSeqWithArray(PETSC_COMM_SELF,1,win->map->n+w->nghost,array,&vw->localrep);
129:     PetscMemcpy(vw->localrep->ops,w->localrep->ops,sizeof(struct _VecOps));
130:     VecRestoreArray(*v,&array);
131:     PetscLogObjectParent((PetscObject)*v,(PetscObject)vw->localrep);
132:     vw->localupdate = w->localupdate;
133:     if (vw->localupdate) {
134:       PetscObjectReference((PetscObject)vw->localupdate);
135:     }
136:   }

138:   /* New vector should inherit stashing property of parent */
139:   (*v)->stash.donotstash   = win->stash.donotstash;
140:   (*v)->stash.ignorenegidx = win->stash.ignorenegidx;

142:   /* change type_name appropriately */
143:   PetscObjectChangeTypeName((PetscObject)(*v),VECMPIVIENNACL);

145:   PetscObjectListDuplicate(((PetscObject)win)->olist,&((PetscObject)(*v))->olist);
146:   PetscFunctionListDuplicate(((PetscObject)win)->qlist,&((PetscObject)(*v))->qlist);
147:   (*v)->map->bs   = PetscAbs(win->map->bs);
148:   (*v)->bstash.bs = win->bstash.bs;
149:   return(0);
150: }

152: PetscErrorCode VecDotNorm2_MPIViennaCL(Vec s,Vec t,PetscScalar *dp,PetscScalar *nm)
153: {
155:   PetscScalar    work[2],sum[2];

158:   VecDotNorm2_SeqViennaCL(s,t,work,work+1);
159:   MPIU_Allreduce((void*)&work,(void*)&sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
160:   *dp  = sum[0];
161:   *nm  = sum[1];
162:   return(0);
163: }

165: PETSC_EXTERN PetscErrorCode VecCreate_MPIViennaCL(Vec vv)
166: {

170:   VecCreate_MPI_Private(vv,PETSC_FALSE,0,0);
171:   PetscObjectChangeTypeName((PetscObject)vv,VECMPIVIENNACL);

173:   vv->ops->dotnorm2        = VecDotNorm2_MPIViennaCL;
174:   vv->ops->waxpy           = VecWAXPY_SeqViennaCL;
175:   vv->ops->duplicate       = VecDuplicate_MPIViennaCL;
176:   vv->ops->dot             = VecDot_MPIViennaCL;
177:   vv->ops->mdot            = VecMDot_MPIViennaCL;
178:   vv->ops->tdot            = VecTDot_MPIViennaCL;
179:   vv->ops->norm            = VecNorm_MPIViennaCL;
180:   vv->ops->scale           = VecScale_SeqViennaCL;
181:   vv->ops->copy            = VecCopy_SeqViennaCL;
182:   vv->ops->set             = VecSet_SeqViennaCL;
183:   vv->ops->swap            = VecSwap_SeqViennaCL;
184:   vv->ops->axpy            = VecAXPY_SeqViennaCL;
185:   vv->ops->axpby           = VecAXPBY_SeqViennaCL;
186:   vv->ops->maxpy           = VecMAXPY_SeqViennaCL;
187:   vv->ops->aypx            = VecAYPX_SeqViennaCL;
188:   vv->ops->axpbypcz        = VecAXPBYPCZ_SeqViennaCL;
189:   vv->ops->pointwisemult   = VecPointwiseMult_SeqViennaCL;
190:   vv->ops->setrandom       = VecSetRandom_SeqViennaCL;
191:   vv->ops->dot_local       = VecDot_SeqViennaCL;
192:   vv->ops->tdot_local      = VecTDot_SeqViennaCL;
193:   vv->ops->norm_local      = VecNorm_SeqViennaCL;
194:   vv->ops->mdot_local      = VecMDot_SeqViennaCL;
195:   vv->ops->destroy         = VecDestroy_MPIViennaCL;
196:   vv->ops->pointwisedivide = VecPointwiseDivide_SeqViennaCL;
197:   vv->ops->placearray      = VecPlaceArray_SeqViennaCL;
198:   vv->ops->replacearray    = VecReplaceArray_SeqViennaCL;
199:   vv->ops->resetarray      = VecResetArray_SeqViennaCL;
200:   /*
201:      get values?
202:   */
203:   VecViennaCLAllocateCheck(vv);
204:   vv->valid_GPU_array      = PETSC_OFFLOAD_GPU;
205:   VecSet(vv,0.0);
206:   return(0);
207: }


210: PETSC_EXTERN PetscErrorCode VecCreate_ViennaCL(Vec v)
211: {
213:   PetscMPIInt    size;

216:   MPI_Comm_size(PetscObjectComm((PetscObject)v),&size);
217:   if (size == 1) {
218:     VecSetType(v,VECSEQVIENNACL);
219:   } else {
220:     VecSetType(v,VECMPIVIENNACL);
221:   }
222:   return(0);
223: }