Actual source code: veccomp.c

slepc-3.19.2 2023-09-05
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 <slepc/private/vecimplslepc.h>

 13: /* Private MPI datatypes and operators */
 14: static MPI_Datatype MPIU_NORM2=0, MPIU_NORM1_AND_2=0;
 15: static PetscBool VecCompInitialized = PETSC_FALSE;
 16: MPI_Op MPIU_NORM2_SUM=0;

 18: /* Private functions */
 19: static inline void SumNorm2(PetscReal*,PetscReal*,PetscReal*,PetscReal*);
 20: static inline PetscReal GetNorm2(PetscReal,PetscReal);
 21: static inline void AddNorm2(PetscReal*,PetscReal*,PetscReal);
 22: static PetscErrorCode VecCompSetSubVecs_Comp(Vec,PetscInt,Vec*);
 23: static PetscErrorCode VecCompGetSubVecs_Comp(Vec,PetscInt*,const Vec**);

 25: #include "veccomp0.h"

 28: #include "veccomp0.h"

 30: static inline void SumNorm2(PetscReal *ssq0,PetscReal *scale0,PetscReal *ssq1,PetscReal *scale1)
 31: {
 32:   PetscReal q;
 33:   if (*scale0 > *scale1) {
 34:     q = *scale1/(*scale0);
 35:     *ssq1 = *ssq0 + q*q*(*ssq1);
 36:     *scale1 = *scale0;
 37:   } else {
 38:     q = *scale0/(*scale1);
 39:     *ssq1 += q*q*(*ssq0);
 40:   }
 41: }

 43: static inline PetscReal GetNorm2(PetscReal ssq,PetscReal scale)
 44: {
 45:   return scale*PetscSqrtReal(ssq);
 46: }

 48: static inline void AddNorm2(PetscReal *ssq,PetscReal *scale,PetscReal x)
 49: {
 50:   PetscReal absx,q;
 51:   if (x != 0.0) {
 52:     absx = PetscAbs(x);
 53:     if (*scale < absx) {
 54:       q = *scale/absx;
 55:       *ssq = 1.0 + *ssq*q*q;
 56:       *scale = absx;
 57:     } else {
 58:       q = absx/(*scale);
 59:       *ssq += q*q;
 60:     }
 61:   }
 62: }

 64: SLEPC_EXTERN void MPIAPI SlepcSumNorm2_Local(void *in,void *out,PetscMPIInt *cnt,MPI_Datatype *datatype)
 65: {
 66:   PetscInt  i,count = *cnt;
 67:   PetscReal *xin = (PetscReal*)in,*xout = (PetscReal*)out;

 69:   PetscFunctionBegin;
 70:   if (*datatype == MPIU_NORM2) {
 71:     for (i=0;i<count;i++) {
 72:       SumNorm2(&xin[i*2],&xin[i*2+1],&xout[i*2],&xout[i*2+1]);
 73:     }
 74:   } else if (*datatype == MPIU_NORM1_AND_2) {
 75:     for (i=0;i<count;i++) {
 76:       xout[i*3] += xin[i*3];
 77:       SumNorm2(&xin[i*3+1],&xin[i*3+2],&xout[i*3+1],&xout[i*3+2]);
 78:     }
 79:   } else {
 80:     (void)(*PetscErrorPrintf)("Can only handle MPIU_NORM* data types");
 81:     MPI_Abort(PETSC_COMM_WORLD,1);
 82:   }
 83:   PetscFunctionReturnVoid();
 84: }

 86: static PetscErrorCode VecCompNormEnd(void)
 87: {
 88:   PetscFunctionBegin;
 89:   PetscCallMPI(MPI_Type_free(&MPIU_NORM2));
 90:   PetscCallMPI(MPI_Type_free(&MPIU_NORM1_AND_2));
 91:   PetscCallMPI(MPI_Op_free(&MPIU_NORM2_SUM));
 92:   VecCompInitialized = PETSC_FALSE;
 93:   PetscFunctionReturn(PETSC_SUCCESS);
 94: }

 96: static PetscErrorCode VecCompNormInit(void)
 97: {
 98:   PetscFunctionBegin;
 99:   PetscCallMPI(MPI_Type_contiguous(2,MPIU_REAL,&MPIU_NORM2));
100:   PetscCallMPI(MPI_Type_commit(&MPIU_NORM2));
101:   PetscCallMPI(MPI_Type_contiguous(3,MPIU_REAL,&MPIU_NORM1_AND_2));
102:   PetscCallMPI(MPI_Type_commit(&MPIU_NORM1_AND_2));
103:   PetscCallMPI(MPI_Op_create(SlepcSumNorm2_Local,PETSC_TRUE,&MPIU_NORM2_SUM));
104:   PetscCall(PetscRegisterFinalize(VecCompNormEnd));
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }

108: PetscErrorCode VecDestroy_Comp(Vec v)
109: {
110:   Vec_Comp       *vs = (Vec_Comp*)v->data;
111:   PetscInt       i;

113:   PetscFunctionBegin;
114: #if defined(PETSC_USE_LOG)
115:   PetscCall(PetscLogObjectState((PetscObject)v,"Length=%" PetscInt_FMT,v->map->n));
116: #endif
117:   for (i=0;i<vs->nx;i++) PetscCall(VecDestroy(&vs->x[i]));
118:   if (--vs->n->friends <= 0) PetscCall(PetscFree(vs->n));
119:   PetscCall(PetscFree(vs->x));
120:   PetscCall(PetscFree(vs));
121:   PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",NULL));
122:   PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",NULL));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: static struct _VecOps DvOps = {
127:   PetscDesignatedInitializer(duplicate,VecDuplicate_Comp),
128:   PetscDesignatedInitializer(duplicatevecs,VecDuplicateVecs_Comp),
129:   PetscDesignatedInitializer(destroyvecs,VecDestroyVecs_Comp),
130:   PetscDesignatedInitializer(dot,VecDot_Comp_MPI),
131:   PetscDesignatedInitializer(mdot,VecMDot_Comp_MPI),
132:   PetscDesignatedInitializer(norm,VecNorm_Comp_MPI),
133:   PetscDesignatedInitializer(tdot,VecTDot_Comp_MPI),
134:   PetscDesignatedInitializer(mtdot,VecMTDot_Comp_MPI),
135:   PetscDesignatedInitializer(scale,VecScale_Comp),
136:   PetscDesignatedInitializer(copy,VecCopy_Comp),
137:   PetscDesignatedInitializer(set,VecSet_Comp),
138:   PetscDesignatedInitializer(swap,VecSwap_Comp),
139:   PetscDesignatedInitializer(axpy,VecAXPY_Comp),
140:   PetscDesignatedInitializer(axpby,VecAXPBY_Comp),
141:   PetscDesignatedInitializer(maxpy,VecMAXPY_Comp),
142:   PetscDesignatedInitializer(aypx,VecAYPX_Comp),
143:   PetscDesignatedInitializer(waxpy,VecWAXPY_Comp),
144:   PetscDesignatedInitializer(axpbypcz,VecAXPBYPCZ_Comp),
145:   PetscDesignatedInitializer(pointwisemult,VecPointwiseMult_Comp),
146:   PetscDesignatedInitializer(pointwisedivide,VecPointwiseDivide_Comp),
147:   PetscDesignatedInitializer(setvalues,NULL),
148:   PetscDesignatedInitializer(assemblybegin,NULL),
149:   PetscDesignatedInitializer(assemblyend,NULL),
150:   PetscDesignatedInitializer(getarray,NULL),
151:   PetscDesignatedInitializer(getsize,VecGetSize_Comp),
152:   PetscDesignatedInitializer(getlocalsize,VecGetLocalSize_Comp),
153:   PetscDesignatedInitializer(restorearray,NULL),
154:   PetscDesignatedInitializer(max,VecMax_Comp),
155:   PetscDesignatedInitializer(min,VecMin_Comp),
156:   PetscDesignatedInitializer(setrandom,VecSetRandom_Comp),
157:   PetscDesignatedInitializer(setoption,NULL),
158:   PetscDesignatedInitializer(setvaluesblocked,NULL),
159:   PetscDesignatedInitializer(destroy,VecDestroy_Comp),
160:   PetscDesignatedInitializer(view,VecView_Comp),
161:   PetscDesignatedInitializer(placearray,NULL),
162:   PetscDesignatedInitializer(replacearray,NULL),
163:   PetscDesignatedInitializer(dot_local,VecDot_Comp_Seq),
164:   PetscDesignatedInitializer(tdot_local,VecTDot_Comp_Seq),
165:   PetscDesignatedInitializer(norm_local,VecNorm_Comp_Seq),
166:   PetscDesignatedInitializer(mdot_local,VecMDot_Comp_Seq),
167:   PetscDesignatedInitializer(mtdot_local,VecMTDot_Comp_Seq),
168:   PetscDesignatedInitializer(load,NULL),
169:   PetscDesignatedInitializer(reciprocal,VecReciprocal_Comp),
170:   PetscDesignatedInitializer(conjugate,VecConjugate_Comp),
171:   PetscDesignatedInitializer(setlocaltoglobalmapping,NULL),
172:   PetscDesignatedInitializer(setvalueslocal,NULL),
173:   PetscDesignatedInitializer(resetarray,NULL),
174:   PetscDesignatedInitializer(setfromoptions,NULL),
175:   PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Comp),
176:   PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Comp),
177:   PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Comp),
178:   PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Comp),
179:   PetscDesignatedInitializer(getvalues,NULL),
180:   PetscDesignatedInitializer(sqrt,VecSqrtAbs_Comp),
181:   PetscDesignatedInitializer(abs,VecAbs_Comp),
182:   PetscDesignatedInitializer(exp,VecExp_Comp),
183:   PetscDesignatedInitializer(log,VecLog_Comp),
184:   PetscDesignatedInitializer(shift,VecShift_Comp),
185:   PetscDesignatedInitializer(create,NULL),
186:   PetscDesignatedInitializer(stridegather,NULL),
187:   PetscDesignatedInitializer(stridescatter,NULL),
188:   PetscDesignatedInitializer(dotnorm2,VecDotNorm2_Comp_MPI),
189:   PetscDesignatedInitializer(getsubvector,NULL),
190:   PetscDesignatedInitializer(restoresubvector,NULL),
191:   PetscDesignatedInitializer(getarrayread,NULL),
192:   PetscDesignatedInitializer(restorearrayread,NULL),
193:   PetscDesignatedInitializer(stridesubsetgather,NULL),
194:   PetscDesignatedInitializer(stridesubsetscatter,NULL),
195:   PetscDesignatedInitializer(viewnative,NULL),
196:   PetscDesignatedInitializer(loadnative,NULL),
197:   PetscDesignatedInitializer(getlocalvector,NULL)
198: };

