Actual source code: veccomp.c

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 <slepc/private/vecimplslepc.h>     /*I "slepcvec.h" I*/

 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: PETSC_STATIC_INLINE void SumNorm2(PetscReal*,PetscReal*,PetscReal*,PetscReal*);
 20: PETSC_STATIC_INLINE PetscReal GetNorm2(PetscReal,PetscReal);
 21: PETSC_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: PETSC_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: PETSC_STATIC_INLINE PetscReal GetNorm2(PetscReal ssq,PetscReal scale)
 44: {
 45:   return scale*PetscSqrtReal(ssq);
 46: }

 48: PETSC_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 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;

 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:     (*PetscErrorPrintf)("Can only handle MPIU_NORM* data types");
 81:     MPI_Abort(MPI_COMM_WORLD,1);
 82:   }
 83:   PetscFunctionReturnVoid();
 84: }

 86: static PetscErrorCode VecCompNormEnd(void)
 87: {

 91:   MPI_Type_free(&MPIU_NORM2);
 92:   MPI_Type_free(&MPIU_NORM1_AND_2);
 93:   MPI_Op_free(&MPIU_NORM2_SUM);
 94:   VecCompInitialized = PETSC_FALSE;
 95:   return(0);
 96: }

 98: static PetscErrorCode VecCompNormInit(void)
 99: {

103:   MPI_Type_contiguous(2,MPIU_REAL,&MPIU_NORM2);
104:   MPI_Type_commit(&MPIU_NORM2);
105:   MPI_Type_contiguous(3,MPIU_REAL,&MPIU_NORM1_AND_2);
106:   MPI_Type_commit(&MPIU_NORM1_AND_2);
107:   MPI_Op_create(SlepcSumNorm2_Local,PETSC_TRUE,&MPIU_NORM2_SUM);
108:   PetscRegisterFinalize(VecCompNormEnd);
109:   return(0);
110: }

112: PetscErrorCode VecDestroy_Comp(Vec v)
113: {
114:   Vec_Comp       *vs = (Vec_Comp*)v->data;
115:   PetscInt       i;

119: #if defined(PETSC_USE_LOG)
120:   PetscLogObjectState((PetscObject)v,"Length=%D",v->map->n);
121: #endif
122:   for (i=0;i<vs->nx;i++) {
123:     VecDestroy(&vs->x[i]);
124:   }
125:   if (--vs->n->friends <= 0) {
126:     PetscFree(vs->n);
127:   }
128:   PetscFree(vs->x);
129:   PetscFree(vs);
130:   PetscObjectComposeFunction((PetscObject)v,"VecCompSetSubVecs_C",NULL);
131:   PetscObjectComposeFunction((PetscObject)v,"VecCompGetSubVecs_C",NULL);
132:   return(0);
133: }

135: static struct _VecOps DvOps = {VecDuplicate_Comp, /* 1 */
136:             VecDuplicateVecs_Comp,
137:             VecDestroyVecs_Comp,
138:             VecDot_Comp_MPI,
139:             VecMDot_Comp_MPI,
140:             VecNorm_Comp_MPI,
141:             VecTDot_Comp_MPI,
142:             VecMTDot_Comp_MPI,
143:             VecScale_Comp,
144:             VecCopy_Comp, /* 10 */
145:             VecSet_Comp,
146:             VecSwap_Comp,
147:             VecAXPY_Comp,
148:             VecAXPBY_Comp,
149:             VecMAXPY_Comp,
150:             VecAYPX_Comp,
151:             VecWAXPY_Comp,
152:             VecAXPBYPCZ_Comp,
153:             VecPointwiseMult_Comp,
154:             VecPointwiseDivide_Comp,
155:             0, /* 20 */
156:             0,0,
157:             0 /*VecGetArray_Seq*/,
158:             VecGetSize_Comp,
159:             VecGetLocalSize_Comp,
160:             0/*VecRestoreArray_Seq*/,
161:             VecMax_Comp,
162:             VecMin_Comp,
163:             VecSetRandom_Comp,
164:             0, /* 30 */
165:             0,
166:             VecDestroy_Comp,
167:             VecView_Comp,
168:             0/*VecPlaceArray_Seq*/,
169:             0/*VecReplaceArray_Seq*/,
170:             VecDot_Comp_Seq,
171:             VecTDot_Comp_Seq,
172:             VecNorm_Comp_Seq,
173:             VecMDot_Comp_Seq,
174:             VecMTDot_Comp_Seq, /* 40 */
175:             0,
176:             VecReciprocal_Comp,
177:             VecConjugate_Comp,
178:             0,0,
179:             0/*VecResetArray_Seq*/,
180:             0,
181:             VecMaxPointwiseDivide_Comp,
182:             VecPointwiseMax_Comp,
183:             VecPointwiseMaxAbs_Comp,
184:             VecPointwiseMin_Comp,
185:             0,
186:             VecSqrtAbs_Comp,
187:             VecAbs_Comp,
188:             VecExp_Comp,
189:             VecLog_Comp,
190:             VecShift_Comp,
191:             0,
192:             0,
193:             0,
194:             VecDotNorm2_Comp_MPI
195:           };

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

