Actual source code: vector.c

petsc-3.13.1 2020-05-02
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:   VecSetErrorIfLocked(w,1);
213:   (*w->ops->pointwisemax)(w,x,y);
214:   PetscObjectStateIncrease((PetscObject)w);
215:   return(0);
216: }

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:   VecSetErrorIfLocked(w,1);
253:   (*w->ops->pointwisemin)(w,x,y);
254:   PetscObjectStateIncrease((PetscObject)w);
255:   return(0);
256: }

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

261:    Logically Collective on Vec

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

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

269:    Level: advanced

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

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

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

297: /*@
298:    VecPointwiseDivide - Computes the componentwise division w = x/y.

300:    Logically Collective on Vec

302:    Input Parameters:
303: .  x, y  - the vectors

305:    Output Parameter:
306: .  w - the result

308:    Level: advanced

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

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

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


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

340:    Collective on Vec

342:    Input Parameters:
343: .  v - a vector to mimic

345:    Output Parameter:
346: .  newv - location to put new vector

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

352:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
353:    vectors.

355:    Level: beginner

357: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
358: @*/
359: PetscErrorCode  VecDuplicate(Vec v,Vec *newv)
360: {

367:   (*v->ops->duplicate)(v,newv);
368:   PetscObjectStateIncrease((PetscObject)*newv);
369:   return(0);
370: }

372: /*@
373:    VecDestroy - Destroys a vector.

375:    Collective on Vec

377:    Input Parameters:
378: .  v  - the vector

380:    Level: beginner

382: .seealso: VecDuplicate(), VecDestroyVecs()
383: @*/
384: PetscErrorCode  VecDestroy(Vec *v)
385: {

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

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

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

407:    Collective on Vec

409:    Input Parameters:
410: +  m - the number of vectors to obtain
411: -  v - a vector to mimic

413:    Output Parameter:
414: .  V - location to put pointer to array of vectors

416:    Notes:
417:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
418:    vector.

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

425:    Level: intermediate

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

437:   (*v->ops->duplicatevecs)(v,m,V);
438:   return(0);
439: }

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

444:    Collective on Vec

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

450:    Fortran Note:
451:    The Fortran interface is slightly different from that given below.
452:    See the Fortran chapter of the users manual

454:    Level: intermediate

456: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
457: @*/
458: PetscErrorCode  VecDestroyVecs(PetscInt m,Vec *vv[])
459: {

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

473: /*@C
474:    VecViewFromOptions - View from Options

476:    Collective on Vec

478:    Input Parameters:
479: +  A - the vector
480: .  obj - Optional object
481: -  name - command line option

483:    Level: intermediate
484: .seealso:  Vec, VecView, PetscObjectViewFromOptions(), VecCreate()
485: @*/
486: PetscErrorCode  VecViewFromOptions(Vec A,PetscObject obj,const char name[])
487: {

492:   PetscObjectViewFromOptions((PetscObject)A,obj,name);
493:   return(0);
494: }

496: /*@C
497:    VecView - Views a vector object.

499:    Collective on Vec

501:    Input Parameters:
502: +  vec - the vector
503: -  viewer - an optional visualization context

505:    Notes:
506:    The available visualization contexts include
507: +     PETSC_VIEWER_STDOUT_SELF - for sequential vectors
508: .     PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
509: -     PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm

511:    You can change the format the vector is printed using the
512:    option PetscViewerPushFormat().

514:    The user can open alternative viewers with
515: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
516: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
517:          specified file; corresponding input uses VecLoad()
518: .    PetscViewerDrawOpen() - Outputs vector to an X window display
519: .    PetscViewerSocketOpen() - Outputs vector to Socket viewer
520: -    PetscViewerHDF5Open() - Outputs vector to HDF5 file viewer

522:    The user can call PetscViewerPushFormat() to specify the output
523:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
524:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
525: +    PETSC_VIEWER_DEFAULT - default, prints vector contents
526: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
527: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
528: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
529:          format common among all vector types

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

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

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

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


548:    Notes for HDF5 Viewer:
549:      The name of the Vec (given with PetscObjectSetName() is the name that is used
550:      for the object in the HDF5 file. If you wish to store the same Vec into multiple
551:      datasets in the same file (typically with different values), you must change its
552:      name each time before calling the VecView(). To load the same vector,
553:      the name of the Vec object passed to VecLoad() must be the same.

555:      If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
556:      If the function PetscViewerHDF5SetBaseDimension2()is called then even if the block size is one it will
557:      be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
558:      See also PetscViewerHDF5SetTimestep() which adds an additional complication to reading and writing Vecs
559:      with the HDF5 viewer.

561:    Level: beginner


564: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
565:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
566:           PetscRealView(), PetscScalarView(), PetscIntView(), PetscViewerHDF5SetTimestep()
567: @*/
568: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
569: {
570:   PetscErrorCode    ierr;
571:   PetscBool         iascii;
572:   PetscViewerFormat format;
573:   PetscMPIInt       size;

578:   if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);}
580:   PetscViewerGetFormat(viewer,&format);
581:   MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
582:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);

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

586:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
587:   if (iascii) {
588:     PetscInt rows,bs;

590:     PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
591:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
592:       PetscViewerASCIIPushTab(viewer);
593:       VecGetSize(vec,&rows);
594:       VecGetBlockSize(vec,&bs);
595:       if (bs != 1) {
596:         PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
597:       } else {
598:         PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
599:       }
600:       PetscViewerASCIIPopTab(viewer);
601:     }
602:   }
603:   VecLockReadPush(vec);
604:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
605:   if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
606:     (*vec->ops->viewnative)(vec,viewer);
607:   } else {
608:     (*vec->ops->view)(vec,viewer);
609:   }
610:   VecLockReadPop(vec);
611:   PetscLogEventEnd(VEC_View,vec,viewer,0,0);
612:   return(0);
613: }

615: #if defined(PETSC_USE_DEBUG)
616:  #include <../src/sys/totalview/tv_data_display.h>
617: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
618: {
619:   const PetscScalar *values;
620:   char              type[32];
621:   PetscErrorCode    ierr;


624:   TV_add_row("Local rows", "int", &v->map->n);
625:   TV_add_row("Global rows", "int", &v->map->N);
626:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
627:   VecGetArrayRead((Vec)v,&values);
628:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
629:   TV_add_row("values",type, values);
630:   VecRestoreArrayRead((Vec)v,&values);
631:   return TV_format_OK;
632: }
633: #endif

635: /*@
636:    VecGetSize - Returns the global number of elements of the vector.

638:    Not Collective

640:    Input Parameter:
641: .  x - the vector

643:    Output Parameters:
644: .  size - the global length of the vector

646:    Level: beginner

648: .seealso: VecGetLocalSize()
649: @*/
650: PetscErrorCode  VecGetSize(Vec x,PetscInt *size)
651: {

658:   (*x->ops->getsize)(x,size);
659:   return(0);
660: }

662: /*@
663:    VecGetLocalSize - Returns the number of elements of the vector stored
664:    in local memory. This routine may be implementation dependent, so use
665:    with care.

667:    Not Collective

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

672:    Output Parameter:
673: .  size - the length of the local piece of the vector

675:    Level: beginner

677: .seealso: VecGetSize()
678: @*/
679: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
680: {

687:   (*x->ops->getlocalsize)(x,size);
688:   return(0);
689: }

691: /*@C
692:    VecGetOwnershipRange - Returns the range of indices owned by
693:    this processor, assuming that the vectors are laid out with the
694:    first n1 elements on the first processor, next n2 elements on the
695:    second, etc.  For certain parallel layouts this range may not be
696:    well defined.

698:    Not Collective

700:    Input Parameter:
701: .  x - the vector

703:    Output Parameters:
704: +  low - the first local element, pass in NULL if not interested
705: -  high - one more than the last local element, pass in NULL if not interested

707:    Note:
708:    The high argument is one more than the last element stored locally.

710:    Fortran: PETSC_NULL_INTEGER should be used instead of NULL

712:    Level: beginner


715: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
716: @*/
717: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
718: {
724:   if (low)  *low  = x->map->rstart;
725:   if (high) *high = x->map->rend;
726:   return(0);
727: }

