Actual source code: vector.c

petsc-3.12.3 2020-01-03
Report Typos and Errors
  1: /*
  2:      Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
  3:    These are the vector functions the user calls.
  4: */
  5:  #include <petsc/private/vecimpl.h>

  7: /* Logging support */
  8: PetscClassId  VEC_CLASSID;
  9: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_Dot, VEC_MDot, VEC_TDot;
 10: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
 11: PetscLogEvent VEC_MTDot, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
 12: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load;
 13: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication,VEC_ReduceBegin,VEC_ReduceEnd,VEC_Ops;
 14: PetscLogEvent VEC_DotNorm2, VEC_AXPBYPCZ;
 15: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
 16: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
 17: PetscLogEvent VEC_CUDACopyFromGPUSome, VEC_CUDACopyToGPUSome;

 19: /*@
 20:    VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
 21:        to be communicated to other processors during the VecAssemblyBegin/End() process

 23:     Not collective

 25:    Input Parameter:
 26: .   vec - the vector

 28:    Output Parameters:
 29: +   nstash   - the size of the stash
 30: .   reallocs - the number of additional mallocs incurred.
 31: .   bnstash   - the size of the block stash
 32: -   breallocs - the number of additional mallocs incurred.in the block stash

 34:    Level: advanced

 36: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()

 38: @*/
 39: PetscErrorCode  VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
 40: {

 44:   VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
 45:   VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
 46:   return(0);
 47: }

 49: /*@
 50:    VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
 51:    by the routine VecSetValuesLocal() to allow users to insert vector entries
 52:    using a local (per-processor) numbering.

 54:    Logically Collective on Vec

 56:    Input Parameters:
 57: +  x - vector
 58: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

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

 63:    Level: intermediate

 65: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
 66:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
 67: @*/
 68: PetscErrorCode  VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
 69: {


 76:   if (x->ops->setlocaltoglobalmapping) {
 77:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
 78:   } else {
 79:     PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
 80:   }
 81:   return(0);
 82: }

 84: /*@
 85:    VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()

 87:    Not Collective

 89:    Input Parameter:
 90: .  X - the vector

 92:    Output Parameter:
 93: .  mapping - the mapping

 95:    Level: advanced


 98: .seealso:  VecSetValuesLocal()
 99: @*/
100: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
101: {
106:   *mapping = X->map->mapping;
107:   return(0);
108: }

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

114:    Collective on Vec

116:    Input Parameter:
117: .  vec - the vector

119:    Level: beginner

121: .seealso: VecAssemblyEnd(), VecSetValues()
122: @*/
123: PetscErrorCode  VecAssemblyBegin(Vec vec)
124: {

130:   VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
131:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
132:   if (vec->ops->assemblybegin) {
133:     (*vec->ops->assemblybegin)(vec);
134:   }
135:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
136:   PetscObjectStateIncrease((PetscObject)vec);
137:   return(0);
138: }

140: /*@
141:    VecAssemblyEnd - Completes assembling the vector.  This routine should
142:    be called after VecAssemblyBegin().

144:    Collective on Vec

146:    Input Parameter:
147: .  vec - the vector

149:    Options Database Keys:
150: +  -vec_view - Prints vector in ASCII format
151: .  -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
152: .  -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
153: .  -vec_view draw - Activates vector viewing using drawing tools
154: .  -display <name> - Sets display name (default is host)
155: .  -draw_pause <sec> - Sets number of seconds to pause after display
156: -  -vec_view socket - Activates vector viewing using a socket

158:    Level: beginner

160: .seealso: VecAssemblyBegin(), VecSetValues()
161: @*/
162: PetscErrorCode  VecAssemblyEnd(Vec vec)
163: {

168:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
170:   if (vec->ops->assemblyend) {
171:     (*vec->ops->assemblyend)(vec);
172:   }
173:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
174:   VecViewFromOptions(vec,NULL,"-vec_view");
175:   return(0);
176: }

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

181:    Logically Collective on Vec

183:    Input Parameters:
184: .  x, y  - the vectors

186:    Output Parameter:
187: .  w - the result

189:    Level: advanced

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

195: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
196: @*/
197: PetscErrorCode  VecPointwiseMax(Vec w,Vec x,Vec y)
198: {

210:   VecCheckSameSize(w,1,x,2);
211:   VecCheckSameSize(w,1,y,3);
212:   (*w->ops->pointwisemax)(w,x,y);
213:   PetscObjectStateIncrease((PetscObject)w);
214:   return(0);
215: }


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

221:    Logically Collective on Vec

223:    Input Parameters:
224: .  x, y  - the vectors

226:    Output Parameter:
227: .  w - the result

229:    Level: advanced

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

235: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
236: @*/
237: PetscErrorCode  VecPointwiseMin(Vec w,Vec x,Vec y)
238: {

250:   VecCheckSameSize(w,1,x,2);
251:   VecCheckSameSize(w,1,y,3);
252:   (*w->ops->pointwisemin)(w,x,y);
253:   PetscObjectStateIncrease((PetscObject)w);
254:   return(0);
255: }

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

