Actual source code: vector.c

petsc-3.11.2 2019-05-18
Report Typos and Errors

  2: /*
  3:      Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
  4:    These are the vector functions the user calls.
  5: */
  6:  #include <petsc/private/vecimpl.h>

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

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

 24:     Not collective

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

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

 35:    Level: advanced

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

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

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

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

 55:    Logically Collective on Vec

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

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

 64:    Level: intermediate

 66:    Concepts: vector^setting values with local numbering

 68: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
 69:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
 70: @*/
 71: PetscErrorCode  VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
 72: {


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

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

 90:    Not Collective

 92:    Input Parameter:
 93: .  X - the vector

 95:    Output Parameter:
 96: .  mapping - the mapping

 98:    Level: advanced

100:    Concepts: vectors^local to global mapping
101:    Concepts: local to global mapping^for vectors

103: .seealso:  VecSetValuesLocal()
104: @*/
105: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
106: {
111:   *mapping = X->map->mapping;
112:   return(0);
113: }

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

119:    Collective on Vec

121:    Input Parameter:
122: .  vec - the vector

124:    Level: beginner

126:    Concepts: assembly^vectors

128: .seealso: VecAssemblyEnd(), VecSetValues()
129: @*/
130: PetscErrorCode  VecAssemblyBegin(Vec vec)
131: {

137:   VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
138:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
139:   if (vec->ops->assemblybegin) {
140:     (*vec->ops->assemblybegin)(vec);
141:   }
142:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
143:   PetscObjectStateIncrease((PetscObject)vec);
144:   return(0);
145: }

147: /*@
148:    VecAssemblyEnd - Completes assembling the vector.  This routine should
149:    be called after VecAssemblyBegin().

151:    Collective on Vec

153:    Input Parameter:
154: .  vec - the vector

156:    Options Database Keys:
157: +  -vec_view - Prints vector in ASCII format
158: .  -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
159: .  -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
160: .  -vec_view draw - Activates vector viewing using drawing tools
161: .  -display <name> - Sets display name (default is host)
162: .  -draw_pause <sec> - Sets number of seconds to pause after display
163: -  -vec_view socket - Activates vector viewing using a socket

165:    Level: beginner

167: .seealso: VecAssemblyBegin(), VecSetValues()
168: @*/
169: PetscErrorCode  VecAssemblyEnd(Vec vec)
170: {

175:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
177:   if (vec->ops->assemblyend) {
178:     (*vec->ops->assemblyend)(vec);
179:   }
180:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
181:   VecViewFromOptions(vec,NULL,"-vec_view");
182:   return(0);
183: }

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

188:    Logically Collective on Vec

190:    Input Parameters:
191: .  x, y  - the vectors

193:    Output Parameter:
194: .  w - the result

196:    Level: advanced

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

202:    Concepts: vector^pointwise multiply

204: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
205: @*/
206: PetscErrorCode  VecPointwiseMax(Vec w,Vec x,Vec y)
207: {

219:   VecCheckSameSize(w,1,x,2);
220:   VecCheckSameSize(w,1,y,3);
221:   (*w->ops->pointwisemax)(w,x,y);
222:   PetscObjectStateIncrease((PetscObject)w);
223:   return(0);
224: }


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

230:    Logically Collective on Vec

232:    Input Parameters:
233: .  x, y  - the vectors

235:    Output Parameter:
236: .  w - the result

238:    Level: advanced

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

244:    Concepts: vector^pointwise multiply

246: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
247: @*/
248: PetscErrorCode  VecPointwiseMin(Vec w,Vec x,Vec y)
249: {

261:   VecCheckSameSize(w,1,x,2);
262:   VecCheckSameSize(w,1,y,3);
263:   (*w->ops->pointwisemin)(w,x,y);
264:   PetscObjectStateIncrease((PetscObject)w);
265:   return(0);
266: }

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

271:    Logically Collective on Vec

273:    Input Parameters:
274: .  x, y  - the vectors

276:    Output Parameter:
277: .  w - the result

279:    Level: advanced

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

284:    Concepts: vector^pointwise multiply

286: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
287: @*/
288: PetscErrorCode  VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
289: {

301:   VecCheckSameSize(w,1,x,2);
302:   VecCheckSameSize(w,1,y,3);
303:   (*w->ops->pointwisemaxabs)(w,x,y);
304:   PetscObjectStateIncrease((PetscObject)w);
305:   return(0);
306: }

308: /*@
309:    VecPointwiseDivide - Computes the componentwise division w = x/y.

311:    Logically Collective on Vec

313:    Input Parameters:
314: .  x, y  - the vectors

316:    Output Parameter:
317: .  w - the result

319:    Level: advanced

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

324:    Concepts: vector^pointwise divide

326: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
327: @*/
328: PetscErrorCode  VecPointwiseDivide(Vec w,Vec x,Vec y)
329: {

341:   VecCheckSameSize(w,1,x,2);
342:   VecCheckSameSize(w,1,y,3);
343:   (*w->ops->pointwisedivide)(w,x,y);
344:   PetscObjectStateIncrease((PetscObject)w);
345:   return(0);
346: }


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

352:    Collective on Vec

354:    Input Parameters:
355: .  v - a vector to mimic

357:    Output Parameter:
358: .  newv - location to put new vector

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

364:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
365:    vectors.

367:    Level: beginner

369: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
370: @*/
371: PetscErrorCode  VecDuplicate(Vec v,Vec *newv)
372: {

379:   (*v->ops->duplicate)(v,newv);
380:   PetscObjectStateIncrease((PetscObject)*newv);
381:   return(0);
382: }

384: /*@
385:    VecDestroy - Destroys a vector.

387:    Collective on Vec

389:    Input Parameters:
390: .  v  - the vector

392:    Level: beginner

394: .seealso: VecDuplicate(), VecDestroyVecs()
395: @*/
396: PetscErrorCode  VecDestroy(Vec *v)
397: {

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

405:   PetscObjectSAWsViewOff((PetscObject)*v);
406:   /* destroy the internal part */
407:   if ((*v)->ops->destroy) {
408:     (*(*v)->ops->destroy)(*v);
409:   }
410:   /* destroy the external/common part */
411:   PetscLayoutDestroy(&(*v)->map);
412:   PetscHeaderDestroy(v);
413:   return(0);
414: }

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

419:    Collective on Vec

421:    Input Parameters:
422: +  m - the number of vectors to obtain
423: -  v - a vector to mimic

425:    Output Parameter:
426: .  V - location to put pointer to array of vectors

428:    Notes:
429:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
430:    vector.

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

437:    Level: intermediate

439: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
440: @*/
441: PetscErrorCode  VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
442: {

449:   (*v->ops->duplicatevecs)(v,m,V);
450:   return(0);
451: }

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

456:    Collective on Vec

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

462:    Fortran Note:
463:    The Fortran interface is slightly different from that given below.
464:    See the Fortran chapter of the users manual

466:    Level: intermediate

468: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
469: @*/
470: PetscErrorCode  VecDestroyVecs(PetscInt m,Vec *vv[])
471: {

476:   if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
477:   if (!m || !*vv) {*vv  = NULL; return(0);}
480:   (*(**vv)->ops->destroyvecs)(m,*vv);
481:   *vv  = NULL;
482:   return(0);
483: }

485: /*@C
486:    VecView - Views a vector object.

488:    Collective on Vec

490:    Input Parameters:
491: +  vec - the vector
492: -  viewer - an optional visualization context

494:    Notes:
495:    The available visualization contexts include
496: +     PETSC_VIEWER_STDOUT_SELF - for sequential vectors
497: .     PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
498: -     PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm

500:    You can change the format the vector is printed using the
501:    option PetscViewerPushFormat().

503:    The user can open alternative visualization contexts with
504: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
505: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
506:          specified file; corresponding input uses VecLoad()
507: .    PetscViewerDrawOpen() - Outputs vector to an X window display
508: -    PetscViewerSocketOpen() - Outputs vector to Socket viewer

510:    The user can call PetscViewerPushFormat() to specify the output
511:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
512:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
513: +    PETSC_VIEWER_DEFAULT - default, prints vector contents
514: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
515: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
516: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
517:          format common among all vector types

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

522:    Notes for binary viewer: If you pass multiply vectors to a binary viewer you can read them back in in the same order
523: $     with VecLoad().
524: $
525: $    If the blocksize of the vector is greater than one then you must provide a unique prefix to
526: $    the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
527: $    vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
528: $    information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
529: $    filename. If you copy the binary file, make sure you copy the associated .info file with it.

531:    Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
532: $    for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
533: $    obviously) several times, you must change its name each time before calling the VecView(). The name you use
534: $    here should equal the name that you use in the Vec object that you use with VecLoad().

536:    See the manual page for VecLoad() on the exact format the binary viewer stores
537:    the values in the file.

539:    Level: beginner

541:    Concepts: vector^printing
542:    Concepts: vector^saving to disk

544: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
545:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
546:           PetscRealView(), PetscScalarView(), PetscIntView()
547: @*/
548: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
549: {
550:   PetscErrorCode    ierr;
551:   PetscBool         iascii;
552:   PetscViewerFormat format;
553:   PetscMPIInt       size;

558:   if (!viewer) {
559:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
560:   }
562:   PetscViewerGetFormat(viewer,&format);
563:   MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
564:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);

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

568:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
569:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
570:   if (iascii) {
571:     PetscInt rows,bs;

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

597: #if defined(PETSC_USE_DEBUG)
598: #include <../src/sys/totalview/tv_data_display.h>
599: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
600: {
601:   const PetscScalar *values;
602:   char              type[32];
603:   PetscErrorCode    ierr;


606:   TV_add_row("Local rows", "int", &v->map->n);
607:   TV_add_row("Global rows", "int", &v->map->N);
608:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
609:   VecGetArrayRead((Vec)v,&values);
610:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
611:   TV_add_row("values",type, values);
612:   VecRestoreArrayRead((Vec)v,&values);
613:   return TV_format_OK;
614: }
615: #endif

617: /*@
618:    VecGetSize - Returns the global number of elements of the vector.

620:    Not Collective

622:    Input Parameter:
623: .  x - the vector

625:    Output Parameters:
626: .  size - the global length of the vector

628:    Level: beginner

630:    Concepts: vector^local size

632: .seealso: VecGetLocalSize()
633: @*/
634: PetscErrorCode  VecGetSize(Vec x,PetscInt *size)
635: {

642:   (*x->ops->getsize)(x,size);
643:   return(0);
644: }

646: /*@
647:    VecGetLocalSize - Returns the number of elements of the vector stored
648:    in local memory. This routine may be implementation dependent, so use
649:    with care.

651:    Not Collective

653:    Input Parameter:
654: .  x - the vector

656:    Output Parameter:
657: .  size - the length of the local piece of the vector

659:    Level: beginner

661:    Concepts: vector^size

663: .seealso: VecGetSize()
664: @*/
665: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
666: {

673:   (*x->ops->getlocalsize)(x,size);
674:   return(0);
675: }

677: /*@C
678:    VecGetOwnershipRange - Returns the range of indices owned by
679:    this processor, assuming that the vectors are laid out with the
680:    first n1 elements on the first processor, next n2 elements on the
681:    second, etc.  For certain parallel layouts this range may not be
682:    well defined.

684:    Not Collective

686:    Input Parameter:
687: .  x - the vector

689:    Output Parameters:
690: +  low - the first local element, pass in NULL if not interested
691: -  high - one more than the last local element, pass in NULL if not interested

693:    Note:
694:    The high argument is one more than the last element stored locally.

696:    Fortran: PETSC_NULL_INTEGER should be used instead of NULL

698:    Level: beginner

700:    Concepts: ownership^of vectors
701:    Concepts: vector^ownership of elements

703: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
704: @*/
705: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
706: {
712:   if (low)  *low  = x->map->rstart;
713:   if (high) *high = x->map->rend;
714:   return(0);
715: }

717: /*@C
718:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
719:    assuming that the vectors are laid out with the
720:    first n1 elements on the first processor, next n2 elements on the
721:    second, etc.  For certain parallel layouts this range may not be
722:    well defined.

724:    Not Collective

726:    Input Parameter:
727: .  x - the vector

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

732:    Note:
733:    The high argument is one more than the last element stored locally.

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

737:    Level: beginner

739:    Concepts: ownership^of vectors
740:    Concepts: vector^ownership of elements

742: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
743: @*/
744: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
745: {

751:   PetscLayoutGetRanges(x->map,ranges);
752:   return(0);
753: }

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

758:    Collective on Vec

760:    Input Parameter:
761: +  x - the vector
762: .  op - the option
763: -  flag - turn the option on or off

765:    Supported Options:
766: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
767:           entries destined to be stored on a separate processor. This can be used
768:           to eliminate the global reduction in the VecAssemblyXXXX() if you know
769:           that you have only used VecSetValues() to set local elements
770: .     VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
771:           in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
772:           ignored.
773: -     VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
774:           entries will always be a subset (possibly equal) of the off-process entries set on the
775:           first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
776:           changed this flag afterwards. If this assembly is not such first assembly, then this
777:           assembly can reuse the communication pattern setup in that first assembly, thus avoiding
778:           a global reduction. Subsequent assemblies setting off-process values should use the same
779:           InsertMode as the first assembly.

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

784:    Level: intermediate

786: @*/
787: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool flag)
788: {

794:   if (x->ops->setoption) {
795:     (*x->ops->setoption)(x,op,flag);
796:   }
797:   return(0);
798: }

800: /* Default routines for obtaining and releasing; */
801: /* may be used by any implementation */
802: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
803: {
805:   PetscInt       i;

810:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
811:   PetscMalloc1(m,V);
812:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
813:   return(0);
814: }

816: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
817: {
819:   PetscInt       i;

823:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
824:   PetscFree(v);
825:   return(0);
826: }

828: /*@
829:    VecResetArray - Resets a vector to use its default memory. Call this
830:    after the use of VecPlaceArray().

832:    Not Collective

834:    Input Parameters:
835: .  vec - the vector

837:    Level: developer

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

841: @*/
842: PetscErrorCode  VecResetArray(Vec vec)
843: {

849:   if (vec->ops->resetarray) {
850:     (*vec->ops->resetarray)(vec);
851:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
852:   PetscObjectStateIncrease((PetscObject)vec);
853:   return(0);
854: }

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

860:   Collective on PetscViewer

862:   Input Parameters:
863: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
864:            some related function before a call to VecLoad().
865: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
866:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

868:    Level: intermediate

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

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

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

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

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

891:   Notes for advanced users:
892:   Most users should not need to know the details of the binary storage
893:   format, since VecLoad() and VecView() completely hide these details.
894:   But for anyone who's interested, the standard binary vector storage
895:   format is
896: .vb
897:      int    VEC_FILE_CLASSID
898:      int    number of rows
899:      PetscScalar *values of all entries
900: .ve

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

908:   Concepts: vector^loading from file

910: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
911: @*/
912: PetscErrorCode  VecLoad(Vec newvec, PetscViewer viewer)
913: {
914:   PetscErrorCode    ierr;
915:   PetscBool         isbinary,ishdf5,isadios,isadios2;
916:   PetscViewerFormat format;

921:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
922:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
923:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
924:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
925:   if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

927:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
928:   if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
929:     VecSetType(newvec, VECSTANDARD);
930:   }
931:   PetscViewerGetFormat(viewer,&format);
932:   if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
933:     (*newvec->ops->loadnative)(newvec,viewer);
934:   } else {
935:     (*newvec->ops->load)(newvec,viewer);
936:   }
937:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
938:   return(0);
939: }


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

