Actual source code: veccomp0.h

slepc-3.11.2 2019-07-30
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2019, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
  8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  9: */

 11: #include <petsc/private/vecimpl.h>

 13: #if defined(__WITH_MPI__)
 15: #else
 17: #endif

 22: PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
 23: {
 24:   PetscScalar    sum = 0.0,work;
 25:   PetscInt       i;
 27:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

 30:   SlepcValidVecComp(a,1);
 31:   SlepcValidVecComp(b,2);
 32:   if (as->x[0]->ops->dot_local) {
 33:     for (i=0,sum=0.0;i<as->n->n;i++) {
 34:       as->x[i]->ops->dot_local(as->x[i],bs->x[i],&work);
 35:       sum += work;
 36:     }
 37: #if defined(__WITH_MPI__)
 38:     work = sum;
 39:     MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
 40: #endif
 41:   } else {
 42:     for (i=0,sum=0.0;i<as->n->n;i++) {
 43:       VecDot(as->x[i],bs->x[i],&work);
 44:       sum += work;
 45:     }
 46:   }
 47:   *z = sum;
 48:   return(0);
 49: }

 51: PetscErrorCode __SUF__(VecMDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
 52: {
 53:   PetscScalar    *work,*work0,*r;
 55:   Vec_Comp       *as = (Vec_Comp*)a->data;
 56:   Vec            *bx;
 57:   PetscInt       i,j;

 60:   SlepcValidVecComp(a,1);
 61:   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);

 63:   if (as->n->n == 0) {
 64:     *z = 0;
 65:     return(0);
 66:   }

 68:   PetscMalloc2(n,&work0,n,&bx);

 70: #if defined(__WITH_MPI__)
 71:   if (as->x[0]->ops->mdot_local) {
 72:     r = work0; work = z;
 73:   } else
 74: #endif
 75:   {
 76:     r = z; work = work0;
 77:   }

 79:   /* z[i] <- a.x' * b[i].x */
 80:   for (i=0;i<n;i++) r[i] = 0.0;
 81:   for (j=0;j<as->n->n;j++) {
 82:     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
 83:     if (as->x[0]->ops->mdot_local) {
 84:       as->x[j]->ops->mdot_local(as->x[j],n,bx,work);
 85:     } else {
 86:       VecMDot(as->x[j],n,bx,work);
 87:     }
 88:     for (i=0;i<n;i++) r[i] += work[i];
 89:   }

 91:   /* If def(__WITH_MPI__) and exists mdot_local */
 92: #if defined(__WITH_MPI__)
 93:   if (as->x[0]->ops->mdot_local) {
 94:     /* z[i] <- Allreduce(work[i]) */
 95:     MPI_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
 96:   }
 97: #endif

 99:   PetscFree2(work0,bx);
100:   return(0);
101: }

103: PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
104: {
105:   PetscScalar    sum = 0.0,work;
106:   PetscInt       i;
108:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

111:   SlepcValidVecComp(a,1);
112:   SlepcValidVecComp(b,2);
113:   if (as->x[0]->ops->tdot_local) {
114:     for (i=0,sum=0.0;i<as->n->n;i++) {
115:       as->x[i]->ops->tdot_local(as->x[i],bs->x[i],&work);
116:       sum += work;
117:     }
118: #if defined(__WITH_MPI__)
119:     work = sum;
120:     MPI_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
121: #endif
122:   } else {
123:     for (i=0,sum=0.0;i<as->n->n;i++) {
124:       VecTDot(as->x[i],bs->x[i],&work);
125:       sum += work;
126:     }
127:   }
128:   *z = sum;
129:   return(0);
130: }