205:   if (m<=0) SETERRQ1(PetscObjectComm((PetscObject)w),PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
206:   PetscMalloc1(m,V);
207:   for (i=0;i<m;i++) { VecDuplicate(w,*V+i); }
208:   return(0);
209: }

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

218:   if (m<=0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
219:   for (i=0;i<m;i++) { VecDestroy(&v[i]); }
220:   PetscFree(v);
221:   return(0);
222: }

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

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

237:   /* Allocate a new Vec_Comp */
238:   if (v->data) { PetscFree(v->data); }
239:   PetscNewLog(v,&s);
240:   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) { PetscMalloc1(nx,&s->x); }
247:     if (x) { PetscMemcpy(s->x,x,sizeof(Vec)*nx); }
248:   } else s->x = x;

250:   s->nx = nx;

252:   if (nx) {
253:     /* Allocate the shared structure, if it is not given */
254:     if (!n) {
255:       for (i=0;i<nx;i++) {
256:         VecGetSize(x[i],&k);
257:         N+= k;
258:         VecGetLocalSize(x[i],&k);
259:         lN+= k;
260:       }
261:       PetscNewLog(v,&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:     VecSetSizes(v,s->n->lN,s->n->N);
274:   }

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

282: SLEPC_EXTERN PetscErrorCode VecCreate_Comp(Vec V)
283: {

287:   VecCreate_Comp_Private(V,NULL,0,PETSC_FALSE,NULL);
288:   return(0);
289: }

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

295:    Collective on Vec

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

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

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

311:    Level: developer

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

322:   VecCreate(comm,V);
323:   PetscMalloc1(n,&x);
324:   PetscLogObjectMemory((PetscObject)*V,n*sizeof(Vec));
325:   for (i=0;i<n;i++) {
326:     VecCreate(comm,&x[i]);
327:     VecSetSizes(x[i],PETSC_DECIDE,Nx[i]);
328:     VecSetType(x[i],t);
329:   }
330:   VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL);
331:   return(0);
332: }

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

338:    Collective on Vec

340:    Input Parameters:
341: +  x - array of Vecs
342: .  n - number of child vectors
343: -  Vparent - (optional) template vector

345:    Output Parameter:
346: .  V - new vector

348:    Level: developer

350: .seealso: VecCreateComp()
351: @*/
352: PetscErrorCode VecCreateCompWithVecs(Vec *x,PetscInt n,Vec Vparent,Vec *V)
353: {
355:   PetscInt       i;

358:   VecCreate(PetscObjectComm((PetscObject)x[0]),V);
359:   for (i=0;i<n;i++) {
360:     PetscObjectReference((PetscObject)x[i]);
361:   }
362:   VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL);
363:   return(0);
364: }