200: PetscErrorCode VecDuplicateVecs_Comp(Vec w,PetscInt m,Vec *V[])
201: {
202:   PetscInt       i;

204:   PetscFunctionBegin;
207:   PetscCheck(m>0,PetscObjectComm((PetscObject)w),PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
208:   PetscCall(PetscMalloc1(m,V));
209:   for (i=0;i<m;i++) PetscCall(VecDuplicate(w,*V+i));
210:   PetscFunctionReturn(PETSC_SUCCESS);
211: }

213: PetscErrorCode VecDestroyVecs_Comp(PetscInt m,Vec v[])
214: {
215:   PetscInt       i;

217:   PetscFunctionBegin;
219:   PetscCheck(m>0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %" PetscInt_FMT,m);
220:   for (i=0;i<m;i++) PetscCall(VecDestroy(&v[i]));
221:   PetscCall(PetscFree(v));
222:   PetscFunctionReturn(PETSC_SUCCESS);
223: }

225: static PetscErrorCode VecCreate_Comp_Private(Vec v,Vec *x,PetscInt nx,PetscBool x_to_me,Vec_Comp_N *n)
226: {
227:   Vec_Comp       *s;
228:   PetscInt       N=0,lN=0,i,k;

230:   PetscFunctionBegin;
231:   if (!VecCompInitialized) {
232:     VecCompInitialized = PETSC_TRUE;
233:     PetscCall(VecRegister(VECCOMP,VecCreate_Comp));
234:     PetscCall(VecCompNormInit());
235:   }

237:   /* Allocate a new Vec_Comp */
238:   if (v->data) PetscCall(PetscFree(v->data));
239:   PetscCall(PetscNew(&s));
240:   PetscCall(PetscMemcpy(v->ops,&DvOps,sizeof(DvOps)));
241:   v->data        = (void*)s;
242:   v->petscnative = PETSC_FALSE;

244:   /* Allocate the array of Vec, if it is needed to be done */
245:   if (!x_to_me) {
246:     if (nx) PetscCall(PetscMalloc1(nx,&s->x));
247:     if (x) PetscCall(PetscArraycpy(s->x,x,nx));
248:   } else s->x = x;

250:   s->nx = nx;

252:   if (nx && x) {
253:     /* Allocate the shared structure, if it is not given */
254:     if (!n) {
255:       for (i=0;i<nx;i++) {
256:         PetscCall(VecGetSize(x[i],&k));
257:         N+= k;
258:         PetscCall(VecGetLocalSize(x[i],&k));
259:         lN+= k;
260:       }
261:       PetscCall(PetscNew(&n));
262:       s->n = n;
263:       n->n = nx;
264:       n->N = N;
265:       n->lN = lN;
266:       n->friends = 1;
267:     } else { /* If not, check in the vector in the shared structure */
268:       s->n = n;
269:       s->n->friends++;
270:     }

272:     /* Set the virtual sizes as the real sizes of the vector */
273:     PetscCall(VecSetSizes(v,s->n->lN,s->n->N));
274:   }

276:   PetscCall(PetscObjectChangeTypeName((PetscObject)v,VECCOMP));
277:   PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",VecCompSetSubVecs_Comp));
278:   PetscCall(PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",VecCompGetSubVecs_Comp));
279:   PetscFunctionReturn(PETSC_SUCCESS);
280: }

282: SLEPC_EXTERN PetscErrorCode VecCreate_Comp(Vec V)
283: {
284:   PetscFunctionBegin;
285:   PetscCall(VecCreate_Comp_Private(V,NULL,0,PETSC_FALSE,NULL));
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: /*@C
290:    VecCreateComp - Creates a new vector containing several subvectors,
291:    each stored separately.

293:    Collective

295:    Input Parameters:
296: +  comm - communicator for the new Vec
297: .  Nx   - array of (initial) global sizes of child vectors
298: .  n    - number of child vectors
299: .  t    - type of the child vectors
300: -  Vparent - (optional) template vector

302:    Output Parameter:
303: .  V - new vector

305:    Notes:
306:    This is similar to PETSc's VecNest but customized for SLEPc's needs. In particular,
307:    the number of child vectors can be modified dynamically, with VecCompSetSubVecs().

309:    Level: developer

311: .seealso: VecCreateCompWithVecs(), VecCompSetSubVecs()
312: @*/
313: PetscErrorCode VecCreateComp(MPI_Comm comm,PetscInt *Nx,PetscInt n,VecType t,Vec Vparent,Vec *V)
314: {
315:   Vec            *x;
316:   PetscInt       i;

318:   PetscFunctionBegin;
319:   PetscCall(VecCreate(comm,V));
320:   PetscCall(PetscMalloc1(n,&x));
321:   for (i=0;i<n;i++) {
322:     PetscCall(VecCreate(comm,&x[i]));
323:     PetscCall(VecSetSizes(x[i],PETSC_DECIDE,Nx[i]));
324:     PetscCall(VecSetType(x[i],t));
325:   }
326:   PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
327:   PetscFunctionReturn(PETSC_SUCCESS);
328: }

330: /*@C
331:    VecCreateCompWithVecs - Creates a new vector containing several subvectors,
332:    each stored separately, from an array of Vecs.

334:    Collective

336:    Input Parameters:
337: +  x - array of Vecs
338: .  n - number of child vectors
339: -  Vparent - (optional) template vector

341:    Output Parameter:
342: .  V - new vector

344:    Level: developer

346: .seealso: VecCreateComp()
347: @*/
348: PetscErrorCode VecCreateCompWithVecs(Vec *x,PetscInt n,Vec Vparent,Vec *V)
349: {
350:   PetscInt       i;

352:   PetscFunctionBegin;
356:   PetscCall(VecCreate(PetscObjectComm((PetscObject)x[0]),V));
357:   for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)x[i]));
358:   PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
359:   PetscFunctionReturn(PETSC_SUCCESS);
360: }

362: PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
363: {
364:   Vec            *x;
365:   PetscInt       i;
366:   Vec_Comp       *s = (Vec_Comp*)win->data;

368:   PetscFunctionBegin;
369:   SlepcValidVecComp(win,1);
370:   PetscCall(VecCreate(PetscObjectComm((PetscObject)win),V));
371:   PetscCall(PetscMalloc1(s->nx,&x));
372:   for (i=0;i<s->nx;i++) {
373:     if (s->x[i]) PetscCall(VecDuplicate(s->x[i],&x[i]));
374:     else x[i] = NULL;
375:   }
376:   PetscCall(VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n));
377:   PetscFunctionReturn(PETSC_SUCCESS);
378: }