260:    Logically Collective on Vec

262:    Input Parameters:
263: .  x, y  - the vectors

265:    Output Parameter:
266: .  w - the result

268:    Level: advanced

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

273: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
274: @*/
275: PetscErrorCode  VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
276: {

288:   VecCheckSameSize(w,1,x,2);
289:   VecCheckSameSize(w,1,y,3);
290:   (*w->ops->pointwisemaxabs)(w,x,y);
291:   PetscObjectStateIncrease((PetscObject)w);
292:   return(0);
293: }

295: /*@
296:    VecPointwiseDivide - Computes the componentwise division w = x/y.

298:    Logically Collective on Vec

300:    Input Parameters:
301: .  x, y  - the vectors

303:    Output Parameter:
304: .  w - the result

306:    Level: advanced

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

311: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
312: @*/
313: PetscErrorCode  VecPointwiseDivide(Vec w,Vec x,Vec y)
314: {

326:   VecCheckSameSize(w,1,x,2);
327:   VecCheckSameSize(w,1,y,3);
328:   (*w->ops->pointwisedivide)(w,x,y);
329:   PetscObjectStateIncrease((PetscObject)w);
330:   return(0);
331: }


334: /*@
335:    VecDuplicate - Creates a new vector of the same type as an existing vector.

337:    Collective on Vec

339:    Input Parameters:
340: .  v - a vector to mimic

342:    Output Parameter:
343: .  newv - location to put new vector

345:    Notes:
346:    VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
347:    for the new vector.  Use VecCopy() to copy a vector.

349:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
350:    vectors.

352:    Level: beginner

354: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
355: @*/
356: PetscErrorCode  VecDuplicate(Vec v,Vec *newv)
357: {

364:   (*v->ops->duplicate)(v,newv);
365:   PetscObjectStateIncrease((PetscObject)*newv);
366:   return(0);
367: }

369: /*@
370:    VecDestroy - Destroys a vector.

372:    Collective on Vec

374:    Input Parameters:
375: .  v  - the vector

377:    Level: beginner

379: .seealso: VecDuplicate(), VecDestroyVecs()
380: @*/
381: PetscErrorCode  VecDestroy(Vec *v)
382: {

386:   if (!*v) return(0);
388:   if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}

390:   PetscObjectSAWsViewOff((PetscObject)*v);
391:   /* destroy the internal part */
392:   if ((*v)->ops->destroy) {
393:     (*(*v)->ops->destroy)(*v);
394:   }
395:   /* destroy the external/common part */
396:   PetscLayoutDestroy(&(*v)->map);
397:   PetscHeaderDestroy(v);
398:   return(0);
399: }

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

404:    Collective on Vec

406:    Input Parameters:
407: +  m - the number of vectors to obtain
408: -  v - a vector to mimic

410:    Output Parameter:
411: .  V - location to put pointer to array of vectors

413:    Notes:
414:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
415:    vector.

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

422:    Level: intermediate

424: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
425: @*/
426: PetscErrorCode  VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
427: {

434:   (*v->ops->duplicatevecs)(v,m,V);
435:   return(0);
436: }

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

441:    Collective on Vec

443:    Input Parameters:
444: +  vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
445: -  m - the number of vectors previously obtained, if zero no vectors are destroyed

447:    Fortran Note:
448:    The Fortran interface is slightly different from that given below.
449:    See the Fortran chapter of the users manual

451:    Level: intermediate

453: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
454: @*/
455: PetscErrorCode  VecDestroyVecs(PetscInt m,Vec *vv[])
456: {

461:   if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
462:   if (!m || !*vv) {*vv  = NULL; return(0);}
465:   (*(**vv)->ops->destroyvecs)(m,*vv);
466:   *vv  = NULL;
467:   return(0);
468: }