729: /*@C
730:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
731:    assuming that the vectors are laid out with the
732:    first n1 elements on the first processor, next n2 elements on the
733:    second, etc.  For certain parallel layouts this range may not be
734:    well defined.

736:    Not Collective

738:    Input Parameter:
739: .  x - the vector

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

744:    Note:
745:    The high argument is one more than the last element stored locally.

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

749:    Level: beginner


752: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
753: @*/
754: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
755: {

761:   PetscLayoutGetRanges(x->map,ranges);
762:   return(0);
763: }

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

768:    Collective on Vec

770:    Input Parameter:
771: +  x - the vector
772: .  op - the option
773: -  flag - turn the option on or off

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

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

794:    Level: intermediate

796: @*/
797: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool flag)
798: {

804:   if (x->ops->setoption) {
805:     (*x->ops->setoption)(x,op,flag);
806:   }
807:   return(0);
808: }

810: /* Default routines for obtaining and releasing; */
811: /* may be used by any implementation */
812: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
813: {
815:   PetscInt       i;

820:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
821:   PetscMalloc1(m,V);
822:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
823:   return(0);
824: }

826: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
827: {
829:   PetscInt       i;

833:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
834:   PetscFree(v);
835:   return(0);
836: }

838: /*@
839:    VecResetArray - Resets a vector to use its default memory. Call this
840:    after the use of VecPlaceArray().

842:    Not Collective

844:    Input Parameters:
845: .  vec - the vector

847:    Level: developer

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

851: @*/
852: PetscErrorCode  VecResetArray(Vec vec)
853: {

859:   if (vec->ops->resetarray) {
860:     (*vec->ops->resetarray)(vec);
861:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
862:   PetscObjectStateIncrease((PetscObject)vec);
863:   return(0);
864: }

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

870:   Collective on PetscViewer

872:   Input Parameters:
873: + vec - the newly loaded vector, this needs to have been created with VecCreate() or
874:            some related function before a call to VecLoad().
875: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
876:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

878:    Level: intermediate

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

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

887:   If the type or size of vec is not set before a call to VecLoad, PETSc
888:   sets the type and the local and global sizes. If type and/or
889:   sizes are already set, then the same are used.

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

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

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

906:   Notes for advanced users when using the binary viewer:
907:   Most users should not need to know the details of the binary storage
908:   format, since VecLoad() and VecView() completely hide these details.
909:   But for anyone who's interested, the standard binary vector storage
910:   format is
911: .vb
912:      PetscInt    VEC_FILE_CLASSID
913:      PetscInt    number of rows
914:      PetscScalar *values of all entries
915: .ve

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

922: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
923: @*/
924: PetscErrorCode  VecLoad(Vec vec, PetscViewer viewer)
925: {
926:   PetscErrorCode    ierr;
927:   PetscBool         isbinary,ishdf5,isadios,isadios2;
928:   PetscViewerFormat format;

934:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
935:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
936:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
937:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
938:   if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

940:   VecSetErrorIfLocked(vec,1);
941:   if (!((PetscObject)vec)->type_name && !vec->ops->create) {
942:     VecSetType(vec, VECSTANDARD);
943:   }
944:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
945:   PetscViewerGetFormat(viewer,&format);
946:   if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
947:     (*vec->ops->loadnative)(vec,viewer);
948:   } else {
949:     (*vec->ops->load)(vec,viewer);
950:   }
951:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
952:   return(0);
953: }


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

959:    Logically Collective on Vec

961:    Input Parameter:
962: .  vec - the vector

964:    Output Parameter:
965: .  vec - the vector reciprocal

967:    Level: intermediate

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

971: @*/
972: PetscErrorCode  VecReciprocal(Vec vec)
973: {

979:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
980:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
981:   VecSetErrorIfLocked(vec,1);
982:   (*vec->ops->reciprocal)(vec);
983:   PetscObjectStateIncrease((PetscObject)vec);
984:   return(0);
985: }