380: static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
381: {
382:   Vec_Comp *s = (Vec_Comp*)win->data;

384:   PetscFunctionBegin;
385:   if (x) *x = s->x;
386:   if (n) *n = s->n->n;
387:   PetscFunctionReturn(PETSC_SUCCESS);
388: }

390: /*@C
391:    VecCompGetSubVecs - Returns the entire array of vectors defining a
392:    compound vector.

394:    Collective

396:    Input Parameter:
397: .  win - compound vector

399:    Output Parameters:
400: +  n - number of child vectors
401: -  x - array of child vectors

403:    Level: developer

405: .seealso: VecCreateComp()
406: @*/
407: PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec **x)
408: {
409:   PetscFunctionBegin;
411:   PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
412:   PetscFunctionReturn(PETSC_SUCCESS);
413: }

415: static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
416: {
417:   Vec_Comp       *s = (Vec_Comp*)win->data;
418:   PetscInt       i,N,nlocal;
419:   Vec_Comp_N     *nn;

421:   PetscFunctionBegin;
422:   PetscCheck(s,PetscObjectComm((PetscObject)win),PETSC_ERR_ORDER,"Must call VecSetSizes first");
423:   if (!s->nx) {
424:     /* vector has been created via VecCreate+VecSetType+VecSetSizes, so allocate data structures */
425:     PetscCall(PetscMalloc1(n,&s->x));
426:     PetscCall(VecGetSize(win,&N));
427:     PetscCheck(N%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Global dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,N,n);
428:     PetscCall(VecGetLocalSize(win,&nlocal));
429:     PetscCheck(nlocal%n==0,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Local dimension %" PetscInt_FMT " is not divisible by %" PetscInt_FMT,nlocal,n);
430:     s->nx = n;
431:     for (i=0;i<n;i++) {
432:       PetscCall(VecCreate(PetscObjectComm((PetscObject)win),&s->x[i]));
433:       PetscCall(VecSetSizes(s->x[i],nlocal/n,N/n));
434:       PetscCall(VecSetFromOptions(s->x[i]));
435:     }
436:     if (!s->n) {
437:       PetscCall(PetscNew(&nn));
438:       s->n = nn;
439:       nn->N = N;
440:       nn->lN = nlocal;
441:       nn->friends = 1;
442:     }
443:   } else PetscCheck(n<=s->nx,PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Number of child vectors cannot be larger than %" PetscInt_FMT,s->nx);
444:   if (x) PetscCall(PetscArraycpy(s->x,x,n));
445:   s->n->n = n;
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: /*@C
450:    VecCompSetSubVecs - Resets the number of subvectors defining a compound vector,
451:    or replaces the subvectors.

453:    Collective

455:    Input Parameters:
456: +  win - compound vector
457: .  n - number of child vectors
458: -  x - array of child vectors

460:    Note:
461:    It is not possible to increase the number of subvectors with respect to the
462:    number set at its creation.

464:    Level: developer

466: .seealso: VecCreateComp(), VecCompGetSubVecs()
467: @*/
468: PetscErrorCode VecCompSetSubVecs(Vec win,PetscInt n,Vec *x)
469: {
470:   PetscFunctionBegin;
473:   PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
474:   PetscFunctionReturn(PETSC_SUCCESS);
475: }

477: PetscErrorCode VecAXPY_Comp(Vec v,PetscScalar alpha,Vec w)
478: {
479:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
480:   PetscInt       i;

482:   PetscFunctionBegin;
483:   SlepcValidVecComp(v,1);
484:   SlepcValidVecComp(w,3);
485:   for (i=0;i<vs->n->n;i++) PetscCall(VecAXPY(vs->x[i],alpha,ws->x[i]));
486:   PetscFunctionReturn(PETSC_SUCCESS);
487: }

489: PetscErrorCode VecAYPX_Comp(Vec v,PetscScalar alpha,Vec w)
490: {
491:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
492:   PetscInt       i;

494:   PetscFunctionBegin;
495:   SlepcValidVecComp(v,1);
496:   SlepcValidVecComp(w,3);
497:   for (i=0;i<vs->n->n;i++) PetscCall(VecAYPX(vs->x[i],alpha,ws->x[i]));
498:   PetscFunctionReturn(PETSC_SUCCESS);
499: }

501: PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,Vec w)
502: {
503:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
504:   PetscInt       i;

506:   PetscFunctionBegin;
507:   SlepcValidVecComp(v,1);
508:   SlepcValidVecComp(w,4);
509:   for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBY(vs->x[i],alpha,beta,ws->x[i]));
510:   PetscFunctionReturn(PETSC_SUCCESS);
511: }

513: PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
514: {
515:   Vec_Comp       *vs = (Vec_Comp*)v->data;
516:   Vec            *wx;
517:   PetscInt       i,j;

519:   PetscFunctionBegin;
520:   SlepcValidVecComp(v,1);
521:   for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);

523:   PetscCall(PetscMalloc1(n,&wx));

525:   for (j=0;j<vs->n->n;j++) {
526:     for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
527:     PetscCall(VecMAXPY(vs->x[j],n,alpha,wx));
528:   }

530:   PetscCall(PetscFree(wx));
531:   PetscFunctionReturn(PETSC_SUCCESS);
532: }