470: /*@C
471:    VecView - Views a vector object.

473:    Collective on Vec

475:    Input Parameters:
476: +  vec - the vector
477: -  viewer - an optional visualization context

479:    Notes:
480:    The available visualization contexts include
481: +     PETSC_VIEWER_STDOUT_SELF - for sequential vectors
482: .     PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
483: -     PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm

485:    You can change the format the vector is printed using the
486:    option PetscViewerPushFormat().

488:    The user can open alternative viewers with
489: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
490: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
491:          specified file; corresponding input uses VecLoad()
492: .    PetscViewerDrawOpen() - Outputs vector to an X window display
493: .    PetscViewerSocketOpen() - Outputs vector to Socket viewer
494: -    PetscViewerHDF5Open() - Outputs vector to HDF5 file viewer 

496:    The user can call PetscViewerPushFormat() to specify the output
497:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
498:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
499: +    PETSC_VIEWER_DEFAULT - default, prints vector contents
500: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
501: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
502: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
503:          format common among all vector types

505:    Notes:
506:     You can pass any number of vector objects, or other PETSc objects to the same viewer.

508:    Notes for binary viewer: 
509:      If you pass multiple vectors to a binary viewer you can read them back in in the same order
510:      with VecLoad().

512:      If the blocksize of the vector is greater than one then you must provide a unique prefix to
513:      the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
514:      vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
515:      information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
516:      filename. If you copy the binary file, make sure you copy the associated .info file with it.
517:      
518:      See the manual page for VecLoad() on the exact format the binary viewer stores
519:      the values in the file.


522:    Notes for HDF5 Viewer: 
523:      The name of the Vec (given with PetscObjectSetName() is the name that is used
524:      for the object in the HDF5 file. If you wish to store the same Vec into multiple
525:      datasets in the same file (typically with different values), you must change its
526:      name each time before calling the VecView(). To load the same vector,
527:      the name of the Vec object passed to VecLoad() must be the same.
528:  
529:      If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
530:      If the function PetscViewerHDF5SetBaseDimension2()is called then even if the block size is one it will
531:      be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
532:      See also PetscViewerHDF5SetTimestep() which adds an additional complication to reading and writing Vecs
533:      with the HDF5 viewer.
534:      
535:    Level: beginner


538: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
539:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
540:           PetscRealView(), PetscScalarView(), PetscIntView(), PetscViewerHDF5SetTimestep()
541: @*/
542: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
543: {
544:   PetscErrorCode    ierr;
545:   PetscBool         iascii;
546:   PetscViewerFormat format;
547:   PetscMPIInt       size;

552:   if (!viewer) {
553:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
554:   }
556:   PetscViewerGetFormat(viewer,&format);
557:   MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
558:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);

560:   if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");

562:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
563:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
564:   if (iascii) {
565:     PetscInt rows,bs;

567:     PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
568:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
569:       PetscViewerASCIIPushTab(viewer);
570:       VecGetSize(vec,&rows);
571:       VecGetBlockSize(vec,&bs);
572:       if (bs != 1) {
573:         PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
574:       } else {
575:         PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
576:       }
577:       PetscViewerASCIIPopTab(viewer);
578:     }
579:   }
580:   VecLockReadPush(vec);
581:   if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
582:     (*vec->ops->viewnative)(vec,viewer);
583:   } else {
584:     (*vec->ops->view)(vec,viewer);
585:   }
586:   VecLockReadPop(vec);
587:   PetscLogEventEnd(VEC_View,vec,viewer,0,0);
588:   return(0);
589: }

591: #if defined(PETSC_USE_DEBUG)
592: #include <../src/sys/totalview/tv_data_display.h>
593: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
594: {
595:   const PetscScalar *values;
596:   char              type[32];
597:   PetscErrorCode    ierr;


600:   TV_add_row("Local rows", "int", &v->map->n);
601:   TV_add_row("Global rows", "int", &v->map->N);
602:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
603:   VecGetArrayRead((Vec)v,&values);
604:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
605:   TV_add_row("values",type, values);
606:   VecRestoreArrayRead((Vec)v,&values);
607:   return TV_format_OK;
608: }
609: #endif

611: /*@
612:    VecGetSize - Returns the global number of elements of the vector.

614:    Not Collective

616:    Input Parameter:
617: .  x - the vector

619:    Output Parameters:
620: .  size - the global length of the vector

622:    Level: beginner

624: .seealso: VecGetLocalSize()
625: @*/
626: PetscErrorCode  VecGetSize(Vec x,PetscInt *size)
627: {

634:   (*x->ops->getsize)(x,size);
635:   return(0);
636: }

638: /*@
639:    VecGetLocalSize - Returns the number of elements of the vector stored
640:    in local memory. This routine may be implementation dependent, so use
641:    with care.

643:    Not Collective

645:    Input Parameter:
646: .  x - the vector

648:    Output Parameter:
649: .  size - the length of the local piece of the vector

651:    Level: beginner

653: .seealso: VecGetSize()
654: @*/
655: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
656: {

663:   (*x->ops->getlocalsize)(x,size);
664:   return(0);
665: }

667: /*@C
668:    VecGetOwnershipRange - Returns the range of indices owned by
669:    this processor, assuming that the vectors are laid out with the
670:    first n1 elements on the first processor, next n2 elements on the
671:    second, etc.  For certain parallel layouts this range may not be
672:    well defined.

674:    Not Collective

676:    Input Parameter:
677: .  x - the vector

679:    Output Parameters:
680: +  low - the first local element, pass in NULL if not interested
681: -  high - one more than the last local element, pass in NULL if not interested

683:    Note:
684:    The high argument is one more than the last element stored locally.

686:    Fortran: PETSC_NULL_INTEGER should be used instead of NULL

688:    Level: beginner


691: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
692: @*/
693: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
694: {
700:   if (low)  *low  = x->map->rstart;
701:   if (high) *high = x->map->rend;
702:   return(0);
703: }

