Actual source code: veccomp0.h

slepc-3.17.2 2022-08-09
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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;
 26:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

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

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

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

 59:   if (as->n->n == 0) {
 60:     *z = 0;
 61:     PetscFunctionReturn(0);
 62:   }

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

 66: #if defined(__WITH_MPI__)
 67:   if (as->x[0]->ops->mdot_local) {
 68:     r = work0; work = z;
 69:   } else
 70: #endif
 71:   {
 72:     r = z; work = work0;
 73:   }

 75:   /* z[i] <- a.x' * b[i].x */
 76:   for (i=0;i<n;i++) r[i] = 0.0;
 77:   for (j=0;j<as->n->n;j++) {
 78:     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
 79:     if (as->x[0]->ops->mdot_local) as->x[j]->ops->mdot_local(as->x[j],n,bx,work);
 80:     else VecMDot(as->x[j],n,bx,work);
 81:     for (i=0;i<n;i++) r[i] += work[i];
 82:   }

 84:   /* If def(__WITH_MPI__) and exists mdot_local */
 85: #if defined(__WITH_MPI__)
 86:   if (as->x[0]->ops->mdot_local) {
 87:     /* z[i] <- Allreduce(work[i]) */
 88:     MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
 89:   }
 90: #endif

 92:   PetscFree2(work0,bx);
 93:   PetscFunctionReturn(0);
 94: }

 96: PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
 97: {
 98:   PetscScalar    sum = 0.0,work;
 99:   PetscInt       i;
100:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

102:   SlepcValidVecComp(a,1);
103:   SlepcValidVecComp(b,2);
104:   if (as->x[0]->ops->tdot_local) {
105:     for (i=0,sum=0.0;i<as->n->n;i++) {
106:       as->x[i]->ops->tdot_local(as->x[i],bs->x[i],&work);
107:       sum += work;
108:     }
109: #if defined(__WITH_MPI__)
110:     work = sum;
111:     MPIU_Allreduce(&work,&sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
112: #endif
113:   } else {
114:     for (i=0,sum=0.0;i<as->n->n;i++) {
115:       VecTDot(as->x[i],bs->x[i],&work);
116:       sum += work;
117:     }
118:   }
119:   *z = sum;
120:   PetscFunctionReturn(0);
121: }

123: PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
124: {
125:   PetscScalar    *work,*work0,*r;
126:   Vec_Comp       *as = (Vec_Comp*)a->data;
127:   Vec            *bx;
128:   PetscInt       i,j;

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

133:   if (as->n->n == 0) {
134:     *z = 0;
135:     PetscFunctionReturn(0);
136:   }

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

140: #if defined(__WITH_MPI__)
141:   if (as->x[0]->ops->mtdot_local) {
142:     r = work0; work = z;
143:   } else
144: #endif
145:   {
146:     r = z; work = work0;
147:   }

149:   /* z[i] <- a.x' * b[i].x */
150:   for (i=0;i<n;i++) r[i] = 0.0;
151:   for (j=0;j<as->n->n;j++) {
152:     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
153:     if (as->x[0]->ops->mtdot_local) as->x[j]->ops->mtdot_local(as->x[j],n,bx,work);
154:     else VecMTDot(as->x[j],n,bx,work);
155:     for (i=0;i<n;i++) r[i] += work[i];
156:   }

158:   /* If def(__WITH_MPI__) and exists mtdot_local */
159: #if defined(__WITH_MPI__)
160:   if (as->x[0]->ops->mtdot_local) {
161:     /* z[i] <- Allreduce(work[i]) */
162:     MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a));
163:   }
164: #endif

166:   PetscFree2(work0,bx);
167:   PetscFunctionReturn(0);
168: }