534: PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,Vec w,Vec z)
535: {
536:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
537:   PetscInt       i;

539:   PetscFunctionBegin;
540:   SlepcValidVecComp(v,1);
541:   SlepcValidVecComp(w,3);
542:   SlepcValidVecComp(z,4);
543:   for (i=0;i<vs->n->n;i++) PetscCall(VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]));
544:   PetscFunctionReturn(PETSC_SUCCESS);
545: }

547: PetscErrorCode VecAXPBYPCZ_Comp(Vec v,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec w,Vec z)
548: {
549:   Vec_Comp        *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
550:   PetscInt        i;

552:   PetscFunctionBegin;
553:   SlepcValidVecComp(v,1);
554:   SlepcValidVecComp(w,5);
555:   SlepcValidVecComp(z,6);
556:   for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]));
557:   PetscFunctionReturn(PETSC_SUCCESS);
558: }

560: PetscErrorCode VecGetSize_Comp(Vec v,PetscInt *size)
561: {
562:   Vec_Comp *vs = (Vec_Comp*)v->data;

564:   PetscFunctionBegin;
566:   if (vs->n) {
567:     SlepcValidVecComp(v,1);
568:     *size = vs->n->N;
569:   } else *size = v->map->N;
570:   PetscFunctionReturn(PETSC_SUCCESS);
571: }