945:    Logically Collective on Vec

947:    Input Parameter:
948: .  vec - the vector

950:    Output Parameter:
951: .  vec - the vector reciprocal

953:    Level: intermediate

955:    Concepts: vector^reciprocal

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

959: @*/
960: PetscErrorCode  VecReciprocal(Vec vec)
961: {

967:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
968:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
969:   (*vec->ops->reciprocal)(vec);
970:   PetscObjectStateIncrease((PetscObject)vec);
971:   return(0);
972: }

974: /*@C
975:     VecSetOperation - Allows user to set a vector operation.

977:    Logically Collective on Vec

979:     Input Parameters:
980: +   vec - the vector
981: .   op - the name of the operation
982: -   f - the function that provides the operation.

984:    Level: advanced

986:     Usage:
987: $      PetscErrorCode userview(Vec,PetscViewer);
988: $      VecCreateMPI(comm,m,M,&x);
989: $      VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);

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

997:     This function is not currently available from Fortran.

999: .keywords: vector, set, operation

1001: .seealso: VecCreate(), MatShellSetOperation()
1002: @*/
1003: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1004: {
1007:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
1008:     vec->ops->viewnative = vec->ops->view;
1009:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1010:     vec->ops->loadnative = vec->ops->load;
1011:   }
1012:   (((void(**)(void))vec->ops)[(int)op]) = f;
1013:   return(0);
1014: }