987: /*@C
988:     VecSetOperation - Allows user to set a vector operation.

990:    Logically Collective on Vec

992:     Input Parameters:
993: +   vec - the vector
994: .   op - the name of the operation
995: -   f - the function that provides the operation.

997:    Level: advanced

999:     Usage:
1000: $      PetscErrorCode userview(Vec,PetscViewer);
1001: $      VecCreateMPI(comm,m,M,&x);
1002: $      VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);

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

1010:     This function is not currently available from Fortran.

1012: .seealso: VecCreate(), MatShellSetOperation()
1013: @*/
1014: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1015: {
1018:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
1019:     vec->ops->viewnative = vec->ops->view;
1020:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1021:     vec->ops->loadnative = vec->ops->load;
1022:   }
1023:   (((void(**)(void))vec->ops)[(int)op]) = f;
1024:   return(0);
1025: }


1028: /*@
1029:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1030:    used during the assembly process to store values that belong to
1031:    other processors.

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

1035:    Input Parameters:
1036: +  vec   - the vector
1037: .  size  - the initial size of the stash.
1038: -  bsize - the initial size of the block-stash(if used).

1040:    Options Database Keys:
1041: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1042: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1044:    Level: intermediate

1046:    Notes:
1047:      The block-stash is used for values set with VecSetValuesBlocked() while
1048:      the stash is used for values set with VecSetValues()

1050:      Run with the option -info and look for output of the form
1051:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1052:      to determine the appropriate value, MM, to use for size and
1053:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1054:      to determine the value, BMM to use for bsize


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

1059: @*/
1060: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1061: {

1066:   VecStashSetInitialSize_Private(&vec->stash,size);
1067:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1068:   return(0);
1069: }

1071: /*@
1072:    VecConjugate - Conjugates a vector.

1074:    Logically Collective on Vec

1076:    Input Parameters:
1077: .  x - the vector

1079:    Level: intermediate

1081: @*/
1082: PetscErrorCode  VecConjugate(Vec x)
1083: {
1084: #if defined(PETSC_USE_COMPLEX)

1090:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1091:   VecSetErrorIfLocked(x,1);
1092:   (*x->ops->conjugate)(x);
1093:   /* we need to copy norms here */
1094:   PetscObjectStateIncrease((PetscObject)x);
1095:   return(0);
1096: #else
1097:   return(0);
1098: #endif
1099: }

1101: /*@
1102:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1104:    Logically Collective on Vec

1106:    Input Parameters:
1107: .  x, y  - the vectors

1109:    Output Parameter:
1110: .  w - the result

1112:    Level: advanced

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

1117: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1118: @*/
1119: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1120: {

1132:   VecCheckSameSize(w,1,x,2);
1133:   VecCheckSameSize(w,2,y,3);
1134:   VecSetErrorIfLocked(w,1);
1135:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1136:   (*w->ops->pointwisemult)(w,x,y);
1137:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1138:   PetscObjectStateIncrease((PetscObject)w);
1139:   return(0);
1140: }

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

1145:    Logically Collective on Vec

1147:    Input Parameters:
1148: +  x  - the vector
1149: -  rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1150:           it will create one internally.

1152:    Output Parameter:
1153: .  x  - the vector

1155:    Example of Usage:
1156: .vb
1157:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1158:      VecSetRandom(x,rctx);
1159:      PetscRandomDestroy(rctx);
1160: .ve

1162:    Level: intermediate


1165: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1166: @*/
1167: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1168: {
1170:   PetscRandom    randObj = NULL;

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

1179:   if (!rctx) {
1180:     MPI_Comm comm;
1181:     PetscObjectGetComm((PetscObject)x,&comm);
1182:     PetscRandomCreate(comm,&randObj);
1183:     PetscRandomSetFromOptions(randObj);
1184:     rctx = randObj;
1185:   }

1187:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1188:   (*x->ops->setrandom)(x,rctx);
1189:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1191:   PetscRandomDestroy(&randObj);
1192:   PetscObjectStateIncrease((PetscObject)x);
1193:   return(0);
1194: }