132: PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
133: {
134:   PetscScalar    *work,*work0,*r;
136:   Vec_Comp       *as = (Vec_Comp*)a->data;
137:   Vec            *bx;
138:   PetscInt       i,j;

141:   SlepcValidVecComp(a,1);
142:   for (i=0;i<n;i++) SlepcValidVecComp(b[i],3);

144:   if (as->n->n == 0) {
145:     *z = 0;
146:     return(0);
147:   }

149:   PetscMalloc2(n,&work0,n,&bx);

151: #if defined(__WITH_MPI__)
152:   if (as->x[0]->ops->mtdot_local) {
153:     r = work0; work = z;
154:   } else
155: #endif
156:   {
157:     r = z; work = work0;
158:   }

160:   /* z[i] <- a.x' * b[i].x */
161:   for (i=0;i<n;i++) r[i] = 0.0;
162:   for (j=0;j<as->n->n;j++) {
163:     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
164:     if (as->x[0]->ops->mtdot_local) {
165:       as->x[j]->ops->mtdot_local(as->x[j],n,bx,work);
166:     } else {
167:       VecMTDot(as->x[j],n,bx,work);
168:     }
169:     for (i=0;i<n;i++) r[i] += work[i];
170:   }

172:   /* If def(__WITH_MPI__) and exists mtdot_local */
173: #if defined(__WITH_MPI__)
174:   if (as->x[0]->ops->mtdot_local) {
175:     /* z[i] <- Allreduce(work[i]) */
176:     MPI_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
177:   }
178: #endif

180:   PetscFree2(work0,bx);
181:   return(0);
182: }

184: PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
185: {
186:   PetscReal      work[3],s=0.0;
188:   Vec_Comp       *as = (Vec_Comp*)a->data;
189:   PetscInt       i;

192:   SlepcValidVecComp(a,1);
193:   /* Initialize norm */
194:   switch (t) {
195:     case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
196:     case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
197:     case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
198:   }
199:   for (i=0;i<as->n->n;i++) {
200:     if (as->x[0]->ops->norm_local) {
201:       as->x[0]->ops->norm_local(as->x[i],t,work);
202:     } else {
203:       VecNorm(as->x[i],t,work);
204:     }
205:     /* norm+= work */
206:     switch (t) {
207:       case NORM_1: *norm+= *work; break;
208:       case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
209:       case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
210:       case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
211:     }
212:   }

214:   /* If def(__WITH_MPI__) and exists norm_local */
215: #if defined(__WITH_MPI__)
216:   if (as->x[0]->ops->norm_local) {
217:     PetscReal work0[3];
218:     /* norm <- Allreduce(work) */
219:     switch (t) {
220:     case NORM_1:
221:       work[0] = *norm;
222:       MPI_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a));
223:       break;
224:     case NORM_2: case NORM_FROBENIUS:
225:       work[0] = *norm; work[1] = s;
226:       MPI_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
227:       *norm = GetNorm2(work0[0],work0[1]);
228:       break;
229:     case NORM_1_AND_2:
230:       work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
231:       MPI_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
232:       norm[0] = work0[0];
233:       norm[1] = GetNorm2(work0[1],work0[2]);
234:       break;
235:     case NORM_INFINITY:
236:       work[0] = *norm;
237:       MPI_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a));
238:       break;
239:     }
240:   }
241: #else
242:   /* Norm correction */
243:   switch (t) {
244:     case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
245:     case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
246:     default: ;
247:   }
248: #endif
249:   return(0);
250: }

252: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
253: {
254:   PetscErrorCode    ierr;
255:   PetscScalar       dp0,nm0,dp1,nm1;
256:   const PetscScalar *vx,*wx;
257:   Vec_Comp          *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
258:   PetscInt          i,n;
259:   PetscBool         t0,t1;
260: #if defined(__WITH_MPI__)
261:   PetscScalar       work[4];
262: #endif

265:   /* Compute recursively the local part */
266:   dp0 = nm0 = 0.0;
267:   PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0);
268:   PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1);
269:   if (t0 && t1) {
270:     SlepcValidVecComp(v,1);
271:     SlepcValidVecComp(w,2);
272:     for (i=0;i<vs->n->n;i++) {
273:       VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1);
274:       dp0 += dp1;
275:       nm0 += nm1;
276:     }
277:   } else if (!t0 && !t1) {
278:     VecGetLocalSize(v,&n);
279:     VecGetArrayRead(v,&vx);
280:     VecGetArrayRead(w,&wx);
281:     for (i=0;i<n;i++) {
282:       dp0 += vx[i]*PetscConj(wx[i]);
283:       nm0 += wx[i]*PetscConj(wx[i]);
284:     }
285:     VecRestoreArrayRead(v,&vx);
286:     VecRestoreArrayRead(w,&wx);
287:   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");

289: #if defined(__WITH_MPI__)
290:     /* [dp, nm] <- Allreduce([dp0, nm0]) */
291:     work[0] = dp0; work[1] = nm0;
292:     MPI_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
293:     *dp = work[2]; *nm = work[3];
294: #else
295:     *dp = dp0; *nm = nm0;
296: #endif
297:   return(0);
298: }