1017: /*@
1018:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1019:    used during the assembly process to store values that belong to
1020:    other processors.

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

1024:    Input Parameters:
1025: +  vec   - the vector
1026: .  size  - the initial size of the stash.
1027: -  bsize - the initial size of the block-stash(if used).

1029:    Options Database Keys:
1030: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1031: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1033:    Level: intermediate

1035:    Notes:
1036:      The block-stash is used for values set with VecSetValuesBlocked() while
1037:      the stash is used for values set with VecSetValues()

1039:      Run with the option -info and look for output of the form
1040:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1041:      to determine the appropriate value, MM, to use for size and
1042:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1043:      to determine the value, BMM to use for bsize

1045:    Concepts: vector^stash
1046:    Concepts: stash^vector

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

1050: @*/
1051: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1052: {

1057:   VecStashSetInitialSize_Private(&vec->stash,size);
1058:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1059:   return(0);
1060: }

1062: /*@
1063:    VecConjugate - Conjugates a vector.

1065:    Logically Collective on Vec

1067:    Input Parameters:
1068: .  x - the vector

1070:    Level: intermediate

1072:    Concepts: vector^conjugate

1074: @*/
1075: PetscErrorCode  VecConjugate(Vec x)
1076: {
1077: #if defined(PETSC_USE_COMPLEX)

1083:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1084:   (*x->ops->conjugate)(x);
1085:   /* we need to copy norms here */
1086:   PetscObjectStateIncrease((PetscObject)x);
1087:   return(0);
1088: #else
1089:   return(0);
1090: #endif
1091: }

