Actual source code: veccomp.c

      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(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;
    206:   PetscAssertPointer(V,3);
    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;
    218:   PetscAssertPointer(v,2);
    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:   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: /*MC
    290:    VECCOMP - VECCOMP = "comp" - Vector type consisting of several subvectors,
    291:    each stored separately.

    293:    Level: developer

    295:    Notes:
    296:    This is similar to PETSc's `VECNEST` but customized for SLEPc's needs. In particular,
    297:    the number of child vectors can be modified dynamically, with `VecCompSetSubVecs()`.

    299: .seealso: `Vec`, `VecType`, `VecCreateComp()`, `VecCreateCompWithVecs()`
    300: M*/

    302: /*@
    303:    VecCreateComp - Creates a new vector containing several subvectors,
    304:    each stored separately.

    306:    Collective

    308:    Input Parameters:
    309: +  comm - communicator for the new `Vec`
    310: .  Nx   - array of (initial) global sizes of child vectors
    311: .  n    - number of child vectors
    312: .  t    - type of the child vectors
    313: -  Vparent - (optional) template vector

    315:    Output Parameter:
    316: .  V - new vector

    318:    Notes:
    319:    This is similar to PETSc's `VECNEST` but customized for SLEPc's needs. In particular,
    320:    the number of child vectors can be modified dynamically, with `VecCompSetSubVecs()`.

    322:    Level: developer

    324: .seealso: `VecCreateCompWithVecs()`, `VecCompSetSubVecs()`
    325: @*/
    326: PetscErrorCode VecCreateComp(MPI_Comm comm,PetscInt Nx[],PetscInt n,VecType t,Vec Vparent,Vec *V)
    327: {
    328:   Vec            *x;
    329:   PetscInt       i;

    331:   PetscFunctionBegin;
    332:   PetscCall(VecCreate(comm,V));
    333:   PetscCall(PetscMalloc1(n,&x));
    334:   for (i=0;i<n;i++) {
    335:     PetscCall(VecCreate(comm,&x[i]));
    336:     PetscCall(VecSetSizes(x[i],PETSC_DECIDE,Nx[i]));
    337:     PetscCall(VecSetType(x[i],t));
    338:   }
    339:   PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_TRUE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
    340:   PetscFunctionReturn(PETSC_SUCCESS);
    341: }

    343: /*@
    344:    VecCreateCompWithVecs - Creates a new vector containing several subvectors,
    345:    each stored separately, from an array of `Vec`s.

    347:    Collective

    349:    Input Parameters:
    350: +  x - array of `Vec`s
    351: .  n - number of child vectors
    352: -  Vparent - (optional) template vector

    354:    Output Parameter:
    355: .  V - new vector

    357:    Level: developer

    359: .seealso: `VecCreateComp()`
    360: @*/
    361: PetscErrorCode VecCreateCompWithVecs(Vec x[],PetscInt n,Vec Vparent,Vec *V)
    362: {
    363:   PetscInt       i;

    365:   PetscFunctionBegin;
    366:   PetscAssertPointer(x,1);
    369:   PetscCall(VecCreate(PetscObjectComm((PetscObject)x[0]),V));
    370:   for (i=0;i<n;i++) PetscCall(PetscObjectReference((PetscObject)x[i]));
    371:   PetscCall(VecCreate_Comp_Private(*V,x,n,PETSC_FALSE,Vparent?((Vec_Comp*)Vparent->data)->n:NULL));
    372:   PetscFunctionReturn(PETSC_SUCCESS);
    373: }

    375: PetscErrorCode VecDuplicate_Comp(Vec win,Vec *V)
    376: {
    377:   Vec            *x;
    378:   PetscInt       i;
    379:   Vec_Comp       *s = (Vec_Comp*)win->data;

    381:   PetscFunctionBegin;
    382:   SlepcValidVecComp(win,1);
    383:   PetscCall(VecCreate(PetscObjectComm((PetscObject)win),V));
    384:   PetscCall(PetscMalloc1(s->nx,&x));
    385:   for (i=0;i<s->nx;i++) {
    386:     if (s->x[i]) PetscCall(VecDuplicate(s->x[i],&x[i]));
    387:     else x[i] = NULL;
    388:   }
    389:   PetscCall(VecCreate_Comp_Private(*V,x,s->nx,PETSC_TRUE,s->n));
    390:   PetscFunctionReturn(PETSC_SUCCESS);
    391: }

    393: static PetscErrorCode VecCompGetSubVecs_Comp(Vec win,PetscInt *n,const Vec **x)
    394: {
    395:   Vec_Comp *s = (Vec_Comp*)win->data;

    397:   PetscFunctionBegin;
    398:   if (x) *x = s->x;
    399:   if (n) *n = s->n->n;
    400:   PetscFunctionReturn(PETSC_SUCCESS);
    401: }

    403: /*@C
    404:    VecCompGetSubVecs - Returns the entire array of vectors defining a
    405:    compound vector.

    407:    Collective

    409:    Input Parameter:
    410: .  win - compound vector

    412:    Output Parameters:
    413: +  n - number of child vectors
    414: -  x - array of child vectors

    416:    Level: developer

    418: .seealso: `VecCreateComp()`
    419: @*/
    420: PetscErrorCode VecCompGetSubVecs(Vec win,PetscInt *n,const Vec *x[])
    421: {
    422:   PetscFunctionBegin;
    424:   PetscUseMethod(win,"VecCompGetSubVecs_C",(Vec,PetscInt*,const Vec**),(win,n,x));
    425:   PetscFunctionReturn(PETSC_SUCCESS);
    426: }

    428: static PetscErrorCode VecCompSetSubVecs_Comp(Vec win,PetscInt n,Vec *x)
    429: {
    430:   Vec_Comp       *s = (Vec_Comp*)win->data;
    431:   PetscInt       i,N,nlocal;
    432:   Vec_Comp_N     *nn;

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

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

    466:    Collective

    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: {
    483:   PetscFunctionBegin;
    486:   PetscTryMethod(win,"VecCompSetSubVecs_C",(Vec,PetscInt,Vec*),(win,n,x));
    487:   PetscFunctionReturn(PETSC_SUCCESS);
    488: }

    490: PetscErrorCode VecAXPY_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(VecAXPY(vs->x[i],alpha,ws->x[i]));
    499:   PetscFunctionReturn(PETSC_SUCCESS);
    500: }

    502: PetscErrorCode VecAYPX_Comp(Vec v,PetscScalar alpha,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,3);
    510:   for (i=0;i<vs->n->n;i++) PetscCall(VecAYPX(vs->x[i],alpha,ws->x[i]));
    511:   PetscFunctionReturn(PETSC_SUCCESS);
    512: }

    514: PetscErrorCode VecAXPBY_Comp(Vec v,PetscScalar alpha,PetscScalar beta,Vec w)
    515: {
    516:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
    517:   PetscInt       i;

    519:   PetscFunctionBegin;
    520:   SlepcValidVecComp(v,1);
    521:   SlepcValidVecComp(w,4);
    522:   for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBY(vs->x[i],alpha,beta,ws->x[i]));
    523:   PetscFunctionReturn(PETSC_SUCCESS);
    524: }

    526: PetscErrorCode VecMAXPY_Comp(Vec v,PetscInt n,const PetscScalar *alpha,Vec *w)
    527: {
    528:   Vec_Comp       *vs = (Vec_Comp*)v->data;
    529:   Vec            *wx;
    530:   PetscInt       i,j;

    532:   PetscFunctionBegin;
    533:   SlepcValidVecComp(v,1);
    534:   for (i=0;i<n;i++) SlepcValidVecComp(w[i],4);

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

    538:   for (j=0;j<vs->n->n;j++) {
    539:     for (i=0;i<n;i++) wx[i] = ((Vec_Comp*)w[i]->data)->x[j];
    540:     PetscCall(VecMAXPY(vs->x[j],n,alpha,wx));
    541:   }

    543:   PetscCall(PetscFree(wx));
    544:   PetscFunctionReturn(PETSC_SUCCESS);
    545: }

    547: PetscErrorCode VecWAXPY_Comp(Vec v,PetscScalar alpha,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,3);
    555:   SlepcValidVecComp(z,4);
    556:   for (i=0;i<vs->n->n;i++) PetscCall(VecWAXPY(vs->x[i],alpha,ws->x[i],zs->x[i]));
    557:   PetscFunctionReturn(PETSC_SUCCESS);
    558: }

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

    565:   PetscFunctionBegin;
    566:   SlepcValidVecComp(v,1);
    567:   SlepcValidVecComp(w,5);
    568:   SlepcValidVecComp(z,6);
    569:   for (i=0;i<vs->n->n;i++) PetscCall(VecAXPBYPCZ(vs->x[i],alpha,beta,gamma,ws->x[i],zs->x[i]));
    570:   PetscFunctionReturn(PETSC_SUCCESS);
    571: }

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

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

    586: PetscErrorCode VecGetLocalSize_Comp(Vec v,PetscInt *size)
    587: {
    588:   Vec_Comp *vs = (Vec_Comp*)v->data;

    590:   PetscFunctionBegin;
    591:   PetscAssertPointer(size,2);
    592:   if (vs->n) {
    593:     SlepcValidVecComp(v,1);
    594:     *size = vs->n->lN;
    595:   } else *size = v->map->n;
    596:   PetscFunctionReturn(PETSC_SUCCESS);
    597: }

    599: PetscErrorCode VecMax_Comp(Vec v,PetscInt *idx,PetscReal *z)
    600: {
    601:   Vec_Comp       *vs = (Vec_Comp*)v->data;
    602:   PetscInt       idxp,s=0,s0;
    603:   PetscReal      zp,z0;
    604:   PetscInt       i;

    606:   PetscFunctionBegin;
    607:   SlepcValidVecComp(v,1);
    608:   if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);

    610:   if (vs->n->n > 0) PetscCall(VecMax(vs->x[0],idx?&idxp:NULL,&zp));
    611:   else {
    612:     zp = PETSC_MIN_REAL;
    613:     if (idx) idxp = -1;
    614:   }
    615:   for (i=1;i<vs->n->n;i++) {
    616:     PetscCall(VecGetSize(vs->x[i-1],&s0));
    617:     s += s0;
    618:     PetscCall(VecMax(vs->x[i],idx?&idxp:NULL,&z0));
    619:     if (zp < z0) {
    620:       if (idx) *idx = s+idxp;
    621:       zp = z0;
    622:     }
    623:   }
    624:   if (z) *z = zp;
    625:   PetscFunctionReturn(PETSC_SUCCESS);
    626: }

    628: PetscErrorCode VecMin_Comp(Vec v,PetscInt *idx,PetscReal *z)
    629: {
    630:   Vec_Comp       *vs = (Vec_Comp*)v->data;
    631:   PetscInt       idxp,s=0,s0;
    632:   PetscReal      zp,z0;
    633:   PetscInt       i;

    635:   PetscFunctionBegin;
    636:   SlepcValidVecComp(v,1);
    637:   if (!idx && !z) PetscFunctionReturn(PETSC_SUCCESS);

    639:   if (vs->n->n > 0) PetscCall(VecMin(vs->x[0],idx?&idxp:NULL,&zp));
    640:   else {
    641:     zp = PETSC_MAX_REAL;
    642:     if (idx) idxp = -1;
    643:   }
    644:   for (i=1;i<vs->n->n;i++) {
    645:     PetscCall(VecGetSize(vs->x[i-1],&s0));
    646:     s += s0;
    647:     PetscCall(VecMin(vs->x[i],idx?&idxp:NULL,&z0));
    648:     if (zp > z0) {
    649:       if (idx) *idx = s+idxp;
    650:       zp = z0;
    651:     }
    652:   }
    653:   if (z) *z = zp;
    654:   PetscFunctionReturn(PETSC_SUCCESS);
    655: }

    657: PetscErrorCode VecMaxPointwiseDivide_Comp(Vec v,Vec w,PetscReal *m)
    658: {
    659:   Vec_Comp       *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
    660:   PetscReal      work;
    661:   PetscInt       i;

    663:   PetscFunctionBegin;
    664:   SlepcValidVecComp(v,1);
    665:   SlepcValidVecComp(w,2);
    666:   if (!m || vs->n->n == 0) PetscFunctionReturn(PETSC_SUCCESS);
    667:   PetscCall(VecMaxPointwiseDivide(vs->x[0],ws->x[0],m));
    668:   for (i=1;i<vs->n->n;i++) {
    669:     PetscCall(VecMaxPointwiseDivide(vs->x[i],ws->x[i],&work));
    670:     *m = PetscMax(*m,work);
    671:   }
    672:   PetscFunctionReturn(PETSC_SUCCESS);
    673: }

    680: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v) \
    681: { \
    682:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
    683:   PetscInt        i; \
    684: \
    685:   PetscFunctionBegin; \
    686:   SlepcValidVecComp(v,1); \
    687:   for (i=0;i<vs->n->n;i++) { \
    688:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i])); \
    689:   } \
    690:   PetscFunctionReturn(PETSC_SUCCESS);\
    691: }

    693: __FUNC_TEMPLATE1__(Conjugate)
    694: __FUNC_TEMPLATE1__(Reciprocal)
    695: __FUNC_TEMPLATE1__(SqrtAbs)
    696: __FUNC_TEMPLATE1__(Abs)
    697: __FUNC_TEMPLATE1__(Exp)
    698: __FUNC_TEMPLATE1__(Log)

    701: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,T0 __a) \
    702: { \
    703:   Vec_Comp        *vs = (Vec_Comp*)v->data; \
    704:   PetscInt        i; \
    705: \
    706:   PetscFunctionBegin; \
    707:   SlepcValidVecComp(v,1); \
    708:   for (i=0;i<vs->n->n;i++) { \
    709:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],__a)); \
    710:   } \
    711:   PetscFunctionReturn(PETSC_SUCCESS);\
    712: }

    714: __FUNC_TEMPLATE2__(Set,PetscScalar)
    715: __FUNC_TEMPLATE2__(View,PetscViewer)
    716: __FUNC_TEMPLATE2__(Scale,PetscScalar)
    717: __FUNC_TEMPLATE2__(SetRandom,PetscRandom)
    718: __FUNC_TEMPLATE2__(Shift,PetscScalar)

    721: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w) \
    722: { \
    723:   Vec_Comp        *vs = (Vec_Comp*)v->data,\
    724:                   *ws = (Vec_Comp*)w->data; \
    725:   PetscInt        i; \
    726: \
    727:   PetscFunctionBegin; \
    728:   SlepcValidVecComp(v,1); \
    729:   SlepcValidVecComp(w,2); \
    730:   for (i=0;i<vs->n->n;i++) { \
    731:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i])); \
    732:   } \
    733:   PetscFunctionReturn(PETSC_SUCCESS);\
    734: }

    736: __FUNC_TEMPLATE3__(Copy)
    737: __FUNC_TEMPLATE3__(Swap)

    740: PetscErrorCode __COMPOSE3__(Vec,NAME,_Comp)(Vec v,Vec w,Vec z) \
    741: { \
    742:   Vec_Comp        *vs = (Vec_Comp*)v->data, \
    743:                   *ws = (Vec_Comp*)w->data, \
    744:                   *zs = (Vec_Comp*)z->data; \
    745:   PetscInt        i; \
    746: \
    747:   PetscFunctionBegin; \
    748:   SlepcValidVecComp(v,1); \
    749:   SlepcValidVecComp(w,2); \
    750:   SlepcValidVecComp(z,3); \
    751:   for (i=0;i<vs->n->n;i++) { \
    752:     PetscCall(__COMPOSE2__(Vec,NAME)(vs->x[i],ws->x[i],zs->x[i])); \
    753:   } \
    754:   PetscFunctionReturn(PETSC_SUCCESS);\
    755: }

    757: __FUNC_TEMPLATE4__(PointwiseMax)
    758: __FUNC_TEMPLATE4__(PointwiseMaxAbs)
    759: __FUNC_TEMPLATE4__(PointwiseMin)
    760: __FUNC_TEMPLATE4__(PointwiseMult)
    761: __FUNC_TEMPLATE4__(PointwiseDivide)