705: /*@C
706:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
707:    assuming that the vectors are laid out with the
708:    first n1 elements on the first processor, next n2 elements on the
709:    second, etc.  For certain parallel layouts this range may not be
710:    well defined.

712:    Not Collective

714:    Input Parameter:
715: .  x - the vector

717:    Output Parameters:
718: .  range - array of length size+1 with the start and end+1 for each process

720:    Note:
721:    The high argument is one more than the last element stored locally.

723:    Fortran: You must PASS in an array of length size+1

725:    Level: beginner


728: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
729: @*/
730: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
731: {

737:   PetscLayoutGetRanges(x->map,ranges);
738:   return(0);
739: }

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

744:    Collective on Vec

746:    Input Parameter:
747: +  x - the vector
748: .  op - the option
749: -  flag - turn the option on or off

751:    Supported Options:
752: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
753:           entries destined to be stored on a separate processor. This can be used
754:           to eliminate the global reduction in the VecAssemblyXXXX() if you know
755:           that you have only used VecSetValues() to set local elements
756: .     VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
757:           in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
758:           ignored.
759: -     VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
760:           entries will always be a subset (possibly equal) of the off-process entries set on the
761:           first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
762:           changed this flag afterwards. If this assembly is not such first assembly, then this
763:           assembly can reuse the communication pattern setup in that first assembly, thus avoiding
764:           a global reduction. Subsequent assemblies setting off-process values should use the same
765:           InsertMode as the first assembly.

767:    Developer Note:
768:    The InsertMode restriction could be removed by packing the stash messages out of place.

770:    Level: intermediate

772: @*/
773: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool flag)
774: {

780:   if (x->ops->setoption) {
781:     (*x->ops->setoption)(x,op,flag);
782:   }
783:   return(0);
784: }

786: /* Default routines for obtaining and releasing; */
787: /* may be used by any implementation */
788: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
789: {
791:   PetscInt       i;

796:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
797:   PetscMalloc1(m,V);
798:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
799:   return(0);
800: }

802: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
803: {
805:   PetscInt       i;

809:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
810:   PetscFree(v);
811:   return(0);
812: }

814: /*@
815:    VecResetArray - Resets a vector to use its default memory. Call this
816:    after the use of VecPlaceArray().

818:    Not Collective

820:    Input Parameters:
821: .  vec - the vector

823:    Level: developer

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

827: @*/
828: PetscErrorCode  VecResetArray(Vec vec)
829: {

835:   if (vec->ops->resetarray) {
836:     (*vec->ops->resetarray)(vec);
837:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
838:   PetscObjectStateIncrease((PetscObject)vec);
839:   return(0);
840: }

842: /*@C
843:   VecLoad - Loads a vector that has been stored in binary or HDF5 format
844:   with VecView().

846:   Collective on PetscViewer

848:   Input Parameters:
849: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
850:            some related function before a call to VecLoad().
851: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
852:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

854:    Level: intermediate

856:   Notes:
857:   Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
858:   before calling this.

860:   The input file must contain the full global vector, as
861:   written by the routine VecView().

863:   If the type or size of newvec is not set before a call to VecLoad, PETSc
864:   sets the type and the local and global sizes. If type and/or
865:   sizes are already set, then the same are used.

867:   If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
868:   the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
869:   vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
870:   information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
871:   filename. If you copy the binary file, make sure you copy the associated .info file with it.

873:   If using HDF5, you must assign the Vec the same name as was used in the Vec
874:   that was stored in the file using PetscObjectSetName(). Otherwise you will
875:   get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT".

877:   If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
878:   in loading the vector. Hence, for example, using Matlab notation h5create('vector.dat','/Test_Vec',[27 1]);
879:   will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
880:   vectors communicator and the rest of the processes having zero entries

882:   Notes for advanced users when using the binary viewer:
883:   Most users should not need to know the details of the binary storage
884:   format, since VecLoad() and VecView() completely hide these details.
885:   But for anyone who's interested, the standard binary vector storage
886:   format is
887: .vb
888:      PetscInt    VEC_FILE_CLASSID
889:      PetscInt    number of rows
890:      PetscScalar *values of all entries
891: .ve

893:    In addition, PETSc automatically uses byte swapping to work on all machines; the files
894:    are written ALWAYS using big-endian ordering. On small-endian machines the numbers
895:    are converted to the small-endian format when they are read in from the file.
896:    See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.

898: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
899: @*/
900: PetscErrorCode  VecLoad(Vec newvec, PetscViewer viewer)
901: {
902:   PetscErrorCode    ierr;
903:   PetscBool         isbinary,ishdf5,isadios,isadios2;
904:   PetscViewerFormat format;

909:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
910:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
911:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
912:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
913:   if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

915:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
916:   if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
917:     VecSetType(newvec, VECSTANDARD);
918:   }
919:   PetscViewerGetFormat(viewer,&format);
920:   if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
921:     (*newvec->ops->loadnative)(newvec,viewer);
922:   } else {
923:     (*newvec->ops->load)(newvec,viewer);
924:   }
925:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
926:   return(0);
927: }


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

933:    Logically Collective on Vec

935:    Input Parameter:
936: .  vec - the vector

938:    Output Parameter:
939: .  vec - the vector reciprocal