1093: /*@
1094:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1096:    Logically Collective on Vec

1098:    Input Parameters:
1099: .  x, y  - the vectors

1101:    Output Parameter:
1102: .  w - the result

1104:    Level: advanced

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

1109:    Concepts: vector^pointwise multiply

1111: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1112: @*/
1113: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1114: {

1126:   VecCheckSameSize(w,1,x,2);
1127:   VecCheckSameSize(w,2,y,3);
1128:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1129:   (*w->ops->pointwisemult)(w,x,y);
1130:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1131:   PetscObjectStateIncrease((PetscObject)w);
1132:   return(0);
1133: }

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

1138:    Logically Collective on Vec

1140:    Input Parameters:
1141: +  x  - the vector
1142: -  rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1143:           it will create one internally.

1145:    Output Parameter:
1146: .  x  - the vector

1148:    Example of Usage:
1149: .vb
1150:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1151:      VecSetRandom(x,rctx);
1152:      PetscRandomDestroy(rctx);
1153: .ve

1155:    Level: intermediate

1157:    Concepts: vector^setting to random
1158:    Concepts: random^vector

1160: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1161: @*/
1162: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1163: {
1165:   PetscRandom    randObj = NULL;

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

1173:   if (!rctx) {
1174:     MPI_Comm comm;
1175:     PetscObjectGetComm((PetscObject)x,&comm);
1176:     PetscRandomCreate(comm,&randObj);
1177:     PetscRandomSetFromOptions(randObj);
1178:     rctx = randObj;
1179:   }

1181:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1182:   (*x->ops->setrandom)(x,rctx);
1183:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1185:   PetscRandomDestroy(&randObj);
1186:   PetscObjectStateIncrease((PetscObject)x);
1187:   return(0);
1188: }