366: PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
367: {
369:   Vec            *x;
370:   PetscInt       i;
371:   Vec_Comp       *s = (Vec_Comp*)win->data;

374:   SlepcValidVecComp(win,1);
375:   VecCreate(PetscObjectComm((PetscObject)win),V);
376:   PetscMalloc1(s->nx,&x);
377:   PetscLogObjectMemory((PetscObject)*V,s->nx*sizeof(Vec));
378:   for (i=0;i<s->nx;i++) {
379:     if (s->x[i]) {
380:       VecDuplicate(s->x[i],&x[i]);
381:     } else x[i] = NULL;
382:   }
383:   VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n);
384:   return(0);
385: }

387: static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
388: {
389:   Vec_Comp *s = (Vec_Comp*)win->data;

392:   if (x) *x = s->x;
393:   if (n) *n = s->n->n;
394:   return(0);
395: }

397: /*@C
398:    VecCompGetSubVecs - Returns the entire array of vectors defining a
399:    compound vector.

401:    Collective on Vec

403:    Input Parameter:
404: .  win - compound vector

406:    Output Parameters:
407: +  n - number of child vectors
408: -  x - array of child vectors

410:    Level: developer

412: .seealso: VecCreateComp()
413: @*/
414: PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec **x)
415: {

420:   PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
421:   return(0);
422: }

424: static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
425: {
426:   Vec_Comp       *s = (Vec_Comp*)win->data;
427:   PetscInt       i,N,nlocal;
428:   Vec_Comp_N     *nn;

432:   if (!s) SETERRQ(PetscObjectComm((PetscObject)win),PETSC_ERR_ORDER,"Must call VecSetSizes first");
433:   if (!s->nx) {
434:     /* vector has been created via VecCreate+VecSetType+VecSetSizes, so allocate data structures */
435:     PetscMalloc1(n,&s->x);
436:     PetscLogObjectMemory((PetscObject)win,n*sizeof(Vec));
437:     VecGetSize(win,&N);
438:     if (N%n) SETERRQ2(PetscObjectComm((PetscObject)win),1,"Global dimension %D is not divisible by %D",N,n);
439:     VecGetLocalSize(win,&nlocal);
440:     if (nlocal%n) SETERRQ2(PetscObjectComm((PetscObject)win),1,"Local dimension %D is not divisible by %D",nlocal,n);
441:     s->nx = n;
442:     for (i=0;i<n;i++) {
443:       VecCreate(PetscObjectComm((PetscObject)win),&s->x[i]);
444:       VecSetSizes(s->x[i],nlocal/n,N/n);
445:       VecSetFromOptions(s->x[i]);
446:     }
447:     if (!s->n) {
448:       PetscNewLog(win,&nn);
449:       s->n = nn;
450:       nn->N = N;
451:       nn->lN = nlocal;
452:       nn->friends = 1;
453:     }
454:   } else if (n > s->nx) SETERRQ1(PetscObjectComm((PetscObject)win),PETSC_ERR_SUP,"Number of child vectors cannot be larger than %D",s->nx);
455:   if (x) {
456:     PetscMemcpy(s->x,x,sizeof(Vec)*n);
457:   }
458:   s->n->n = n;
459:   return(0);
460: }

462: /*@C
463:    VecCompSetSubVecs - Resets the number of subvectors defining a compound vector,
464:    or replaces the subvectors.

466:    Collective on Vec

468:    Input Parameters:
469: +  win - compound vector
470: .  n - number of child vectors
471: -  x - array of child vectors

473:    Note:
474:    It is not possible to increase the number of subvectors with respect to the
475:    number set at its creation.

477:    Level: developer

479: .seealso: VecCreateComp(), VecCompGetSubVecs()
480: @*/
481: PetscErrorCode VecCompSetSubVecs(Vec win,PetscInt n,Vec *x)
482: {

488:   PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
489:   return(0);
490: }

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

499:   SlepcValidVecComp(v,1);
500:   SlepcValidVecComp(w,3);
501:   for (i=0;i<vs->n->n;i++) {
502:     VecAXPY(vs->x[i],alpha,ws->x[i]);
503:   }
504:   return(0);
505: }