941:    Level: intermediate

943: .seealso: VecLog(), VecExp(), VecSqrtAbs()

945: @*/
946: PetscErrorCode  VecReciprocal(Vec vec)
947: {

953:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
954:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
955:   (*vec->ops->reciprocal)(vec);
956:   PetscObjectStateIncrease((PetscObject)vec);
957:   return(0);
958: }

960: /*@C
961:     VecSetOperation - Allows user to set a vector operation.

963:    Logically Collective on Vec

965:     Input Parameters:
966: +   vec - the vector
967: .   op - the name of the operation
968: -   f - the function that provides the operation.

970:    Level: advanced

972:     Usage:
973: $      PetscErrorCode userview(Vec,PetscViewer);
974: $      VecCreateMPI(comm,m,M,&x);
975: $      VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);

977:     Notes:
978:     See the file include/petscvec.h for a complete list of matrix
979:     operations, which all have the form VECOP_<OPERATION>, where
980:     <OPERATION> is the name (in all capital letters) of the
981:     user interface routine (e.g., VecView() -> VECOP_VIEW).

983:     This function is not currently available from Fortran.

985: .seealso: VecCreate(), MatShellSetOperation()
986: @*/
987: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
988: {
991:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
992:     vec->ops->viewnative = vec->ops->view;
993:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
994:     vec->ops->loadnative = vec->ops->load;
995:   }
996:   (((void(**)(void))vec->ops)[(int)op]) = f;
997:   return(0);
998: }


1001: /*@
1002:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1003:    used during the assembly process to store values that belong to
1004:    other processors.

1006:    Not Collective, different processes can have different size stashes

1008:    Input Parameters:
1009: +  vec   - the vector
1010: .  size  - the initial size of the stash.
1011: -  bsize - the initial size of the block-stash(if used).

1013:    Options Database Keys:
1014: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1015: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1017:    Level: intermediate

1019:    Notes:
1020:      The block-stash is used for values set with VecSetValuesBlocked() while
1021:      the stash is used for values set with VecSetValues()

1023:      Run with the option -info and look for output of the form
1024:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1025:      to determine the appropriate value, MM, to use for size and
1026:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1027:      to determine the value, BMM to use for bsize


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

1032: @*/
1033: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1034: {

1039:   VecStashSetInitialSize_Private(&vec->stash,size);
1040:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1041:   return(0);
1042: }

1044: /*@
1045:    VecConjugate - Conjugates a vector.

1047:    Logically Collective on Vec

1049:    Input Parameters:
1050: .  x - the vector

1052:    Level: intermediate

1054: @*/
1055: PetscErrorCode  VecConjugate(Vec x)
1056: {
1057: #if defined(PETSC_USE_COMPLEX)

1063:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1064:   (*x->ops->conjugate)(x);
1065:   /* we need to copy norms here */
1066:   PetscObjectStateIncrease((PetscObject)x);
1067:   return(0);
1068: #else
1069:   return(0);
1070: #endif
1071: }

1073: /*@
1074:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1076:    Logically Collective on Vec

1078:    Input Parameters:
1079: .  x, y  - the vectors

1081:    Output Parameter:
1082: .  w - the result

1084:    Level: advanced

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

1089: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1090: @*/
1091: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1092: {

1104:   VecCheckSameSize(w,1,x,2);
1105:   VecCheckSameSize(w,2,y,3);
1106:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1107:   (*w->ops->pointwisemult)(w,x,y);
1108:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1109:   PetscObjectStateIncrease((PetscObject)w);
1110:   return(0);
1111: }

1113: /*@
1114:    VecSetRandom - Sets all components of a vector to random numbers.

1116:    Logically Collective on Vec

1118:    Input Parameters:
1119: +  x  - the vector
1120: -  rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1121:           it will create one internally.

1123:    Output Parameter:
1124: .  x  - the vector

1126:    Example of Usage:
1127: .vb
1128:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1129:      VecSetRandom(x,rctx);
1130:      PetscRandomDestroy(rctx);
1131: .ve

1133:    Level: intermediate


1136: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1137: @*/
1138: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1139: {
1141:   PetscRandom    randObj = NULL;

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

1149:   if (!rctx) {
1150:     MPI_Comm comm;
1151:     PetscObjectGetComm((PetscObject)x,&comm);
1152:     PetscRandomCreate(comm,&randObj);
1153:     PetscRandomSetFromOptions(randObj);
1154:     rctx = randObj;
1155:   }

1157:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1158:   (*x->ops->setrandom)(x,rctx);
1159:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1161:   PetscRandomDestroy(&randObj);
1162:   PetscObjectStateIncrease((PetscObject)x);
1163:   return(0);
1164: }

1166: /*@
1167:   VecZeroEntries - puts a 0.0 in each element of a vector

1169:   Logically Collective on Vec

1171:   Input Parameter:
1172: . vec - The vector

1174:   Level: beginner

1176:   Developer Note: This routine does not need to exist since the exact functionality is obtained with
1177:      VecSet(vec,0);  I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1178:      like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1179:      this routine should not exist.

1181: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1182: @*/
1183: PetscErrorCode  VecZeroEntries(Vec vec)
1184: {

1188:   VecSet(vec,0);
1189:   return(0);
1190: }

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

