Actual source code: vector.c

      1: #define PETSCVEC_DLL

      3: /*
      4:      Provides the interface functions for all vector operations.
      5:    These are the vector functions the user calls.
      6: */
     7:  #include vecimpl.h

      9: /* Logging support */
     10: PetscCookie PETSCVEC_DLLEXPORT VEC_COOKIE = 0;
     11: PetscEvent  VEC_View = 0, VEC_Max = 0, VEC_Min = 0, VEC_DotBarrier = 0, VEC_Dot = 0, VEC_MDotBarrier = 0, VEC_MDot = 0, VEC_TDot = 0;
     12: PetscEvent  VEC_Norm = 0, VEC_Normalize = 0, VEC_Scale = 0, VEC_Copy = 0, VEC_Set = 0, VEC_AXPY = 0, VEC_AYPX = 0, VEC_WAXPY = 0;
     13: PetscEvent  VEC_MTDot = 0, VEC_NormBarrier = 0, VEC_MAXPY = 0, VEC_Swap = 0, VEC_AssemblyBegin = 0, VEC_ScatterBegin = 0, VEC_ScatterEnd = 0;
     14: PetscEvent  VEC_AssemblyEnd = 0, VEC_PointwiseMult = 0, VEC_SetValues = 0, VEC_Load = 0, VEC_ScatterBarrier = 0;
     15: PetscEvent  VEC_SetRandom = 0, VEC_ReduceArithmetic = 0, VEC_ReduceBarrier = 0, VEC_ReduceCommunication = 0;

     17: /* ugly globals for VecSetValue() and VecSetValueLocal() */
     18: PetscInt    PETSCVEC_DLLEXPORT VecSetValue_Row = 0;
     19: PetscScalar PETSCVEC_DLLEXPORT VecSetValue_Value = 0.0;

     21: /*@
     22:   VecZeroEntries - puts a 0.0 in each element of a vector

     24:   Collective on Vec

     26:   Input Parameter:
     27: . vec - The vector

     29:   Level: beginner

     31: .keywords: Vec, set, options, database
     32: .seealso: VecCreate(), VecPrintHelp(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
     33: @*/
     34: PetscErrorCode PETSCVEC_DLLEXPORT VecZeroEntries (Vec vec)
     35: {
     38:   VecSet(vec,0.0);
     39:   return(0);
     40: }

     44: /*
     45:   VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
     46:   processor and a PETSc MPI vector on more than one processor.

     48:   Collective on Vec

     50:   Input Parameter:
     51: . vec - The vector

     53:   Level: intermediate

     55: .keywords: Vec, set, options, database, type
     56: .seealso: VecSetFromOptions(), VecSetType()
     57: */
     58: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec)
     59: {
     60:   PetscTruth     opt;
     61:   const char     *defaultType;
     62:   char           typeName[256];
     63:   PetscMPIInt    size;

     67:   if (vec->type_name) {
     68:     defaultType = vec->type_name;
     69:   } else {
     70:     MPI_Comm_size(vec->comm, &size);
     71:     if (size > 1) {
     72:       defaultType = VECMPI;
     73:     } else {
     74:       defaultType = VECSEQ;
     75:     }
     76:   }

     78:   if (!VecRegisterAllCalled) {
     79:     VecRegisterAll(PETSC_NULL);
     80:   }
     81:   PetscOptionsList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
     82:   if (opt) {
     83:     VecSetType(vec, typeName);
     84:   } else {
     85:     VecSetType(vec, defaultType);
     86:   }
     87:   return(0);
     88: }

     92: /*@C
     93:   VecSetFromOptions - Configures the vector from the options database.

     95:   Collective on Vec

     97:   Input Parameter:
     98: . vec - The vector

    100:   Notes:  To see all options, run your program with the -help option, or consult the users manual.
    101:           Must be called after VecCreate() but before the vector is used.

    103:   Level: beginner

    105:   Concepts: vectors^setting options
    106:   Concepts: vectors^setting type

    108: .keywords: Vec, set, options, database
    109: .seealso: VecCreate(), VecPrintHelp(), VecSetOptionsPrefix()
    110: @*/
    111: PetscErrorCode PETSCVEC_DLLEXPORT VecSetFromOptions(Vec vec)
    112: {
    113:   PetscTruth     opt;

    119:   PetscOptionsBegin(vec->comm, vec->prefix, "Vector options", "Vec");

    121:   /* Handle generic vector options */
    122:   PetscOptionsHasName(PETSC_NULL, "-help", &opt);
    123:   if (opt) {
    124:     VecPrintHelp(vec);
    125:   }

    127:   /* Handle vector type options */
    128:   VecSetTypeFromOptions_Private(vec);

    130:   /* Handle specific vector options */
    131:   if (vec->ops->setfromoptions) {
    132:     (*vec->ops->setfromoptions)(vec);
    133:   }
    134:   PetscOptionsEnd();

    136:   VecViewFromOptions(vec, vec->name);
    137:   return(0);
    138: }

    142: /*@
    143:   VecPrintHelp - Prints some options for the Vec.

    145:   Input Parameter:
    146: . vec - The vector

    148:   Options Database Keys:
    149: $  -help, -h

    151:   Level: intermediate

    153: .keywords: Vec, help
    154: .seealso: VecSetFromOptions()
    155: @*/
    156: PetscErrorCode PETSCVEC_DLLEXPORT VecPrintHelp(Vec vec)
    157: {
    160:   return(0);
    161: }

    165: /*@
    166:   VecSetSizes - Sets the local and global sizes, and checks to determine compatibility

    168:   Collective on Vec

    170:   Input Parameters:
    171: + v - the vector
    172: . n - the local size (or PETSC_DECIDE to have it set)
    173: - N - the global size (or PETSC_DECIDE)

    175:   Notes:
    176:   n and N cannot be both PETSC_DECIDE
    177:   If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.

    179:   Level: intermediate

    181: .seealso: VecGetSize(), PetscSplitOwnership()
    182: @*/
    183: PetscErrorCode PETSCVEC_DLLEXPORT VecSetSizes(Vec v, PetscInt n, PetscInt N)
    184: {
    187:   if (N > 0 && n > N) SETERRQ2(PETSC_ERR_ARG_INCOMP,"Local size %D cannot be larger than global size %D",n,N);
    188:   if ((v->n >= 0 || v->N >= 0) && (v->n != n || v->N != N)) SETERRQ4(PETSC_ERR_SUP,"Cannot change/reset vector sizes to %D local %D global after previously setting them to %D local %D global",n,N,v->n,v->N);
    189:   v->n = n;
    190:   v->N = N;
    191:   return(0);
    192: }

    196: /*@
    197:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
    198:    and VecSetValuesBlockedLocal().

    200:    Collective on Vec

    202:    Input Parameter:
    203: +  v - the vector
    204: -  bs - the blocksize

    206:    Notes:
    207:    All vectors obtained by VecDuplicate() inherit the same blocksize.

    209:    Level: advanced

    211: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlocked(), VecGetBlockSize()

    213:   Concepts: block size^vectors
    214: @*/
    215: PetscErrorCode PETSCVEC_DLLEXPORT VecSetBlockSize(Vec v,PetscInt bs)
    216: {
    219:   if (bs <= 0) bs = 1;
    220:   if (bs == v->bs) return(0);
    221:   if (v->bs != -1) SETERRQ2(PETSC_ERR_ARG_WRONGSTATE,"Cannot reset blocksize. Current size %D new %D",v->bs,bs);
    222:   if (v->N != -1 && v->N % bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Vector length not divisible by blocksize %D %D",v->N,bs);
    223:   if (v->n != -1 && v->n % bs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local vector length not divisible by blocksize %D %D\n\
    224:    Try setting blocksize before setting the vector type",v->n,bs);
    225:
    226:   v->bs        = bs;
    227:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
    228:   return(0);
    229: }

    233: /*@
    234:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
    235:    and VecSetValuesBlockedLocal().

    237:    Collective on Vec

    239:    Input Parameter:
    240: .  v - the vector

    242:    Output Parameter:
    243: .  bs - the blocksize

    245:    Notes:
    246:    All vectors obtained by VecDuplicate() inherit the same blocksize.

    248:    Level: advanced

    250: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlocked(), VecSetBlockSize()

    252:    Concepts: vector^block size
    253:    Concepts: block^vector

    255: @*/
    256: PetscErrorCode PETSCVEC_DLLEXPORT VecGetBlockSize(Vec v,PetscInt *bs)
    257: {
    261:   *bs = v->bs;
    262:   return(0);
    263: }

    267: /*@
    268:    VecValid - Checks whether a vector object is valid.

    270:    Not Collective

    272:    Input Parameter:
    273: .  v - the object to check

    275:    Output Parameter:
    276:    flg - flag indicating vector status, either
    277:    PETSC_TRUE if vector is valid, or PETSC_FALSE otherwise.

    279:    Level: developer

    281: @*/
    282: PetscErrorCode PETSCVEC_DLLEXPORT VecValid(Vec v,PetscTruth *flg)
    283: {
    286:   if (!v)                           *flg = PETSC_FALSE;
    287:   else if (v->cookie != VEC_COOKIE) *flg = PETSC_FALSE;
    288:   else                              *flg = PETSC_TRUE;
    289:   return(0);
    290: }

    294: /*@C
    295:    VecSetOptionsPrefix - Sets the prefix used for searching for all
    296:    Vec options in the database.

    298:    Collective on Vec

    300:    Input Parameter:
    301: +  v - the Vec context
    302: -  prefix - the prefix to prepend to all option names

    304:    Notes:
    305:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    306:    The first character of all runtime options is AUTOVECICALLY the hyphen.

    308:    Level: advanced

    310: .keywords: Vec, set, options, prefix, database

    312: .seealso: VecSetFromOptions()
    313: @*/
    314: PetscErrorCode PETSCVEC_DLLEXPORT VecSetOptionsPrefix(Vec v,const char prefix[])
    315: {

    320:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
    321:   return(0);
    322: }

    326: /*@C
    327:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
    328:    Vec options in the database.

    330:    Collective on Vec

    332:    Input Parameters:
    333: +  v - the Vec context
    334: -  prefix - the prefix to prepend to all option names

    336:    Notes:
    337:    A hyphen (-) must NOT be given at the beginning of the prefix name.
    338:    The first character of all runtime options is AUTOVECICALLY the hyphen.

    340:    Level: advanced

    342: .keywords: Vec, append, options, prefix, database

    344: .seealso: VecGetOptionsPrefix()
    345: @*/
    346: PetscErrorCode PETSCVEC_DLLEXPORT VecAppendOptionsPrefix(Vec v,const char prefix[])
    347: {
    349:
    352:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
    353:   return(0);
    354: }

    358: /*@C
    359:    VecGetOptionsPrefix - Sets the prefix used for searching for all
    360:    Vec options in the database.

    362:    Not Collective

    364:    Input Parameter:
    365: .  A - the Vec context

    367:    Output Parameter:
    368: .  prefix - pointer to the prefix string used

    370:    Notes: On the fortran side, the user should pass in a string 'prefix' of
    371:    sufficient length to hold the prefix.

    373:    Level: advanced

    375: .keywords: Vec, get, options, prefix, database

    377: .seealso: VecAppendOptionsPrefix()
    378: @*/
    379: PetscErrorCode PETSCVEC_DLLEXPORT VecGetOptionsPrefix(Vec v,const char *prefix[])
    380: {

    385:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
    386:   return(0);
    387: }

    391: /*@
    392:    VecSetUp - Sets up the internal vector data structures for the later use.

    394:    Collective on Vec

    396:    Input Parameters:
    397: .  v - the Vec context

    399:    Notes:
    400:    For basic use of the Vec classes the user need not explicitly call
    401:    VecSetUp(), since these actions will happen automatically.

    403:    Level: advanced

    405: .keywords: Vec, setup

    407: .seealso: VecCreate(), VecDestroy()
    408: @*/
    409: PetscErrorCode PETSCVEC_DLLEXPORT VecSetUp(Vec v)
    410: {

    415:   VecSetFromOptions(v);
    416:   return(0);
    417: }

    421: /*@
    422:    VecDot - Computes the vector dot product.

    424:    Collective on Vec

    426:    Input Parameters:
    427: .  x, y - the vectors

    429:    Output Parameter:
    430: .  alpha - the dot product

    432:    Performance Issues:
    433: +    per-processor memory bandwidth
    434: .    interprocessor latency
    435: -    work load inbalance that causes certain processes to arrive much earlier than
    436:      others

    438:    Notes for Users of Complex Numbers:
    439:    For complex vectors, VecDot() computes
    440: $     val = (x,y) = y^H x,
    441:    where y^H denotes the conjugate transpose of y.

    443:    Use VecTDot() for the indefinite form
    444: $     val = (x,y) = y^T x,
    445:    where y^T denotes the transpose of y.

    447:    Level: intermediate

    449:    Concepts: inner product
    450:    Concepts: vector^inner product

    452: .seealso: VecMDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd()
    453: @*/
    454: PetscErrorCode PETSCVEC_DLLEXPORT VecDot(Vec x,Vec y,PetscScalar *val)
    455: {

    465:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    466:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    468:   PetscLogEventBarrierBegin(VEC_DotBarrier,x,y,0,0,x->comm);
    469:   (*x->ops->dot)(x,y,val);
    470:   PetscLogEventBarrierEnd(VEC_DotBarrier,x,y,0,0,x->comm);
    471:   return(0);
    472: }

    476: /*@
    477:    VecNorm  - Computes the vector norm.

    479:    Collective on Vec

    481:    Input Parameters:
    482: +  x - the vector
    483: -  type - one of NORM_1, NORM_2, NORM_INFINITY.  Also available
    484:           NORM_1_AND_2, which computes both norms and stores them
    485:           in a two element array.

    487:    Output Parameter:
    488: .  val - the norm

    490:    Notes:
    491: $     NORM_1 denotes sum_i |x_i|
    492: $     NORM_2 denotes sqrt(sum_i (x_i)^2)
    493: $     NORM_INFINITY denotes max_i |x_i|

    495:    Level: intermediate

    497:    Performance Issues:
    498: +    per-processor memory bandwidth
    499: .    interprocessor latency
    500: -    work load inbalance that causes certain processes to arrive much earlier than
    501:      others

    503:    Compile Option:
    504:    PETSC_HAVE_SLOW_BLAS_NORM2 will cause a C (loop unrolled) version of the norm to be used, rather
    505:  than the BLAS. This should probably only be used when one is using the FORTRAN BLAS routines
    506:  (as opposed to vendor provided) because the FORTRAN BLAS NRM2() routine is very slow.

    508:    Concepts: norm
    509:    Concepts: vector^norm

    511: .seealso: VecDot(), VecTDot(), VecNorm(), VecDotBegin(), VecDotEnd(),
    512:           VecNormBegin(), VecNormEnd()

    514: @*/
    515: PetscErrorCode PETSCVEC_DLLEXPORT VecNorm(Vec x,NormType type,PetscReal *val)
    516: {
    517:   PetscTruth     flg;
    519:   PetscInt       type_id;

    526:   /*
    527:    * Cached data?
    528:    */
    529:   if (type!=NORM_1_AND_2) {
    530:     VecNormComposedDataID(type,&type_id);
    531:     PetscObjectComposedDataGetReal((PetscObject)x,type_id,*val,flg);
    532:     if (flg) return(0);
    533:   }
    534:

    536:   PetscLogEventBarrierBegin(VEC_NormBarrier,x,0,0,0,x->comm);
    537:   (*x->ops->norm)(x,type,val);
    538:   PetscLogEventBarrierEnd(VEC_NormBarrier,x,0,0,0,x->comm);

    540:   if (type!=NORM_1_AND_2) {
    541:     PetscObjectComposedDataSetReal((PetscObject)x,type_id,*val);
    542:   }

    544:   return(0);
    545: }

    549: PetscErrorCode PETSCVEC_DLLEXPORT VecNormComposedDataID(NormType type,PetscInt *type_id)
    550: {
    551:   static PetscInt id_norm1=-1,id_norm2=-1,id_normInf=-1,id_normF=-1,id_norm12=-1;
    552:   PetscErrorCode  ierr;

    555:   switch (type) {
    556:   case NORM_1 :
    557:     if (id_norm1==-1) {
    558:       PetscObjectComposedDataRegister(&id_norm1);}
    559:     *type_id = id_norm1; break;
    560:   case NORM_2 :
    561:     if (id_norm2==-1) {
    562:       PetscObjectComposedDataRegister(&id_norm2);}
    563:     *type_id = id_norm2; break;
    564:   case NORM_1_AND_2 :
    565:     /* we don't handle this one yet */
    566:     if (id_norm1==-1) {
    567:       PetscObjectComposedDataRegister(&id_norm1);}
    568:     if (id_norm2==-1) {
    569:       PetscObjectComposedDataRegister(&id_norm2);}
    570:     *type_id = id_norm12; break;
    571:   case NORM_INFINITY :
    572:     if (id_normInf==-1) {
    573:       PetscObjectComposedDataRegister(&id_normInf);}
    574:     *type_id = id_normInf; break;
    575:   case NORM_FROBENIUS :
    576:     if (id_normF==-1) {
    577:       PetscObjectComposedDataRegister(&id_normF);}
    578:     *type_id = id_normF; break;
    579:   }
    580:   return(0);
    581: }

    585: /*@
    586:    VecNormalize - Normalizes a vector by 2-norm.

    588:    Collective on Vec

    590:    Input Parameters:
    591: +  x - the vector

    593:    Output Parameter:
    594: .  x - the normalized vector
    595: -  val - the vector norm before normalization

    597:    Level: intermediate

    599:    Concepts: vector^normalizing
    600:    Concepts: normalizing^vector

    602: @*/
    603: PetscErrorCode PETSCVEC_DLLEXPORT VecNormalize(Vec x,PetscReal *val)
    604: {

    611:   PetscLogEventBegin(VEC_Normalize,x,0,0,0);
    612:   VecNorm(x,NORM_2,val);
    613:   if (!*val) {
    614:     PetscLogInfo((x,"VecNormalize:Vector of zero norm can not be normalized; Returning only the zero norm\n"));
    615:   } else {
    616:     PetscScalar tmp = 1.0/(*val);
    617:     VecScale(x,tmp);
    618:   }

    620:   PetscLogEventEnd(VEC_Normalize,x,0,0,0);
    621:   return(0);
    622: }

    626: /*@C
    627:    VecMax - Determines the maximum vector component and its location.

    629:    Collective on Vec

    631:    Input Parameter:
    632: .  x - the vector

    634:    Output Parameters:
    635: +  val - the maximum component
    636: -  p - the location of val

    638:    Notes:
    639:    Returns the value PETSC_MIN and p = -1 if the vector is of length 0.

    641:    Level: intermediate

    643:    Concepts: maximum^of vector
    644:    Concepts: vector^maximum value

    646: .seealso: VecNorm(), VecMin()
    647: @*/
    648: PetscErrorCode PETSCVEC_DLLEXPORT VecMax(Vec x,PetscInt *p,PetscReal *val)
    649: {

    656:   PetscLogEventBegin(VEC_Max,x,0,0,0);
    657:   (*x->ops->max)(x,p,val);
    658:   PetscLogEventEnd(VEC_Max,x,0,0,0);
    659:   return(0);
    660: }

    664: /*@
    665:    VecMin - Determines the minimum vector component and its location.

    667:    Collective on Vec

    669:    Input Parameters:
    670: .  x - the vector

    672:    Output Parameter:
    673: +  val - the minimum component
    674: -  p - the location of val

    676:    Level: intermediate

    678:    Notes:
    679:    Returns the value PETSC_MAX and p = -1 if the vector is of length 0.

    681:    Concepts: minimum^of vector
    682:    Concepts: vector^minimum entry

    684: .seealso: VecMax()
    685: @*/
    686: PetscErrorCode PETSCVEC_DLLEXPORT VecMin(Vec x,PetscInt *p,PetscReal *val)
    687: {

    694:   PetscLogEventBegin(VEC_Min,x,0,0,0);
    695:   (*x->ops->min)(x,p,val);
    696:   PetscLogEventEnd(VEC_Min,x,0,0,0);
    697:   return(0);
    698: }

    702: /*@
    703:    VecTDot - Computes an indefinite vector dot product. That is, this
    704:    routine does NOT use the complex conjugate.

    706:    Collective on Vec

    708:    Input Parameters:
    709: .  x, y - the vectors

    711:    Output Parameter:
    712: .  val - the dot product

    714:    Notes for Users of Complex Numbers:
    715:    For complex vectors, VecTDot() computes the indefinite form
    716: $     val = (x,y) = y^T x,
    717:    where y^T denotes the transpose of y.

    719:    Use VecDot() for the inner product
    720: $     val = (x,y) = y^H x,
    721:    where y^H denotes the conjugate transpose of y.

    723:    Level: intermediate

    725:    Concepts: inner product^non-Hermitian
    726:    Concepts: vector^inner product
    727:    Concepts: non-Hermitian inner product

    729: .seealso: VecDot(), VecMTDot()
    730: @*/
    731: PetscErrorCode PETSCVEC_DLLEXPORT VecTDot(Vec x,Vec y,PetscScalar *val)
    732: {

    742:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    743:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    745:   PetscLogEventBegin(VEC_TDot,x,y,0,0);
    746:   (*x->ops->tdot)(x,y,val);
    747:   PetscLogEventEnd(VEC_TDot,x,y,0,0);
    748:   return(0);
    749: }

    753: /*@
    754:    VecScale - Scales a vector.

    756:    Collective on Vec

    758:    Input Parameters:
    759: +  x - the vector
    760: -  alpha - the scalar

    762:    Output Parameter:
    763: .  x - the scaled vector

    765:    Note:
    766:    For a vector with n components, VecScale() computes
    767: $      x[i] = alpha * x[i], for i=1,...,n.

    769:    Level: intermediate

    771:    Concepts: vector^scaling
    772:    Concepts: scaling^vector

    774: @*/
    775: PetscErrorCode PETSCVEC_DLLEXPORT VecScale (Vec x, PetscScalar alpha)
    776: {
    777:   PetscReal      scale,norm1=0.0,norm2=0.0,normInf=0.0,normF=0.0;
    778:   PetscTruth     flg1,flg2,flgInf,flgF;
    780:   PetscInt       type_id1,type_id2,type_idInf,type_idF;

    785:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
    786:   PetscLogEventBegin(VEC_Scale,x,0,0,0);
    787:   (*x->ops->scale)(x,alpha);

    789:   /*
    790:    * Update cached data
    791:    */
    792:   /* see if we have cached norms */
    793:   /* 1 */
    794:   VecNormComposedDataID(NORM_1,&type_id1);
    795:   PetscObjectComposedDataGetReal((PetscObject)x,type_id1,norm1,flg1);
    796:   /* 2 */
    797:   VecNormComposedDataID(NORM_2,&type_id2);
    798:   PetscObjectComposedDataGetReal((PetscObject)x,type_id2,norm2,flg2);
    799:   /* inf */
    800:   VecNormComposedDataID(NORM_INFINITY,&type_idInf);
    801:   PetscObjectComposedDataGetReal((PetscObject)x,type_idInf,normInf,flgInf);
    802:   /* frobenius */
    803:   VecNormComposedDataID(NORM_FROBENIUS,&type_idF);
    804:   PetscObjectComposedDataGetReal((PetscObject)x,type_idF,normF,flgF);

    806:   /* in general we consider this object touched */
    807:   PetscObjectStateIncrease((PetscObject)x);

    809:   /* however, norms can be simply updated */
    810:   scale = PetscAbsScalar(alpha);
    811:   /* 1 */
    812:   if (flg1) {
    813:     PetscObjectComposedDataSetReal((PetscObject)x,type_id1,scale*norm1);
    814:   }
    815:   /* 2 */
    816:   if (flg2) {
    817:     PetscObjectComposedDataSetReal((PetscObject)x,type_id2,scale*norm2);
    818:   }
    819:   /* inf */
    820:   if (flgInf) {
    821:     PetscObjectComposedDataSetReal((PetscObject)x,type_idInf,scale*normInf);
    822:   }
    823:   /* frobenius */
    824:   if (flgF) {
    825:     PetscObjectComposedDataSetReal((PetscObject)x,type_idF,scale*normF);
    826:   }

    828:   PetscLogEventEnd(VEC_Scale,x,0,0,0);
    829:   return(0);
    830: }

    834: /*@
    835:    VecCopy - Copies a vector.

    837:    Collective on Vec

    839:    Input Parameter:
    840: .  x - the vector

    842:    Output Parameter:
    843: .  y - the copy

    845:    Notes:
    846:    For default parallel PETSc vectors, both x and y must be distributed in
    847:    the same manner; local copies are done.

    849:    Level: beginner

    851: .seealso: VecDuplicate()
    852: @*/
    853: PetscErrorCode PETSCVEC_DLLEXPORT VecCopy(Vec x,Vec y)
    854: {
    855:   PetscTruth     flg;
    856:   PetscReal      norm=0.0;
    858:   PetscInt       type_id;

    866:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    867:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    869:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
    870:   (*x->ops->copy)(x,y);

    872:   /*
    873:    * Update cached data
    874:    */
    875:   /* in general we consider this object touched */
    876:   PetscObjectStateIncrease((PetscObject)y);
    877:   /* however, norms can be simply copied over */
    878:   /* 2 */
    879:   VecNormComposedDataID(NORM_2,&type_id);
    880:   PetscObjectComposedDataGetReal((PetscObject)x,type_id,norm,flg);
    881:   if (flg) {
    882:     PetscObjectComposedDataSetReal((PetscObject)y,type_id,norm);
    883:   }
    884:   /* 1 */
    885:   VecNormComposedDataID(NORM_1,&type_id);
    886:   PetscObjectComposedDataGetReal((PetscObject)x,type_id,norm,flg);
    887:   if (flg) {
    888:     PetscObjectComposedDataSetReal((PetscObject)y,type_id,norm);
    889:   }
    890:   /* inf */
    891:   VecNormComposedDataID(NORM_INFINITY,&type_id);
    892:   PetscObjectComposedDataGetReal((PetscObject)x,type_id,norm,flg);
    893:   if (flg) {
    894:     PetscObjectComposedDataSetReal((PetscObject)y,type_id,norm);
    895:   }
    896:   /* frobenius */
    897:   VecNormComposedDataID(NORM_FROBENIUS,&type_id);
    898:   PetscObjectComposedDataGetReal((PetscObject)x,type_id,norm,flg);
    899:   if (flg) {
    900:     PetscObjectComposedDataSetReal((PetscObject)y,type_id,norm);
    901:   }

    903:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
    904:   return(0);
    905: }

    909: /*@
    910:    VecSet - Sets all components of a vector to a single scalar value.

    912:    Collective on Vec

    914:    Input Parameters:
    915: +  x  - the vector
    916: -  alpha - the scalar

    918:    Output Parameter:
    919: .  x  - the vector

    921:    Note:
    922:    For a vector of dimension n, VecSet() computes
    923: $     x[i] = alpha, for i=1,...,n,
    924:    so that all vector entries then equal the identical
    925:    scalar value, alpha.  Use the more general routine
    926:    VecSetValues() to set different vector entries.

    928:    You CANNOT call this after you have called VecSetValues() but before you call
    929:    VecAssemblyBegin/End().

    931:    Level: beginner

    933: .seealso VecSetValues(), VecSetValuesBlocked(), VecSetRandom()

    935:    Concepts: vector^setting to constant

    937: @*/
    938: PetscErrorCode PETSCVEC_DLLEXPORT VecSet(Vec x,PetscScalar alpha)
    939: {
    940:   PetscReal      val;
    942:   PetscInt       type_id;

    947:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"You cannot call this after you have called VecSetValues() but\n before you have called VecAssemblyBegin/End()");
    948: #if defined (PETSC_USE_DEBUG)
    949:  {
    950:    PetscReal alpha_local,alpha_max;
    951:    alpha_local = PetscAbsScalar(alpha);
    952:    MPI_Allreduce(&alpha_local,&alpha_max,1,MPIU_REAL,MPI_MAX,x->comm);
    953:    if (alpha_local != alpha_max) SETERRQ(PETSC_ERR_ARG_WRONG,"Same value should be used across all processors");
    954:  }
    955: #endif

    957:   PetscLogEventBegin(VEC_Set,x,0,0,0);
    958:   (*x->ops->set)(x,alpha);
    959:   PetscLogEventEnd(VEC_Set,x,0,0,0);

    961:   /*
    962:    * Update cached data
    963:    */
    964:   /* in general we consider this object touched */
    965:   PetscObjectStateIncrease((PetscObject)x);
    966:   /* however, norms can be simply set */
    967:   /* 1 */
    968:   val = PetscAbsScalar(alpha);
    969:   VecNormComposedDataID(NORM_1,&type_id);
    970:   PetscObjectComposedDataSetReal((PetscObject)x,type_id,x->N * val);
    971:   /* inf */
    972:   VecNormComposedDataID(NORM_INFINITY,&type_id);
    973:   PetscObjectComposedDataSetReal((PetscObject)x,type_id,val);
    974:   /* 2 */
    975:   val = sqrt((double)x->N) * val;
    976:   VecNormComposedDataID(NORM_2,&type_id);
    977:   PetscObjectComposedDataSetReal((PetscObject)x,type_id,val);
    978:   /* frobenius */
    979:   VecNormComposedDataID(NORM_FROBENIUS,&type_id);
    980:   PetscObjectComposedDataSetReal((PetscObject)x,type_id,val);

    982:   return(0);
    983: }

    987: /*@C
    988:    VecSetRandom - Sets all components of a vector to random numbers.

    990:    Collective on Vec

    992:    Input Parameters:
    993: +  x  - the vector
    994: -  rctx - the random number context, formed by PetscRandomCreate(), or PETSC_NULL and
    995:           it will create one internally.

    997:    Output Parameter:
    998: .  x  - the vector

    1000:    Example of Usage:
    1001: .vb
    1002:      PetscRandomCreate(PETSC_COMM_WORLD,RANDOM_DEFAULT,&rctx);
    1003:      VecSetRandom(x,rctx);
    1004:      PetscRandomDestroy(rctx);
    1005: .ve

    1007:    Level: intermediate

    1009:    Concepts: vector^setting to random
    1010:    Concepts: random^vector

    1012: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
    1013: @*/
    1014: PetscErrorCode PETSCVEC_DLLEXPORT VecSetRandom(Vec x,PetscRandom rctx)
    1015: {
    1017:   PetscRandom    randObj = PETSC_NULL;

    1023:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");

    1025:   if (!rctx) {
    1026:     MPI_Comm    comm;
    1027:     PetscObjectGetComm((PetscObject)x,&comm);
    1028:     PetscRandomCreate(comm,RANDOM_DEFAULT,&randObj);
    1029:     rctx = randObj;
    1030:   }

    1032:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
    1033:   (*x->ops->setrandom)(x,rctx);
    1034:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
    1035:
    1036:   if (randObj) {
    1037:     PetscRandomDestroy(randObj);
    1038:   }
    1039:   PetscObjectStateIncrease((PetscObject)x);
    1040:   return(0);
    1041: }

    1045: /*@
    1046:    VecAXPY - Computes y = alpha x + y.

    1048:    Collective on Vec

    1050:    Input Parameters:
    1051: +  alpha - the scalar
    1052: -  x, y  - the vectors

    1054:    Output Parameter:
    1055: .  y - output vector

    1057:    Level: intermediate

    1059:    Concepts: vector^BLAS
    1060:    Concepts: BLAS

    1062: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY()
    1063: @*/
    1064: PetscErrorCode PETSCVEC_DLLEXPORT VecAXPY(Vec y,PetscScalar alpha,Vec x)
    1065: {

    1074:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1075:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1077:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
    1078:   (*y->ops->axpy)(y,alpha,x);
    1079:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
    1080:   PetscObjectStateIncrease((PetscObject)y);
    1081:   return(0);
    1082: }

    1086: /*@
    1087:    VecAXPBY - Computes y = alpha x + beta y.

    1089:    Collective on Vec

    1091:    Input Parameters:
    1092: +  alpha,beta - the scalars
    1093: -  x, y  - the vectors

    1095:    Output Parameter:
    1096: .  y - output vector

    1098:    Level: intermediate

    1100:    Concepts: BLAS
    1101:    Concepts: vector^BLAS

    1103: .seealso: VecAYPX(), VecMAXPY(), VecWAXPY(), VecAXPY()
    1104: @*/
    1105: PetscErrorCode PETSCVEC_DLLEXPORT VecAXPBY(Vec y,PetscScalar alpha,PetscScalar beta,Vec x)
    1106: {

    1115:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1116:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1118:   PetscLogEventBegin(VEC_AXPY,x,y,0,0);
    1119:   (*y->ops->axpby)(y,alpha,beta,x);
    1120:   PetscLogEventEnd(VEC_AXPY,x,y,0,0);
    1121:   PetscObjectStateIncrease((PetscObject)y);
    1122:   return(0);
    1123: }

    1127: /*@
    1128:    VecAYPX - Computes y = x + alpha y.

    1130:    Collective on Vec

    1132:    Input Parameters:
    1133: +  alpha - the scalar
    1134: -  x, y  - the vectors

    1136:    Output Parameter:
    1137: .  y - output vector

    1139:    Level: intermediate

    1141:    Concepts: vector^BLAS
    1142:    Concepts: BLAS

    1144: .seealso: VecAXPY(), VecWAXPY()
    1145: @*/
    1146: PetscErrorCode PETSCVEC_DLLEXPORT VecAYPX(Vec y,PetscScalar alpha,Vec x)
    1147: {

    1156:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1157:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1159:   PetscLogEventBegin(VEC_AYPX,x,y,0,0);
    1160:    (*y->ops->aypx)(y,alpha,x);
    1161:   PetscLogEventEnd(VEC_AYPX,x,y,0,0);
    1162:   PetscObjectStateIncrease((PetscObject)y);
    1163:   return(0);
    1164: }

    1168: /*@
    1169:    VecSwap - Swaps the vectors x and y.

    1171:    Collective on Vec

    1173:    Input Parameters:
    1174: .  x, y  - the vectors

    1176:    Level: advanced

    1178:    Concepts: vector^swapping values

    1180: @*/
    1181: PetscErrorCode PETSCVEC_DLLEXPORT VecSwap(Vec x,Vec y)
    1182: {
    1183:   PetscReal      norm1x=0.0,norm2x=0.0,normInfx=0.0,normFx=0.0;
    1184:   PetscReal      norm1y=0.0,norm2y=0.0,normInfy=0.0,normFy=0.0;
    1185:   PetscTruth     flg1x,flg2x,flgInfx,flgFx;
    1186:   PetscTruth     flg1y,flg2y,flgInfy,flgFy;
    1187:   PetscInt       type_id1,type_id2,type_idInf,type_idF;

    1196:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
    1197:   if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
    1198:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1199:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1201:   PetscLogEventBegin(VEC_Swap,x,y,0,0);

    1203:   /* See if we have cached norms */
    1204:   /* 1 */
    1205:   VecNormComposedDataID(NORM_1,&type_id1);
    1206:   PetscObjectComposedDataGetReal((PetscObject)x,type_id1,norm1x,flg1x);
    1207:   PetscObjectComposedDataGetReal((PetscObject)y,type_id1,norm1y,flg1y);
    1208:   /* 2 */
    1209:   VecNormComposedDataID(NORM_2,&type_id2);
    1210:   PetscObjectComposedDataGetReal((PetscObject)x,type_id2,norm2x,flg2x);
    1211:   PetscObjectComposedDataGetReal((PetscObject)y,type_id2,norm2y,flg2y);
    1212:   /* inf */
    1213:   VecNormComposedDataID(NORM_INFINITY,&type_idInf);
    1214:   PetscObjectComposedDataGetReal((PetscObject)x,type_idInf,normInfx,flgInfx);
    1215:   PetscObjectComposedDataGetReal((PetscObject)y,type_idInf,normInfy,flgInfy);
    1216:   /* frobenius */
    1217:   VecNormComposedDataID(NORM_FROBENIUS,&type_idF);
    1218:   PetscObjectComposedDataGetReal((PetscObject)x,type_idF,normFx,flgFx);
    1219:   PetscObjectComposedDataGetReal((PetscObject)y,type_idF,normFy,flgFy);

    1221:   /* Do the actual swap */
    1222:   (*x->ops->swap)(x,y);
    1223:   PetscObjectStateIncrease((PetscObject)x);
    1224:   PetscObjectStateIncrease((PetscObject)y);

    1226:   /* Swap any cached norms */
    1227:   /* 1 */
    1228:   if (flg1x) {
    1229:     PetscObjectComposedDataSetReal((PetscObject)y,type_id1,norm1x);
    1230:   }
    1231:   if (flg1y) {
    1232:     PetscObjectComposedDataSetReal((PetscObject)x,type_id1,norm1y);
    1233:   }
    1234:   /* 2 */
    1235:   if (flg2x) {
    1236:     PetscObjectComposedDataSetReal((PetscObject)y,type_id2,norm2x);
    1237:   }
    1238:   if (flg2y) {
    1239:     PetscObjectComposedDataSetReal((PetscObject)x,type_id2,norm2y);
    1240:   }
    1241:   /* inf */
    1242:   if (flgInfx) {
    1243:     PetscObjectComposedDataSetReal((PetscObject)y,type_idInf,normInfx);
    1244:   }
    1245:   if (flgInfy) {
    1246:     PetscObjectComposedDataSetReal((PetscObject)x,type_idInf,normInfy);
    1247:   }
    1248:   /* frobenius */
    1249:   if (flgFx) {
    1250:     PetscObjectComposedDataSetReal((PetscObject)y,type_idF,normFx);
    1251:   }
    1252:   if (flgFy) {
    1253:     PetscObjectComposedDataSetReal((PetscObject)x,type_idF,normFy);
    1254:   }

    1256:   PetscLogEventEnd(VEC_Swap,x,y,0,0);

    1258:   return(0);
    1259: }

    1263: /*@
    1264:    VecWAXPY - Computes w = alpha x + y.

    1266:    Collective on Vec

    1268:    Input Parameters:
    1269: +  alpha - the scalar
    1270: -  x, y  - the vectors

    1272:    Output Parameter:
    1273: .  w - the result

    1275:    Level: intermediate

    1277:    Concepts: vector^BLAS
    1278:    Concepts: BLAS

    1280: .seealso: VecAXPY(), VecAYPX()
    1281: @*/
    1282: PetscErrorCode PETSCVEC_DLLEXPORT VecWAXPY(Vec w,PetscScalar alpha,Vec x,Vec y)
    1283: {

    1295:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1296:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1298:   PetscLogEventBegin(VEC_WAXPY,x,y,w,0);
    1299:    (*w->ops->waxpy)(w,alpha,x,y);
    1300:   PetscLogEventEnd(VEC_WAXPY,x,y,w,0);
    1301:   PetscObjectStateIncrease((PetscObject)w);
    1302:   return(0);
    1303: }

    1307: /*@
    1308:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

    1310:    Collective on Vec

    1312:    Input Parameters:
    1313: .  x, y  - the vectors

    1315:    Output Parameter:
    1316: .  w - the result

    1318:    Level: advanced

    1320:    Notes: any subset of the x, y, and w may be the same vector.

    1322:    Concepts: vector^pointwise multiply

    1324: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
    1325: @*/
    1326: PetscErrorCode PETSCVEC_DLLEXPORT VecPointwiseMult(Vec w, Vec x,Vec y)
    1327: {

    1339:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1340:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1342:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
    1343:   (*w->ops->pointwisemult)(w,x,y);
    1344:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
    1345:   PetscObjectStateIncrease((PetscObject)w);
    1346:   return(0);
    1347: }

    1351: /*@
    1352:    VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).

    1354:    Collective on Vec

    1356:    Input Parameters:
    1357: .  x, y  - the vectors

    1359:    Output Parameter:
    1360: .  w - the result

    1362:    Level: advanced

    1364:    Notes: any subset of the x, y, and w may be the same vector.
    1365:           For complex numbers compares only the real part

    1367:    Concepts: vector^pointwise multiply

    1369: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
    1370: @*/
    1371: PetscErrorCode PETSCVEC_DLLEXPORT VecPointwiseMax(Vec w,Vec x,Vec y)
    1372: {

    1384:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1385:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1387:   (*w->ops->pointwisemax)(w,x,y);
    1388:   PetscObjectStateIncrease((PetscObject)w);
    1389:   return(0);
    1390: }

    1395: /*@
    1396:    VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).

    1398:    Collective on Vec

    1400:    Input Parameters:
    1401: .  x, y  - the vectors

    1403:    Output Parameter:
    1404: .  w - the result

    1406:    Level: advanced

    1408:    Notes: any subset of the x, y, and w may be the same vector.
    1409:           For complex numbers compares only the real part

    1411:    Concepts: vector^pointwise multiply

    1413: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
    1414: @*/
    1415: PetscErrorCode PETSCVEC_DLLEXPORT VecPointwiseMin(Vec w,Vec x,Vec y)
    1416: {

    1428:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1429:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1431:   (*w->ops->pointwisemin)(w,x,y);
    1432:   PetscObjectStateIncrease((PetscObject)w);
    1433:   return(0);
    1434: }

    1438: /*@
    1439:    VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).

    1441:    Collective on Vec

    1443:    Input Parameters:
    1444: .  x, y  - the vectors

    1446:    Output Parameter:
    1447: .  w - the result

    1449:    Level: advanced

    1451:    Notes: any subset of the x, y, and w may be the same vector.

    1453:    Concepts: vector^pointwise multiply

    1455: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
    1456: @*/
    1457: PetscErrorCode PETSCVEC_DLLEXPORT VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
    1458: {

    1470:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1471:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1473:   (*w->ops->pointwisemaxabs)(w,x,y);
    1474:   PetscObjectStateIncrease((PetscObject)w);
    1475:   return(0);
    1476: }

    1480: /*@
    1481:    VecPointwiseDivide - Computes the componentwise division w = x/y.

    1483:    Collective on Vec

    1485:    Input Parameters:
    1486: .  x, y  - the vectors

    1488:    Output Parameter:
    1489: .  w - the result

    1491:    Level: advanced

    1493:    Notes: any subset of the x, y, and w may be the same vector.

    1495:    Concepts: vector^pointwise divide

    1497: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
    1498: @*/
    1499: PetscErrorCode PETSCVEC_DLLEXPORT VecPointwiseDivide(Vec w,Vec x,Vec y)
    1500: {

    1512:   if (x->N != y->N || x->N != w->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1513:   if (x->n != y->n || x->n != w->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1515:   (*w->ops->pointwisedivide)(w,x,y);
    1516:   PetscObjectStateIncrease((PetscObject)w);
    1517:   return(0);
    1518: }

    1522: /*@
    1523:    VecMaxPointwiseDivide - Computes the maximum of the componentwise division w = abs(x/y).

    1525:    Collective on Vec

    1527:    Input Parameters:
    1528: .  x, y  - the vectors

    1530:    Output Parameter:
    1531: .  max - the result

    1533:    Level: advanced

    1535:    Notes: any subset of the x, y, and w may be the same vector.

    1537: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs()
    1538: @*/
    1539: PetscErrorCode PETSCVEC_DLLEXPORT VecMaxPointwiseDivide(Vec x,Vec y,PetscReal *max)
    1540: {

    1550:   if (x->N != y->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    1551:   if (x->n != y->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    1553:   (*x->ops->maxpointwisedivide)(x,y,max);
    1554:   return(0);
    1555: }

    1559: /*@C
    1560:    VecDuplicate - Creates a new vector of the same type as an existing vector.

    1562:    Collective on Vec

    1564:    Input Parameters:
    1565: .  v - a vector to mimic

    1567:    Output Parameter:
    1568: .  newv - location to put new vector

    1570:    Notes:
    1571:    VecDuplicate() does not copy the vector, but rather allocates storage
    1572:    for the new vector.  Use VecCopy() to copy a vector.

    1574:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
    1575:    vectors.

    1577:    Level: beginner

    1579: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
    1580: @*/
    1581: PetscErrorCode PETSCVEC_DLLEXPORT VecDuplicate(Vec x,Vec *newv)
    1582: {

    1589:   (*x->ops->duplicate)(x,newv);
    1590:   PetscObjectStateIncrease((PetscObject)*newv);
    1591:   return(0);
    1592: }

    1596: /*@C
    1597:    VecDestroy - Destroys a vector.

    1599:    Collective on Vec

    1601:    Input Parameters:
    1602: .  v  - the vector

    1604:    Level: beginner

    1606: .seealso: VecDuplicate(), VecDestroyVecs()
    1607: @*/
    1608: PetscErrorCode PETSCVEC_DLLEXPORT VecDestroy(Vec v)
    1609: {

    1614:   if (--v->refct > 0) return(0);
    1615:   /* destroy the internal part */
    1616:   if (v->ops->destroy) {
    1617:     (*v->ops->destroy)(v);
    1618:   }
    1619:   /* destroy the external/common part */
    1620:   if (v->mapping) {
    1621:     ISLocalToGlobalMappingDestroy(v->mapping);
    1622:   }
    1623:   if (v->bmapping) {
    1624:     ISLocalToGlobalMappingDestroy(v->bmapping);
    1625:   }
    1626:   if (v->map) {
    1627:     PetscMapDestroy(v->map);
    1628:   }
    1629:   PetscHeaderDestroy(v);
    1630:   return(0);
    1631: }

    1635: /*@C
    1636:    VecDuplicateVecs - Creates several vectors of the same type as an existing vector.

    1638:    Collective on Vec

    1640:    Input Parameters:
    1641: +  m - the number of vectors to obtain
    1642: -  v - a vector to mimic

    1644:    Output Parameter:
    1645: .  V - location to put pointer to array of vectors

    1647:    Notes:
    1648:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
    1649:    vector.

    1651:    Fortran Note:
    1652:    The Fortran interface is slightly different from that given below, it
    1653:    requires one to pass in V a Vec (integer) array of size at least m.
    1654:    See the Fortran chapter of the users manual and petsc/src/vec/examples for details.

    1656:    Level: intermediate

    1658: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
    1659: @*/
    1660: PetscErrorCode PETSCVEC_DLLEXPORT VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
    1661: {

    1668:   (*v->ops->duplicatevecs)(v, m,V);
    1669:   return(0);
    1670: }

    1674: /*@C
    1675:    VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().

    1677:    Collective on Vec

    1679:    Input Parameters:
    1680: +  vv - pointer to array of vector pointers
    1681: -  m - the number of vectors previously obtained

    1683:    Fortran Note:
    1684:    The Fortran interface is slightly different from that given below.
    1685:    See the Fortran chapter of the users manual and
    1686:    petsc/src/vec/examples for details.

    1688:    Level: intermediate

    1690: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
    1691: @*/
    1692: PetscErrorCode PETSCVEC_DLLEXPORT VecDestroyVecs(Vec vv[],PetscInt m)
    1693: {

    1700:   (*(*vv)->ops->destroyvecs)(vv,m);
    1701:   return(0);
    1702: }

    1706: /*@
    1707:    VecSetValues - Inserts or adds values into certain locations of a vector.

    1709:    Not Collective

    1711:    Input Parameters:
    1712: +  x - vector to insert in
    1713: .  ni - number of elements to add
    1714: .  ix - indices where to add
    1715: .  y - array of values
    1716: -  iora - either INSERT_VALUES or ADD_VALUES, where
    1717:    ADD_VALUES adds values to any existing entries, and
    1718:    INSERT_VALUES replaces existing entries with new values

    1720:    Notes:
    1721:    VecSetValues() sets x[ix[i]] = y[i], for i=0,...,ni-1.

    1723:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
    1724:    options cannot be mixed without intervening calls to the assembly
    1725:    routines.

    1727:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
    1728:    MUST be called after all calls to VecSetValues() have been completed.

    1730:    VecSetValues() uses 0-based indices in Fortran as well as in C.

    1732:    Negative indices may be passed in ix, these rows are
    1733:    simply ignored. This allows easily inserting element load matrices
    1734:    with homogeneous Dirchlet boundary conditions that you don't want represented
    1735:    in the vector.

    1737:    Level: beginner

    1739:    Concepts: vector^setting values

    1741: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesLocal(),
    1742:            VecSetValue(), VecSetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecGetValues()
    1743: @*/
    1744: PetscErrorCode PETSCVEC_DLLEXPORT VecSetValues(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
    1745: {

    1753:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
    1754:   (*x->ops->setvalues)(x,ni,ix,y,iora);
    1755:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
    1756:   PetscObjectStateIncrease((PetscObject)x);
    1757:   return(0);
    1758: }

    1762: /*@
    1763:    VecGetValues - Gets values from certain locations of a vector. Currently
    1764:           can only get values on the same processor

    1766:     Collective on Vec
    1767:
    1768:    Input Parameters:
    1769: +  x - vector to get values from
    1770: .  ni - number of elements to get
    1771: -  ix - indices where to get them from (in global 1d numbering)

    1773:    Output Parameter:
    1774: .   y - array of values

    1776:    Notes:
    1777:    The user provides the allocated array y; it is NOT allocated in this routine

    1779:    VecGetValues() gets y[i] = x[ix[i]], for i=0,...,ni-1.

    1781:    VecAssemblyBegin() and VecAssemblyEnd()  MUST be called before calling this

    1783:    VecGetValues() uses 0-based indices in Fortran as well as in C.

    1785:    Level: beginner

    1787:    Concepts: vector^getting values

    1789: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecGetValuesLocal(),
    1790:            VecGetValue(), VecGetValuesBlocked(), InsertMode, INSERT_VALUES, ADD_VALUES, VecSetValues()
    1791: @*/
    1792: PetscErrorCode PETSCVEC_DLLEXPORT VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
    1793: {

    1801:   (*x->ops->getvalues)(x,ni,ix,y);
    1802:   return(0);
    1803: }

    1807: /*@
    1808:    VecSetValuesBlocked - Inserts or adds blocks of values into certain locations of a vector.

    1810:    Not Collective

    1812:    Input Parameters:
    1813: +  x - vector to insert in
    1814: .  ni - number of blocks to add
    1815: .  ix - indices where to add in block count, rather than element count
    1816: .  y - array of values
    1817: -  iora - either INSERT_VALUES or ADD_VALUES, where
    1818:    ADD_VALUES adds values to any existing entries, and
    1819:    INSERT_VALUES replaces existing entries with new values

    1821:    Notes:
    1822:    VecSetValuesBlocked() sets x[bs*ix[i]+j] = y[bs*i+j],
    1823:    for j=0,...,bs, for i=0,...,ni-1. where bs was set with VecSetBlockSize().

    1825:    Calls to VecSetValuesBlocked() with the INSERT_VALUES and ADD_VALUES
    1826:    options cannot be mixed without intervening calls to the assembly
    1827:    routines.

    1829:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
    1830:    MUST be called after all calls to VecSetValuesBlocked() have been completed.

    1832:    VecSetValuesBlocked() uses 0-based indices in Fortran as well as in C.

    1834:    Negative indices may be passed in ix, these rows are
    1835:    simply ignored. This allows easily inserting element load matrices
    1836:    with homogeneous Dirchlet boundary conditions that you don't want represented
    1837:    in the vector.

    1839:    Level: intermediate

    1841:    Concepts: vector^setting values blocked

    1843: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValuesBlockedLocal(),
    1844:            VecSetValues()
    1845: @*/
    1846: PetscErrorCode PETSCVEC_DLLEXPORT VecSetValuesBlocked(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
    1847: {

    1855:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
    1856:   (*x->ops->setvaluesblocked)(x,ni,ix,y,iora);
    1857:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
    1858:   PetscObjectStateIncrease((PetscObject)x);
    1859:   return(0);
    1860: }

    1864: /*@
    1865:    VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
    1866:    by the routine VecSetValuesLocal() to allow users to insert vector entries
    1867:    using a local (per-processor) numbering.

    1869:    Collective on Vec

    1871:    Input Parameters:
    1872: +  x - vector
    1873: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

    1875:    Notes:
    1876:    All vectors obtained with VecDuplicate() from this vector inherit the same mapping.

    1878:    Level: intermediate

    1880:    Concepts: vector^setting values with local numbering

    1882: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
    1883:            VecSetLocalToGlobalMappingBlocked(), VecSetValuesBlockedLocal()
    1884: @*/
    1885: PetscErrorCode PETSCVEC_DLLEXPORT VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
    1886: {

    1893:   if (x->mapping) {
    1894:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for vector");
    1895:   }

    1897:   if (x->ops->setlocaltoglobalmapping) {
    1898:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
    1899:   } else {
    1900:     x->mapping = mapping;
    1901:     PetscObjectReference((PetscObject)mapping);
    1902:   }
    1903:   return(0);
    1904: }

    1908: /*@
    1909:    VecSetLocalToGlobalMappingBlock - Sets a local numbering to global numbering used
    1910:    by the routine VecSetValuesBlockedLocal() to allow users to insert vector entries
    1911:    using a local (per-processor) numbering.

    1913:    Collective on Vec

    1915:    Input Parameters:
    1916: +  x - vector
    1917: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

    1919:    Notes:
    1920:    All vectors obtained with VecDuplicate() from this vector inherit the same mapping.

    1922:    Level: intermediate

    1924:    Concepts: vector^setting values blocked with local numbering

    1926: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
    1927:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
    1928: @*/
    1929: PetscErrorCode PETSCVEC_DLLEXPORT VecSetLocalToGlobalMappingBlock(Vec x,ISLocalToGlobalMapping mapping)
    1930: {

    1937:   if (x->bmapping) {
    1938:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Mapping already set for vector");
    1939:   }
    1940:   x->bmapping = mapping;
    1941:   PetscObjectReference((PetscObject)mapping);
    1942:   return(0);
    1943: }

    1947: /*@
    1948:    VecSetValuesLocal - Inserts or adds values into certain locations of a vector,
    1949:    using a local ordering of the nodes.

    1951:    Not Collective

    1953:    Input Parameters:
    1954: +  x - vector to insert in
    1955: .  ni - number of elements to add
    1956: .  ix - indices where to add
    1957: .  y - array of values
    1958: -  iora - either INSERT_VALUES or ADD_VALUES, where
    1959:    ADD_VALUES adds values to any existing entries, and
    1960:    INSERT_VALUES replaces existing entries with new values

    1962:    Level: intermediate

    1964:    Notes:
    1965:    VecSetValuesLocal() sets x[ix[i]] = y[i], for i=0,...,ni-1.

    1967:    Calls to VecSetValues() with the INSERT_VALUES and ADD_VALUES
    1968:    options cannot be mixed without intervening calls to the assembly
    1969:    routines.

    1971:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
    1972:    MUST be called after all calls to VecSetValuesLocal() have been completed.

    1974:    VecSetValuesLocal() uses 0-based indices in Fortran as well as in C.

    1976:    Concepts: vector^setting values with local numbering

    1978: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetLocalToGlobalMapping(),
    1979:            VecSetValuesBlockedLocal()
    1980: @*/
    1981: PetscErrorCode PETSCVEC_DLLEXPORT VecSetValuesLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
    1982: {
    1984:   PetscInt       lixp[128],*lix = lixp;

    1992:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
    1993:   if (!x->ops->setvalueslocal) {
    1994:     if (!x->mapping) {
    1995:       SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMapping()");
    1996:     }
    1997:     if (ni > 128) {
    1998:       PetscMalloc(ni*sizeof(PetscInt),&lix);
    1999:     }
    2000:     ISLocalToGlobalMappingApply(x->mapping,ni,(PetscInt*)ix,lix);
    2001:     (*x->ops->setvalues)(x,ni,lix,y,iora);
    2002:     if (ni > 128) {
    2003:       PetscFree(lix);
    2004:     }
    2005:   } else {
    2006:     (*x->ops->setvalueslocal)(x,ni,ix,y,iora);
    2007:   }
    2008:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
    2009:   PetscObjectStateIncrease((PetscObject)x);
    2010:   return(0);
    2011: }

    2015: /*@
    2016:    VecSetValuesBlockedLocal - Inserts or adds values into certain locations of a vector,
    2017:    using a local ordering of the nodes.

    2019:    Not Collective

    2021:    Input Parameters:
    2022: +  x - vector to insert in
    2023: .  ni - number of blocks to add
    2024: .  ix - indices where to add in block count, not element count
    2025: .  y - array of values
    2026: -  iora - either INSERT_VALUES or ADD_VALUES, where
    2027:    ADD_VALUES adds values to any existing entries, and
    2028:    INSERT_VALUES replaces existing entries with new values

    2030:    Level: intermediate

    2032:    Notes:
    2033:    VecSetValuesBlockedLocal() sets x[bs*ix[i]+j] = y[bs*i+j],
    2034:    for j=0,..bs-1, for i=0,...,ni-1, where bs has been set with VecSetBlockSize().

    2036:    Calls to VecSetValuesBlockedLocal() with the INSERT_VALUES and ADD_VALUES
    2037:    options cannot be mixed without intervening calls to the assembly
    2038:    routines.

    2040:    These values may be cached, so VecAssemblyBegin() and VecAssemblyEnd()
    2041:    MUST be called after all calls to VecSetValuesBlockedLocal() have been completed.

    2043:    VecSetValuesBlockedLocal() uses 0-based indices in Fortran as well as in C.

    2046:    Concepts: vector^setting values blocked with local numbering

    2048: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesBlocked(),
    2049:            VecSetLocalToGlobalMappingBlocked()
    2050: @*/
    2051: PetscErrorCode PETSCVEC_DLLEXPORT VecSetValuesBlockedLocal(Vec x,PetscInt ni,const PetscInt ix[],const PetscScalar y[],InsertMode iora)
    2052: {
    2054:   PetscInt       lixp[128],*lix = lixp;

    2061:   if (!x->bmapping) {
    2062:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Local to global never set with VecSetLocalToGlobalMappingBlocked()");
    2063:   }
    2064:   if (ni > 128) {
    2065:     PetscMalloc(ni*sizeof(PetscInt),&lix);
    2066:   }

    2068:   PetscLogEventBegin(VEC_SetValues,x,0,0,0);
    2069:   ISLocalToGlobalMappingApply(x->bmapping,ni,(PetscInt*)ix,lix);
    2070:   (*x->ops->setvaluesblocked)(x,ni,lix,y,iora);
    2071:   PetscLogEventEnd(VEC_SetValues,x,0,0,0);
    2072:   if (ni > 128) {
    2073:     PetscFree(lix);
    2074:   }
    2075:   PetscObjectStateIncrease((PetscObject)x);
    2076:   return(0);
    2077: }

    2081: /*@
    2082:    VecAssemblyBegin - Begins assembling the vector.  This routine should
    2083:    be called after completing all calls to VecSetValues().

    2085:    Collective on Vec

    2087:    Input Parameter:
    2088: .  vec - the vector

    2090:    Level: beginner

    2092:    Concepts: assembly^vectors

    2094: .seealso: VecAssemblyEnd(), VecSetValues()
    2095: @*/
    2096: PetscErrorCode PETSCVEC_DLLEXPORT VecAssemblyBegin(Vec vec)
    2097: {
    2099:   PetscTruth     flg;

    2105:   PetscOptionsHasName(vec->prefix,"-vec_view_stash",&flg);
    2106:   if (flg) {
    2107:     VecStashView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
    2108:   }

    2110:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
    2111:   if (vec->ops->assemblybegin) {
    2112:     (*vec->ops->assemblybegin)(vec);
    2113:   }
    2114:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
    2115:   PetscObjectStateIncrease((PetscObject)vec);
    2116:   return(0);
    2117: }

    2121: /*@
    2122:    VecAssemblyEnd - Completes assembling the vector.  This routine should
    2123:    be called after VecAssemblyBegin().

    2125:    Collective on Vec

    2127:    Input Parameter:
    2128: .  vec - the vector

    2130:    Options Database Keys:
    2131: +  -vec_view - Prints vector in ASCII format
    2132: .  -vec_view_matlab - Prints vector in ASCII Matlab format to stdout
    2133: .  -vec_view_matlab_file - Prints vector in Matlab format to matlaboutput.mat
    2134: .  -vec_view_draw - Activates vector viewing using drawing tools
    2135: .  -display <name> - Sets display name (default is host)
    2136: .  -draw_pause <sec> - Sets number of seconds to pause after display
    2137: -  -vec_view_socket - Activates vector viewing using a socket
    2138:
    2139:    Level: beginner

    2141: .seealso: VecAssemblyBegin(), VecSetValues()
    2142: @*/
    2143: PetscErrorCode PETSCVEC_DLLEXPORT VecAssemblyEnd(Vec vec)
    2144: {
    2146:   PetscTruth     flg;

    2150:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
    2152:   if (vec->ops->assemblyend) {
    2153:     (*vec->ops->assemblyend)(vec);
    2154:   }
    2155:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
    2156:   PetscOptionsBegin(vec->comm,vec->prefix,"Vector Options","Vec");
    2157:     PetscOptionsName("-vec_view","Print vector to stdout","VecView",&flg);
    2158:     if (flg) {
    2159:       VecView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
    2160:     }
    2161:     PetscOptionsName("-vec_view_matlab","Print vector to stdout in a format Matlab can read","VecView",&flg);
    2162:     if (flg) {
    2163:       PetscViewerPushFormat(PETSC_VIEWER_STDOUT_(vec->comm),PETSC_VIEWER_ASCII_MATLAB);
    2164:       VecView(vec,PETSC_VIEWER_STDOUT_(vec->comm));
    2165:       PetscViewerPopFormat(PETSC_VIEWER_STDOUT_(vec->comm));
    2166:     }
    2167: #if defined(PETSC_HAVE_MATLAB)
    2168:     PetscOptionsName("-vec_view_matlab_file","Print vector to matlaboutput.mat format Matlab can read","VecView",&flg);
    2169:     if (flg) {
    2170:       VecView(vec,PETSC_VIEWER_MATLAB_(vec->comm));
    2171:     }
    2172: #endif
    2173: #if defined(PETSC_USE_SOCKET_VIEWER)
    2174:     PetscOptionsName("-vec_view_socket","Send vector to socket (can be read from matlab)","VecView",&flg);
    2175:     if (flg) {
    2176:       VecView(vec,PETSC_VIEWER_SOCKET_(vec->comm));
    2177:       PetscViewerFlush(PETSC_VIEWER_SOCKET_(vec->comm));
    2178:     }
    2179: #endif
    2180:     PetscOptionsName("-vec_view_binary","Save vector to file in binary format","VecView",&flg);
    2181:     if (flg) {
    2182:       VecView(vec,PETSC_VIEWER_BINARY_(vec->comm));
    2183:       PetscViewerFlush(PETSC_VIEWER_BINARY_(vec->comm));
    2184:     }
    2185:   PetscOptionsEnd();
    2186:   /* These invoke PetscDrawGetDraw which invokes PetscOptionsBegin/End, */
    2187:   /* hence they should not be inside the above PetscOptionsBegin/End block. */
    2188:   PetscOptionsHasName(vec->prefix,"-vec_view_draw",&flg);
    2189:   if (flg) {
    2190:     VecView(vec,PETSC_VIEWER_DRAW_(vec->comm));
    2191:     PetscViewerFlush(PETSC_VIEWER_DRAW_(vec->comm));
    2192:   }
    2193:   PetscOptionsHasName(vec->prefix,"-vec_view_draw_lg",&flg);
    2194:   if (flg) {
    2195:     PetscViewerSetFormat(PETSC_VIEWER_DRAW_(vec->comm),PETSC_VIEWER_DRAW_LG);
    2196:     VecView(vec,PETSC_VIEWER_DRAW_(vec->comm));
    2197:     PetscViewerFlush(PETSC_VIEWER_DRAW_(vec->comm));
    2198:   }
    2199:   return(0);
    2200: }

    2205: /*@C
    2206:    VecMTDot - Computes indefinite vector multiple dot products.
    2207:    That is, it does NOT use the complex conjugate.

    2209:    Collective on Vec

    2211:    Input Parameters:
    2212: +  nv - number of vectors
    2213: .  x - one vector
    2214: -  y - array of vectors.  Note that vectors are pointers

    2216:    Output Parameter:
    2217: .  val - array of the dot products

    2219:    Notes for Users of Complex Numbers:
    2220:    For complex vectors, VecMTDot() computes the indefinite form
    2221: $      val = (x,y) = y^T x,
    2222:    where y^T denotes the transpose of y.

    2224:    Use VecMDot() for the inner product
    2225: $      val = (x,y) = y^H x,
    2226:    where y^H denotes the conjugate transpose of y.

    2228:    Level: intermediate

    2230:    Concepts: inner product^multiple
    2231:    Concepts: vector^multiple inner products

    2233: .seealso: VecMDot(), VecTDot()
    2234: @*/
    2235: PetscErrorCode PETSCVEC_DLLEXPORT VecMTDot(PetscInt nv,Vec x,const Vec y[],PetscScalar *val)
    2236: {

    2247:   if (x->N != (*y)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    2248:   if (x->n != (*y)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    2250:   PetscLogEventBegin(VEC_MTDot,x,*y,0,0);
    2251:   (*x->ops->mtdot)(nv,x,y,val);
    2252:   PetscLogEventEnd(VEC_MTDot,x,*y,0,0);
    2253:   return(0);
    2254: }

    2258: /*@C
    2259:    VecMDot - Computes vector multiple dot products.

    2261:    Collective on Vec

    2263:    Input Parameters:
    2264: +  nv - number of vectors
    2265: .  x - one vector
    2266: -  y - array of vectors.

    2268:    Output Parameter:
    2269: .  val - array of the dot products

    2271:    Notes for Users of Complex Numbers:
    2272:    For complex vectors, VecMDot() computes
    2273: $     val = (x,y) = y^H x,
    2274:    where y^H denotes the conjugate transpose of y.

    2276:    Use VecMTDot() for the indefinite form
    2277: $     val = (x,y) = y^T x,
    2278:    where y^T denotes the transpose of y.

    2280:    Level: intermediate

    2282:    Concepts: inner product^multiple
    2283:    Concepts: vector^multiple inner products

    2285: .seealso: VecMTDot(), VecDot()
    2286: @*/
    2287: PetscErrorCode PETSCVEC_DLLEXPORT VecMDot(PetscInt nv,Vec x,const Vec y[],PetscScalar *val)
    2288: {

    2299:   if (x->N != (*y)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    2300:   if (x->n != (*y)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    2302:   PetscLogEventBarrierBegin(VEC_MDotBarrier,x,*y,0,0,x->comm);
    2303:   (*x->ops->mdot)(nv,x,y,val);
    2304:   PetscLogEventBarrierEnd(VEC_MDotBarrier,x,*y,0,0,x->comm);
    2305:   return(0);
    2306: }

    2310: /*@C
    2311:    VecMAXPY - Computes y = y + sum alpha[j] x[j]

    2313:    Collective on Vec

    2315:    Input Parameters:
    2316: +  nv - number of scalars and x-vectors
    2317: .  alpha - array of scalars
    2318: .  y - one vector
    2319: -  x - array of vectors

    2321:    Level: intermediate

    2323:    Concepts: BLAS

    2325: .seealso: VecAXPY(), VecWAXPY(), VecAYPX()
    2326: @*/
    2327: PetscErrorCode PETSCVEC_DLLEXPORT VecMAXPY(Vec y,PetscInt nv,const PetscScalar alpha[],Vec *x)
    2328: {

    2339:   if (y->N != (*x)->N) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
    2340:   if (y->n != (*x)->n) SETERRQ(PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

    2342:   PetscLogEventBegin(VEC_MAXPY,*x,y,0,0);
    2343:   (*y->ops->maxpy)(y,nv,alpha,x);
    2344:   PetscLogEventEnd(VEC_MAXPY,*x,y,0,0);
    2345:   PetscObjectStateIncrease((PetscObject)y);
    2346:   return(0);
    2347: }

    2349: /*MC
    2350:    VecGetArray - Returns a pointer to a contiguous array that contains this
    2351:    processor's portion of the vector data. For the standard PETSc
    2352:    vectors, VecGetArray() returns a pointer to the local data array and
    2353:    does not use any copies. If the underlying vector data is not stored
    2354:    in a contiquous array this routine will copy the data to a contiquous
    2355:    array and return a pointer to that. You MUST call VecRestoreArray()
    2356:    when you no longer need access to the array.

    2358:    Synopsis:
    2359:    PetscErrorCode VecGetArray(Vec x,PetscScalar *a[])

    2361:    Not Collective

    2363:    Input Parameter:
    2364: .  x - the vector

    2366:    Output Parameter:
    2367: .  a - location to put pointer to the array

    2369:    Fortran Note:
    2370:    This routine is used differently from Fortran 77
    2371: $    Vec         x
    2372: $    PetscScalar x_array(1)
    2373: $    PetscOffset i_x
    2374: $    PetscErrorCode ierr
    2375: $       call VecGetArray(x,x_array,i_x,ierr)
    2376: $
    2377: $   Access first local entry in vector with
    2378: $      value = x_array(i_x + 1)
    2379: $
    2380: $      ...... other code
    2381: $       call VecRestoreArray(x,x_array,i_x,ierr)
    2382:    For Fortran 90 see VecGetArrayF90()

    2384:    See the Fortran chapter of the users manual and
    2385:    petsc/src/snes/examples/tutorials/ex5f.F for details.

    2387:    Level: beginner

    2389:    Concepts: vector^accessing local values

    2391: .seealso: VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(), VecGetArray2d()
    2392: M*/
    2395: PetscErrorCode VecGetArray_Private(Vec x,PetscScalar *a[])
    2396: {

    2403:   (*x->ops->getarray)(x,a);
    2404:   return(0);
    2405: }

    2410: /*@C
    2411:    VecGetArrays - Returns a pointer to the arrays in a set of vectors
    2412:    that were created by a call to VecDuplicateVecs().  You MUST call
    2413:    VecRestoreArrays() when you no longer need access to the array.

    2415:    Not Collective

    2417:    Input Parameter:
    2418: +  x - the vectors
    2419: -  n - the number of vectors

    2421:    Output Parameter:
    2422: .  a - location to put pointer to the array

    2424:    Fortran Note:
    2425:    This routine is not supported in Fortran.

    2427:    Level: intermediate

    2429: .seealso: VecGetArray(), VecRestoreArrays()
    2430: @*/
    2431: PetscErrorCode PETSCVEC_DLLEXPORT VecGetArrays(const Vec x[],PetscInt n,PetscScalar **a[])
    2432: {
    2434:   PetscInt       i;
    2435:   PetscScalar    **q;

    2441:   if (n <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Must get at least one array n = %D",n);
    2442:   PetscMalloc(n*sizeof(PetscScalar*),&q);
    2443:   for (i=0; i<n; ++i) {
    2444:     VecGetArray(x[i],&q[i]);
    2445:   }
    2446:   *a = q;
    2447:   return(0);
    2448: }

    2452: /*@C
    2453:    VecRestoreArrays - Restores a group of vectors after VecGetArrays()
    2454:    has been called.

    2456:    Not Collective

    2458:    Input Parameters:
    2459: +  x - the vector
    2460: .  n - the number of vectors
    2461: -  a - location of pointer to arrays obtained from VecGetArrays()

    2463:    Notes:
    2464:    For regular PETSc vectors this routine does not involve any copies. For
    2465:    any special vectors that do not store local vector data in a contiguous
    2466:    array, this routine will copy the data back into the underlying
    2467:    vector data structure from the arrays obtained with VecGetArrays().

    2469:    Fortran Note:
    2470:    This routine is not supported in Fortran.

    2472:    Level: intermediate

    2474: .seealso: VecGetArrays(), VecRestoreArray()
    2475: @*/
    2476: PetscErrorCode PETSCVEC_DLLEXPORT VecRestoreArrays(const Vec x[],PetscInt n,PetscScalar **a[])
    2477: {
    2479:   PetscInt       i;
    2480:   PetscScalar    **q = *a;

    2487:   for(i=0;i<n;++i) {
    2488:     VecRestoreArray(x[i],&q[i]);
    2489:  }
    2490:   PetscFree(q);
    2491:   return(0);
    2492: }

    2494: /*MC
    2495:    VecRestoreArray - Restores a vector after VecGetArray() has been called.

    2497:    Not Collective

    2499:    Synopsis:
    2500:    PetscErrorCode VecRestoreArray(Vec x,PetscScalar *a[])

    2502:    Input Parameters:
    2503: +  x - the vector
    2504: -  a - location of pointer to array obtained from VecGetArray()

    2506:    Level: beginner

    2508:    Notes:
    2509:    For regular PETSc vectors this routine does not involve any copies. For
    2510:    any special vectors that do not store local vector data in a contiguous
    2511:    array, this routine will copy the data back into the underlying
    2512:    vector data structure from the array obtained with VecGetArray().

    2514:    This routine actually zeros out the a pointer. This is to prevent accidental
    2515:    us of the array after it has been restored. If you pass null for a it will
    2516:    not zero the array pointer a.

    2518:    Fortran Note:
    2519:    This routine is used differently from Fortran 77
    2520: $    Vec         x
    2521: $    PetscScalar x_array(1)
    2522: $    PetscOffset i_x
    2523: $    PetscErrorCode ierr
    2524: $       call VecGetArray(x,x_array,i_x,ierr)
    2525: $
    2526: $   Access first local entry in vector with
    2527: $      value = x_array(i_x + 1)
    2528: $
    2529: $      ...... other code
    2530: $       call VecRestoreArray(x,x_array,i_x,ierr)

    2532:    See the Fortran chapter of the users manual and
    2533:    petsc/src/snes/examples/tutorials/ex5f.F for details.
    2534:    For Fortran 90 see VecRestoreArrayF90()

    2536: .seealso: VecGetArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(), VecRestoreArray2d()
    2537: M*/
    2540: PetscErrorCode VecRestoreArray_Private(Vec x,PetscScalar *a[])
    2541: {

    2548: #if defined(PETSC_USE_DEBUG)
    2549:   CHKMEMQ;
    2550: #endif
    2551:   if (x->ops->restorearray) {
    2552:     (*x->ops->restorearray)(x,a);
    2553:   }
    2554:   PetscObjectStateIncrease((PetscObject)x);
    2555:   return(0);
    2556: }

    2558: #undef  __FUNCT__
    2560: /*@
    2561:   VecViewFromOptions - This function visualizes the vector based upon user options.

    2563:   Collective on Vec

    2565:   Input Parameters:
    2566: . vec   - The vector
    2567: . title - The title

    2569:   Level: intermediate

    2571: .keywords: Vec, view, options, database
    2572: .seealso: VecSetFromOptions(), VecView()
    2573: @*/
    2574: PetscErrorCode PETSCVEC_DLLEXPORT VecViewFromOptions(Vec vec, char *title)
    2575: {
    2576:   PetscViewer    viewer;
    2577:   PetscDraw      draw;
    2578:   PetscTruth     opt;
    2579:   char           *titleStr;
    2580:   char           typeName[1024];
    2581:   char           fileName[PETSC_MAX_PATH_LEN];
    2582:   size_t         len;

    2586:   PetscOptionsHasName(vec->prefix, "-vec_view", &opt);
    2587:   if (opt) {
    2588:     PetscOptionsGetString(vec->prefix, "-vec_view", typeName, 1024, &opt);
    2589:     PetscStrlen(typeName, &len);
    2590:     if (len > 0) {
    2591:       PetscViewerCreate(vec->comm, &viewer);
    2592:       PetscViewerSetType(viewer, typeName);
    2593:       PetscOptionsGetString(vec->prefix, "-vec_view_file", fileName, 1024, &opt);
    2594:       if (opt) {
    2595:         PetscViewerSetFilename(viewer, fileName);
    2596:       } else {
    2597:         PetscViewerSetFilename(viewer, vec->name);
    2598:       }
    2599:       VecView(vec, viewer);
    2600:       PetscViewerFlush(viewer);
    2601:       PetscViewerDestroy(viewer);
    2602:     } else {
    2603:       VecView(vec, PETSC_VIEWER_STDOUT_(vec->comm));
    2604:     }
    2605:   }
    2606:   PetscOptionsHasName(vec->prefix, "-vec_view_draw", &opt);
    2607:   if (opt) {
    2608:     PetscViewerDrawOpen(vec->comm, 0, 0, 0, 0, 300, 300, &viewer);
    2609:     PetscViewerDrawGetDraw(viewer, 0, &draw);
    2610:     if (title) {
    2611:       titleStr = title;
    2612:     } else {
    2613:       PetscObjectName((PetscObject) vec);
    2614:       titleStr = vec->name;
    2615:     }
    2616:     PetscDrawSetTitle(draw, titleStr);
    2617:     VecView(vec, viewer);
    2618:     PetscViewerFlush(viewer);
    2619:     PetscDrawPause(draw);
    2620:     PetscViewerDestroy(viewer);
    2621:   }
    2622:   return(0);
    2623: }

    2627: /*@C
    2628:    VecView - Views a vector object.

    2630:    Collective on Vec

    2632:    Input Parameters:
    2633: +  v - the vector
    2634: -  viewer - an optional visualization context

    2636:    Notes:
    2637:    The available visualization contexts include
    2638: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
    2639: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
    2640:          output where only the first processor opens
    2641:          the file.  All other processors send their
    2642:          data to the first processor to print.

    2644:    You can change the format the vector is printed using the
    2645:    option PetscViewerSetFormat().

    2647:    The user can open alternative visualization contexts with
    2648: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
    2649: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
    2650:          specified file; corresponding input uses VecLoad()
    2651: .    PetscViewerDrawOpen() - Outputs vector to an X window display
    2652: -    PetscViewerSocketOpen() - Outputs vector to Socket viewer

    2654:    The user can call PetscViewerSetFormat() to specify the output
    2655:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
    2656:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
    2657: +    PETSC_VIEWER_ASCII_DEFAULT - default, prints vector contents
    2658: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in Matlab format
    2659: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
    2660: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
    2661:          format common among all vector types

    2663:    Level: beginner

    2665:    Concepts: vector^printing
    2666:    Concepts: vector^saving to disk

    2668: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
    2669:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
    2670:           PetscRealView(), PetscScalarView(), PetscIntView()
    2671: @*/
    2672: PetscErrorCode PETSCVEC_DLLEXPORT VecView(Vec vec,PetscViewer viewer)
    2673: {
    2674:   PetscErrorCode    ierr;
    2675:   PetscViewerFormat format;

    2680:   if (!viewer) viewer = PETSC_VIEWER_STDOUT_(vec->comm);
    2683:   if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");

    2685:   /*
    2686:      Check if default viewer has been overridden, but user request it anyways
    2687:   */
    2688:   PetscViewerGetFormat(viewer,&format);
    2689:   if (vec->ops->viewnative && format == PETSC_VIEWER_NATIVE) {
    2690:     PetscViewerPopFormat(viewer);
    2691:     (*vec->ops->viewnative)(vec,viewer);
    2692:     PetscViewerPushFormat(viewer,PETSC_VIEWER_NATIVE);
    2693:   } else {
    2694:     (*vec->ops->view)(vec,viewer);
    2695:   }
    2696:   return(0);
    2697: }

    2701: /*@
    2702:    VecGetSize - Returns the global number of elements of the vector.

    2704:    Not Collective

    2706:    Input Parameter:
    2707: .  x - the vector

    2709:    Output Parameters:
    2710: .  size - the global length of the vector

    2712:    Level: beginner

    2714:    Concepts: vector^local size

    2716: .seealso: VecGetLocalSize()
    2717: @*/
    2718: PetscErrorCode PETSCVEC_DLLEXPORT VecGetSize(Vec x,PetscInt *size)
    2719: {

    2726:   (*x->ops->getsize)(x,size);
    2727:   return(0);
    2728: }

    2732: /*@
    2733:    VecGetLocalSize - Returns the number of elements of the vector stored
    2734:    in local memory. This routine may be implementation dependent, so use
    2735:    with care.

    2737:    Not Collective

    2739:    Input Parameter:
    2740: .  x - the vector

    2742:    Output Parameter:
    2743: .  size - the length of the local piece of the vector

    2745:    Level: beginner

    2747:    Concepts: vector^size

    2749: .seealso: VecGetSize()
    2750: @*/
    2751: PetscErrorCode PETSCVEC_DLLEXPORT VecGetLocalSize(Vec x,PetscInt *size)
    2752: {

    2759:   (*x->ops->getlocalsize)(x,size);
    2760:   return(0);
    2761: }

    2765: /*@C
    2766:    VecGetOwnershipRange - Returns the range of indices owned by
    2767:    this processor, assuming that the vectors are laid out with the
    2768:    first n1 elements on the first processor, next n2 elements on the
    2769:    second, etc.  For certain parallel layouts this range may not be
    2770:    well defined.

    2772:    Not Collective

    2774:    Input Parameter:
    2775: .  x - the vector

    2777:    Output Parameters:
    2778: +  low - the first local element, pass in PETSC_NULL if not interested
    2779: -  high - one more than the last local element, pass in PETSC_NULL if not interested

    2781:    Note:
    2782:    The high argument is one more than the last element stored locally.

    2784:    Fortran: PETSC_NULL_INTEGER should be used instead of PETSC_NULL

    2786:    Level: beginner

    2788:    Concepts: ownership^of vectors
    2789:    Concepts: vector^ownership of elements

    2791: @*/
    2792: PetscErrorCode PETSCVEC_DLLEXPORT VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
    2793: {

    2801:   PetscMapGetLocalRange(x->map,low,high);
    2802:   return(0);
    2803: }

    2807: /*@C
    2808:    VecGetPetscMap - Returns the map associated with the vector

    2810:    Not Collective

    2812:    Input Parameter:
    2813: .  x - the vector

    2815:    Output Parameters:
    2816: .  map - the map

    2818:    Level: developer

    2820: @*/
    2821: PetscErrorCode PETSCVEC_DLLEXPORT VecGetPetscMap(Vec x,PetscMap *map)
    2822: {
    2827:   *map = x->map;
    2828:   return(0);
    2829: }

    2833: /*@
    2834:    VecSetOption - Sets an option for controling a vector's behavior.

    2836:    Collective on Vec

    2838:    Input Parameter:
    2839: +  x - the vector
    2840: -  op - the option

    2842:    Supported Options:
    2843: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
    2844:       entries destined to be stored on a seperate processor. This can be used
    2845:       to eliminate the global reduction in the VecAssemblyXXXX() if you know
    2846:       that you have only used VecSetValues() to set local elements
    2847: -   VEC_TREAT_OFF_PROC_ENTRIES restores the treatment of off processor entries.

    2849:    Level: intermediate

    2851: @*/
    2852: PetscErrorCode PETSCVEC_DLLEXPORT VecSetOption(Vec x,VecOption op)
    2853: {

    2859:   if (x->ops->setoption) {
    2860:     (*x->ops->setoption)(x,op);
    2861:   }
    2862:   return(0);
    2863: }

    2867: /* Default routines for obtaining and releasing; */
    2868: /* may be used by any implementation */
    2869: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
    2870: {
    2872:   PetscInt       i;

    2877:   if (m <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
    2878:   PetscMalloc(m*sizeof(Vec*),V);
    2879:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
    2880:   return(0);
    2881: }

    2885: PetscErrorCode VecDestroyVecs_Default(Vec v[], PetscInt m)
    2886: {
    2888:   PetscInt       i;

    2892:   if (m <= 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
    2893:   for (i=0; i<m; i++) {VecDestroy(v[i]);}
    2894:   PetscFree(v);
    2895:   return(0);
    2896: }

    2900: /*@
    2901:    VecPlaceArray - Allows one to replace the array in a vector with an
    2902:    array provided by the user. This is useful to avoid copying an array
    2903:    into a vector.

    2905:    Not Collective

    2907:    Input Parameters:
    2908: +  vec - the vector
    2909: -  array - the array

    2911:    Notes:
    2912:    You can return to the original array with a call to VecResetArray()

    2914:    Level: developer

    2916: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecResetArray()

    2918: @*/
    2919: PetscErrorCode PETSCVEC_DLLEXPORT VecPlaceArray(Vec vec,const PetscScalar array[])
    2920: {

    2927:   if (vec->ops->placearray) {
    2928:     (*vec->ops->placearray)(vec,array);
    2929:   } else {
    2930:     SETERRQ(PETSC_ERR_SUP,"Cannot place array in this type of vector");
    2931:   }
    2932:   PetscObjectStateIncrease((PetscObject)vec);
    2933:   return(0);
    2934: }

    2938: /*@
    2939:    VecResetArray - Resets a vector to use its default memory. Call this
    2940:    after the use of VecPlaceArray().

    2942:    Not Collective

    2944:    Input Parameters:
    2945: .  vec - the vector

    2947:    Level: developer

    2949: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()

    2951: @*/
    2952: PetscErrorCode PETSCVEC_DLLEXPORT VecResetArray(Vec vec)
    2953: {

    2959:   if (vec->ops->resetarray) {
    2960:     (*vec->ops->resetarray)(vec);
    2961:   } else {
    2962:     SETERRQ(PETSC_ERR_SUP,"Cannot reset array in this type of vector");
    2963:   }
    2964:   PetscObjectStateIncrease((PetscObject)vec);
    2965:   return(0);
    2966: }

    2970: /*@C
    2971:    VecReplaceArray - Allows one to replace the array in a vector with an
    2972:    array provided by the user. This is useful to avoid copying an array
    2973:    into a vector.

    2975:    Not Collective

    2977:    Input Parameters:
    2978: +  vec - the vector
    2979: -  array - the array

    2981:    Notes:
    2982:    This permanently replaces the array and frees the memory associated
    2983:    with the old array.

    2985:    The memory passed in MUST be obtained with PetscMalloc() and CANNOT be
    2986:    freed by the user. It will be freed when the vector is destroy.

    2988:    Not supported from Fortran

    2990:    Level: developer

    2992: .seealso: VecGetArray(), VecRestoreArray(), VecPlaceArray(), VecResetArray()

    2994: @*/
    2995: PetscErrorCode PETSCVEC_DLLEXPORT VecReplaceArray(Vec vec,const PetscScalar array[])
    2996: {

    3002:   if (vec->ops->replacearray) {
    3003:     (*vec->ops->replacearray)(vec,array);
    3004:  } else {
    3005:     SETERRQ(PETSC_ERR_SUP,"Cannot replace array in this type of vector");
    3006:   }
    3007:   PetscObjectStateIncrease((PetscObject)vec);
    3008:   return(0);
    3009: }

    3011: /*MC
    3012:     VecDuplicateVecsF90 - Creates several vectors of the same type as an existing vector
    3013:     and makes them accessible via a Fortran90 pointer.

    3015:     Synopsis:
    3016:     VecDuplicateVecsF90(Vec x,int n,{Vec, pointer :: y(:)},integer ierr)

    3018:     Collective on Vec

    3020:     Input Parameters:
    3021: +   x - a vector to mimic
    3022: -   n - the number of vectors to obtain

    3024:     Output Parameters:
    3025: +   y - Fortran90 pointer to the array of vectors
    3026: -   ierr - error code

    3028:     Example of Usage:
    3029: .vb
    3030:     Vec x
    3031:     Vec, pointer :: y(:)
    3032:     ....
    3033:     call VecDuplicateVecsF90(x,2,y,ierr)
    3034:     call VecSet(y(2),alpha,ierr)
    3035:     call VecSet(y(2),alpha,ierr)
    3036:     ....
    3037:     call VecDestroyVecsF90(y,2,ierr)
    3038: .ve

    3040:     Notes:
    3041:     Not yet supported for all F90 compilers

    3043:     Use VecDestroyVecsF90() to free the space.

    3045:     Level: beginner

    3047: .seealso:  VecDestroyVecsF90(), VecDuplicateVecs()

    3049: M*/

    3051: /*MC
    3052:     VecRestoreArrayF90 - Restores a vector to a usable state after a call to
    3053:     VecGetArrayF90().

    3055:     Synopsis:
    3056:     VecRestoreArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

    3058:     Not collective

    3060:     Input Parameters:
    3061: +   x - vector
    3062: -   xx_v - the Fortran90 pointer to the array

    3064:     Output Parameter:
    3065: .   ierr - error code

    3067:     Example of Usage:
    3068: .vb
    3069:     PetscScalar, pointer :: xx_v(:)
    3070:     ....
    3071:     call VecGetArrayF90(x,xx_v,ierr)
    3072:     a = xx_v(3)
    3073:     call VecRestoreArrayF90(x,xx_v,ierr)
    3074: .ve
    3075:
    3076:     Notes:
    3077:     Not yet supported for all F90 compilers

    3079:     Level: beginner

    3081: .seealso:  VecGetArrayF90(), VecGetArray(), VecRestoreArray()

    3083: M*/

    3085: /*MC
    3086:     VecDestroyVecsF90 - Frees a block of vectors obtained with VecDuplicateVecsF90().

    3088:     Synopsis:
    3089:     VecDestroyVecsF90({Vec, pointer :: x(:)},integer n,integer ierr)

    3091:     Input Parameters:
    3092: +   x - pointer to array of vector pointers
    3093: -   n - the number of vectors previously obtained

    3095:     Output Parameter:
    3096: .   ierr - error code

    3098:     Notes:
    3099:     Not yet supported for all F90 compilers

    3101:     Level: beginner

    3103: .seealso:  VecDestroyVecs(), VecDuplicateVecsF90()

    3105: M*/

    3107: /*MC
    3108:     VecGetArrayF90 - Accesses a vector array from Fortran90. For default PETSc
    3109:     vectors, VecGetArrayF90() returns a pointer to the local data array. Otherwise,
    3110:     this routine is implementation dependent. You MUST call VecRestoreArrayF90()
    3111:     when you no longer need access to the array.

    3113:     Synopsis:
    3114:     VecGetArrayF90(Vec x,{Scalar, pointer :: xx_v(:)},integer ierr)

    3116:     Not Collective

    3118:     Input Parameter:
    3119: .   x - vector

    3121:     Output Parameters:
    3122: +   xx_v - the Fortran90 pointer to the array
    3123: -   ierr - error code

    3125:     Example of Usage:
    3126: .vb
    3127:     PetscScalar, pointer :: xx_v(:)
    3128:     ....
    3129:     call VecGetArrayF90(x,xx_v,ierr)
    3130:     a = xx_v(3)
    3131:     call VecRestoreArrayF90(x,xx_v,ierr)
    3132: .ve

    3134:     Notes:
    3135:     Not yet supported for all F90 compilers

    3137:     Level: beginner

    3139: .seealso:  VecRestoreArrayF90(), VecGetArray(), VecRestoreArray()

    3141: M*/

    3145: /*@C
    3146:   VecLoadIntoVector - Loads a vector that has been stored in binary format
    3147:   with VecView().

    3149:   Collective on PetscViewer

    3151:   Input Parameters:
    3152: + viewer - binary file viewer, obtained from PetscViewerBinaryOpen()
    3153: - vec - vector to contain files values (must be of correct length)

    3155:   Level: intermediate

    3157:   Notes:
    3158:   The input file must contain the full global vector, as
    3159:   written by the routine VecView().

    3161:   Use VecLoad() to create the vector as the values are read in

    3163:   Notes for advanced users:
    3164:   Most users should not need to know the details of the binary storage
    3165:   format, since VecLoad() and VecView() completely hide these details.
    3166:   But for anyone who's interested, the standard binary matrix storage
    3167:   format is
    3168: .vb
    3169:      int    VEC_FILE_COOKIE
    3170:      int    number of rows
    3171:      PetscScalar *values of all nonzeros
    3172: .ve

    3174:    Note for Cray users, the int's stored in the binary file are 32 bit
    3175: integers; not 64 as they are represented in the memory, so if you
    3176: write your own routines to read/write these binary files from the Cray
    3177: you need to adjust the integer sizes that you read in, see
    3178: PetscBinaryRead() and PetscBinaryWrite() to see how this may be
    3179: done.

    3181:    In addition, PETSc automatically does the byte swapping for
    3182: machines that store the bytes reversed, e.g.  DEC alpha, freebsd,
    3183: linux, Windows and the paragon; thus if you write your own binary
    3184: read/write routines you have to swap the bytes; see PetscBinaryRead()
    3185: and PetscBinaryWrite() to see how this may be done.

    3187:    Concepts: vector^loading from file

    3189: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
    3190: @*/
    3191: PetscErrorCode PETSCVEC_DLLEXPORT VecLoadIntoVector(PetscViewer viewer,Vec vec)
    3192: {

    3199:   if (!vec->ops->loadintovector) {
    3200:     SETERRQ(PETSC_ERR_SUP,"Vector does not support load");
    3201:   }
    3202:   (*vec->ops->loadintovector)(viewer,vec);
    3203:   PetscObjectStateIncrease((PetscObject)vec);
    3204:   return(0);
    3205: }

    3209: /*@
    3210:    VecReciprocal - Replaces each component of a vector by its reciprocal.

    3212:    Collective on Vec

    3214:    Input Parameter:
    3215: .  v - the vector

    3217:    Output Parameter:
    3218: .  v - the vector reciprocal

    3220:    Level: intermediate

    3222:    Concepts: vector^reciprocal

    3224: @*/
    3225: PetscErrorCode PETSCVEC_DLLEXPORT VecReciprocal(Vec vec)
    3226: {

    3232:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
    3233:   if (!vec->ops->reciprocal) {
    3234:     SETERRQ(PETSC_ERR_SUP,"Vector does not support reciprocal operation");
    3235:   }
    3236:   (*vec->ops->reciprocal)(vec);
    3237:   PetscObjectStateIncrease((PetscObject)vec);
    3238:   return(0);
    3239: }

    3243: PetscErrorCode PETSCVEC_DLLEXPORT VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
    3244: {
    3247:   /* save the native version of the viewer */
    3248:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
    3249:     vec->ops->viewnative = vec->ops->view;
    3250:   }
    3251:   (((void(**)(void))vec->ops)[(int)op]) = f;
    3252:   return(0);
    3253: }

    3257: /*@
    3258:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
    3259:    used during the assembly process to store values that belong to
    3260:    other processors.

    3262:    Collective on Vec

    3264:    Input Parameters:
    3265: +  vec   - the vector
    3266: .  size  - the initial size of the stash.
    3267: -  bsize - the initial size of the block-stash(if used).

    3269:    Options Database Keys:
    3270: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
    3271: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

    3273:    Level: intermediate

    3275:    Notes:
    3276:      The block-stash is used for values set with VecSetValuesBlocked() while
    3277:      the stash is used for values set with VecSetValues()

    3279:      Run with the option -log_info and look for output of the form
    3280:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
    3281:      to determine the appropriate value, MM, to use for size and
    3282:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
    3283:      to determine the value, BMM to use for bsize

    3285:    Concepts: vector^stash
    3286:    Concepts: stash^vector

    3288: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()

    3290: @*/
    3291: PetscErrorCode PETSCVEC_DLLEXPORT VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
    3292: {

    3297:   VecStashSetInitialSize_Private(&vec->stash,size);
    3298:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
    3299:   return(0);
    3300: }

    3304: /*@
    3305:    VecStashView - Prints the entries in the vector stash and block stash.

    3307:    Collective on Vec

    3309:    Input Parameters:
    3310: +  vec   - the vector
    3311: -  viewer - the viewer

    3313:    Level: advanced

    3315:    Concepts: vector^stash
    3316:    Concepts: stash^vector

    3318: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()

    3320: @*/
    3321: PetscErrorCode PETSCVEC_DLLEXPORT VecStashView(Vec v,PetscViewer viewer)
    3322: {
    3324:   PetscMPIInt    rank;
    3325:   PetscInt       i,j;
    3326:   PetscTruth     match;
    3327:   VecStash       *s;
    3328:   PetscScalar    val;

    3335:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&match);
    3336:   if (!match) SETERRQ1(PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
    3337:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
    3338:   MPI_Comm_rank(v->comm,&rank);
    3339:   s = &v->bstash;

    3341:   /* print block stash */
    3342:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
    3343:   for (i=0; i<s->n; i++) {
    3344:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
    3345:     for (j=0; j<s->bs; j++) {
    3346:       val = s->array[i*s->bs+j];
    3347: #if defined(PETSC_USE_COMPLEX)
    3348:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
    3349: #else
    3350:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
    3351: #endif
    3352:     }
    3353:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
    3354:   }
    3355:   PetscViewerFlush(viewer);

    3357:   s = &v->stash;

    3359:   /* print basic stash */
    3360:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
    3361:   for (i=0; i<s->n; i++) {
    3362:     val = s->array[i];
    3363: #if defined(PETSC_USE_COMPLEX)
    3364:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
    3365: #else
    3366:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
    3367: #endif
    3368:   }
    3369:   PetscViewerFlush(viewer);

    3371:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
    3372:   return(0);
    3373: }

    3377: /*@C
    3378:    VecGetArray2d - Returns a pointer to a 2d contiguous array that contains this
    3379:    processor's portion of the vector data.  You MUST call VecRestoreArray2d()
    3380:    when you no longer need access to the array.

    3382:    Not Collective

    3384:    Input Parameter:
    3385: +  x - the vector
    3386: .  m - first dimension of two dimensional array
    3387: .  n - second dimension of two dimensional array
    3388: .  mstart - first index you will use in first coordinate direction (often 0)
    3389: -  nstart - first index in the second coordinate direction (often 0)

    3391:    Output Parameter:
    3392: .  a - location to put pointer to the array

    3394:    Level: beginner

    3396:   Notes:
    3397:    For a vector obtained from DACreateLocalVector() mstart and nstart are likely
    3398:    obtained from the corner indices obtained from DAGetGhostCorners() while for
    3399:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
    3400:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray2d().
    3401:
    3402:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3404:    Concepts: vector^accessing local values as 2d array

    3406: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
    3407:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
    3408:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
    3409: @*/
    3410: PetscErrorCode PETSCVEC_DLLEXPORT VecGetArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
    3411: {
    3413:   PetscInt       i,N;
    3414:   PetscScalar    *aa;

    3420:   VecGetLocalSize(x,&N);
    3421:   if (m*n != N) SETERRQ3(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 2d array dimensions %D by %D",N,m,n);
    3422:   VecGetArray(x,&aa);

    3424:   PetscMalloc(m*sizeof(PetscScalar*),a);
    3425:   for (i=0; i<m; i++) (*a)[i] = aa + i*n - nstart;
    3426:   *a -= mstart;
    3427:   return(0);
    3428: }

    3432: /*@C
    3433:    VecRestoreArray2d - Restores a vector after VecGetArray2d() has been called.

    3435:    Not Collective

    3437:    Input Parameters:
    3438: +  x - the vector
    3439: .  m - first dimension of two dimensional array
    3440: .  n - second dimension of the two dimensional array
    3441: .  mstart - first index you will use in first coordinate direction (often 0)
    3442: .  nstart - first index in the second coordinate direction (often 0)
    3443: -  a - location of pointer to array obtained from VecGetArray2d()

    3445:    Level: beginner

    3447:    Notes:
    3448:    For regular PETSc vectors this routine does not involve any copies. For
    3449:    any special vectors that do not store local vector data in a contiguous
    3450:    array, this routine will copy the data back into the underlying
    3451:    vector data structure from the array obtained with VecGetArray().

    3453:    This routine actually zeros out the a pointer.

    3455: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
    3456:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
    3457:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
    3458: @*/
    3459: PetscErrorCode PETSCVEC_DLLEXPORT VecRestoreArray2d(Vec x,PetscInt m,PetscInt n,PetscInt mstart,PetscInt nstart,PetscScalar **a[])
    3460: {
    3462:   void           *dummy;

    3468:   dummy = (void*)(*a + mstart);
    3469:   PetscFree(dummy);
    3470:   VecRestoreArray(x,PETSC_NULL);
    3471:   return(0);
    3472: }

    3476: /*@C
    3477:    VecGetArray1d - Returns a pointer to a 1d contiguous array that contains this
    3478:    processor's portion of the vector data.  You MUST call VecRestoreArray1d()
    3479:    when you no longer need access to the array.

    3481:    Not Collective

    3483:    Input Parameter:
    3484: +  x - the vector
    3485: .  m - first dimension of two dimensional array
    3486: -  mstart - first index you will use in first coordinate direction (often 0)

    3488:    Output Parameter:
    3489: .  a - location to put pointer to the array

    3491:    Level: beginner

    3493:   Notes:
    3494:    For a vector obtained from DACreateLocalVector() mstart are likely
    3495:    obtained from the corner indices obtained from DAGetGhostCorners() while for
    3496:    DACreateGlobalVector() they are the corner indices from DAGetCorners().
    3497:
    3498:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3500: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
    3501:           VecRestoreArray2d(), DAVecGetArray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
    3502:           VecGetArray2d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
    3503: @*/
    3504: PetscErrorCode PETSCVEC_DLLEXPORT VecGetArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
    3505: {
    3507:   PetscInt       N;

    3513:   VecGetLocalSize(x,&N);
    3514:   if (m != N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Local array size %D does not match 1d array dimensions %D",N,m);
    3515:   VecGetArray(x,a);
    3516:   *a  -= mstart;
    3517:   return(0);
    3518: }

    3522: /*@C
    3523:    VecRestoreArray1d - Restores a vector after VecGetArray1d() has been called.

    3525:    Not Collective

    3527:    Input Parameters:
    3528: +  x - the vector
    3529: .  m - first dimension of two dimensional array
    3530: .  mstart - first index you will use in first coordinate direction (often 0)
    3531: -  a - location of pointer to array obtained from VecGetArray21()

    3533:    Level: beginner

    3535:    Notes:
    3536:    For regular PETSc vectors this routine does not involve any copies. For
    3537:    any special vectors that do not store local vector data in a contiguous
    3538:    array, this routine will copy the data back into the underlying
    3539:    vector data structure from the array obtained with VecGetArray1d().

    3541:    This routine actually zeros out the a pointer.

    3543:    Concepts: vector^accessing local values as 1d array

    3545: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
    3546:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
    3547:           VecGetArray1d(), VecRestoreArray2d(), VecGetArray4d(), VecRestoreArray4d()
    3548: @*/
    3549: PetscErrorCode PETSCVEC_DLLEXPORT VecRestoreArray1d(Vec x,PetscInt m,PetscInt mstart,PetscScalar *a[])
    3550: {

    3556:   VecRestoreArray(x,PETSC_NULL);
    3557:   return(0);
    3558: }

    3562: /*@
    3563:    VecConjugate - Conjugates a vector.

    3565:    Collective on Vec

    3567:    Input Parameters:
    3568: .  x - the vector

    3570:    Level: intermediate

    3572:    Concepts: vector^conjugate

    3574: @*/
    3575: PetscErrorCode PETSCVEC_DLLEXPORT VecConjugate(Vec x)
    3576: {
    3577: #ifdef PETSC_USE_COMPLEX

    3583:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
    3584:   (*x->ops->conjugate)(x);
    3585:   /* we need to copy norms here */
    3586:   PetscObjectStateIncrease((PetscObject)x);
    3587:   return(0);
    3588: #else
    3589:   return(0);
    3590: #endif
    3591: }

    3595: /*@C
    3596:    VecGetArray3d - Returns a pointer to a 3d contiguous array that contains this
    3597:    processor's portion of the vector data.  You MUST call VecRestoreArray3d()
    3598:    when you no longer need access to the array.

    3600:    Not Collective

    3602:    Input Parameter:
    3603: +  x - the vector
    3604: .  m - first dimension of three dimensional array
    3605: .  n - second dimension of three dimensional array
    3606: .  p - third dimension of three dimensional array
    3607: .  mstart - first index you will use in first coordinate direction (often 0)
    3608: .  nstart - first index in the second coordinate direction (often 0)
    3609: -  pstart - first index in the third coordinate direction (often 0)

    3611:    Output Parameter:
    3612: .  a - location to put pointer to the array

    3614:    Level: beginner

    3616:   Notes:
    3617:    For a vector obtained from DACreateLocalVector() mstart, nstart, and pstart are likely
    3618:    obtained from the corner indices obtained from DAGetGhostCorners() while for
    3619:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
    3620:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray3d().
    3621:
    3622:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3624:    Concepts: vector^accessing local values as 3d array

    3626: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
    3627:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
    3628:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
    3629: @*/
    3630: PetscErrorCode PETSCVEC_DLLEXPORT VecGetArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
    3631: {
    3633:   PetscInt       i,N,j;
    3634:   PetscScalar    *aa,**b;

    3640:   VecGetLocalSize(x,&N);
    3641:   if (m*n*p != N) SETERRQ4(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 3d array dimensions %D by %D by %D",N,m,n,p);
    3642:   VecGetArray(x,&aa);

    3644:   PetscMalloc(m*sizeof(PetscScalar**)+m*n*sizeof(PetscScalar*),a);
    3645:   b    = (PetscScalar **)((*a) + m);
    3646:   for (i=0; i<m; i++)   (*a)[i] = b + i*n - nstart;
    3647:   for (i=0; i<m; i++) {
    3648:     for (j=0; j<n; j++) {
    3649:       b[i*n+j] = aa + i*n*p + j*p - pstart;
    3650:     }
    3651:   }
    3652:   *a -= mstart;
    3653:   return(0);
    3654: }

    3658: /*@C
    3659:    VecRestoreArray3d - Restores a vector after VecGetArray3d() has been called.

    3661:    Not Collective

    3663:    Input Parameters:
    3664: +  x - the vector
    3665: .  m - first dimension of three dimensional array
    3666: .  n - second dimension of the three dimensional array
    3667: .  p - third dimension of the three dimensional array
    3668: .  mstart - first index you will use in first coordinate direction (often 0)
    3669: .  nstart - first index in the second coordinate direction (often 0)
    3670: .  pstart - first index in the third coordinate direction (often 0)
    3671: -  a - location of pointer to array obtained from VecGetArray3d()

    3673:    Level: beginner

    3675:    Notes:
    3676:    For regular PETSc vectors this routine does not involve any copies. For
    3677:    any special vectors that do not store local vector data in a contiguous
    3678:    array, this routine will copy the data back into the underlying
    3679:    vector data structure from the array obtained with VecGetArray().

    3681:    This routine actually zeros out the a pointer.

    3683: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
    3684:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
    3685:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
    3686: @*/
    3687: PetscErrorCode PETSCVEC_DLLEXPORT VecRestoreArray3d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscScalar ***a[])
    3688: {
    3690:   void           *dummy;

    3696:   dummy = (void*)(*a + mstart);
    3697:   PetscFree(dummy);
    3698:   VecRestoreArray(x,PETSC_NULL);
    3699:   return(0);
    3700: }

    3704: /*@C
    3705:    VecGetArray4d - Returns a pointer to a 4d contiguous array that contains this
    3706:    processor's portion of the vector data.  You MUST call VecRestoreArray4d()
    3707:    when you no longer need access to the array.

    3709:    Not Collective

    3711:    Input Parameter:
    3712: +  x - the vector
    3713: .  m - first dimension of four dimensional array
    3714: .  n - second dimension of four dimensional array
    3715: .  p - third dimension of four dimensional array
    3716: .  q - fourth dimension of four dimensional array
    3717: .  mstart - first index you will use in first coordinate direction (often 0)
    3718: .  nstart - first index in the second coordinate direction (often 0)
    3719: .  pstart - first index in the third coordinate direction (often 0)
    3720: -  qstart - first index in the fourth coordinate direction (often 0)

    3722:    Output Parameter:
    3723: .  a - location to put pointer to the array

    3725:    Level: beginner

    3727:   Notes:
    3728:    For a vector obtained from DACreateLocalVector() mstart, nstart, and pstart are likely
    3729:    obtained from the corner indices obtained from DAGetGhostCorners() while for
    3730:    DACreateGlobalVector() they are the corner indices from DAGetCorners(). In both cases
    3731:    the arguments from DAGet[Ghost}Corners() are reversed in the call to VecGetArray3d().
    3732:
    3733:    For standard PETSc vectors this is an inexpensive call; it does not copy the vector values.

    3735:    Concepts: vector^accessing local values as 3d array

    3737: .seealso: VecGetArray(), VecRestoreArray(), VecGetArrays(), VecGetArrayF90(), VecPlaceArray(),
    3738:           VecRestoreArray2d(), DAVecGetarray(), DAVecRestoreArray(), VecGetArray3d(), VecRestoreArray3d(),
    3739:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d()
    3740: @*/
    3741: PetscErrorCode PETSCVEC_DLLEXPORT VecGetArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
    3742: {
    3744:   PetscInt       i,N,j,k;
    3745:   PetscScalar    *aa,***b,**c;

    3751:   VecGetLocalSize(x,&N);
    3752:   if (m*n*p*q != N) SETERRQ5(PETSC_ERR_ARG_INCOMP,"Local array size %D does not match 4d array dimensions %D by %D by %D by %D",N,m,n,p,q);
    3753:   VecGetArray(x,&aa);

    3755:   PetscMalloc(m*sizeof(PetscScalar***)+m*n*sizeof(PetscScalar**)+m*n*p*sizeof(PetscScalar*),a);
    3756:   b    = (PetscScalar ***)((*a) + m);
    3757:   c    = (PetscScalar **)(b + m*n);
    3758:   for (i=0; i<m; i++)   (*a)[i] = b + i*n - nstart;
    3759:   for (i=0; i<m; i++) {
    3760:     for (j=0; j<n; j++) {
    3761:       b[i*n+j] = c + i*n*p + j*p - pstart;
    3762:     }
    3763:   }
    3764:   for (i=0; i<m; i++) {
    3765:     for (j=0; j<n; j++) {
    3766:       for (k=0; k<p; k++) {
    3767:         c[i*n*p+j*p+k] = aa + i*n*p*q + j*p*q + k*q - qstart;
    3768:       }
    3769:     }
    3770:   }
    3771:   *a -= mstart;
    3772:   return(0);
    3773: }

    3777: /*@C
    3778:    VecRestoreArray4d - Restores a vector after VecGetArray3d() has been called.

    3780:    Not Collective

    3782:    Input Parameters:
    3783: +  x - the vector
    3784: .  m - first dimension of four dimensional array
    3785: .  n - second dimension of the four dimensional array
    3786: .  p - third dimension of the four dimensional array
    3787: .  q - fourth dimension of the four dimensional array
    3788: .  mstart - first index you will use in first coordinate direction (often 0)
    3789: .  nstart - first index in the second coordinate direction (often 0)
    3790: .  pstart - first index in the third coordinate direction (often 0)
    3791: .  qstart - first index in the fourth coordinate direction (often 0)
    3792: -  a - location of pointer to array obtained from VecGetArray4d()

    3794:    Level: beginner

    3796:    Notes:
    3797:    For regular PETSc vectors this routine does not involve any copies. For
    3798:    any special vectors that do not store local vector data in a contiguous
    3799:    array, this routine will copy the data back into the underlying
    3800:    vector data structure from the array obtained with VecGetArray().

    3802:    This routine actually zeros out the a pointer.

    3804: .seealso: VecGetArray(), VecRestoreArray(), VecRestoreArrays(), VecRestoreArrayF90(), VecPlaceArray(),
    3805:           VecGetArray2d(), VecGetArray3d(), VecRestoreArray3d(), DAVecGetArray(), DAVecRestoreArray()
    3806:           VecGetArray1d(), VecRestoreArray1d(), VecGetArray4d(), VecRestoreArray4d(), VecGet
    3807: @*/
    3808: PetscErrorCode PETSCVEC_DLLEXPORT VecRestoreArray4d(Vec x,PetscInt m,PetscInt n,PetscInt p,PetscInt q,PetscInt mstart,PetscInt nstart,PetscInt pstart,PetscInt qstart,PetscScalar ****a[])
    3809: {
    3811:   void           *dummy;

    3817:   dummy = (void*)(*a + mstart);
    3818:   PetscFree(dummy);
    3819:   VecRestoreArray(x,PETSC_NULL);
    3820:   return(0);
    3821: }

    3823: EXTERN PetscErrorCode VecStashGetInfo_Private(VecStash*,PetscInt*,PetscInt*);
    3826: /*@
    3827:    VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
    3828:        to be communicated to other processors during the VecAssemblyBegin/End() process

    3830:     Not collective

    3832:    Input Parameter:
    3833: .   vec - the vector

    3835:    Output Parameters:
    3836: +   nstash   - the size of the stash
    3837: .   reallocs - the number of additional mallocs incurred.
    3838: .   bnstash   - the size of the block stash
    3839: -   breallocs - the number of additional mallocs incurred.in the block stash
    3840:
    3841:    Level: advanced

    3843: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
    3844:
    3845: @*/
    3846: PetscErrorCode PETSCVEC_DLLEXPORT VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *brealloc)
    3847: {
    3850:   VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
    3851:   VecStashGetInfo_Private(&vec->bstash,nstash,reallocs);
    3852:   return(0);
    3853: }