Actual source code: vinv.c

      1: #define PETSCVEC_DLL

      3: /*
      4:      Some useful vector utility functions.
      5: */
     6:  #include vecimpl.h

     13: /*@
     14:    VecStrideScale - Scales a subvector of a vector defined
     15:    by a starting point and a stride.

     17:    Collective on Vec

     19:    Input Parameter:
     20: +  v - the vector
     21: .  start - starting point of the subvector (defined by a stride)
     22: -  scale - value to multiply each subvector entry by

     24:    Notes:
     25:    One must call VecSetBlockSize() before this routine to set the stride
     26:    information, or use a vector created from a multicomponent DA.

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

     30:    Level: advanced

     32:    Concepts: scale^on stride of vector
     33:    Concepts: stride^scale

     35: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideScale()
     36: @*/
     37: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScale(Vec v,PetscInt start,PetscScalar *scale)
     38: {
     40:   PetscInt       i,n,bs;
     41:   PetscScalar    *x,xscale = *scale;

     46:   VecGetLocalSize(v,&n);
     47:   VecGetArray(v,&x);

     49:   bs   = v->bs;
     50:   if (start >= bs) {
     51:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
     52:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
     53:   }
     54:   x += start;

     56:   for (i=0; i<n; i+=bs) {
     57:     x[i] *= xscale;
     58:   }
     59:   x -= start;
     60:   VecRestoreArray(v,&x);
     61:   return(0);
     62: }

     66: /*@C
     67:    VecStrideNorm - Computes the norm of subvector of a vector defined
     68:    by a starting point and a stride.

     70:    Collective on Vec

     72:    Input Parameter:
     73: +  v - the vector
     74: .  start - starting point of the subvector (defined by a stride)
     75: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

     77:    Output Parameter:
     78: .  norm - the norm

     80:    Notes:
     81:    One must call VecSetBlockSize() before this routine to set the stride
     82:    information, or use a vector created from a multicomponent DA.

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

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

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

     92:    Level: advanced

     94:    Concepts: norm^on stride of vector
     95:    Concepts: stride^norm

     97: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
     98: @*/
     99: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideNorm(Vec v,PetscInt start,NormType ntype,PetscReal *nrm)
    100: {
    102:   PetscInt       i,n,bs;
    103:   PetscScalar    *x;
    104:   PetscReal      tnorm;
    105:   MPI_Comm       comm;

    110:   VecGetLocalSize(v,&n);
    111:   VecGetArray(v,&x);
    112:   PetscObjectGetComm((PetscObject)v,&comm);

    114:   bs   = v->bs;
    115:   if (start >= bs) {
    116:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
    117:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    118:   }
    119:   x += start;

    121:   if (ntype == NORM_2) {
    122:     PetscScalar sum = 0.0;
    123:     for (i=0; i<n; i+=bs) {
    124:       sum += x[i]*(PetscConj(x[i]));
    125:     }
    126:     tnorm  = PetscRealPart(sum);
    127:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);
    128:     *nrm = sqrt(*nrm);
    129:   } else if (ntype == NORM_1) {
    130:     tnorm = 0.0;
    131:     for (i=0; i<n; i+=bs) {
    132:       tnorm += PetscAbsScalar(x[i]);
    133:     }
    134:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_SUM,comm);
    135:   } else if (ntype == NORM_INFINITY) {
    136:     PetscReal tmp;
    137:     tnorm = 0.0;

    139:     for (i=0; i<n; i+=bs) {
    140:       if ((tmp = PetscAbsScalar(x[i])) > tnorm) tnorm = tmp;
    141:       /* check special case of tmp == NaN */
    142:       if (tmp != tmp) {tnorm = tmp; break;}
    143:     }
    144:     MPI_Allreduce(&tnorm,nrm,1,MPIU_REAL,MPI_MAX,comm);
    145:   } else {
    146:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
    147:   }

    149:   VecRestoreArray(v,&x);
    150:   return(0);
    151: }

    155: /*@
    156:    VecStrideMax - Computes the maximum of subvector of a vector defined
    157:    by a starting point and a stride and optionally its location.

    159:    Collective on Vec

    161:    Input Parameter:
    162: +  v - the vector
    163: -  start - starting point of the subvector (defined by a stride)

    165:    Output Parameter:
    166: +  index - the location where the maximum occurred  (pass PETSC_NULL if not required)
    167: -  nrm - the max

    169:    Notes:
    170:    One must call VecSetBlockSize() before this routine to set the stride
    171:    information, or use a vector created from a multicomponent DA.

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

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

    180:    Level: advanced

    182:    Concepts: maximum^on stride of vector
    183:    Concepts: stride^maximum

    185: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
    186: @*/
    187: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMax(Vec v,PetscInt start,PetscInt *idex,PetscReal *nrm)
    188: {
    190:   PetscInt       i,n,bs,id;
    191:   PetscScalar    *x;
    192:   PetscReal      max,tmp;
    193:   MPI_Comm       comm;

    199:   VecGetLocalSize(v,&n);
    200:   VecGetArray(v,&x);
    201:   PetscObjectGetComm((PetscObject)v,&comm);

    203:   bs   = v->bs;
    204:   if (start >= bs) {
    205:     SETERRQ2(PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n\
    206:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    207:   }
    208:   x += start;

    210:   id = -1;
    211:   if (!n) {
    212:     max = PETSC_MIN;
    213:   } else {
    214: #if defined(PETSC_USE_COMPLEX)
    215:     max = PetscRealPart(x[0]);
    216: #else
    217:     max = x[0];
    218: #endif
    219:     for (i=bs; i<n; i+=bs) {
    220: #if defined(PETSC_USE_COMPLEX)
    221:       if ((tmp = PetscRealPart(x[i])) > max) { max = tmp; id = i;}
    222: #else
    223:       if ((tmp = x[i]) > max) { max = tmp; id = i;}
    224: #endif
    225:     }
    226:   }
    227:   VecRestoreArray(v,&x);

    229:   if (!idex) {
    230:     MPI_Allreduce(&max,nrm,1,MPIU_REAL,MPI_MAX,comm);
    231:   } else {
    232:     PetscReal in[2],out[2];
    233:     PetscInt  rstart;

    235:     VecGetOwnershipRange(v,&rstart,PETSC_NULL);
    236:     in[0] = max;
    237:     in[1] = rstart+id;
    238:     MPI_Allreduce(in,out,2,MPIU_REAL,VecMax_Local_Op,v->comm);
    239:     *nrm  = out[0];
    240:     *idex = (PetscInt)out[1];
    241:   }

    243:   return(0);
    244: }

    248: /*@
    249:    VecStrideMin - Computes the minimum of subvector of a vector defined
    250:    by a starting point and a stride and optionally its location.

    252:    Collective on Vec

    254:    Input Parameter:
    255: +  v - the vector
    256: -  start - starting point of the subvector (defined by a stride)

    258:    Output Parameter:
    259: +  idex - the location where the minimum occurred. (pass PETSC_NULL if not required)
    260: -  nrm - the min

    262:    Level: advanced

    264:    Notes:
    265:    One must call VecSetBlockSize() before this routine to set the stride
    266:    information, or use a vector created from a multicomponent DA.

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

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

    275:    Concepts: minimum^on stride of vector
    276:    Concepts: stride^minimum

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

    292:   VecGetLocalSize(v,&n);
    293:   VecGetArray(v,&x);
    294:   PetscObjectGetComm((PetscObject)v,&comm);

    296:   bs   = v->bs;
    297:   if (start >= bs) {
    298:     SETERRQ2(PETSC_ERR_ARG_WRONG,"Start of stride subvector (%D) is too large for stride\n\
    299:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    300:   }
    301:   x += start;

    303:   id = -1;
    304:   if (!n) {
    305:     min = PETSC_MAX;
    306:   } else {
    307: #if defined(PETSC_USE_COMPLEX)
    308:     min = PetscRealPart(x[0]);
    309: #else
    310:     min = x[0];
    311: #endif
    312:     for (i=bs; i<n; i+=bs) {
    313: #if defined(PETSC_USE_COMPLEX)
    314:       if ((tmp = PetscRealPart(x[i])) < min) { min = tmp; id = i;}
    315: #else
    316:       if ((tmp = x[i]) < min) { min = tmp; id = i;}
    317: #endif
    318:     }
    319:   }
    320:   VecRestoreArray(v,&x);

    322:   if (!idex) {
    323:     MPI_Allreduce(&min,nrm,1,MPIU_REAL,MPI_MIN,comm);
    324:   } else {
    325:     PetscReal in[2],out[2];
    326:     PetscInt  rstart;

    328:     VecGetOwnershipRange(v,&rstart,PETSC_NULL);
    329:     in[0] = min;
    330:     in[1] = rstart+id;
    331:     MPI_Allreduce(in,out,2,MPIU_REAL,VecMin_Local_Op,v->comm);
    332:     *nrm  = out[0];
    333:     *idex = (PetscInt)out[1];
    334:   }

    336:   return(0);
    337: }

    341: /*@
    342:    VecStrideScaleAll - Scales the subvectors of a vector defined
    343:    by a starting point and a stride.

    345:    Collective on Vec

    347:    Input Parameter:
    348: +  v - the vector
    349: -  scales - values to multiply each subvector entry by

    351:    Notes:
    352:    One must call VecSetBlockSize() before this routine to set the stride
    353:    information, or use a vector created from a multicomponent DA.

    356:    Level: advanced

    358:    Concepts: scale^on stride of vector
    359:    Concepts: stride^scale

    361: .seealso: VecNorm(), VecStrideScale(), VecScale(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
    362: @*/
    363: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScaleAll(Vec v,PetscScalar *scales)
    364: {
    366:   PetscInt         i,j,n,bs;
    367:   PetscScalar *x;

    372:   VecGetLocalSize(v,&n);
    373:   VecGetArray(v,&x);

    375:   bs   = v->bs;

    377:   /* need to provide optimized code for each bs */
    378:   for (i=0; i<n; i+=bs) {
    379:     for (j=0; j<bs; j++) {
    380:       x[i+j] *= scales[j];
    381:     }
    382:   }
    383:   VecRestoreArray(v,&x);
    384:   return(0);
    385: }

    389: /*@
    390:    VecStrideNormAll - Computes the norms  subvectors of a vector defined
    391:    by a starting point and a stride.

    393:    Collective on Vec

    395:    Input Parameter:
    396: +  v - the vector
    397: -  ntype - type of norm, one of NORM_1, NORM_2, NORM_INFINITY

    399:    Output Parameter:
    400: .  nrm - the norms

    402:    Notes:
    403:    One must call VecSetBlockSize() before this routine to set the stride
    404:    information, or use a vector created from a multicomponent DA.

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

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

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

    414:    Level: advanced

    416:    Concepts: norm^on stride of vector
    417:    Concepts: stride^norm

    419: .seealso: VecNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin(), VecStrideMax()
    420: @*/
    421: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideNormAll(Vec v,NormType ntype,PetscReal *nrm)
    422: {
    424:   PetscInt         i,j,n,bs;
    425:   PetscScalar *x;
    426:   PetscReal   tnorm[128];
    427:   MPI_Comm    comm;

    432:   VecGetLocalSize(v,&n);
    433:   VecGetArray(v,&x);
    434:   PetscObjectGetComm((PetscObject)v,&comm);

    436:   bs   = v->bs;
    437:   if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    439:   if (ntype == NORM_2) {
    440:     PetscScalar sum[128];
    441:     for (j=0; j<bs; j++) sum[j] = 0.0;
    442:     for (i=0; i<n; i+=bs) {
    443:       for (j=0; j<bs; j++) {
    444:         sum[j] += x[i+j]*(PetscConj(x[i+j]));
    445:       }
    446:     }
    447:     for (j=0; j<bs; j++) {
    448:       tnorm[j]  = PetscRealPart(sum[j]);
    449:     }
    450:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_SUM,comm);
    451:     for (j=0; j<bs; j++) {
    452:       nrm[j] = sqrt(nrm[j]);
    453:     }
    454:   } else if (ntype == NORM_1) {
    455:     for (j=0; j<bs; j++) {
    456:       tnorm[j] = 0.0;
    457:     }
    458:     for (i=0; i<n; i+=bs) {
    459:       for (j=0; j<bs; j++) {
    460:         tnorm[j] += PetscAbsScalar(x[i+j]);
    461:       }
    462:     }
    463:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_SUM,comm);
    464:   } else if (ntype == NORM_INFINITY) {
    465:     PetscReal tmp;
    466:     for (j=0; j<bs; j++) {
    467:       tnorm[j] = 0.0;
    468:     }

    470:     for (i=0; i<n; i+=bs) {
    471:       for (j=0; j<bs; j++) {
    472:         if ((tmp = PetscAbsScalar(x[i+j])) > tnorm[j]) tnorm[j] = tmp;
    473:         /* check special case of tmp == NaN */
    474:         if (tmp != tmp) {tnorm[j] = tmp; break;}
    475:       }
    476:     }
    477:     MPI_Allreduce(tnorm,nrm,bs,MPIU_REAL,MPI_MAX,comm);
    478:   } else {
    479:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown norm type");
    480:   }

    482:   VecRestoreArray(v,&x);
    483:   return(0);
    484: }

    488: /*@
    489:    VecStrideMaxAll - Computes the maximums of subvectors of a vector defined
    490:    by a starting point and a stride and optionally its location.

    492:    Collective on Vec

    494:    Input Parameter:
    495: .  v - the vector

    497:    Output Parameter:
    498: +  index - the location where the maximum occurred (not supported, pass PETSC_NULL,
    499:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
    500: -  nrm - the maximums

    502:    Notes:
    503:    One must call VecSetBlockSize() before this routine to set the stride
    504:    information, or use a vector created from a multicomponent DA.

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

    510:    Level: advanced

    512:    Concepts: maximum^on stride of vector
    513:    Concepts: stride^maximum

    515: .seealso: VecMax(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMin()
    516: @*/
    517: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMaxAll(Vec v,PetscInt *idex,PetscReal *nrm)
    518: {
    520:   PetscInt         i,j,n,bs;
    521:   PetscScalar *x;
    522:   PetscReal   max[128],tmp;
    523:   MPI_Comm    comm;

    528:   if (idex) {
    529:     SETERRQ(PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
    530:   }
    531:   VecGetLocalSize(v,&n);
    532:   VecGetArray(v,&x);
    533:   PetscObjectGetComm((PetscObject)v,&comm);

    535:   bs   = v->bs;
    536:   if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    538:   if (!n) {
    539:     for (j=0; j<bs; j++) {
    540:       max[j] = PETSC_MIN;
    541:     }
    542:   } else {
    543:     for (j=0; j<bs; j++) {
    544: #if defined(PETSC_USE_COMPLEX)
    545:       max[j] = PetscRealPart(x[j]);
    546: #else
    547:       max[j] = x[j];
    548: #endif
    549:     }
    550:     for (i=bs; i<n; i+=bs) {
    551:       for (j=0; j<bs; j++) {
    552: #if defined(PETSC_USE_COMPLEX)
    553:         if ((tmp = PetscRealPart(x[i+j])) > max[j]) { max[j] = tmp;}
    554: #else
    555:         if ((tmp = x[i+j]) > max[j]) { max[j] = tmp; }
    556: #endif
    557:       }
    558:     }
    559:   }
    560:   MPI_Allreduce(max,nrm,bs,MPIU_REAL,MPI_MAX,comm);

    562:   VecRestoreArray(v,&x);
    563:   return(0);
    564: }

    568: /*@
    569:    VecStrideMinAll - Computes the minimum of subvector of a vector defined
    570:    by a starting point and a stride and optionally its location.

    572:    Collective on Vec

    574:    Input Parameter:
    575: .  v - the vector

    577:    Output Parameter:
    578: +  idex - the location where the minimum occurred (not supported, pass PETSC_NULL,
    579:            if you need this, send mail to petsc-maint@mcs.anl.gov to request it)
    580: -  nrm - the minimums

    582:    Level: advanced

    584:    Notes:
    585:    One must call VecSetBlockSize() before this routine to set the stride
    586:    information, or use a vector created from a multicomponent DA.

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

    592:    Concepts: minimum^on stride of vector
    593:    Concepts: stride^minimum

    595: .seealso: VecMin(), VecStrideNorm(), VecStrideGather(), VecStrideScatter(), VecStrideMax()
    596: @*/
    597: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideMinAll(Vec v,PetscInt *idex,PetscReal *nrm)
    598: {
    600:   PetscInt         i,n,bs,j;
    601:   PetscScalar *x;
    602:   PetscReal   min[128],tmp;
    603:   MPI_Comm    comm;

    608:   if (idex) {
    609:     SETERRQ(PETSC_ERR_SUP,"No support yet for returning index; send mail to petsc-maint@mcs.anl.gov asking for it");
    610:   }
    611:   VecGetLocalSize(v,&n);
    612:   VecGetArray(v,&x);
    613:   PetscObjectGetComm((PetscObject)v,&comm);

    615:   bs   = v->bs;
    616:   if (bs > 128) SETERRQ(PETSC_ERR_SUP,"Currently supports only blocksize up to 128");

    618:   if (!n) {
    619:     for (j=0; j<bs; j++) {
    620:       min[j] = PETSC_MAX;
    621:     }
    622:   } else {
    623:     for (j=0; j<bs; j++) {
    624: #if defined(PETSC_USE_COMPLEX)
    625:       min[j] = PetscRealPart(x[j]);
    626: #else
    627:       min[j] = x[j];
    628: #endif
    629:     }
    630:     for (i=bs; i<n; i+=bs) {
    631:       for (j=0; j<bs; j++) {
    632: #if defined(PETSC_USE_COMPLEX)
    633:         if ((tmp = PetscRealPart(x[i+j])) < min[j]) { min[j] = tmp;}
    634: #else
    635:         if ((tmp = x[i+j]) < min[j]) { min[j] = tmp; }
    636: #endif
    637:       }
    638:     }
    639:   }
    640:   MPI_Allreduce(min,nrm,bs,MPIU_REAL,MPI_MIN,comm);

    642:   VecRestoreArray(v,&x);
    643:   return(0);
    644: }

    646: /*----------------------------------------------------------------------------------------------*/
    649: /*@
    650:    VecStrideGatherAll - Gathers all the single components from a multi-component vector into
    651:    seperate vectors.

    653:    Collective on Vec

    655:    Input Parameter:
    656: +  v - the vector
    657: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    659:    Output Parameter:
    660: .  s - the location where the subvectors are stored

    662:    Notes:
    663:    One must call VecSetBlockSize() before this routine to set the stride
    664:    information, or use a vector created from a multicomponent DA.

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

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

    673:    Not optimized; could be easily

    675:    Level: advanced

    677:    Concepts: gather^into strided vector

    679: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
    680:           VecStrideScatterAll()
    681: @*/
    682: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideGatherAll(Vec v,Vec *s,InsertMode addv)
    683: {
    685:   PetscInt       i,n,bs,j,k,*bss,nv,jj,nvc;
    686:   PetscScalar    *x,**y;

    692:   VecGetLocalSize(v,&n);
    693:   VecGetArray(v,&x);
    694:   bs   = v->bs;
    695:   if (bs < 0) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

    697:   PetscMalloc2(bs,PetscReal*,&y,bs,PetscInt,&bss);
    698:   nv   = 0;
    699:   nvc  = 0;
    700:   for (i=0; i<bs; i++) {
    701:     VecGetBlockSize(s[i],&bss[i]);
    702:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    703:     VecGetArray(s[i],&y[i]);
    704:     nvc  += bss[i];
    705:     nv++;
    706:     if (nvc > bs)  SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
    707:     if (nvc == bs) break;
    708:   }

    710:   n =  n/bs;

    712:   jj = 0;
    713:   if (addv == INSERT_VALUES) {
    714:     for (j=0; j<nv; j++) {
    715:       for (k=0; k<bss[j]; k++) {
    716:         for (i=0; i<n; i++) {
    717:           y[j][i*bss[j] + k] = x[bs*i+jj+k];
    718:         }
    719:       }
    720:       jj += bss[j];
    721:     }
    722:   } else if (addv == ADD_VALUES) {
    723:     for (j=0; j<nv; j++) {
    724:       for (k=0; k<bss[j]; k++) {
    725:         for (i=0; i<n; i++) {
    726:           y[j][i*bss[j] + k] += x[bs*i+jj+k];
    727:         }
    728:       }
    729:       jj += bss[j];
    730:     }
    731: #if !defined(PETSC_USE_COMPLEX)
    732:   } else if (addv == MAX_VALUES) {
    733:     for (j=0; j<nv; j++) {
    734:       for (k=0; k<bss[j]; k++) {
    735:         for (i=0; i<n; i++) {
    736:           y[j][i*bss[j] + k] = PetscMax(y[j][i*bss[j] + k],x[bs*i+jj+k]);
    737:         }
    738:       }
    739:       jj += bss[j];
    740:     }
    741: #endif
    742:   } else {
    743:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
    744:   }

    746:   VecRestoreArray(v,&x);
    747:   for (i=0; i<nv; i++) {
    748:     VecRestoreArray(s[i],&y[i]);
    749:   }
    750:   PetscFree2(y,bss);
    751:   return(0);
    752: }

    756: /*@
    757:    VecStrideScatterAll - Scatters all the single components from seperate vectors into
    758:      a multi-component vector.

    760:    Collective on Vec

    762:    Input Parameter:
    763: +  s - the location where the subvectors are stored
    764: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

    766:    Output Parameter:
    767: .  v - the multicomponent vector

    769:    Notes:
    770:    One must call VecSetBlockSize() before this routine to set the stride
    771:    information, or use a vector created from a multicomponent DA.

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

    776:    Not optimized; could be easily

    778:    Level: advanced

    780:    Concepts:  scatter^into strided vector

    782: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGather(),
    783:           VecStrideScatterAll()
    784: @*/
    785: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScatterAll(Vec *s,Vec v,InsertMode addv)
    786: {
    788:   PetscInt        i,n,bs,j,jj,k,*bss,nv,nvc;
    789:   PetscScalar     *x,**y;

    795:   VecGetLocalSize(v,&n);
    796:   VecGetArray(v,&x);
    797:   bs   = v->bs;
    798:   if (bs < 0) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Input vector does not have a valid blocksize set");

    800:   PetscMalloc2(bs,PetscScalar**,&y,bs,PetscInt,&bss);
    801:   nv  = 0;
    802:   nvc = 0;
    803:   for (i=0; i<bs; i++) {
    804:     VecGetBlockSize(s[i],&bss[i]);
    805:     if (bss[i] < 1) bss[i] = 1; /* if user never set it then assume 1  Re: [PETSC #8241] VecStrideGatherAll */
    806:     VecGetArray(s[i],&y[i]);
    807:     nvc  += bss[i];
    808:     nv++;
    809:     if (nvc > bs)  SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of subvectors in subvectors > number of vectors in main vector");
    810:     if (nvc == bs) break;
    811:   }

    813:   n =  n/bs;

    815:   jj = 0;
    816:   if (addv == INSERT_VALUES) {
    817:     for (j=0; j<nv; j++) {
    818:       for (k=0; k<bss[j]; k++) {
    819:         for (i=0; i<n; i++) {
    820:           x[bs*i+jj+k] = y[j][i*bss[j] + k];
    821:         }
    822:       }
    823:       jj += bss[j];
    824:     }
    825:   } else if (addv == ADD_VALUES) {
    826:     for (j=0; j<nv; j++) {
    827:       for (k=0; k<bss[j]; k++) {
    828:         for (i=0; i<n; i++) {
    829:           x[bs*i+jj+k] += y[j][i*bss[j] + k];
    830:         }
    831:       }
    832:       jj += bss[j];
    833:     }
    834: #if !defined(PETSC_USE_COMPLEX)
    835:   } else if (addv == MAX_VALUES) {
    836:     for (j=0; j<nv; j++) {
    837:       for (k=0; k<bss[j]; k++) {
    838:         for (i=0; i<n; i++) {
    839:           x[bs*i+jj+k] = PetscMax(x[bs*i+jj+k],y[j][i*bss[j] + k]);
    840:         }
    841:       }
    842:       jj += bss[j];
    843:     }
    844: #endif
    845:   } else {
    846:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
    847:   }

    849:   VecRestoreArray(v,&x);
    850:   for (i=0; i<nv; i++) {
    851:     VecRestoreArray(s[i],&y[i]);
    852:   }
    853:   PetscFree2(y,bss);
    854:   return(0);
    855: }

    859: /*@
    860:    VecStrideGather - Gathers a single component from a multi-component vector into
    861:    another vector.

    863:    Collective on Vec

    865:    Input Parameter:
    866: +  v - the vector
    867: .  start - starting point of the subvector (defined by a stride)
    868: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

    873:    Notes:
    874:    One must call VecSetBlockSize() before this routine to set the stride
    875:    information, or use a vector created from a multicomponent DA.

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

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

    883:    Not optimized; could be easily

    885:    Level: advanced

    887:    Concepts: gather^into strided vector

    889: .seealso: VecStrideNorm(), VecStrideScatter(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    890:           VecStrideScatterAll()
    891: @*/
    892: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideGather(Vec v,PetscInt start,Vec s,InsertMode addv)
    893: {
    895:   PetscInt       i,n,bs,ns;
    896:   PetscScalar    *x,*y;

    901:   VecGetLocalSize(v,&n);
    902:   VecGetLocalSize(s,&ns);
    903:   VecGetArray(v,&x);
    904:   VecGetArray(s,&y);

    906:   bs   = v->bs;
    907:   if (start >= bs) {
    908:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
    909:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    910:   }
    911:   if (n != ns*bs) {
    912:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for gather from original vector %D",ns*bs,n);
    913:   }
    914:   x += start;
    915:   n =  n/bs;

    917:   if (addv == INSERT_VALUES) {
    918:     for (i=0; i<n; i++) {
    919:       y[i] = x[bs*i];
    920:     }
    921:   } else if (addv == ADD_VALUES) {
    922:     for (i=0; i<n; i++) {
    923:       y[i] += x[bs*i];
    924:     }
    925: #if !defined(PETSC_USE_COMPLEX)
    926:   } else if (addv == MAX_VALUES) {
    927:     for (i=0; i<n; i++) {
    928:       y[i] = PetscMax(y[i],x[bs*i]);
    929:     }
    930: #endif
    931:   } else {
    932:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
    933:   }

    935:   VecRestoreArray(v,&x);
    936:   VecRestoreArray(s,&y);
    937:   return(0);
    938: }

    942: /*@
    943:    VecStrideScatter - Scatters a single component from a vector into a multi-component vector.

    945:    Collective on Vec

    947:    Input Parameter:
    948: +  s - the single-component vector
    949: .  start - starting point of the subvector (defined by a stride)
    950: -  addv - one of ADD_VALUES,INSERT_VALUES,MAX_VALUES

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

    955:    Notes:
    956:    One must call VecSetBlockSize() on the multi-component vector before this
    957:    routine to set the stride  information, or use a vector created from a multicomponent DA.

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

    962:    Not optimized; could be easily

    964:    Level: advanced

    966:    Concepts: scatter^into strided vector

    968: .seealso: VecStrideNorm(), VecStrideGather(), VecStrideMin(), VecStrideMax(), VecStrideGatherAll(),
    969:           VecStrideScatterAll()
    970: @*/
    971: PetscErrorCode PETSCVEC_DLLEXPORT VecStrideScatter(Vec s,PetscInt start,Vec v,InsertMode addv)
    972: {
    974:   PetscInt       i,n,bs,ns;
    975:   PetscScalar    *x,*y;

    980:   VecGetLocalSize(v,&n);
    981:   VecGetLocalSize(s,&ns);
    982:   VecGetArray(v,&x);
    983:   VecGetArray(s,&y);

    985:   bs   = v->bs;
    986:   if (start >= bs) {
    987:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Start of stride subvector (%D) is too large for stride\n\
    988:             Have you set the vector blocksize (%D) correctly with VecSetBlockSize()?",start,bs);
    989:   }
    990:   if (n != ns*bs) {
    991:     SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Subvector length * blocksize %D not correct for scatter to multicomponent vector %D",ns*bs,n);
    992:   }
    993:   x += start;
    994:   n =  n/bs;

    997:   if (addv == INSERT_VALUES) {
    998:     for (i=0; i<n; i++) {
    999:       x[bs*i] = y[i];
    1000:     }
    1001:   } else if (addv == ADD_VALUES) {
    1002:     for (i=0; i<n; i++) {
    1003:       x[bs*i] += y[i];
    1004:     }
    1005: #if !defined(PETSC_USE_COMPLEX)
    1006:   } else if (addv == MAX_VALUES) {
    1007:     for (i=0; i<n; i++) {
    1008:       x[bs*i] = PetscMax(y[i],x[bs*i]);
    1009:     }
    1010: #endif
    1011:   } else {
    1012:     SETERRQ(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown insert type");
    1013:   }

    1016:   VecRestoreArray(v,&x);
    1017:   VecRestoreArray(s,&y);
    1018:   return(0);
    1019: }

    1023: PetscErrorCode VecReciprocal_Default(Vec v)
    1024: {
    1026:   PetscInt         i,n;
    1027:   PetscScalar *x;

    1030:   VecGetLocalSize(v,&n);
    1031:   VecGetArray(v,&x);
    1032:   for (i=0; i<n; i++) {
    1033:     if (x[i] != 0.0) x[i] = 1.0/x[i];
    1034:   }
    1035:   VecRestoreArray(v,&x);
    1036:   return(0);
    1037: }

    1041: /*@
    1042:   VecSqrt - Replaces each component of a vector by the square root of its magnitude.

    1044:   Not collective

    1046:   Input Parameter:
    1047: . v - The vector

    1049:   Output Parameter:
    1050: . v - The vector square root

    1052:   Level: beginner

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

    1056: .keywords: vector, sqrt, square root
    1057: @*/
    1058: PetscErrorCode PETSCVEC_DLLEXPORT VecSqrt(Vec v)
    1059: {
    1060:   PetscScalar *x;
    1061:   PetscInt         i, n;

    1066:   VecGetLocalSize(v, &n);
    1067:   VecGetArray(v, &x);
    1068:   for(i = 0; i < n; i++) {
    1069:     x[i] = sqrt(PetscAbsScalar(x[i]));
    1070:   }
    1071:   VecRestoreArray(v, &x);
    1072:   return(0);
    1073: }

    1077: /*@
    1078:    VecSum - Computes the sum of all the components of a vector.

    1080:    Collective on Vec

    1082:    Input Parameter:
    1083: .  v - the vector

    1085:    Output Parameter:
    1086: .  sum - the result

    1088:    Level: beginner

    1090:    Concepts: sum^of vector entries

    1092: .seealso: VecNorm()
    1093: @*/
    1094: PetscErrorCode PETSCVEC_DLLEXPORT VecSum(Vec v,PetscScalar *sum)
    1095: {
    1097:   PetscInt         i,n;
    1098:   PetscScalar *x,lsum = 0.0;

    1103:   VecGetLocalSize(v,&n);
    1104:   VecGetArray(v,&x);
    1105:   for (i=0; i<n; i++) {
    1106:     lsum += x[i];
    1107:   }
    1108:   MPI_Allreduce(&lsum,sum,1,MPIU_SCALAR,PetscSum_Op,v->comm);
    1109:   VecRestoreArray(v,&x);
    1110:   return(0);
    1111: }

    1115: /*@
    1116:    VecShift - Shifts all of the components of a vector by computing
    1117:    x[i] = x[i] + shift.

    1119:    Collective on Vec

    1121:    Input Parameters:
    1122: +  v - the vector
    1123: -  shift - the shift

    1125:    Output Parameter:
    1126: .  v - the shifted vector

    1128:    Level: intermediate

    1130:    Concepts: vector^adding constant

    1132: @*/
    1133: PetscErrorCode PETSCVEC_DLLEXPORT VecShift(Vec v,PetscScalar shift)
    1134: {
    1136:   PetscInt i,n;
    1137:   PetscScalar *x;

    1141:   VecGetLocalSize(v,&n);
    1142:   VecGetArray(v,&x);
    1143:   for (i=0; i<n; i++) {
    1144:     x[i] += shift;
    1145:   }
    1146:   VecRestoreArray(v,&x);
    1147:   return(0);
    1148: }

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

    1155:    Collective on Vec

    1157:    Input Parameters:
    1158: .  v - the vector

    1160:    Level: intermediate

    1162:    Concepts: vector^absolute value

    1164: @*/
    1165: PetscErrorCode PETSCVEC_DLLEXPORT VecAbs(Vec v)
    1166: {
    1168:   PetscInt       i,n;
    1169:   PetscScalar    *x;

    1173:   VecGetLocalSize(v,&n);
    1174:   VecGetArray(v,&x);
    1175:   for (i=0; i<n; i++) {
    1176:     x[i] = PetscAbsScalar(x[i]);
    1177:   }
    1178:   VecRestoreArray(v,&x);
    1179:   return(0);
    1180: }

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

    1187:   Input Parameters:
    1188: + vec   - The vector
    1189: . order - The ordering
    1190: - inv   - The flag for inverting the permutation

    1192:   Level: beginner

    1194:   Note: This function does not yet support parallel Index Sets

    1196: .seealso: MatPermute()
    1197: .keywords: vec, permute
    1198: @*/
    1199: PetscErrorCode PETSCVEC_DLLEXPORT VecPermute(Vec x, IS row, PetscTruth inv)
    1200: {
    1201:   PetscScalar    *array, *newArray;
    1202:   PetscInt       *idx;
    1203:   PetscInt       i;

    1207:   ISGetIndices(row, &idx);
    1208:   VecGetArray(x, &array);
    1209:   PetscMalloc(x->n*sizeof(PetscScalar), &newArray);
    1210: #ifdef PETSC_USE_DEBUG
    1211:   for(i = 0; i < x->n; i++) {
    1212:     if ((idx[i] < 0) || (idx[i] >= x->n)) {
    1213:       SETERRQ2(PETSC_ERR_ARG_CORRUPT, "Permutation index %D is out of bounds: %D", i, idx[i]);
    1214:     }
    1215:   }
    1216: #endif
    1217:   if (!inv) {
    1218:     for(i = 0; i < x->n; i++) newArray[i]      = array[idx[i]];
    1219:   } else {
    1220:     for(i = 0; i < x->n; i++) newArray[idx[i]] = array[i];
    1221:   }
    1222:   VecRestoreArray(x, &array);
    1223:   ISRestoreIndices(row, &idx);
    1224:   VecReplaceArray(x, newArray);
    1225:   return(0);
    1226: }

    1230: /*@
    1231:    VecEqual - Compares two vectors.

    1233:    Collective on Vec

    1235:    Input Parameters:
    1236: +  vec1 - the first matrix
    1237: -  vec2 - the second matrix

    1239:    Output Parameter:
    1240: .  flg - PETSC_TRUE if the vectors are equal; PETSC_FALSE otherwise.

    1242:    Level: intermediate

    1244:    Concepts: equal^two vectors
    1245:    Concepts: vector^equality

    1247: @*/
    1248: PetscErrorCode PETSCVEC_DLLEXPORT VecEqual(Vec vec1,Vec vec2,PetscTruth *flg)
    1249: {
    1250:   PetscScalar    *v1,*v2;
    1252:   PetscInt       n1,n2;
    1253:   PetscTruth     flg1;

    1256:   if (vec1 == vec2) {
    1257:     *flg = PETSC_TRUE;
    1258:   } else {
    1259:     VecGetSize(vec1,&n1);
    1260:     VecGetSize(vec2,&n2);
    1261:     if (n1 != n2) {
    1262:       flg1 = PETSC_FALSE;
    1263:     } else {
    1264:       VecGetArray(vec1,&v1);
    1265:       VecGetArray(vec2,&v2);
    1266:       PetscMemcmp(v1,v2,n1*sizeof(PetscScalar),&flg1);
    1267:       VecRestoreArray(vec1,&v1);
    1268:       VecRestoreArray(vec2,&v2);
    1269:     }
    1270:     /* combine results from all processors */
    1271:     MPI_Allreduce(&flg1,flg,1,MPI_INT,MPI_MIN,vec1->comm);
    1272:   }
    1273:   return(0);
    1274: }