507: PetscErrorCode VecAYPX_Comp(Vec v,PetscScalar alpha,Vec w)
508: {
510:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
511:   PetscInt       i;

514:   SlepcValidVecComp(v,1);
515:   SlepcValidVecComp(w,3);
516:   for (i=0;i<vs->n->n;i++) {
517:     VecAYPX(vs->x[i],alpha,ws->x[i]);
518:   }
519:   return(0);
520: }

522: PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,Vec w)
523: {
525:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
526:   PetscInt       i;

529:   SlepcValidVecComp(v,1);
530:   SlepcValidVecComp(w,4);
531:   for (i=0;i<vs->n->n;i++) {
532:     VecAXPBY(vs->x[i],alpha,beta,ws->x[i]);
533:   }
534:   return(0);
535: }

537: PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
538: {
540:   Vec_Comp       *vs = (Vec_Comp*)v->data;
541:   Vec            *wx;
542:   PetscInt       i,j;

545:   SlepcValidVecComp(v,1);
546:   for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);

548:   PetscMalloc1(n,&wx);

550:   for (j=0;j<vs->n->n;j++) {
551:     for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
552:     VecMAXPY(vs->x[j],n,alpha,wx);
553:   }

555:   PetscFree(wx);
556:   return(0);
557: }

559: PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,Vec w,Vec z)
560: {
562:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
563:   PetscInt       i;

566:   SlepcValidVecComp(v,1);
567:   SlepcValidVecComp(w,3);
568:   SlepcValidVecComp(z,4);
569:   for (i=0;i<vs->n->n;i++) {
570:     VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]);
571:   }
572:   return(0);
573: }

575: PetscErrorCode VecAXPBYPCZ_Comp(Vec v,PetscScalar alpha,PetscScalar beta,PetscScalar gamma,Vec w,Vec z)
576: {
577:   PetscErrorCode  ierr;
578:   Vec_Comp        *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data,*zs = (Vec_Comp*)z->data;
579:   PetscInt        i;

582:   SlepcValidVecComp(v,1);
583:   SlepcValidVecComp(w,5);
584:   SlepcValidVecComp(z,6);
585:   for (i=0;i<vs->n->n;i++) {
586:     VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]);
587:   }
588:   return(0);
589: }

591: PetscErrorCode VecGetSize_Comp(Vec v,PetscInt *size)
592: {
593:   Vec_Comp *vs = (Vec_Comp*)v->data;

597:   if (vs->n) {
598:     SlepcValidVecComp(v,1);
599:     *size = vs->n->N;
600:   } else *size = v->map->N;
601:   return(0);
602: }

604: PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
605: {
606:   Vec_Comp *vs = (Vec_Comp*)v->data;

610:   if (vs->n) {
611:     SlepcValidVecComp(v,1);
612:     *size = vs->n->lN;
613:   } else *size = v->map->n;
614:   return(0);
615: }

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

626:   SlepcValidVecComp(v,1);
627:   if (!idx && !z) return(0);

629:   if (vs->n->n > 0) {
630:     VecMax(vs->x[0],idx?&idxp:NULL,&zp);
631:   } else {
632:     zp = PETSC_MIN_REAL;
633:     if (idx) idxp = -1;
634:   }
635:   for (i=1;i<vs->n->n;i++) {
636:     VecGetSize(vs->x[i-1],&s0);
637:     s += s0;
638:     VecMax(vs->x[i],idx?&idxp:NULL,&z0);
639:     if (zp < z0) {
640:       if (idx) *idx = s+idxp;
641:       zp = z0;
642:     }
643:   }
644:   if (z) *z = zp;
645:   return(0);
646: }

