Actual source code: veccomp.c

slepc-main 2024-11-15
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(getlocaltoglobalmapping,NULL),
173:   PetscDesignatedInitializer(setvalueslocal,NULL),
174:   PetscDesignatedInitializer(resetarray,NULL),
175:   PetscDesignatedInitializer(setfromoptions,NULL),
176:   PetscDesignatedInitializer(maxpointwisedivide,VecMaxPointwiseDivide_Comp),
177:   PetscDesignatedInitializer(pointwisemax,VecPointwiseMax_Comp),
178:   PetscDesignatedInitializer(pointwisemaxabs,VecPointwiseMaxAbs_Comp),
179:   PetscDesignatedInitializer(pointwisemin,VecPointwiseMin_Comp),
180:   PetscDesignatedInitializer(getvalues,NULL),
181:   PetscDesignatedInitializer(sqrt,VecSqrtAbs_Comp),
182:   PetscDesignatedInitializer(abs,VecAbs_Comp),
183:   PetscDesignatedInitializer(exp,VecExp_Comp),
184:   PetscDesignatedInitializer(log,VecLog_Comp),
185:   PetscDesignatedInitializer(shift,VecShift_Comp),
186:   PetscDesignatedInitializer(create,NULL),
187:   PetscDesignatedInitializer(stridegather,NULL),
188:   PetscDesignatedInitializer(stridescatter,NULL),
189:   PetscDesignatedInitializer(dotnorm2,VecDotNorm2_Comp_MPI),
190:   PetscDesignatedInitializer(getsubvector,NULL),
191:   PetscDesignatedInitializer(restoresubvector,NULL),
192:   PetscDesignatedInitializer(getarrayread,NULL),
193:   PetscDesignatedInitializer(restorearrayread,NULL),
194:   PetscDesignatedInitializer(stridesubsetgather,NULL),
195:   PetscDesignatedInitializer(stridesubsetscatter,NULL),
196:   PetscDesignatedInitializer(viewnative,NULL),
197:   PetscDesignatedInitializer(loadnative,NULL),
198:   PetscDesignatedInitializer(getlocalvector,NULL)
199: };

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

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

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

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

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

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

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

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

251:   s->nx = nx;

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

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

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

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

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

294:    Collective

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

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

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

310:    Level: developer

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

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

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

335:    Collective

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

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

345:    Level: developer

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

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

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

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

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

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

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

395:    Collective

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

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

404:    Level: developer

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

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

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

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

454:    Collective

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

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

465:    Level: developer

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

724: __FUNC_TEMPLATE3__(Copy)
725: __FUNC_TEMPLATE3__(Swap)

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

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