1190: /*@
1191:   VecZeroEntries - puts a 0.0 in each element of a vector

1193:   Logically Collective on Vec

1195:   Input Parameter:
1196: . vec - The vector

1198:   Level: beginner

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

1205: .keywords: Vec, set, options, database
1206: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1207: @*/
1208: PetscErrorCode  VecZeroEntries(Vec vec)
1209: {

1213:   VecSet(vec,0);
1214:   return(0);
1215: }

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

1221:   Collective on Vec

1223:   Input Parameter:
1224: . vec - The vector

1226:   Level: intermediate

1228: .keywords: Vec, set, options, database, type
1229: .seealso: VecSetFromOptions(), VecSetType()
1230: */
1231: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1232: {
1233:   PetscBool      opt;
1234:   VecType        defaultType;
1235:   char           typeName[256];
1236:   PetscMPIInt    size;

1240:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1241:   else {
1242:     MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1243:     if (size > 1) defaultType = VECMPI;
1244:     else defaultType = VECSEQ;
1245:   }

1247:   VecRegisterAll();
1248:   PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1249:   if (opt) {
1250:     VecSetType(vec, typeName);
1251:   } else {
1252:     VecSetType(vec, defaultType);
1253:   }
1254:   return(0);
1255: }

1257: /*@
1258:   VecSetFromOptions - Configures the vector from the options database.

1260:   Collective on Vec

1262:   Input Parameter:
1263: . vec - The vector

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

1269:   Level: beginner

1271:   Concepts: vectors^setting options
1272:   Concepts: vectors^setting type

1274: .keywords: Vec, set, options, database
1275: .seealso: VecCreate(), VecSetOptionsPrefix()
1276: @*/
1277: PetscErrorCode  VecSetFromOptions(Vec vec)
1278: {


1284:   PetscObjectOptionsBegin((PetscObject)vec);
1285:   /* Handle vector type options */
1286:   VecSetTypeFromOptions_Private(PetscOptionsObject,vec);

1288:   /* Handle specific vector options */
1289:   if (vec->ops->setfromoptions) {
1290:     (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1291:   }

1293:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1294:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1295:   PetscOptionsEnd();
1296:   return(0);
1297: }

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

1302:   Collective on Vec

1304:   Input Parameters:
1305: + v - the vector
1306: . n - the local size (or PETSC_DECIDE to have it set)
1307: - N - the global size (or PETSC_DECIDE)

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

1313:   Level: intermediate

1315: .seealso: VecGetSize(), PetscSplitOwnership()
1316: @*/
1317: PetscErrorCode  VecSetSizes(Vec v, PetscInt n, PetscInt N)
1318: {

1324:   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);
1325:   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);
1326:   v->map->n = n;
1327:   v->map->N = N;
1328:   if (v->ops->create) {
1329:     (*v->ops->create)(v);
1330:     v->ops->create = 0;
1331:   }
1332:   return(0);
1333: }