1196: /*@
1197:   VecZeroEntries - puts a 0.0 in each element of a vector

1199:   Logically Collective on Vec

1201:   Input Parameter:
1202: . vec - The vector

1204:   Level: beginner

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

1211: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1212: @*/
1213: PetscErrorCode  VecZeroEntries(Vec vec)
1214: {

1218:   VecSet(vec,0);
1219:   return(0);
1220: }

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

1226:   Collective on Vec

1228:   Input Parameter:
1229: . vec - The vector

1231:   Level: intermediate

1233: .seealso: VecSetFromOptions(), VecSetType()
1234: */
1235: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1236: {
1237:   PetscBool      opt;
1238:   VecType        defaultType;
1239:   char           typeName[256];
1240:   PetscMPIInt    size;

1244:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1245:   else {
1246:     MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1247:     if (size > 1) defaultType = VECMPI;
1248:     else defaultType = VECSEQ;
1249:   }

1251:   VecRegisterAll();
1252:   PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1253:   if (opt) {
1254:     VecSetType(vec, typeName);
1255:   } else {
1256:     VecSetType(vec, defaultType);
1257:   }
1258:   return(0);
1259: }

1261: /*@
1262:   VecSetFromOptions - Configures the vector from the options database.

1264:   Collective on Vec

1266:   Input Parameter:
1267: . vec - The vector

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

1273:   Level: beginner


1276: .seealso: VecCreate(), VecSetOptionsPrefix()
1277: @*/
1278: PetscErrorCode  VecSetFromOptions(Vec vec)
1279: {


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

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

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

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

1303:   Collective on Vec

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

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

1314:   Level: intermediate

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

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

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

1340:    Logically Collective on Vec

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

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

1349:    Level: advanced

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

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()


1387: @*/
1388: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1389: {

1395:   PetscLayoutGetBlockSize(v->map,bs);
1396:   return(0);
1397: }

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

1403:    Logically Collective on Vec

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

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

1413:    Level: advanced

1415: .seealso: VecSetFromOptions()
1416: @*/
1417: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1418: {

1423:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1424:   return(0);
1425: }

1427: /*@C
1428:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1429:    Vec options in the database.

1431:    Logically Collective on Vec

1433:    Input Parameters:
1434: +  v - the Vec context
1435: -  prefix - the prefix to prepend to all option names

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

1441:    Level: advanced

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

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

1455: /*@C
1456:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1457:    Vec options in the database.

1459:    Not Collective

1461:    Input Parameter:
1462: .  v - the Vec context

1464:    Output Parameter:
1465: .  prefix - pointer to the prefix string used

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

1471:    Level: advanced

1473: .seealso: VecAppendOptionsPrefix()
1474: @*/
1475: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1476: {

1481:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1482:   return(0);
1483: }

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

1488:    Collective on Vec

1490:    Input Parameters:
1491: .  v - the Vec context

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

1497:    Level: advanced

1499: .seealso: VecCreate(), VecDestroy()
1500: @*/
1501: PetscErrorCode  VecSetUp(Vec v)
1502: {
1503:   PetscMPIInt    size;

1508:   if (v->map->n < 0 && v->map->N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1509:   if (!((PetscObject)v)->type_name) {
1510:     MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1511:     if (size == 1) {
1512:       VecSetType(v, VECSEQ);
1513:     } else {
1514:       VecSetType(v, VECMPI);
1515:     }
1516:   }
1517:   return(0);
1518: }

1520: /*
1521:     These currently expose the PetscScalar/PetscReal in updating the
1522:     cached norm. If we push those down into the implementation these
1523:     will become independent of PetscScalar/PetscReal
1524: */

1526: /*@
1527:    VecCopy - Copies a vector. y <- x

1529:    Logically Collective on Vec

1531:    Input Parameter:
1532: .  x - the vector

1534:    Output Parameter:
1535: .  y - the copy

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

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

1546:    Level: beginner

1548: .seealso: VecDuplicate()
1549: @*/
1550: PetscErrorCode  VecCopy(Vec x,Vec y)
1551: {
1552:   PetscBool      flgs[4];
1553:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1555:   PetscInt       i;

1562:   if (x == y) return(0);
1563:   VecCheckSameLocalSize(x,1,y,2);
1564:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1565:   VecSetErrorIfLocked(y,2);

1567: #if !defined(PETSC_USE_MIXED_PRECISION)
1568:   for (i=0; i<4; i++) {
1569:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1570:   }
1571: #endif

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

1610:   PetscObjectStateIncrease((PetscObject)y);
1611: #if !defined(PETSC_USE_MIXED_PRECISION)
1612:   for (i=0; i<4; i++) {
1613:     if (flgs[i]) {
1614:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1615:     }
1616:   }
1617: #endif

1619:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1620:   return(0);
1621: }

1623: /*@
1624:    VecSwap - Swaps the vectors x and y.

1626:    Logically Collective on Vec

1628:    Input Parameters:
1629: .  x, y  - the vectors

1631:    Level: advanced

1633: @*/
1634: PetscErrorCode  VecSwap(Vec x,Vec y)
1635: {
1636:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1637:   PetscBool      flgxs[4],flgys[4];
1639:   PetscInt       i;

1647:   VecCheckSameSize(x,1,y,2);
1648:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1649:   if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1650:   VecSetErrorIfLocked(x,1);
1651:   VecSetErrorIfLocked(y,2);

1653:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1654:   for (i=0; i<4; i++) {
1655:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1656:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1657:   }
1658:   (*x->ops->swap)(x,y);
1659:   PetscObjectStateIncrease((PetscObject)x);
1660:   PetscObjectStateIncrease((PetscObject)y);
1661:   for (i=0; i<4; i++) {
1662:     if (flgxs[i]) {
1663:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1664:     }
1665:     if (flgys[i]) {
1666:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1667:     }
1668:   }
1669:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1670:   return(0);
1671: }

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

1676:   Collective on VecStash

1678:   Input Parameters:
1679: + obj   - the VecStash object
1680: . bobj - optional other object that provides the prefix
1681: - optionname - option to activate viewing

1683:   Level: intermediate

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

1687: */
1688: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1689: {
1690:   PetscErrorCode    ierr;
1691:   PetscViewer       viewer;
1692:   PetscBool         flg;
1693:   PetscViewerFormat format;
1694:   char              *prefix;

1697:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1698:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1699:   if (flg) {
1700:     PetscViewerPushFormat(viewer,format);
1701:     VecStashView(obj,viewer);
1702:     PetscViewerPopFormat(viewer);
1703:     PetscViewerDestroy(&viewer);
1704:   }
1705:   return(0);
1706: }

1708: /*@
1709:    VecStashView - Prints the entries in the vector stash and block stash.

1711:    Collective on Vec

1713:    Input Parameters:
1714: +  v - the vector
1715: -  viewer - the viewer

1717:    Level: advanced


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

1722: @*/
1723: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1724: {
1726:   PetscMPIInt    rank;
1727:   PetscInt       i,j;
1728:   PetscBool      match;
1729:   VecStash       *s;
1730:   PetscScalar    val;


1737:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1738:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1739:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1740:   MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1741:   s    = &v->bstash;

1743:   /* print block stash */
1744:   PetscViewerASCIIPushSynchronized(viewer);
1745:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1746:   for (i=0; i<s->n; i++) {
1747:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1748:     for (j=0; j<s->bs; j++) {
1749:       val = s->array[i*s->bs+j];
1750: #if defined(PETSC_USE_COMPLEX)
1751:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1752: #else
1753:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1754: #endif
1755:     }
1756:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1757:   }
1758:   PetscViewerFlush(viewer);

1760:   s = &v->stash;

1762:   /* print basic stash */
1763:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1764:   for (i=0; i<s->n; i++) {
1765:     val = s->array[i];
1766: #if defined(PETSC_USE_COMPLEX)
1767:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1768: #else
1769:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1770: #endif
1771:   }
1772:   PetscViewerFlush(viewer);
1773:   PetscViewerASCIIPopSynchronized(viewer);
1774:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1775:   return(0);
1776: }

1778: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1779: {
1780:   PetscInt       i,N,rstart,rend;
1782:   PetscScalar    *xx;
1783:   PetscReal      *xreal;
1784:   PetscBool      iset;

1787:   VecGetOwnershipRange(v,&rstart,&rend);
1788:   VecGetSize(v,&N);
1789:   PetscCalloc1(N,&xreal);
1790:   PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1791:   if (iset) {
1792:     VecGetArray(v,&xx);
1793:     for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1794:     VecRestoreArray(v,&xx);
1795:   }
1796:   PetscFree(xreal);
1797:   if (set) *set = iset;
1798:   return(0);
1799: }

1801: /*@
1802:    VecGetLayout - get PetscLayout describing vector layout

1804:    Not Collective

1806:    Input Arguments:
1807: .  x - the vector

1809:    Output Arguments:
1810: .  map - the layout

1812:    Level: developer

1814: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1815: @*/
1816: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1817: {

1821:   *map = x->map;
1822:   return(0);
1823: }

1825: /*@
1826:    VecSetLayout - set PetscLayout describing vector layout

1828:    Not Collective

1830:    Input Arguments:
1831: +  x - the vector
1832: -  map - the layout

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

1837:    Level: developer

1839: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1840: @*/
1841: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1842: {

1847:   PetscLayoutReference(map,&x->map);
1848:   return(0);
1849: }

1851: PetscErrorCode VecSetInf(Vec xin)
1852: {
1853:   PetscInt       i,n = xin->map->n;
1854:   PetscScalar    *xx;
1855:   PetscScalar    zero=0.0,one=1.0,inf=one/zero;

1859:   VecGetArray(xin,&xx);
1860:   for (i=0; i<n; i++) xx[i] = inf;
1861:   VecRestoreArray(xin,&xx);
1862:   return(0);
1863: }

1865: /*@
1866:      VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU

1868:    Input Parameters:
1869: +   v - the vector
1870: -   flg - bind to the CPU if value of PETSC_TRUE

1872:    Level: intermediate
1873: @*/
1874: PetscErrorCode VecBindToCPU(Vec v,PetscBool flg)
1875: {
1876: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)

1880:   if (v->boundtocpu == flg) return(0);
1881:   v->boundtocpu = flg;
1882:   if (v->ops->bindtocpu) {
1883:     (*v->ops->bindtocpu)(v,flg);
1884:   }
1885:   return(0);
1886: #else
1887:   return 0;
1888: #endif
1889: }

1891: /*@C
1892:   VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.

1894:   Logically Collective on Vec

1896:   Input Parameters:
1897: +  v    - the vector
1898: -  mbytes - minimum data size in bytes

1900:   Options Database Keys:

1902: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.
1903:                                   Note that this takes a PetscScalar, to accommodate large values;
1904:                                   specifying -1 ensures that pinned memory will always be used.

1906:   Level: developer

1908: .seealso: VecGetPinnedMemoryMin()
1909: @*/
1910: PetscErrorCode VecSetPinnedMemoryMin(Vec v,size_t mbytes)
1911: {
1912: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1914:   v->minimum_bytes_pinned_memory = mbytes;
1915:   return(0);
1916: #else
1917:   return 0;
1918: #endif
1919: }

1921: /*@C
1922:   VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.

1924:   Logically Collective on Vec

1926:   Input Parameters:
1927: .  v    - the vector

1929:   Output Parameters:
1930: .  mbytes - minimum data size in bytes

1932:   Level: developer

1934: .seealso: VecSetPinnedMemoryMin()
1935: @*/
1936: PetscErrorCode VecGetPinnedMemoryMin(Vec v,size_t *mbytes)
1937: {
1938: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1940:   *mbytes = v->minimum_bytes_pinned_memory;
1941:   return(0);
1942: #else
1943:   return 0;
1944: #endif
1945: }