Actual source code: vinv.c

      2: /*
      3:      Some useful vector utility functions.
      4: */
      5: #include <../src/vec/vec/impls/mpi/pvecimpl.h>

      7: /*@
      8:    VecStrideSet - Sets a subvector of a vector defined
      9:    by a starting point and a stride with a given value

     11:    Logically Collective on Vec

     13:    Input Parameter:
     14: +  v - the vector
     15: .  start - starting point of the subvector (defined by a stride)
     16: -  s - value to set for each entry in that subvector

     18:    Notes:
     19:    One must call VecSetBlockSize() before this routine to set the stride
     20:    information, or use a vector created from a multicomponent DMDA.

     22:    This will only work if the desire subvector is a stride subvector

     24:    Level: advanced

     27: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
     28: @*/
     29: PetscErrorCode  VecStrideSet(Vec v,PetscInt start,PetscScalar s)
     30: {
     32:   PetscInt       i,n,bs;
     33:   PetscScalar    *x;

     39:   VecSetErrorIfLocked(v,1);

     41:   VecGetLocalSize(v,&n);
     42:   VecGetArray(v,&x);
     43:   VecGetBlockSize(v,&bs);
     44:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
     45:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n  Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
     46:   x += start;
     47:   for (i=0; i<n; i+=bs) x[i] = s;
     48:   x -= start;
     49:   VecRestoreArray(v,&x);
     50:   return(0);
     51: }

     53: /*@
     54:    VecStrideScale - Scales a subvector of a vector defined
     55:    by a starting point and a stride.

     57:    Logically Collective on Vec

     59:    Input Parameter:
     60: +  v - the vector
     61: .  start - starting point of the subvector (defined by a stride)
     62: -  scale - value to multiply each subvector entry by

     64:    Notes:
     65:    One must call VecSetBlockSize() before this routine to set the stride
     66:    information, or use a vector created from a multicomponent DMDA.

     68:    This will only work if the desire subvector is a stride subvector

     70:    Level: advanced

     73: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
     74: @*/
     75: PetscErrorCode  VecStrideScale(Vec v,PetscInt start,PetscScalar scale)
     76: {
     78:   PetscInt       i,n,bs;
     79:   PetscScalar    *x;

     85:   VecSetErrorIfLocked(v,1);

     87:   VecGetLocalSize(v,&n);
     88:   VecGetArray(v,&x);
     89:   VecGetBlockSize(v,&bs);
     90:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
     91:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n  Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
     92:   x += start;
     93:   for (i=0; i<n; i+=bs) x[i] *= scale;
     94:   x -= start;
     95:   VecRestoreArray(v,&x);
     96:   return(0);
     97: }

     99: /*@
    100:    VecStrideNorm - Computes the norm of subvector of a vector defined
    101:    by a starting point and a stride.

    103:    Collective on Vec

    105:    Input Parameter:
    106: +  v - the vector
    107: .  start - starting point of the subvector (defined by a stride)
    108: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

    110:    Output Parameter:
    111: .  norm - the norm

    113:    Notes:
    114:    One must call VecSetBlockSize() before this routine to set the stride
    115:    information, or use a vector created from a multicomponent DMDA.

    117:    If x is the array representing the vector x then this computes the norm
    118:    of the array (x[start],x[start+stride],x[start+2*stride], ....)

    120:    This is useful for computing, say the norm of the pressure variable when
    121:    the pressure is stored (interlaced) with other variables, say density etc.

    123:    This will only work if the desire subvector is a stride subvector

    125:    Level: advanced

    128: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
    129: @*/
    130: PetscErrorCode  VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
    131: {
    132:   PetscErrorCode    ierr;
    133:   PetscInt          i,n,bs;
    134:   const PetscScalar *x;
    135:   PetscReal         tnorm;
    136:   MPI_Comm          comm;

    141:   VecGetLocalSize(v,&n);
    142:   VecGetArrayRead(v,&x);
    143:   PetscObjectGetComm((PetscObject)v,&comm);

    145:   VecGetBlockSize(v,&bs);
    146:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
    147:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    148:   x += start;

    150:   if (ntype == NORM_2) {
    151:     PetscScalar sum = 0.0;
    152:     for (i=0; i<n; i+=bs) sum += x[i]*(PetscConj(x[i]));
    153:     tnorm = PetscRealPart(sum);
    154:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
    155:     *nrm  = PetscSqrtReal(*nrm);
    156:   } else if (ntype == NORM_1) {
    157:     tnorm = 0.0;
    158:     for (i=0; i<n; i+=bs) tnorm += PetscAbsScalar(x[i]);
    159:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_SUM,comm);
    160:   } else if (ntype == NORM_INFINITY) {
    161:     PetscReal tmp;
    162:     tnorm = 0.0;

    164:     for (i=0; i<n; i+=bs) {
    165:       if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
    166:       /* check special case of tmp == NaN */
    167:       if (tmp != tmp) {tnorm = tmp; break;}
    168:     }
    169:     MPIU_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPIU_MAX,comm);
    170:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
    171:   VecRestoreArrayRead(v,&x);
    172:   return(0);
    173: }

    175: /*@
    176:    VecStrideMax - Computes the maximum of subvector of a vector defined
    177:    by a starting point and a stride and optionally its location.

    179:    Collective on Vec

    181:    Input Parameter:
    182: +  v - the vector
    183: -  start - starting point of the subvector (defined by a stride)

    185:    Output Parameter:
    186: +  index - the location where the maximum occurred  (pass NULL if not required)
    187: -  nrm - the maximum value in the subvector

    189:    Notes:
    190:    One must call VecSetBlockSize() before this routine to set the stride
    191:    information, or use a vector created from a multicomponent DMDA.

    193:    If xa is the array representing the vector x, then this computes the max
    194:    of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

    196:    This is useful for computing, say the maximum of the pressure variable when
    197:    the pressure is stored (interlaced) with other variables, e.g., density, etc.
    198:    This will only work if the desire subvector is a stride subvector.

    200:    Level: advanced

    203: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
    204: @*/
    205: PetscErrorCode  VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
    206: {
    207:   PetscErrorCode    ierr;
    208:   PetscInt          i,n,bs,id;
    209:   const PetscScalar *x;
    210:   PetscReal         max,tmp;
    211:   MPI_Comm          comm;

    217:   VecGetLocalSize(v,&n);
    218:   VecGetArrayRead(v,&x);
    219:   PetscObjectGetComm((PetscObject)v,&comm);

    221:   VecGetBlockSize(v,&bs);
    222:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
    223:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    224:   x += start;

    226:   id = -1;
    227:   if (!n) max = PETSC_MIN_REAL;
    228:   else {
    229:     id  = 0;
    230:     max = PetscRealPart(x[0]);
    231:     for (i=bs; i<n; i+=bs) {
    232:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
    233:     }
    234:   }
    235:   VecRestoreArrayRead(v,&x);

    237:   if (!idex) {
    238:     MPIU_Allreduce(&max,nrm,1,MPIU_REAL,MPIU_MAX,comm);
    239:   } else {
    240:     PetscReal in[2],out[2];
    241:     PetscInt  rstart;

    243:     VecGetOwnershipRange(v,&rstart,NULL);
    244:     in[0] = max;
    245:     in[1] = rstart+id+start;
    246:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MAXINDEX_OP,PetscObjectComm((PetscObject)v));
    247:     *nrm  = out[0];
    248:     *idex = (PetscInt)out[1];
    249:   }
    250:   return(0);
    251: }

    253: /*@
    254:    VecStrideMin - Computes the minimum of subvector of a vector defined
    255:    by a starting point and a stride and optionally its location.

    257:    Collective on Vec

    259:    Input Parameter:
    260: +  v - the vector
    261: -  start - starting point of the subvector (defined by a stride)

    263:    Output Parameter:
    264: +  idex - the location where the minimum occurred. (pass NULL if not required)
    265: -  nrm - the minimum value in the subvector

    267:    Level: advanced

    269:    Notes:
    270:    One must call VecSetBlockSize() before this routine to set the stride
    271:    information, or use a vector created from a multicomponent DMDA.

    273:    If xa is the array representing the vector x, then this computes the min
    274:    of the array (xa[start],xa[start+stride],xa[start+2*stride], ....)

    276:    This is useful for computing, say the minimum of the pressure variable when
    277:    the pressure is stored (interlaced) with other variables, e.g., density, etc.
    278:    This will only work if the desire subvector is a stride subvector.

    281: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
    282: @*/
    283: PetscErrorCode  VecStrideMin(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
    284: {
    285:   PetscErrorCode    ierr;
    286:   PetscInt          i,n,bs,id;
    287:   const PetscScalar *x;
    288:   PetscReal         min,tmp;
    289:   MPI_Comm          comm;

    295:   VecGetLocalSize(v,&n);
    296:   VecGetArrayRead(v,&x);
    297:   PetscObjectGetComm((PetscObject)v,&comm);

    299:   VecGetBlockSize(v,&bs);
    300:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
    301:   else if (start >= bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\nHave you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    302:   x += start;

    304:   id = -1;
    305:   if (!n) min = PETSC_MAX_REAL;
    306:   else {
    307:     id = 0;
    308:     min = PetscRealPart(x[0]);
    309:     for (i=bs; i<n; i+=bs) {
    310:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
    311:     }
    312:   }
    313:   VecRestoreArrayRead(v,&x);

    315:   if (!idex) {
    316:     MPIU_Allreduce(&min,nrm,1,MPIU_REAL,MPIU_MIN,comm);
    317:   } else {
    318:     PetscReal in[2],out[2];
    319:     PetscInt  rstart;

    321:     VecGetOwnershipRange(v,&rstart,NULL);
    322:     in[0] = min;
    323:     in[1] = rstart+id;
    324:     MPIU_Allreduce(in,out,2,MPIU_REAL,MPIU_MININDEX_OP,PetscObjectComm((PetscObject)v));
    325:     *nrm  = out[0];
    326:     *idex = (PetscInt)out[1];
    327:   }
    328:   return(0);
    329: }

    331: /*@
    332:    VecStrideScaleAll - Scales the subvectors of a vector defined
    333:    by a starting point and a stride.

    335:    Logically Collective on Vec

    337:    Input Parameter:
    338: +  v - the vector
    339: -  scales - values to multiply each subvector entry by

    341:    Notes:
    342:    One must call VecSetBlockSize() before this routine to set the stride
    343:    information, or use a vector created from a multicomponent DMDA.

    345:    The dimension of scales must be the same as the vector block size

    348:    Level: advanced

    351: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
    352: @*/
    353: PetscErrorCode  VecStrideScaleAll(Vec v,const PetscScalar *scales)
    354: {
    356:   PetscInt       i,j,n,bs;
    357:   PetscScalar    *x;

    362:   VecSetErrorIfLocked(v,1);

    364:   VecGetLocalSize(v,&n);
    365:   VecGetArray(v,&x);
    366:   VecGetBlockSize(v,&bs);

    368:   /* need to provide optimized code for each bs */
    369:   for (i=0; i<n; i+=bs) {
    370:     for (j=0; j<bs; j++) x[i+j] *= scales[j];
    371:   }
    372:   VecRestoreArray(v,&x);
    373:   return(0);
    374: }

    376: /*@
    377:    VecStrideNormAll - Computes the norms of subvectors of a vector defined
    378:    by a starting point and a stride.

    380:    Collective on Vec

    382:    Input Parameter:
    383: +  v - the vector
    384: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

    386:    Output Parameter:
    387: .  nrm - the norms

    389:    Notes:
    390:    One must call VecSetBlockSize() before this routine to set the stride
    391:    information, or use a vector created from a multicomponent DMDA.

    393:    If x is the array representing the vector x then this computes the norm
    394:    of the array (x[start],x[start+stride],x[start+2*stride], ....) for each start < stride

    396:    The dimension of nrm must be the same as the vector block size

    398:    This will only work if the desire subvector is a stride subvector

    400:    Level: advanced

    403: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
    404: @*/
    405: PetscErrorCode  VecStrideNormAll(Vec v,NormType ntype,PetscReal nrm[])
    406: {
    407:   PetscErrorCode    ierr;
    408:   PetscInt          i,j,n,bs;
    409:   const PetscScalar *x;
    410:   PetscReal         tnorm[128];
    411:   MPI_Comm          comm;

    416:   VecGetLocalSize(v,&n);
    417:   VecGetArrayRead(v,&x);
    418:   PetscObjectGetComm((PetscObject)v,&comm);

    420:   VecGetBlockSize(v,&bs);
    421:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    423:   if (ntype == NORM_2) {
    424:     PetscScalar sum[128];
    425:     for (j=0; j<bs; j++) sum[j] = 0.0;
    426:     for (i=0; i<n; i+=bs) {
    427:       for (j=0; j<bs; j++) sum[j] += x[i+j]*(PetscConj(x[i+j]));
    428:     }
    429:     for (j=0; j<bs; j++) tnorm[j]  = PetscRealPart(sum[j]);

    431:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
    432:     for (j=0; j<bs; j++) nrm[j] = PetscSqrtReal(nrm[j]);
    433:   } else if (ntype == NORM_1) {
    434:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

    436:     for (i=0; i<n; i+=bs) {
    437:       for (j=0; j<bs; j++) tnorm[j] += PetscAbsScalar(x[i+j]);
    438:     }

    440:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_SUM,comm);
    441:   } else if (ntype == NORM_INFINITY) {
    442:     PetscReal tmp;
    443:     for (j=0; j<bs; j++) tnorm[j] = 0.0;

    445:     for (i=0; i<n; i+=bs) {
    446:       for (j=0; j<bs; j++) {
    447:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
    448:         /* check special case of tmp == NaN */
    449:         if (tmp != tmp) {tnorm[j] = tmp; break;}
    450:       }
    451:     }
    452:     MPIU_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPIU_MAX,comm);
    453:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
    454:   VecRestoreArrayRead(v,&x);
    455:   return(0);
    456: }

    458: /*@
    459:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
    460:    by a starting point and a stride and optionally its location.

    462:    Collective on Vec

    464:    Input Parameter:
    465: .  v - the vector

    467:    Output Parameter:
    468: +  index - the location where the maximum occurred (not supported, pass NULL,
    469:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
    470: -  nrm - the maximum values of each subvector

    472:    Notes:
    473:    One must call VecSetBlockSize() before this routine to set the stride
    474:    information, or use a vector created from a multicomponent DMDA.

    476:    The dimension of nrm must be the same as the vector block size

    478:    Level: advanced

    481: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
    482: @*/
    483: PetscErrorCode  VecStrideMaxAll(Vec v,PetscInt idex[],PetscReal nrm[])
    484: {
    485:   PetscErrorCode    ierr;
    486:   PetscInt          i,j,n,bs;
    487:   const PetscScalar *x;
    488:   PetscReal         max[128],tmp;
    489:   MPI_Comm          comm;

    494:   if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
    495:   VecGetLocalSize(v,&n);
    496:   VecGetArrayRead(v,&x);
    497:   PetscObjectGetComm((PetscObject)v,&comm);

    499:   VecGetBlockSize(v,&bs);
    500:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    502:   if (!n) {
    503:     for (j=0; j<bs; j++) max[j] = PETSC_MIN_REAL;
    504:   } else {
    505:     for (j=0; j<bs; j++) max[j] = PetscRealPart(x[j]);

    507:     for (i=bs; i<n; i+=bs) {
    508:       for (j=0; j<bs; j++) {
    509:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) max[j] = tmp;
    510:       }
    511:     }
    512:   }
    513:   MPIU_Allreduce(max,nrm,bs,MPIU_REAL,MPIU_MAX,comm);

    515:   VecRestoreArrayRead(v,&x);
    516:   return(0);
    517: }

    519: /*@
    520:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
    521:    by a starting point and a stride and optionally its location.

    523:    Collective on Vec

    525:    Input Parameter:
    526: .  v - the vector

    528:    Output Parameter:
    529: +  idex - the location where the minimum occurred (not supported, pass NULL,
    530:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
    531: -  nrm - the minimums of each subvector

    533:    Level: advanced

    535:    Notes:
    536:    One must call VecSetBlockSize() before this routine to set the stride
    537:    information, or use a vector created from a multicomponent DMDA.

    539:    The dimension of nrm must be the same as the vector block size

    542: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
    543: @*/
    544: PetscErrorCode  VecStrideMinAll(Vec v,PetscInt idex[],PetscReal nrm[])
    545: {
    546:   PetscErrorCode    ierr;
    547:   PetscInt          i,n,bs,j;
    548:   const PetscScalar *x;
    549:   PetscReal         min[128],tmp;
    550:   MPI_Comm          comm;

    555:   if (idex) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
    556:   VecGetLocalSize(v,&n);
    557:   VecGetArrayRead(v,&x);
    558:   PetscObjectGetComm((PetscObject)v,&comm);

    560:   VecGetBlockSize(v,&bs);
    561:   if (bs > 128) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    563:   if (!n) {
    564:     for (j=0; j<bs; j++) min[j] = PETSC_MAX_REAL;
    565:   } else {
    566:     for (j=0; j<bs; j++) min[j] = PetscRealPart(x[j]);

    568:     for (i=bs; i<n; i+=bs) {
    569:       for (j=0; j<bs; j++) {
    570:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) min[j] = tmp;
    571:       }
    572:     }
    573:   }
    574:   MPIU_Allreduce(min,nrm,bs,MPIU_REAL,MPIU_MIN,comm);

    576:   VecRestoreArrayRead(v,&x);
    577:   return(0);
    578: }

    580: /*----------------------------------------------------------------------------------------------*/
    581: /*@
    582:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
    583:    separate vectors.

    585:    Collective on Vec

    587:    Input Parameter:
    588: +  v - the vector
    589: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    591:    Output Parameter:
    592: .  s - the location where the subvectors are stored

    594:    Notes:
    595:    One must call VecSetBlockSize() before this routine to set the stride
    596:    information, or use a vector created from a multicomponent DMDA.

    598:    If x is the array representing the vector x then this gathers
    599:    the arrays (x[start],x[start+stride],x[start+2*stride], ....)
    600:    for start=0,1,2,...bs-1

    602:    The parallel layout of the vector and the subvector must be the same;
    603:    i.e., nlocal of v = stride*(nlocal of s)

    605:    Not optimized; could be easily

    607:    Level: advanced

    609: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
    610:           VecStrideScatterAll()
    611: @*/
    612: PetscErrorCode  VecStrideGatherAll(Vec v,Vec s[],InsertMode addv)
    613: {
    614:   PetscErrorCode    ierr;
    615:   PetscInt          i,n,n2,bs,j,k,*bss = NULL,nv,jj,nvc;
    616:   PetscScalar       **y;
    617:   const PetscScalar *x;

    623:   VecGetLocalSize(v,&n);
    624:   VecGetLocalSize(s[0],&n2);
    625:   VecGetArrayRead(v,&x);
    626:   VecGetBlockSize(v,&bs);
    627:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

    629:   PetscMalloc2(bs,&y,bs,&bss);
    630:   nv   = 0;
    631:   nvc  = 0;
    632:   for (i=0; i<bs; i++) {
    633:     VecGetBlockSize(s[i],&bss[i]);
    634:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    635:     VecGetArray(s[i],&y[i]);
    636:     nvc += bss[i];
    637:     nv++;
    638:     if (nvc > bs)  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
    639:     if (nvc == bs) break;
    640:   }

    642:   n =  n/bs;

    644:   jj = 0;
    645:   if (addv == INSERT_VALUES) {
    646:     for (j=0; j<nv; j++) {
    647:       for (k=0; k<bss[j]; k++) {
    648:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = x[bs*i+jj+k];
    649:       }
    650:       jj += bss[j];
    651:     }
    652:   } else if (addv == ADD_VALUES) {
    653:     for (j=0; j<nv; j++) {
    654:       for (k=0; k<bss[j]; k++) {
    655:         for (i=0; i<n; i++) y[j][i*bss[j] + k] += x[bs*i+jj+k];
    656:       }
    657:       jj += bss[j];
    658:     }
    659: #if !defined(PETSC_USE_COMPLEX)
    660:   } else if (addv == MAX_VALUES) {
    661:     for (j=0; j<nv; j++) {
    662:       for (k=0; k<bss[j]; k++) {
    663:         for (i=0; i<n; i++) y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
    664:       }
    665:       jj += bss[j];
    666:     }
    667: #endif
    668:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    670:   VecRestoreArrayRead(v,&x);
    671:   for (i=0; i<nv; i++) {
    672:     VecRestoreArray(s[i],&y[i]);
    673:   }

    675:   PetscFree2(y,bss);
    676:   return(0);
    677: }

    679: /*@
    680:    VecStrideScatterAll - Scatters all the single components from separate vectors into
    681:      a multi-component vector.

    683:    Collective on Vec

    685:    Input Parameter:
    686: +  s - the location where the subvectors are stored
    687: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    689:    Output Parameter:
    690: .  v - the multicomponent vector

    692:    Notes:
    693:    One must call VecSetBlockSize() before this routine to set the stride
    694:    information, or use a vector created from a multicomponent DMDA.

    696:    The parallel layout of the vector and the subvector must be the same;
    697:    i.e., nlocal of v = stride*(nlocal of s)

    699:    Not optimized; could be easily

    701:    Level: advanced

    703: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
    704:           VecStrideScatterAll()
    705: @*/
    706: PetscErrorCode  VecStrideScatterAll(Vec s[],Vec v,InsertMode addv)
    707: {
    708:   PetscErrorCode    ierr;
    709:   PetscInt          i,n,n2,bs,j,jj,k,*bss = NULL,nv,nvc;
    710:   PetscScalar       *x;
    711:   PetscScalar const **y;

    717:   VecGetLocalSize(v,&n);
    718:   VecGetLocalSize(s[0],&n2);
    719:   VecGetArray(v,&x);
    720:   VecGetBlockSize(v,&bs);
    721:   if (bs <= 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

    723:   PetscMalloc2(bs,(PetscScalar***)&y,bs,&bss);
    724:   nv   = 0;
    725:   nvc  = 0;
    726:   for (i=0; i<bs; i++) {
    727:     VecGetBlockSize(s[i],&bss[i]);
    728:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    729:     VecGetArrayRead(s[i],&y[i]);
    730:     nvc += bss[i];
    731:     nv++;
    732:     if (nvc > bs) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
    733:     if (nvc == bs) break;
    734:   }

    736:   n =  n/bs;

    738:   jj = 0;
    739:   if (addv == INSERT_VALUES) {
    740:     for (j=0; j<nv; j++) {
    741:       for (k=0; k<bss[j]; k++) {
    742:         for (i=0; i<n; i++) x[bs*i+jj+k] = y[j][i*bss[j] + k];
    743:       }
    744:       jj += bss[j];
    745:     }
    746:   } else if (addv == ADD_VALUES) {
    747:     for (j=0; j<nv; j++) {
    748:       for (k=0; k<bss[j]; k++) {
    749:         for (i=0; i<n; i++) x[bs*i+jj+k] += y[j][i*bss[j] + k];
    750:       }
    751:       jj += bss[j];
    752:     }
    753: #if !defined(PETSC_USE_COMPLEX)
    754:   } else if (addv == MAX_VALUES) {
    755:     for (j=0; j<nv; j++) {
    756:       for (k=0; k<bss[j]; k++) {
    757:         for (i=0; i<n; i++) x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
    758:       }
    759:       jj += bss[j];
    760:     }
    761: #endif
    762:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    764:   VecRestoreArray(v,&x);
    765:   for (i=0; i<nv; i++) {
    766:     VecRestoreArrayRead(s[i],&y[i]);
    767:   }
    768:   PetscFree2(*(PetscScalar***)&y,bss);
    769:   return(0);
    770: }

    772: /*@
    773:    VecStrideGather - Gathers a single component from a multi-component vector into
    774:    another vector.

    776:    Collective on Vec

    778:    Input Parameter:
    779: +  v - the vector
    780: .  start - starting point of the subvector (defined by a stride)
    781: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    783:    Output Parameter:
    784: .  s - the location where the subvector is stored

    786:    Notes:
    787:    One must call VecSetBlockSize() before this routine to set the stride
    788:    information, or use a vector created from a multicomponent DMDA.

    790:    If x is the array representing the vector x then this gathers
    791:    the array (x[start],x[start+stride],x[start+2*stride], ....)

    793:    The parallel layout of the vector and the subvector must be the same;
    794:    i.e., nlocal of v = stride*(nlocal of s)

    796:    Not optimized; could be easily

    798:    Level: advanced

    800: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    801:           VecStrideScatterAll()
    802: @*/
    803: PetscErrorCode  VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
    804: {

    810:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
    811:   if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
    812:   if (!v->ops->stridegather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
    813:   (*v->ops->stridegather)(v,start,s,addv);
    814:   return(0);
    815: }

    817: /*@
    818:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

    820:    Collective on Vec

    822:    Input Parameter:
    823: +  s - the single-component vector
    824: .  start - starting point of the subvector (defined by a stride)
    825: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    827:    Output Parameter:
    828: .  v - the location where the subvector is scattered (the multi-component vector)

    830:    Notes:
    831:    One must call VecSetBlockSize() on the multi-component vector before this
    832:    routine to set the stride  information, or use a vector created from a multicomponent DMDA.

    834:    The parallel layout of the vector and the subvector must be the same;
    835:    i.e., nlocal of v = stride*(nlocal of s)

    837:    Not optimized; could be easily

    839:    Level: advanced

    841: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    842:           VecStrideScatterAll(), VecStrideSubSetScatter(), VecStrideSubSetGather()
    843: @*/
    844: PetscErrorCode  VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
    845: {

    851:   if (start < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative start %D",start);
    852:   if (start >= v->map->bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,v->map->bs);
    853:   if (!v->ops->stridescatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
    854:   (*v->ops->stridescatter)(s,start,v,addv);
    855:   return(0);
    856: }

    858: /*@
    859:    VecStrideSubSetGather - Gathers a subset of components from a multi-component vector into
    860:    another vector.

    862:    Collective on Vec

    864:    Input Parameter:
    865: +  v - the vector
    866: .  nidx - the number of indices
    867: .  idxv - the indices of the components 0 <= idxv[0] ...idxv[nidx-1] < bs(v), they need not be sorted
    868: .  idxs - the indices of the components 0 <= idxs[0] ...idxs[nidx-1] < bs(s), they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
    869: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    871:    Output Parameter:
    872: .  s - the location where the subvector is stored

    874:    Notes:
    875:    One must call VecSetBlockSize() on both vectors before this routine to set the stride
    876:    information, or use a vector created from a multicomponent DMDA.

    879:    The parallel layout of the vector and the subvector must be the same;

    881:    Not optimized; could be easily

    883:    Level: advanced

    885: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideGather(), VecStrideSubSetScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    886:           VecStrideScatterAll()
    887: @*/
    888: PetscErrorCode  VecStrideSubSetGather(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
    889: {

    895:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
    896:   if (!v->ops->stridesubsetgather) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
    897:   (*v->ops->stridesubsetgather)(v,nidx,idxv,idxs,s,addv);
    898:   return(0);
    899: }

    901: /*@
    902:    VecStrideSubSetScatter - Scatters components from a vector into a subset of components of a multi-component vector.

    904:    Collective on Vec

    906:    Input Parameter:
    907: +  s - the smaller-component vector
    908: .  nidx - the number of indices in idx
    909: .  idxs - the indices of the components in the smaller-component vector, 0 <= idxs[0] ...idxs[nidx-1] < bs(s) they need not be sorted, may be null if nidx == bs(s) or is PETSC_DETERMINE
    910: .  idxv - the indices of the components in the larger-component vector, 0 <= idx[0] ...idx[nidx-1] < bs(v) they need not be sorted
    911: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    913:    Output Parameter:
    914: .  v - the location where the subvector is into scattered (the multi-component vector)

    916:    Notes:
    917:    One must call VecSetBlockSize() on the vectors before this
    918:    routine to set the stride  information, or use a vector created from a multicomponent DMDA.

    920:    The parallel layout of the vector and the subvector must be the same;

    922:    Not optimized; could be easily

    924:    Level: advanced

    926: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideGather(), VecStrideSubSetGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    927:           VecStrideScatterAll()
    928: @*/
    929: PetscErrorCode  VecStrideSubSetScatter(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
    930: {

    936:   if (nidx == PETSC_DETERMINE) nidx = s->map->bs;
    937:   if (!v->ops->stridesubsetscatter) SETERRQ(PetscObjectComm((PetscObject)s),PETSC_ERR_SUP,"Not implemented for this Vec class");
    938:   (*v->ops->stridesubsetscatter)(s,nidx,idxs,idxv,v,addv);
    939:   return(0);
    940: }

    942: PetscErrorCode  VecStrideGather_Default(Vec v,PetscInt start,Vec s,InsertMode addv)
    943: {
    945:   PetscInt       i,n,bs,ns;
    946:   const PetscScalar *x;
    947:   PetscScalar       *y;

    950:   VecGetLocalSize(v,&n);
    951:   VecGetLocalSize(s,&ns);
    952:   VecGetArrayRead(v,&x);
    953:   VecGetArray(s,&y);

    955:   bs = v->map->bs;
    956:   if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
    957:   x += start;
    958:   n  =  n/bs;

    960:   if (addv == INSERT_VALUES) {
    961:     for (i=0; i<n; i++) y[i] = x[bs*i];
    962:   } else if (addv == ADD_VALUES) {
    963:     for (i=0; i<n; i++) y[i] += x[bs*i];
    964: #if !defined(PETSC_USE_COMPLEX)
    965:   } else if (addv == MAX_VALUES) {
    966:     for (i=0; i<n; i++) y[i] = PetscMax(y[i],x[bs*i]);
    967: #endif
    968:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    970:   VecRestoreArrayRead(v,&x);
    971:   VecRestoreArray(s,&y);
    972:   return(0);
    973: }

    975: PetscErrorCode  VecStrideScatter_Default(Vec s,PetscInt start,Vec v,InsertMode addv)
    976: {
    977:   PetscErrorCode    ierr;
    978:   PetscInt          i,n,bs,ns;
    979:   PetscScalar       *x;
    980:   const PetscScalar *y;

    983:   VecGetLocalSize(v,&n);
    984:   VecGetLocalSize(s,&ns);
    985:   VecGetArray(v,&x);
    986:   VecGetArrayRead(s,&y);

    988:   bs = v->map->bs;
    989:   if (n != ns*bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
    990:   x += start;
    991:   n  =  n/bs;

    993:   if (addv == INSERT_VALUES) {
    994:     for (i=0; i<n; i++) x[bs*i] = y[i];
    995:   } else if (addv == ADD_VALUES) {
    996:     for (i=0; i<n; i++) x[bs*i] += y[i];
    997: #if !defined(PETSC_USE_COMPLEX)
    998:   } else if (addv == MAX_VALUES) {
    999:     for (i=0; i<n; i++) x[bs*i] = PetscMax(y[i],x[bs*i]);
    1000: #endif
    1001:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    1003:   VecRestoreArray(v,&x);
    1004:   VecRestoreArrayRead(s,&y);
    1005:   return(0);
    1006: }

    1008: PetscErrorCode  VecStrideSubSetGather_Default(Vec v,PetscInt nidx,const PetscInt idxv[],const PetscInt idxs[],Vec s,InsertMode addv)
    1009: {
    1010:   PetscErrorCode    ierr;
    1011:   PetscInt          i,j,n,bs,bss,ns;
    1012:   const PetscScalar *x;
    1013:   PetscScalar       *y;

    1016:   VecGetLocalSize(v,&n);
    1017:   VecGetLocalSize(s,&ns);
    1018:   VecGetArrayRead(v,&x);
    1019:   VecGetArray(s,&y);

    1021:   bs  = v->map->bs;
    1022:   bss = s->map->bs;
    1023:   n  =  n/bs;

    1025:   if (PetscDefined(USE_DEBUG)) {
    1026:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
    1027:     for (j=0; j<nidx; j++) {
    1028:       if (idxv[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxv[j]);
    1029:       if (idxv[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxv[j],bs);
    1030:     }
    1031:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not gathering into all locations");
    1032:   }

    1034:   if (addv == INSERT_VALUES) {
    1035:     if (!idxs) {
    1036:       for (i=0; i<n; i++) {
    1037:         for (j=0; j<bss; j++) y[bss*i+j] = x[bs*i+idxv[j]];
    1038:       }
    1039:     } else {
    1040:       for (i=0; i<n; i++) {
    1041:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = x[bs*i+idxv[j]];
    1042:       }
    1043:     }
    1044:   } else if (addv == ADD_VALUES) {
    1045:     if (!idxs) {
    1046:       for (i=0; i<n; i++) {
    1047:         for (j=0; j<bss; j++) y[bss*i+j] += x[bs*i+idxv[j]];
    1048:       }
    1049:     } else {
    1050:       for (i=0; i<n; i++) {
    1051:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] += x[bs*i+idxv[j]];
    1052:       }
    1053:     }
    1054: #if !defined(PETSC_USE_COMPLEX)
    1055:   } else if (addv == MAX_VALUES) {
    1056:     if (!idxs) {
    1057:       for (i=0; i<n; i++) {
    1058:         for (j=0; j<bss; j++) y[bss*i+j] = PetscMax(y[bss*i+j],x[bs*i+idxv[j]]);
    1059:       }
    1060:     } else {
    1061:       for (i=0; i<n; i++) {
    1062:         for (j=0; j<bss; j++) y[bss*i+idxs[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i+idxv[j]]);
    1063:       }
    1064:     }
    1065: #endif
    1066:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    1068:   VecRestoreArrayRead(v,&x);
    1069:   VecRestoreArray(s,&y);
    1070:   return(0);
    1071: }

    1073: PetscErrorCode  VecStrideSubSetScatter_Default(Vec s,PetscInt nidx,const PetscInt idxs[],const PetscInt idxv[],Vec v,InsertMode addv)
    1074: {
    1075:   PetscErrorCode    ierr;
    1076:   PetscInt          j,i,n,bs,ns,bss;
    1077:   PetscScalar       *x;
    1078:   const PetscScalar *y;

    1081:   VecGetLocalSize(v,&n);
    1082:   VecGetLocalSize(s,&ns);
    1083:   VecGetArray(v,&x);
    1084:   VecGetArrayRead(s,&y);

    1086:   bs  = v->map->bs;
    1087:   bss = s->map->bs;
    1088:   n  =  n/bs;

    1090:   if (PetscDefined(USE_DEBUG)) {
    1091:     if (n != ns/bss) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Incompatible layout of vectors");
    1092:     for (j=0; j<bss; j++) {
    1093:       if (idxs) {
    1094:         if (idxs[j] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is negative",j,idxs[j]);
    1095:         if (idxs[j] >= bs) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"idx[%D] %D is greater than or equal to vector blocksize %D",j,idxs[j],bs);
    1096:       }
    1097:     }
    1098:     if (!idxs && bss != nidx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Must provide idxs when not scattering from all locations");
    1099:   }

    1101:   if (addv == INSERT_VALUES) {
    1102:     if (!idxs) {
    1103:       for (i=0; i<n; i++) {
    1104:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+j];
    1105:       }
    1106:     } else {
    1107:       for (i=0; i<n; i++) {
    1108:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = y[bss*i+idxs[j]];
    1109:       }
    1110:     }
    1111:   } else if (addv == ADD_VALUES) {
    1112:     if (!idxs) {
    1113:       for (i=0; i<n; i++) {
    1114:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+j];
    1115:       }
    1116:     } else {
    1117:       for (i=0; i<n; i++) {
    1118:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] += y[bss*i+idxs[j]];
    1119:       }
    1120:     }
    1121: #if !defined(PETSC_USE_COMPLEX)
    1122:   } else if (addv == MAX_VALUES) {
    1123:     if (!idxs) {
    1124:       for (i=0; i<n; i++) {
    1125:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+j],x[bs*i + idxv[j]]);
    1126:       }
    1127:     } else {
    1128:       for (i=0; i<n; i++) {
    1129:         for (j=0; j<bss; j++) x[bs*i + idxv[j]] = PetscMax(y[bss*i+idxs[j]],x[bs*i + idxv[j]]);
    1130:       }
    1131:     }
    1132: #endif
    1133:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");

    1135:   VecRestoreArray(v,&x);
    1136:   VecRestoreArrayRead(s,&y);
    1137:   return(0);
    1138: }

    1140: PetscErrorCode VecReciprocal_Default(Vec v)
    1141: {
    1143:   PetscInt       i,n;
    1144:   PetscScalar    *x;

    1147:   VecGetLocalSize(v,&n);
    1148:   VecGetArray(v,&x);
    1149:   for (i=0; i<n; i++) {
    1150:     if (x[i] != (PetscScalar)0.0) x[i] = (PetscScalar)1.0/x[i];
    1151:   }
    1152:   VecRestoreArray(v,&x);
    1153:   return(0);
    1154: }

    1156: /*@
    1157:   VecExp - Replaces each component of a vector by e^x_i

    1159:   Not collective

    1161:   Input Parameter:
    1162: . v - The vector

    1164:   Output Parameter:
    1165: . v - The vector of exponents

    1167:   Level: beginner

    1169: .seealso:  VecLog(), VecAbs(), VecSqrtAbs(), VecReciprocal()

    1171: @*/
    1172: PetscErrorCode  VecExp(Vec v)
    1173: {
    1174:   PetscScalar    *x;
    1175:   PetscInt       i, n;

    1180:   if (v->ops->exp) {
    1181:     (*v->ops->exp)(v);
    1182:   } else {
    1183:     VecGetLocalSize(v, &n);
    1184:     VecGetArray(v, &x);
    1185:     for (i = 0; i < n; i++) x[i] = PetscExpScalar(x[i]);
    1186:     VecRestoreArray(v, &x);
    1187:   }
    1188:   return(0);
    1189: }

    1191: /*@
    1192:   VecLog - Replaces each component of a vector by log(x_i), the natural logarithm

    1194:   Not collective

    1196:   Input Parameter:
    1197: . v - The vector

    1199:   Output Parameter:
    1200: . v - The vector of logs

    1202:   Level: beginner

    1204: .seealso:  VecExp(), VecAbs(), VecSqrtAbs(), VecReciprocal()

    1206: @*/
    1207: PetscErrorCode  VecLog(Vec v)
    1208: {
    1209:   PetscScalar    *x;
    1210:   PetscInt       i, n;

    1215:   if (v->ops->log) {
    1216:     (*v->ops->log)(v);
    1217:   } else {
    1218:     VecGetLocalSize(v, &n);
    1219:     VecGetArray(v, &x);
    1220:     for (i = 0; i < n; i++) x[i] = PetscLogScalar(x[i]);
    1221:     VecRestoreArray(v, &x);
    1222:   }
    1223:   return(0);
    1224: }

    1226: /*@
    1227:   VecSqrtAbs - Replaces each component of a vector by the square root of its magnitude.

    1229:   Not collective

    1231:   Input Parameter:
    1232: . v - The vector

    1234:   Output Parameter:
    1235: . v - The vector square root

    1237:   Level: beginner

    1239:   Note: The actual function is sqrt(|x_i|)

    1241: .seealso: VecLog(), VecExp(), VecReciprocal(), VecAbs()

    1243: @*/
    1244: PetscErrorCode  VecSqrtAbs(Vec v)
    1245: {
    1246:   PetscScalar    *x;
    1247:   PetscInt       i, n;

    1252:   if (v->ops->sqrt) {
    1253:     (*v->ops->sqrt)(v);
    1254:   } else {
    1255:     VecGetLocalSize(v, &n);
    1256:     VecGetArray(v, &x);
    1257:     for (i = 0; i < n; i++) x[i] = PetscSqrtReal(PetscAbsScalar(x[i]));
    1258:     VecRestoreArray(v, &x);
    1259:   }
    1260:   return(0);
    1261: }

    1263: /*@
    1264:   VecDotNorm2 - computes the inner product of two vectors and the 2-norm squared of the second vector

    1266:   Collective on Vec

    1268:   Input Parameter:
    1269: + s - first vector
    1270: - t - second vector

    1272:   Output Parameter:
    1273: + dp - s'conj(t)
    1274: - nm - t'conj(t)

    1276:   Level: advanced

    1278:   Notes:
    1279:     conj(x) is the complex conjugate of x when x is complex

    1282: .seealso:   VecDot(), VecNorm(), VecDotBegin(), VecNormBegin(), VecDotEnd(), VecNormEnd()

    1284: @*/
    1285: PetscErrorCode  VecDotNorm2(Vec s,Vec t,PetscScalar *dp, PetscReal *nm)
    1286: {
    1287:   const PetscScalar *sx, *tx;
    1288:   PetscScalar       dpx = 0.0, nmx = 0.0,work[2],sum[2];
    1289:   PetscInt          i, n;
    1290:   PetscErrorCode    ierr;

    1300:   if (s->map->N != t->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1301:   if (s->map->n != t->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1303:   PetscLogEventBegin(VEC_DotNorm2,s,t,0,0);
    1304:   if (s->ops->dotnorm2) {
    1305:     (*s->ops->dotnorm2)(s,t,dp,&dpx);
    1306:     *nm  = PetscRealPart(dpx);
    1307:   } else {
    1308:     VecGetLocalSize(s, &n);
    1309:     VecGetArrayRead(s, &sx);
    1310:     VecGetArrayRead(t, &tx);

    1312:     for (i = 0; i<n; i++) {
    1313:       dpx += sx[i]*PetscConj(tx[i]);
    1314:       nmx += tx[i]*PetscConj(tx[i]);
    1315:     }
    1316:     work[0] = dpx;
    1317:     work[1] = nmx;

    1319:     MPIU_Allreduce(work,sum,2,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)s));
    1320:     *dp  = sum[0];
    1321:     *nm  = PetscRealPart(sum[1]);

    1323:     VecRestoreArrayRead(t, &tx);
    1324:     VecRestoreArrayRead(s, &sx);
    1325:     PetscLogFlops(4.0*n);
    1326:   }
    1327:   PetscLogEventEnd(VEC_DotNorm2,s,t,0,0);
    1328:   return(0);
    1329: }

    1331: /*@
    1332:    VecSum - Computes the sum of all the components of a vector.

    1334:    Collective on Vec

    1336:    Input Parameter:
    1337: .  v - the vector

    1339:    Output Parameter:
    1340: .  sum - the result

    1342:    Level: beginner

    1344: .seealso: VecNorm()
    1345: @*/
    1346: PetscErrorCode  VecSum(Vec v,PetscScalar *sum)
    1347: {
    1348:   PetscErrorCode    ierr;
    1349:   PetscInt          i,n;
    1350:   const PetscScalar *x;
    1351:   PetscScalar       lsum = 0.0;

    1356:   VecGetLocalSize(v,&n);
    1357:   VecGetArrayRead(v,&x);
    1358:   for (i=0; i<n; i++) lsum += x[i];
    1359:   MPIU_Allreduce(&lsum,sum,1,MPIU_SCALAR,MPIU_SUM,PetscObjectComm((PetscObject)v));
    1360:   VecRestoreArrayRead(v,&x);
    1361:   return(0);
    1362: }

    1364: /*@
    1365:    VecImaginaryPart - Replaces a complex vector with its imginary part

    1367:    Collective on Vec

    1369:    Input Parameter:
    1370: .  v - the vector

    1372:    Level: beginner

    1374: .seealso: VecNorm(), VecRealPart()
    1375: @*/
    1376: PetscErrorCode  VecImaginaryPart(Vec v)
    1377: {
    1378:   PetscErrorCode    ierr;
    1379:   PetscInt          i,n;
    1380:   PetscScalar       *x;

    1384:   VecGetLocalSize(v,&n);
    1385:   VecGetArray(v,&x);
    1386:   for (i=0; i<n; i++) x[i] = PetscImaginaryPart(x[i]);
    1387:   VecRestoreArray(v,&x);
    1388:   return(0);
    1389: }

    1391: /*@
    1392:    VecRealPart - Replaces a complex vector with its real part

    1394:    Collective on Vec

    1396:    Input Parameter:
    1397: .  v - the vector

    1399:    Level: beginner

    1401: .seealso: VecNorm(), VecImaginaryPart()
    1402: @*/
    1403: PetscErrorCode  VecRealPart(Vec v)
    1404: {
    1405:   PetscErrorCode    ierr;
    1406:   PetscInt          i,n;
    1407:   PetscScalar       *x;

    1411:   VecGetLocalSize(v,&n);
    1412:   VecGetArray(v,&x);
    1413:   for (i=0; i<n; i++) x[i] = PetscRealPart(x[i]);
    1414:   VecRestoreArray(v,&x);
    1415:   return(0);
    1416: }

    1418: /*@
    1419:    VecShift - Shifts all of the components of a vector by computing
    1420:    x[i] = x[i] + shift.

    1422:    Logically Collective on Vec

    1424:    Input Parameters:
    1425: +  v - the vector
    1426: -  shift - the shift

    1428:    Output Parameter:
    1429: .  v - the shifted vector

    1431:    Level: intermediate

    1433: @*/
    1434: PetscErrorCode  VecShift(Vec v,PetscScalar shift)
    1435: {
    1437:   PetscInt       i,n;
    1438:   PetscScalar    *x;

    1443:   VecSetErrorIfLocked(v,1);
    1444:   if (shift == 0.0) return(0);

    1446:   if (v->ops->shift) {
    1447:     (*v->ops->shift)(v,shift);
    1448:   } else {
    1449:     VecGetLocalSize(v,&n);
    1450:     VecGetArray(v,&x);
    1451:     for (i=0; i<n; i++) x[i] += shift;
    1452:     VecRestoreArray(v,&x);
    1453:   }
    1454:   return(0);
    1455: }

    1457: /*@
    1458:    VecAbs - Replaces every element in a vector with its absolute value.

    1460:    Logically Collective on Vec

    1462:    Input Parameters:
    1463: .  v - the vector

    1465:    Level: intermediate

    1467: @*/
    1468: PetscErrorCode  VecAbs(Vec v)
    1469: {
    1471:   PetscInt       i,n;
    1472:   PetscScalar    *x;

    1476:   VecSetErrorIfLocked(v,1);

    1478:   if (v->ops->abs) {
    1479:     (*v->ops->abs)(v);
    1480:   } else {
    1481:     VecGetLocalSize(v,&n);
    1482:     VecGetArray(v,&x);
    1483:     for (i=0; i<n; i++) x[i] = PetscAbsScalar(x[i]);
    1484:     VecRestoreArray(v,&x);
    1485:   }
    1486:   return(0);
    1487: }

    1489: /*@
    1490:   VecPermute - Permutes a vector in place using the given ordering.

    1492:   Input Parameters:
    1493: + vec   - The vector
    1494: . order - The ordering
    1495: - inv   - The flag for inverting the permutation

    1497:   Level: beginner

    1499:   Note: This function does not yet support parallel Index Sets with non-local permutations

    1501: .seealso: MatPermute()
    1502: @*/
    1503: PetscErrorCode  VecPermute(Vec x, IS row, PetscBool inv)
    1504: {
    1505:   const PetscScalar *array;
    1506:   PetscScalar       *newArray;
    1507:   const PetscInt    *idx;
    1508:   PetscInt          i,rstart,rend;
    1509:   PetscErrorCode    ierr;

    1514:   VecSetErrorIfLocked(x,1);
    1515:   VecGetOwnershipRange(x,&rstart,&rend);
    1516:   ISGetIndices(row, &idx);
    1517:   VecGetArrayRead(x, &array);
    1518:   PetscMalloc1(x->map->n, &newArray);
    1519:   if (PetscDefined(USE_DEBUG)) {
    1520:     for (i = 0; i < x->map->n; i++) {
    1521:       if ((idx[i] < rstart) || (idx[i] >= rend)) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
    1522:     }
    1523:   }
    1524:   if (!inv) {
    1525:     for (i = 0; i < x->map->n; i++) newArray[i] = array[idx[i]-rstart];
    1526:   } else {
    1527:     for (i = 0; i < x->map->n; i++) newArray[idx[i]-rstart] = array[i];
    1528:   }
    1529:   VecRestoreArrayRead(x, &array);
    1530:   ISRestoreIndices(row, &idx);
    1531:   VecReplaceArray(x, newArray);
    1532:   return(0);
    1533: }

    1535: /*@
    1536:    VecEqual - Compares two vectors. Returns true if the two vectors are either pointing to the same memory buffer,
    1537:    or if the two vectors have the same local and global layout as well as bitwise equality of all entries.
    1538:    Does NOT take round-off errors into account.

    1540:    Collective on Vec

    1542:    Input Parameters:
    1543: +  vec1 - the first vector
    1544: -  vec2 - the second vector

    1546:    Output Parameter:
    1547: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

    1549:    Level: intermediate

    1552: @*/
    1553: PetscErrorCode  VecEqual(Vec vec1,Vec vec2,PetscBool  *flg)
    1554: {
    1555:   const PetscScalar  *v1,*v2;
    1556:   PetscErrorCode     ierr;
    1557:   PetscInt           n1,n2,N1,N2;
    1558:   PetscBool          flg1;

    1564:   if (vec1 == vec2) *flg = PETSC_TRUE;
    1565:   else {
    1566:     VecGetSize(vec1,&N1);
    1567:     VecGetSize(vec2,&N2);
    1568:     if (N1 != N2) flg1 = PETSC_FALSE;
    1569:     else {
    1570:       VecGetLocalSize(vec1,&n1);
    1571:       VecGetLocalSize(vec2,&n2);
    1572:       if (n1 != n2) flg1 = PETSC_FALSE;
    1573:       else {
    1574:         VecGetArrayRead(vec1,&v1);
    1575:         VecGetArrayRead(vec2,&v2);
    1576:         PetscArraycmp(v1,v2,n1,&flg1);
    1577:         VecRestoreArrayRead(vec1,&v1);
    1578:         VecRestoreArrayRead(vec2,&v2);
    1579:       }
    1580:     }
    1581:     /* combine results from all processors */
    1582:     MPIU_Allreduce(&flg1,flg,1,MPIU_BOOL,MPI_MIN,PetscObjectComm((PetscObject)vec1));
    1583:   }
    1584:   return(0);
    1585: }

    1587: /*@
    1588:    VecUniqueEntries - Compute the number of unique entries, and those entries

    1590:    Collective on Vec

    1592:    Input Parameter:
    1593: .  vec - the vector

    1595:    Output Parameters:
    1596: +  n - The number of unique entries
    1597: -  e - The entries

    1599:    Level: intermediate

    1601: @*/
    1602: PetscErrorCode  VecUniqueEntries(Vec vec, PetscInt *n, PetscScalar **e)
    1603: {
    1604:   const PetscScalar *v;
    1605:   PetscScalar       *tmp, *vals;
    1606:   PetscMPIInt       *N, *displs, l;
    1607:   PetscInt          ng, m, i, j, p;
    1608:   PetscMPIInt       size;
    1609:   PetscErrorCode    ierr;

    1614:   MPI_Comm_size(PetscObjectComm((PetscObject) vec), &size);
    1615:   VecGetLocalSize(vec, &m);
    1616:   VecGetArrayRead(vec, &v);
    1617:   PetscMalloc2(m,&tmp,size,&N);
    1618:   for (i = 0, j = 0, l = 0; i < m; ++i) {
    1619:     /* Can speed this up with sorting */
    1620:     for (j = 0; j < l; ++j) {
    1621:       if (v[i] == tmp[j]) break;
    1622:     }
    1623:     if (j == l) {
    1624:       tmp[j] = v[i];
    1625:       ++l;
    1626:     }
    1627:   }
    1628:   VecRestoreArrayRead(vec, &v);
    1629:   /* Gather serial results */
    1630:   MPI_Allgather(&l, 1, MPI_INT, N, 1, MPI_INT, PetscObjectComm((PetscObject) vec));
    1631:   for (p = 0, ng = 0; p < size; ++p) {
    1632:     ng += N[p];
    1633:   }
    1634:   PetscMalloc2(ng,&vals,size+1,&displs);
    1635:   for (p = 1, displs[0] = 0; p <= size; ++p) {
    1636:     displs[p] = displs[p-1] + N[p-1];
    1637:   }
    1638:   MPI_Allgatherv(tmp, l, MPIU_SCALAR, vals, N, displs, MPIU_SCALAR, PetscObjectComm((PetscObject) vec));
    1639:   /* Find unique entries */
    1640: #ifdef PETSC_USE_COMPLEX
    1641:   SETERRQ(PetscObjectComm((PetscObject) vec), PETSC_ERR_SUP, "Does not work with complex numbers");
    1642: #else
    1643:   *n = displs[size];
    1644:   PetscSortRemoveDupsReal(n, (PetscReal *) vals);
    1645:   if (e) {
    1647:     PetscMalloc1(*n, e);
    1648:     for (i = 0; i < *n; ++i) {
    1649:       (*e)[i] = vals[i];
    1650:     }
    1651:   }
    1652:   PetscFree2(vals,displs);
    1653:   PetscFree2(tmp,N);
    1654:   return(0);
    1655: #endif
    1656: }