1335: /*@
1336:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1337:    and VecSetValuesBlockedLocal().

1339:    Logically Collective on Vec

1341:    Input Parameter:
1342: +  v - the vector
1343: -  bs - the blocksize

1345:    Notes:
1346:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1348:    Level: advanced

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

1352:   Concepts: block size^vectors
1353: @*/
1354: PetscErrorCode  VecSetBlockSize(Vec v,PetscInt bs)
1355: {

1360:   if (bs < 0 || bs == v->map->bs) return(0);
1362:   PetscLayoutSetBlockSize(v->map,bs);
1363:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1364:   return(0);
1365: }

1367: /*@
1368:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1369:    and VecSetValuesBlockedLocal().

1371:    Not Collective

1373:    Input Parameter:
1374: .  v - the vector

1376:    Output Parameter:
1377: .  bs - the blocksize

1379:    Notes:
1380:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1382:    Level: advanced

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

1386:    Concepts: vector^block size
1387:    Concepts: block^vector

1389: @*/
1390: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1391: {

1397:   PetscLayoutGetBlockSize(v->map,bs);
1398:   return(0);
1399: }

1401: /*@C
1402:    VecSetOptionsPrefix - Sets the prefix used for searching for all
1403:    Vec options in the database.

1405:    Logically Collective on Vec

1407:    Input Parameter:
1408: +  v - the Vec context
1409: -  prefix - the prefix to prepend to all option names

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

1415:    Level: advanced

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

1419: .seealso: VecSetFromOptions()
1420: @*/
1421: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1422: {

1427:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1428:   return(0);
1429: }

1431: /*@C
1432:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1433:    Vec options in the database.

1435:    Logically Collective on Vec

1437:    Input Parameters:
1438: +  v - the Vec context
1439: -  prefix - the prefix to prepend to all option names

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

1445:    Level: advanced

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

1449: .seealso: VecGetOptionsPrefix()
1450: @*/
1451: PetscErrorCode  VecAppendOptionsPrefix(Vec v,const char prefix[])
1452: {

1457:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1458:   return(0);
1459: }

1461: /*@C
1462:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1463:    Vec options in the database.

1465:    Not Collective

1467:    Input Parameter:
1468: .  v - the Vec context

1470:    Output Parameter:
1471: .  prefix - pointer to the prefix string used

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

1477:    Level: advanced

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

1481: .seealso: VecAppendOptionsPrefix()
1482: @*/
1483: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1484: {

1489:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1490:   return(0);
1491: }

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