573: PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
574: {
575:   Vec_Comp *vs = (Vec_Comp*)v->data;

577:   PetscFunctionBegin;
579:   if (vs->n) {
580:     SlepcValidVecComp(v,1);
581:     *size = vs->n->lN;
582:   } else *size = v->map->n;
583:   PetscFunctionReturn(PETSC_SUCCESS);
584: }

586: PetscErrorCode VecMax_Comp(Vec v,PetscInt *idx,PetscReal *z)
587: {
588:   Vec_Comp       *vs = (Vec_Comp*)v->data;
589:   PetscInt       idxp,s=0,s0;
590:   PetscReal      zp,z0;
591:   PetscInt       i;

593:   PetscFunctionBegin;
594:   SlepcValidVecComp(v,1);
595:   if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);

597:   if (vs->n->n > 0) PetscCall(VecMax(vs->x[0],idx?&idxp:NULL,&zp));
598:   else {
599:     zp = PETSC_MIN_REAL;
600:     if (idx) idxp = -1;
601:   }
602:   for (i=1;i<vs->n->n;i++) {
603:     PetscCall(VecGetSize(vs->x[i-1],&s0));
604:     s += s0;
605:     PetscCall(VecMax(vs->x[i],idx?&idxp:NULL,&z0));
606:     if (zp < z0) {
607:       if (idx) *idx = s+idxp;
608:       zp = z0;
609:     }
610:   }
611:   if (z) *z = zp;
612:   PetscFunctionReturn(PETSC_SUCCESS);
613: }

615: PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
616: {
617:   Vec_Comp       *vs = (Vec_Comp*)v->data;
618:   PetscInt       idxp,s=0,s0;
619:   PetscReal      zp,z0;
620:   PetscInt       i;

622:   PetscFunctionBegin;
623:   SlepcValidVecComp(v,1);
624:   if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);