1196:   Collective on Vec

1198:   Input Parameter:
1199: . vec - The vector

1201:   Level: intermediate

1203: .seealso: VecSetFromOptions(), VecSetType()
1204: */
1205: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1206: {
1207:   PetscBool      opt;
1208:   VecType        defaultType;
1209:   char           typeName[256];
1210:   PetscMPIInt    size;

1214:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1215:   else {
1216:     MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1217:     if (size > 1) defaultType = VECMPI;
1218:     else defaultType = VECSEQ;
1219:   }

1221:   VecRegisterAll();
1222:   PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1223:   if (opt) {
1224:     VecSetType(vec, typeName);
1225:   } else {
1226:     VecSetType(vec, defaultType);
1227:   }
1228:   return(0);
1229: }

1231: /*@
1232:   VecSetFromOptions - Configures the vector from the options database.

1234:   Collective on Vec

1236:   Input Parameter:
1237: . vec - The vector

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

1243:   Level: beginner


1246: .seealso: VecCreate(), VecSetOptionsPrefix()
1247: @*/
1248: PetscErrorCode  VecSetFromOptions(Vec vec)
1249: {


1255:   PetscObjectOptionsBegin((PetscObject)vec);
1256:   /* Handle vector type options */
1257:   VecSetTypeFromOptions_Private(PetscOptionsObject,vec);

1259:   /* Handle specific vector options */
1260:   if (vec->ops->setfromoptions) {
1261:     (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1262:   }

1264:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1265:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1266:   PetscOptionsEnd();
1267:   return(0);
1268: }

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

1273:   Collective on Vec

1275:   Input Parameters:
1276: + v - the vector
1277: . n - the local size (or PETSC_DECIDE to have it set)
1278: - N - the global size (or PETSC_DECIDE)

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

1284:   Level: intermediate

1286: .seealso: VecGetSize(), PetscSplitOwnership()
1287: @*/
1288: PetscErrorCode  VecSetSizes(Vec v, PetscInt n, PetscInt N)
1289: {

1295:   if (N >= 0 && n > N) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Local size %D cannot be larger than global size %D",n,N);
1296:   if ((v->map->n >= 0 || v->map->N >= 0) && (v->map->n != n || v->map->N != N)) SETERRQ4(PETSC_COMM_SELF,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->map->n,v->map->N);
1297:   v->map->n = n;
1298:   v->map->N = N;
1299:   if (v->ops->create) {
1300:     (*v->ops->create)(v);
1301:     v->ops->create = 0;
1302:   }
1303:   return(0);
1304: }

1306: /*@
1307:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1308:    and VecSetValuesBlockedLocal().

1310:    Logically Collective on Vec

1312:    Input Parameter:
1313: +  v - the vector
1314: -  bs - the blocksize

1316:    Notes:
1317:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1319:    Level: advanced

1321: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()

1323: @*/
1324: PetscErrorCode  VecSetBlockSize(Vec v,PetscInt bs)
1325: {

1330:   if (bs < 0 || bs == v->map->bs) return(0);
1332:   PetscLayoutSetBlockSize(v->map,bs);
1333:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1334:   return(0);
1335: }

1337: /*@
1338:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1339:    and VecSetValuesBlockedLocal().

1341:    Not Collective

1343:    Input Parameter:
1344: .  v - the vector

1346:    Output Parameter:
1347: .  bs - the blocksize

1349:    Notes:
1350:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1352:    Level: advanced

1354: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()


1357: @*/
1358: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1359: {

1365:   PetscLayoutGetBlockSize(v->map,bs);
1366:   return(0);
1367: }

1369: /*@C
1370:    VecSetOptionsPrefix - Sets the prefix used for searching for all
1371:    Vec options in the database.

1373:    Logically Collective on Vec

1375:    Input Parameter:
1376: +  v - the Vec context
1377: -  prefix - the prefix to prepend to all option names

1379:    Notes:
1380:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1381:    The first character of all runtime options is AUTOMATICALLY the hyphen.

1383:    Level: advanced

1385: .seealso: VecSetFromOptions()
1386: @*/
1387: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1388: {

1393:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1394:   return(0);
1395: }

1397: /*@C
1398:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1399:    Vec options in the database.

1401:    Logically Collective on Vec

1403:    Input Parameters:
1404: +  v - the Vec context
1405: -  prefix - the prefix to prepend to all option names

1407:    Notes:
1408:    A hyphen (-) must NOT be given at the beginning of the prefix name.
1409:    The first character of all runtime options is AUTOMATICALLY the hyphen.

1411:    Level: advanced

1413: .seealso: VecGetOptionsPrefix()
1414: @*/
1415: PetscErrorCode  VecAppendOptionsPrefix(Vec v,const char prefix[])
1416: {

1421:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1422:   return(0);
1423: }

1425: /*@C
1426:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1427:    Vec options in the database.

1429:    Not Collective

1431:    Input Parameter:
1432: .  v - the Vec context

1434:    Output Parameter:
1435: .  prefix - pointer to the prefix string used

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

1441:    Level: advanced

1443: .seealso: VecAppendOptionsPrefix()
1444: @*/
1445: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1446: {

1451:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1452:   return(0);
1453: }

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

1458:    Collective on Vec

1460:    Input Parameters:
1461: .  v - the Vec context

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

1467:    Level: advanced

1469: .seealso: VecCreate(), VecDestroy()
1470: @*/
1471: PetscErrorCode  VecSetUp(Vec v)
1472: {
1473:   PetscMPIInt    size;

1478:   if (v->map->n < 0 && v->map->N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1479:   if (!((PetscObject)v)->type_name) {
1480:     MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1481:     if (size == 1) {
1482:       VecSetType(v, VECSEQ);
1483:     } else {
1484:       VecSetType(v, VECMPI);
1485:     }
1486:   }
1487:   return(0);
1488: }

1490: /*
1491:     These currently expose the PetscScalar/PetscReal in updating the
1492:     cached norm. If we push those down into the implementation these
1493:     will become independent of PetscScalar/PetscReal
1494: */

1496: /*@
1497:    VecCopy - Copies a vector. y <- x

1499:    Logically Collective on Vec

1501:    Input Parameter:
1502: .  x - the vector

1504:    Output Parameter:
1505: .  y - the copy

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

1511:    Developer Notes:
1513:    of the vectors to be sequential and one to be parallel so long as both have the same
1514:    local sizes. This is used in some internal functions in PETSc.

1516:    Level: beginner

1518: .seealso: VecDuplicate()
1519: @*/
1520: PetscErrorCode  VecCopy(Vec x,Vec y)
1521: {
1522:   PetscBool      flgs[4];
1523:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1525:   PetscInt       i;

1532:   if (x == y) return(0);
1533:   VecCheckSameLocalSize(x,1,y,2);
1534:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1535:   VecSetErrorIfLocked(y,2);

1537: #if !defined(PETSC_USE_MIXED_PRECISION)
1538:   for (i=0; i<4; i++) {
1539:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1540:   }
1541: #endif

1543:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
1544: #if defined(PETSC_USE_MIXED_PRECISION)
1545:   extern PetscErrorCode VecGetArray(Vec,double**);
1546:   extern PetscErrorCode VecRestoreArray(Vec,double**);
1547:   extern PetscErrorCode VecGetArray(Vec,float**);
1548:   extern PetscErrorCode VecRestoreArray(Vec,float**);
1549:   extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1550:   extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1551:   extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1552:   extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1553:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1554:     PetscInt    i,n;
1555:     const float *xx;
1556:     double      *yy;
1557:     VecGetArrayRead(x,&xx);
1558:     VecGetArray(y,&yy);
1559:     VecGetLocalSize(x,&n);
1560:     for (i=0; i<n; i++) yy[i] = xx[i];
1561:     VecRestoreArrayRead(x,&xx);
1562:     VecRestoreArray(y,&yy);
1563:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1564:     PetscInt     i,n;
1565:     float        *yy;
1566:     const double *xx;
1567:     VecGetArrayRead(x,&xx);
1568:     VecGetArray(y,&yy);
1569:     VecGetLocalSize(x,&n);
1570:     for (i=0; i<n; i++) yy[i] = (float) xx[i];
1571:     VecRestoreArrayRead(x,&xx);
1572:     VecRestoreArray(y,&yy);
1573:   } else {
1574:     (*x->ops->copy)(x,y);
1575:   }
1576: #else
1577:   (*x->ops->copy)(x,y);
1578: #endif

1580:   PetscObjectStateIncrease((PetscObject)y);
1581: #if !defined(PETSC_USE_MIXED_PRECISION)
1582:   for (i=0; i<4; i++) {
1583:     if (flgs[i]) {
1584:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1585:     }
1586:   }
1587: #endif

1589:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1590:   return(0);
1591: }

1593: /*@
1594:    VecSwap - Swaps the vectors x and y.

1596:    Logically Collective on Vec

1598:    Input Parameters:
1599: .  x, y  - the vectors

1601:    Level: advanced

1603: @*/
1604: PetscErrorCode  VecSwap(Vec x,Vec y)
1605: {
1606:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1607:   PetscBool      flgxs[4],flgys[4];
1609:   PetscInt       i;

1617:   VecCheckSameSize(x,1,y,2);
1618:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1619:   if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");

1621:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1622:   for (i=0; i<4; i++) {
1623:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1624:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1625:   }
1626:   (*x->ops->swap)(x,y);
1627:   PetscObjectStateIncrease((PetscObject)x);
1628:   PetscObjectStateIncrease((PetscObject)y);
1629:   for (i=0; i<4; i++) {
1630:     if (flgxs[i]) {
1631:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1632:     }
1633:     if (flgys[i]) {
1634:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1635:     }
1636:   }
1637:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1638:   return(0);
1639: }

1641: /*
1642:   VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.

1644:   Collective on VecStash

1646:   Input Parameters:
1647: + obj   - the VecStash object
1648: . bobj - optional other object that provides the prefix
1649: - optionname - option to activate viewing

1651:   Level: intermediate

1653:   Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView

1655: */
1656: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1657: {
1658:   PetscErrorCode    ierr;
1659:   PetscViewer       viewer;
1660:   PetscBool         flg;
1661:   PetscViewerFormat format;
1662:   char              *prefix;

1665:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1666:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1667:   if (flg) {
1668:     PetscViewerPushFormat(viewer,format);
1669:     VecStashView(obj,viewer);
1670:     PetscViewerPopFormat(viewer);
1671:     PetscViewerDestroy(&viewer);
1672:   }
1673:   return(0);
1674: }

1676: /*@
1677:    VecStashView - Prints the entries in the vector stash and block stash.

1679:    Collective on Vec

1681:    Input Parameters:
1682: +  v - the vector
1683: -  viewer - the viewer

1685:    Level: advanced


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

1690: @*/
1691: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1692: {
1694:   PetscMPIInt    rank;
1695:   PetscInt       i,j;
1696:   PetscBool      match;
1697:   VecStash       *s;
1698:   PetscScalar    val;


1705:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1706:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1707:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1708:   MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1709:   s    = &v->bstash;

1711:   /* print block stash */
1712:   PetscViewerASCIIPushSynchronized(viewer);
1713:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1714:   for (i=0; i<s->n; i++) {
1715:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1716:     for (j=0; j<s->bs; j++) {
1717:       val = s->array[i*s->bs+j];
1718: #if defined(PETSC_USE_COMPLEX)
1719:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1720: #else
1721:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1722: #endif
1723:     }
1724:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1725:   }
1726:   PetscViewerFlush(viewer);

1728:   s = &v->stash;

1730:   /* print basic stash */
1731:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1732:   for (i=0; i<s->n; i++) {
1733:     val = s->array[i];
1734: #if defined(PETSC_USE_COMPLEX)
1735:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1736: #else
1737:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1738: #endif
1739:   }
1740:   PetscViewerFlush(viewer);
1741:   PetscViewerASCIIPopSynchronized(viewer);
1742:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1743:   return(0);
1744: }

1746: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1747: {
1748:   PetscInt       i,N,rstart,rend;
1750:   PetscScalar    *xx;
1751:   PetscReal      *xreal;
1752:   PetscBool      iset;

1755:   VecGetOwnershipRange(v,&rstart,&rend);
1756:   VecGetSize(v,&N);
1757:   PetscCalloc1(N,&xreal);
1758:   PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1759:   if (iset) {
1760:     VecGetArray(v,&xx);
1761:     for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1762:     VecRestoreArray(v,&xx);
1763:   }
1764:   PetscFree(xreal);
1765:   if (set) *set = iset;
1766:   return(0);
1767: }

1769: /*@
1770:    VecGetLayout - get PetscLayout describing vector layout

1772:    Not Collective

1774:    Input Arguments:
1775: .  x - the vector

1777:    Output Arguments:
1778: .  map - the layout

1780:    Level: developer

1782: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1783: @*/
1784: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1785: {

1789:   *map = x->map;
1790:   return(0);
1791: }

1793: /*@
1794:    VecSetLayout - set PetscLayout describing vector layout

1796:    Not Collective

1798:    Input Arguments:
1799: +  x - the vector
1800: -  map - the layout

1802:    Notes:
1803:    It is normally only valid to replace the layout with a layout known to be equivalent.

1805:    Level: developer

1807: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1808: @*/
1809: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1810: {

1815:   PetscLayoutReference(map,&x->map);
1816:   return(0);
1817: }

1819: PetscErrorCode VecSetInf(Vec xin)
1820: {
1821:   PetscInt       i,n = xin->map->n;
1822:   PetscScalar    *xx;
1823:   PetscScalar    zero=0.0,one=1.0,inf=one/zero;

1827:   VecGetArray(xin,&xx);
1828:   for (i=0; i<n; i++) xx[i] = inf;
1829:   VecRestoreArray(xin,&xx);
1830:   return(0);
1831: }

1833: /*@
1834:      VecPinToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU

1836:    Input Parameters:
1837: +   v - the vector
1838: -   flg - pin to the CPU if value of PETSC_TRUE

1840:    Level: intermediate
1841: @*/
1842: PetscErrorCode VecPinToCPU(Vec v,PetscBool flg)
1843: {
1844: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)

1848:   if (v->pinnedtocpu == flg) return(0);
1849:   v->pinnedtocpu = flg;
1850:   if (v->ops->pintocpu) {
1851:     (*v->ops->pintocpu)(v,flg);
1852:   }
1853:   return(0);
1854: #else
1855:   return 0;
1856: #endif
1857: }