1496:    Collective on Vec

1498:    Input Parameters:
1499: .  v - the Vec context

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

1505:    Level: advanced

1507: .keywords: Vec, setup

1509: .seealso: VecCreate(), VecDestroy()
1510: @*/
1511: PetscErrorCode  VecSetUp(Vec v)
1512: {
1513:   PetscMPIInt    size;

1518:   if (!((PetscObject)v)->type_name) {
1519:     MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1520:     if (size == 1) {
1521:       VecSetType(v, VECSEQ);
1522:     } else {
1523:       VecSetType(v, VECMPI);
1524:     }
1525:   }
1526:   return(0);
1527: }

1529: /*
1530:     These currently expose the PetscScalar/PetscReal in updating the
1531:     cached norm. If we push those down into the implementation these
1532:     will become independent of PetscScalar/PetscReal
1533: */

1535: /*@
1536:    VecCopy - Copies a vector. y <- x

1538:    Logically Collective on Vec

1540:    Input Parameter:
1541: .  x - the vector

1543:    Output Parameter:
1544: .  y - the copy

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

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

1555:    Level: beginner

1557: .seealso: VecDuplicate()
1558: @*/
1559: PetscErrorCode  VecCopy(Vec x,Vec y)
1560: {
1561:   PetscBool      flgs[4];
1562:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1564:   PetscInt       i;

1571:   if (x == y) return(0);
1572:   VecCheckSameLocalSize(x,1,y,2);
1573:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1574:   VecSetErrorIfLocked(y,2);

1576: #if !defined(PETSC_USE_MIXED_PRECISION)
1577:   for (i=0; i<4; i++) {
1578:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1579:   }
1580: #endif

1582:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
1583: #if defined(PETSC_USE_MIXED_PRECISION)
1584:   extern PetscErrorCode VecGetArray(Vec,double**);
1585:   extern PetscErrorCode VecRestoreArray(Vec,double**);
1586:   extern PetscErrorCode VecGetArray(Vec,float**);
1587:   extern PetscErrorCode VecRestoreArray(Vec,float**);
1588:   extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1589:   extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1590:   extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1591:   extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1592:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1593:     PetscInt    i,n;
1594:     const float *xx;
1595:     double      *yy;
1596:     VecGetArrayRead(x,&xx);
1597:     VecGetArray(y,&yy);
1598:     VecGetLocalSize(x,&n);
1599:     for (i=0; i<n; i++) yy[i] = xx[i];
1600:     VecRestoreArrayRead(x,&xx);
1601:     VecRestoreArray(y,&yy);
1602:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1603:     PetscInt     i,n;
1604:     float        *yy;
1605:     const double *xx;
1606:     VecGetArrayRead(x,&xx);
1607:     VecGetArray(y,&yy);
1608:     VecGetLocalSize(x,&n);
1609:     for (i=0; i<n; i++) yy[i] = (float) xx[i];
1610:     VecRestoreArrayRead(x,&xx);
1611:     VecRestoreArray(y,&yy);
1612:   } else {
1613:     (*x->ops->copy)(x,y);
1614:   }
1615: #else
1616:   (*x->ops->copy)(x,y);
1617: #endif

1619:   PetscObjectStateIncrease((PetscObject)y);
1620: #if !defined(PETSC_USE_MIXED_PRECISION)
1621:   for (i=0; i<4; i++) {
1622:     if (flgs[i]) {
1623:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1624:     }
1625:   }
1626: #endif

1628:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1629:   return(0);
1630: }

1632: /*@
1633:    VecSwap - Swaps the vectors x and y.

1635:    Logically Collective on Vec

1637:    Input Parameters:
1638: .  x, y  - the vectors

1640:    Level: advanced

1642:    Concepts: vector^swapping values

1644: @*/
1645: PetscErrorCode  VecSwap(Vec x,Vec y)
1646: {
1647:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1648:   PetscBool      flgxs[4],flgys[4];
1650:   PetscInt       i;

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

1662:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1663:   for (i=0; i<4; i++) {
1664:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1665:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1666:   }
1667:   (*x->ops->swap)(x,y);
1668:   PetscObjectStateIncrease((PetscObject)x);
1669:   PetscObjectStateIncrease((PetscObject)y);
1670:   for (i=0; i<4; i++) {
1671:     if (flgxs[i]) {
1672:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1673:     }
1674:     if (flgys[i]) {
1675:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1676:     }
1677:   }
1678:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1679:   return(0);
1680: }

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