626:   if (vs->n->n > 0) PetscCall(VecMin(vs->x[0],idx?&idxp:NULL,&zp));
627:   else {
628:     zp = PETSC_MAX_REAL;
629:     if (idx) idxp = -1;
630:   }
631:   for (i=1;i<vs->n->n;i++) {
632:     PetscCall(VecGetSize(vs->x[i-1],&s0));
633:     s += s0;
634:     PetscCall(VecMin(vs->x[i],idx?&idxp:NULL,&z0));
635:     if (zp > z0) {
636:       if (idx) *idx = s+idxp;
637:       zp = z0;
638:     }
639:   }
640:   if (z) *z = zp;
641:   PetscFunctionReturn(PETSC_SUCCESS);
642: }

644: PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
645: {
646:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
647:   PetscReal      work;
648:   PetscInt       i;

650:   PetscFunctionBegin;
651:   SlepcValidVecComp(v,1);
652:   SlepcValidVecComp(w,2);
653:   if (!m || vs->n->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
654:   PetscCall(VecMaxPointwiseDivide(vs->x[0],ws->x[0],m));
655:   for (i=1;i<vs->n->n;i++) {
656:     PetscCall(VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work));
657:     *m = PetscMax(*m,work);
658:   }
659:   PetscFunctionReturn(PETSC_SUCCESS);
660: }


667: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v) \
668: { \
669:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
670:   PetscInt        i; \
671: \
672:   PetscFunctionBegin; \
673:   SlepcValidVecComp(v,1); \
674:   for (i=0;i<vs->n->n;i++) { \
675:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i])); \
676:   } \
677:   PetscFunctionReturn(PETSC_SUCCESS);\
678: }

680: __FUNC_TEMPLATE1__(Conjugate)
681: __FUNC_TEMPLATE1__(Reciprocal)
682: __FUNC_TEMPLATE1__(SqrtAbs)
683: __FUNC_TEMPLATE1__(Abs)
684: __FUNC_TEMPLATE1__(Exp)
685: __FUNC_TEMPLATE1__(Log)

688: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,T0 __a) \
689: { \
690:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
691:   PetscInt        i; \
692: \
693:   PetscFunctionBegin; \
694:   SlepcValidVecComp(v,1); \
695:   for (i=0;i<vs->n->n;i++) { \
696:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],__a)); \
697:   } \
698:   PetscFunctionReturn(PETSC_SUCCESS);\
699: }

701: __FUNC_TEMPLATE2__(Set,PetscScalar)
702: __FUNC_TEMPLATE2__(View,PetscViewer)
703: __FUNC_TEMPLATE2__(Scale,PetscScalar)
704: __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
705: __FUNC_TEMPLATE2__(Shift,PetscScalar)

708: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w) \
709: { \
710:   Vec_Comp        *vs = (Vec_Comp*)v->data,\
711:                   *ws = (Vec_Comp*)w->data; \
712:   PetscInt        i; \
713: \
714:   PetscFunctionBegin; \
715:   SlepcValidVecComp(v,1); \
716:   SlepcValidVecComp(w,2); \
717:   for (i=0;i<vs->n->n;i++) { \
718:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i])); \
719:   } \
720:   PetscFunctionReturn(PETSC_SUCCESS);\
721: }

723: __FUNC_TEMPLATE3__(Copy)
724: __FUNC_TEMPLATE3__(Swap)

727: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w,Vec z) \
728: { \
729:   Vec_Comp        *vs = (Vec_Comp*)v->data, \
730:                   *ws = (Vec_Comp*)w->data, \
731:                   *zs = (Vec_Comp*)z->data; \
732:   PetscInt        i; \
733: \
734:   PetscFunctionBegin; \
735:   SlepcValidVecComp(v,1); \
736:   SlepcValidVecComp(w,2); \
737:   SlepcValidVecComp(z,3); \
738:   for (i=0;i<vs->n->n;i++) { \
739:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i],zs->x[i])); \
740:   } \
741:   PetscFunctionReturn(PETSC_SUCCESS);\
742: }

744: __FUNC_TEMPLATE4__(PointwiseMax)
745: __FUNC_TEMPLATE4__(PointwiseMaxAbs)
746: __FUNC_TEMPLATE4__(PointwiseMin)
747: __FUNC_TEMPLATE4__(PointwiseMult)
748: __FUNC_TEMPLATE4__(PointwiseDivide)