170: PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
171: {
172:   PetscReal      work[3],s=0.0;
173:   Vec_Comp       *as = (Vec_Comp*)a->data;
174:   PetscInt       i;

176:   SlepcValidVecComp(a,1);
177:   /* Initialize norm */
178:   switch (t) {
179:     case NORM_1: case NORM_INFINITY: *norm = 0.0; break;
180:     case NORM_2: case NORM_FROBENIUS: *norm = 1.0; s = 0.0; break;
181:     case NORM_1_AND_2: norm[0] = 0.0; norm[1] = 1.0; s = 0.0; break;
182:   }
183:   for (i=0;i<as->n->n;i++) {
184:     if (as->x[0]->ops->norm_local) as->x[0]->ops->norm_local(as->x[i],t,work);
185:     else VecNorm(as->x[i],t,work);
186:     /* norm+= work */
187:     switch (t) {
188:       case NORM_1: *norm+= *work; break;
189:       case NORM_2: case NORM_FROBENIUS: AddNorm2(norm,&s,*work); break;
190:       case NORM_1_AND_2: norm[0]+= work[0]; AddNorm2(&norm[1],&s,work[1]); break;
191:       case NORM_INFINITY: *norm = PetscMax(*norm,*work); break;
192:     }
193:   }

195:   /* If def(__WITH_MPI__) and exists norm_local */
196: #if defined(__WITH_MPI__)
197:   if (as->x[0]->ops->norm_local) {
198:     PetscReal work0[3];
199:     /* norm <- Allreduce(work) */
200:     switch (t) {
201:     case NORM_1:
202:       work[0] = *norm;
203:       MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a));
204:       break;
205:     case NORM_2: case NORM_FROBENIUS:
206:       work[0] = *norm; work[1] = s;
207:       MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
208:       *norm = GetNorm2(work0[0],work0[1]);
209:       break;
210:     case NORM_1_AND_2:
211:       work[0] = norm[0]; work[1] = norm[1]; work[2] = s;
212:       MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a));
213:       norm[0] = work0[0];
214:       norm[1] = GetNorm2(work0[1],work0[2]);
215:       break;
216:     case NORM_INFINITY:
217:       work[0] = *norm;
218:       MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a));
219:       break;
220:     }
221:   }
222: #else
223:   /* Norm correction */
224:   switch (t) {
225:     case NORM_2: case NORM_FROBENIUS: *norm = GetNorm2(*norm,s); break;
226:     case NORM_1_AND_2: norm[1] = GetNorm2(norm[1],s); break;
227:     default: ;
228:   }
229: #endif
230:   PetscFunctionReturn(0);
231: }

233: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
234: {
235:   PetscScalar       dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
236:   const PetscScalar *vx,*wx;
237:   Vec_Comp          *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
238:   PetscInt          i,n;
239:   PetscBool         t0,t1;
240: #if defined(__WITH_MPI__)
241:   PetscScalar       work[4];
242: #endif

244:   /* Compute recursively the local part */
245:   PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0);
246:   PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1);
247:   if (t0 && t1) {
248:     SlepcValidVecComp(v,1);
249:     SlepcValidVecComp(w,2);
250:     for (i=0;i<vs->n->n;i++) {
251:       VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1);
252:       dp0 += dp1;
253:       nm0 += nm1;
254:     }
255:   } else if (!t0 && !t1) {
256:     VecGetLocalSize(v,&n);
257:     VecGetArrayRead(v,&vx);
258:     VecGetArrayRead(w,&wx);
259:     for (i=0;i<n;i++) {
260:       dp0 += vx[i]*PetscConj(wx[i]);
261:       nm0 += wx[i]*PetscConj(wx[i]);
262:     }
263:     VecRestoreArrayRead(v,&vx);
264:     VecRestoreArrayRead(w,&wx);
265:   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");

267: #if defined(__WITH_MPI__)
268:     /* [dp, nm] <- Allreduce([dp0, nm0]) */
269:     work[0] = dp0; work[1] = nm0;
270:     MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
271:     *dp = work[2]; *nm = work[3];
272: #else
273:     *dp = dp0; *nm = nm0;
274: #endif
275:   PetscFunctionReturn(0);
276: }