648: PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
649: {
651:   Vec_Comp       *vs = (Vec_Comp*)v->data;
652:   PetscInt       idxp,s=0,s0;
653:   PetscReal      zp,z0;
654:   PetscInt       i;

657:   SlepcValidVecComp(v,1);
658:   if (!idx && !z) return(0);

660:   if (vs->n->n > 0) {
661:     VecMin(vs->x[0],idx?&idxp:NULL,&zp);
662:   } else {
663:     zp = PETSC_MAX_REAL;
664:     if (idx) idxp = -1;
665:   }
666:   for (i=1;i<vs->n->n;i++) {
667:     VecGetSize(vs->x[i-1],&s0);
668:     s += s0;
669:     VecMin(vs->x[i],idx?&idxp:NULL,&z0);
670:     if (zp > z0) {
671:       if (idx) *idx = s+idxp;
672:       zp = z0;
673:     }
674:   }
675:   if (z) *z = zp;
676:   return(0);
677: }

679: PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
680: {
682:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
683:   PetscReal      work;
684:   PetscInt       i;

687:   SlepcValidVecComp(v,1);
688:   SlepcValidVecComp(w,2);
689:   if (!m || vs->n->n == 0) return(0);
690:   VecMaxPointwiseDivide(vs->x[0],ws->x[0],m);
691:   for (i=1;i<vs->n->n;i++) {
692:     VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work);
693:     *m = PetscMax(*m,work);
694:   }
695:   return(0);
696: }



704: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v) \
705: { \
706:   PetscErrorCode  ierr; \
707:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
708:   PetscInt        i; \
709: \
711:   SlepcValidVecComp(v,1); \
712:   for (i=0;i<vs->n->n;i++) { \
713:     __COMPOSE2__(Vec,NAME)(vs->x[i]); \
714:   } \
715:   return(0);\
716: }

718: __FUNC_TEMPLATE1__(Conjugate)
719: __FUNC_TEMPLATE1__(Reciprocal)
720: __FUNC_TEMPLATE1__(SqrtAbs)
721: __FUNC_TEMPLATE1__(Abs)
722: __FUNC_TEMPLATE1__(Exp)
723: __FUNC_TEMPLATE1__(Log)

726: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,T0 __a) \
727: { \
728:   PetscErrorCode  ierr; \
729:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
730:   PetscInt        i; \
731: \
733:   SlepcValidVecComp(v,1); \
734:   for (i=0;i<vs->n->n;i++) { \
735:     __COMPOSE2__(Vec,NAME)(vs->x[i],__a); \
736:   } \
737:   return(0);\
738: }

740: __FUNC_TEMPLATE2__(Set,PetscScalar)
741: __FUNC_TEMPLATE2__(View,PetscViewer)
742: __FUNC_TEMPLATE2__(Scale,PetscScalar)
743: __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
744: __FUNC_TEMPLATE2__(Shift,PetscScalar)

747: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w) \
748: { \
749:   PetscErrorCode  ierr; \
750:   Vec_Comp        *vs = (Vec_Comp*)v->data,\
751:                   *ws = (Vec_Comp*)w->data; \
752:   PetscInt        i; \
753: \
755:   SlepcValidVecComp(v,1); \
756:   SlepcValidVecComp(w,2); \
757:   for (i=0;i<vs->n->n;i++) { \
758:     __COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i]); \
759:   } \
760:   return(0);\
761: }

763: __FUNC_TEMPLATE3__(Copy)
764: __FUNC_TEMPLATE3__(Swap)

767: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w,Vec z) \
768: { \
769:   PetscErrorCode  ierr; \
770:   Vec_Comp        *vs = (Vec_Comp*)v->data, \
771:                   *ws = (Vec_Comp*)w->data, \
772:                   *zs = (Vec_Comp*)z->data; \
773:   PetscInt        i; \
774: \
776:   SlepcValidVecComp(v,1); \
777:   SlepcValidVecComp(w,2); \
778:   SlepcValidVecComp(z,3); \
779:   for (i=0;i<vs->n->n;i++) { \
780:     __COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i],zs->x[i]); \
781:   } \
782:   return(0);\
783: }

785: __FUNC_TEMPLATE4__(PointwiseMax)
786: __FUNC_TEMPLATE4__(PointwiseMaxAbs)
787: __FUNC_TEMPLATE4__(PointwiseMin)
788: __FUNC_TEMPLATE4__(PointwiseMult)
789: __FUNC_TEMPLATE4__(PointwiseDivide)