1685:   Collective on VecStash

1687:   Input Parameters:
1688: + obj   - the VecStash object
1689: . bobj - optional other object that provides the prefix
1690: - optionname - option to activate viewing

1692:   Level: intermediate

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

1696: */
1697: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1698: {
1699:   PetscErrorCode    ierr;
1700:   PetscViewer       viewer;
1701:   PetscBool         flg;
1702:   PetscViewerFormat format;
1703:   char              *prefix;

1706:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1707:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1708:   if (flg) {
1709:     PetscViewerPushFormat(viewer,format);
1710:     VecStashView(obj,viewer);
1711:     PetscViewerPopFormat(viewer);
1712:     PetscViewerDestroy(&viewer);
1713:   }
1714:   return(0);
1715: }

1717: /*@
1718:    VecStashView - Prints the entries in the vector stash and block stash.

1720:    Collective on Vec

1722:    Input Parameters:
1723: +  v - the vector
1724: -  viewer - the viewer

1726:    Level: advanced

1728:    Concepts: vector^stash
1729:    Concepts: stash^vector

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

1733: @*/
1734: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1735: {
1737:   PetscMPIInt    rank;
1738:   PetscInt       i,j;
1739:   PetscBool      match;
1740:   VecStash       *s;
1741:   PetscScalar    val;


1748:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1749:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1750:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1751:   MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1752:   s    = &v->bstash;

1754:   /* print block stash */
1755:   PetscViewerASCIIPushSynchronized(viewer);
1756:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1757:   for (i=0; i<s->n; i++) {
1758:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1759:     for (j=0; j<s->bs; j++) {
1760:       val = s->array[i*s->bs+j];
1761: #if defined(PETSC_USE_COMPLEX)
1762:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1763: #else
1764:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1765: #endif
1766:     }
1767:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1768:   }
1769:   PetscViewerFlush(viewer);

1771:   s = &v->stash;

1773:   /* print basic stash */
1774:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1775:   for (i=0; i<s->n; i++) {
1776:     val = s->array[i];
1777: #if defined(PETSC_USE_COMPLEX)
1778:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1779: #else
1780:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1781: #endif
1782:   }
1783:   PetscViewerFlush(viewer);
1784:   PetscViewerASCIIPopSynchronized(viewer);
1785:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1786:   return(0);
1787: }

1789: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1790: {
1791:   PetscInt       i,N,rstart,rend;
1793:   PetscScalar    *xx;
1794:   PetscReal      *xreal;
1795:   PetscBool      iset;

1798:   VecGetOwnershipRange(v,&rstart,&rend);
1799:   VecGetSize(v,&N);
1800:   PetscCalloc1(N,&xreal);
1801:   PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1802:   if (iset) {
1803:     VecGetArray(v,&xx);
1804:     for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1805:     VecRestoreArray(v,&xx);
1806:   }
1807:   PetscFree(xreal);
1808:   if (set) *set = iset;
1809:   return(0);
1810: }

1812: /*@
1813:    VecGetLayout - get PetscLayout describing vector layout

1815:    Not Collective

1817:    Input Arguments:
1818: .  x - the vector

1820:    Output Arguments:
1821: .  map - the layout

1823:    Level: developer

1825: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1826: @*/
1827: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1828: {

1832:   *map = x->map;
1833:   return(0);
1834: }

1836: /*@
1837:    VecSetLayout - set PetscLayout describing vector layout

1839:    Not Collective

1841:    Input Arguments:
1842: +  x - the vector
1843: -  map - the layout

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

1848:    Level: developer

1850: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1851: @*/
1852: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1853: {

1858:   PetscLayoutReference(map,&x->map);
1859:   return(0);
1860: }

1862: PetscErrorCode VecSetInf(Vec xin)
1863: {
1864:   PetscInt       i,n = xin->map->n;
1865:   PetscScalar    *xx;
1866:   PetscScalar    zero=0.0,one=1.0,inf=one/zero;

1870:   VecGetArray(xin,&xx);
1871:   for (i=0; i<n; i++) xx[i] = inf;
1872:   VecRestoreArray(xin,&xx);
1873:   return(0);
1874: }