Actual source code: veccomp0.h

      1: /*
      2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
      4:    Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain

      6:    This file is part of SLEPc.
      7:    SLEPc is distributed under a 2-clause BSD license (see LICENSE).
      8:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      9: */

     11: #include <petsc/private/vecimpl.h>

     13: #if defined(__WITH_MPI__)
     15: #else
     17: #endif

     22: static PetscErrorCode __SUF__(VecDot_Comp)(Vec a,Vec b,PetscScalar *z)
     23: {
     24:   PetscScalar    sum = 0.0,work;
     25:   PetscInt       i;
     26:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

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

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

     57:   PetscFunctionBegin;
     58:   SlepcValidVecComp(a,1);
     59:   SlepcValidVecsComp(b,n,3);

     61:   if (as->n->n == 0) {
     62:     *z = 0;
     63:     PetscFunctionReturn(PETSC_SUCCESS);
     64:   }

     66:   PetscCall(PetscMalloc2(n,&work0,n,&bx));

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

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

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

     96:   PetscCall(PetscFree2(work0,bx));
     97:   PetscFunctionReturn(PETSC_SUCCESS);
     98: }

    100: static PetscErrorCode __SUF__(VecTDot_Comp)(Vec a,Vec b,PetscScalar *z)
    101: {
    102:   PetscScalar    sum = 0.0,work;
    103:   PetscInt       i;
    104:   Vec_Comp       *as = (Vec_Comp*)a->data,*bs = (Vec_Comp*)b->data;

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

    128: static PetscErrorCode __SUF__(VecMTDot_Comp)(Vec a,PetscInt n,const Vec b[],PetscScalar *z)
    129: {
    130:   PetscScalar    *work,*work0,*r;
    131:   Vec_Comp       *as = (Vec_Comp*)a->data;
    132:   Vec            *bx;
    133:   PetscInt       i,j;

    135:   PetscFunctionBegin;
    136:   SlepcValidVecComp(a,1);
    137:   SlepcValidVecsComp(b,n,3);

    139:   if (as->n->n == 0) {
    140:     *z = 0;
    141:     PetscFunctionReturn(PETSC_SUCCESS);
    142:   }

    144:   PetscCall(PetscMalloc2(n,&work0,n,&bx));

    146: #if defined(__WITH_MPI__)
    147:   if (as->x[0]->ops->mtdot_local) {
    148:     r = work0;
    149:     work = z;
    150:   } else
    151: #endif
    152:   {
    153:     r = z;
    154:     work = work0;
    155:   }

    157:   /* z[i] <- a.x' * b[i].x */
    158:   for (i=0;i<n;i++) r[i] = 0.0;
    159:   for (j=0;j<as->n->n;j++) {
    160:     for (i=0;i<n;i++) bx[i] = ((Vec_Comp*)b[i]->data)->x[j];
    161:     if (as->x[0]->ops->mtdot_local) PetscUseTypeMethod(as->x[j],mtdot_local,n,bx,work);
    162:     else PetscCall(VecMTDot(as->x[j],n,bx,work));
    163:     for (i=0;i<n;i++) r[i] += work[i];
    164:   }

    166:   /* If def(__WITH_MPI__) and exists mtdot_local */
    167: #if defined(__WITH_MPI__)
    168:   if (as->x[0]->ops->mtdot_local) {
    169:     /* z[i] <- Allreduce(work[i]) */
    170:     PetscCallMPI(MPIU_Allreduce(r,z,n,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)a)));
    171:   }
    172: #endif

    174:   PetscCall(PetscFree2(work0,bx));
    175:   PetscFunctionReturn(PETSC_SUCCESS);
    176: }

    178: static PetscErrorCode __SUF__(VecNorm_Comp)(Vec a,NormType t,PetscReal *norm)
    179: {
    180:   PetscReal      work[3],s=0.0;
    181:   Vec_Comp       *as = (Vec_Comp*)a->data;
    182:   PetscInt       i;

    184:   PetscFunctionBegin;
    185:   SlepcValidVecComp(a,1);
    186:   /* Initialize norm */
    187:   switch (t) {
    188:     case NORM_1:
    189:     case NORM_INFINITY:
    190:       *norm = 0.0;
    191:       break;
    192:     case NORM_2:
    193:     case NORM_FROBENIUS:
    194:       *norm = 1.0;
    195:       s = 0.0;
    196:       break;
    197:     case NORM_1_AND_2:
    198:       norm[0] = 0.0;
    199:       norm[1] = 1.0;
    200:       s = 0.0;
    201:       break;
    202:   }
    203:   for (i=0;i<as->n->n;i++) {
    204:     if (as->x[0]->ops->norm_local) PetscUseTypeMethod(as->x[i],norm_local,t,work);
    205:     else PetscCall(VecNorm(as->x[i],t,work));
    206:     /* norm+= work */
    207:     switch (t) {
    208:       case NORM_1:
    209:         *norm += *work;
    210:         break;
    211:       case NORM_2:
    212:       case NORM_FROBENIUS:
    213:         AddNorm2(norm,&s,*work);
    214:         break;
    215:       case NORM_1_AND_2:
    216:         norm[0] += work[0];
    217:         AddNorm2(&norm[1],&s,work[1]);
    218:         break;
    219:       case NORM_INFINITY:
    220:         *norm = PetscMax(*norm,*work);
    221:         break;
    222:     }
    223:   }

    225:   /* If def(__WITH_MPI__) and exists norm_local */
    226: #if defined(__WITH_MPI__)
    227:   if (as->x[0]->ops->norm_local) {
    228:     PetscReal work0[3];
    229:     /* norm <- Allreduce(work) */
    230:     switch (t) {
    231:     case NORM_1:
    232:       work[0] = *norm;
    233:       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_SUM,PetscObjectComm((PetscObject)a)));
    234:       break;
    235:     case NORM_2:
    236:     case NORM_FROBENIUS:
    237:       work[0] = *norm;
    238:       work[1] = s;
    239:       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
    240:       *norm = GetNorm2(work0[0],work0[1]);
    241:       break;
    242:     case NORM_1_AND_2:
    243:       work[0] = norm[0];
    244:       work[1] = norm[1];
    245:       work[2] = s;
    246:       PetscCallMPI(MPIU_Allreduce(work,work0,1,MPIU_NORM1_AND_2,MPIU_NORM2_SUM,PetscObjectComm((PetscObject)a)));
    247:       norm[0] = work0[0];
    248:       norm[1] = GetNorm2(work0[1],work0[2]);
    249:       break;
    250:     case NORM_INFINITY:
    251:       work[0] = *norm;
    252:       PetscCallMPI(MPIU_Allreduce(work,norm,1,MPIU_REAL,MPIU_MAX,PetscObjectComm((PetscObject)a)));
    253:       break;
    254:     }
    255:   }
    256: #else
    257:   /* Norm correction */
    258:   switch (t) {
    259:     case NORM_1:
    260:     case NORM_INFINITY:
    261:       break;
    262:     case NORM_2:
    263:     case NORM_FROBENIUS:
    264:       *norm = GetNorm2(*norm,s);
    265:       break;
    266:     case NORM_1_AND_2:
    267:       norm[1] = GetNorm2(norm[1],s);
    268:       break;
    269:   }
    270: #endif
    271:   PetscFunctionReturn(PETSC_SUCCESS);
    272: }

    274: PetscErrorCode __SUF__(VecDotNorm2_Comp)(Vec v,Vec w,PetscScalar *dp,PetscScalar *nm)
    275: {
    276:   PetscScalar       dp0=0.0,nm0=0.0,dp1=0.0,nm1=0.0;
    277:   const PetscScalar *vx,*wx;
    278:   Vec_Comp          *vs = (Vec_Comp*)v->data,*ws = (Vec_Comp*)w->data;
    279:   PetscInt          i,n;
    280:   PetscBool         t0,t1;
    281: #if defined(__WITH_MPI__)
    282:   PetscScalar       work[4];
    283: #endif

    285:   PetscFunctionBegin;
    286:   /* Compute recursively the local part */
    287:   PetscCall(PetscObjectTypeCompare((PetscObject)v,VECCOMP,&t0));
    288:   PetscCall(PetscObjectTypeCompare((PetscObject)w,VECCOMP,&t1));
    289:   if (t0 && t1) {
    290:     SlepcValidVecComp(v,1);
    291:     SlepcValidVecComp(w,2);
    292:     for (i=0;i<vs->n->n;i++) {
    293:       PetscCall(VecDotNorm2_Comp_Seq(vs->x[i],ws->x[i],&dp1,&nm1));
    294:       dp0 += dp1;
    295:       nm0 += nm1;
    296:     }
    297:   } else if (!t0 && !t1) {
    298:     PetscCall(VecGetLocalSize(v,&n));
    299:     PetscCall(VecGetArrayRead(v,&vx));
    300:     PetscCall(VecGetArrayRead(w,&wx));
    301:     for (i=0;i<n;i++) {
    302:       dp0 += vx[i]*PetscConj(wx[i]);
    303:       nm0 += wx[i]*PetscConj(wx[i]);
    304:     }
    305:     PetscCall(VecRestoreArrayRead(v,&vx));
    306:     PetscCall(VecRestoreArrayRead(w,&wx));
    307:   } else SETERRQ(PetscObjectComm((PetscObject)v),PETSC_ERR_ARG_INCOMP,"Incompatible vector types");

    309: #if defined(__WITH_MPI__)
    310:     /* [dp, nm] <- Allreduce([dp0, nm0]) */
    311:     work[0] = dp0;
    312:     work[1] = nm0;
    313:     PetscCallMPI(MPIU_Allreduce(work,&work[2],2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v)));
    314:     *dp = work[2];
    315:     *nm = work[3];
    316: #else
    317:     *dp = dp0;
    318:     *nm = nm0;
    319: #endif
    320:   PetscFunctionReturn(PETSC_SUCCESS);
    321: }