Actual source code: vector.c

petsc-master 2019-06-17
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: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
 67:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
 68: @*/
 69: PetscErrorCode  VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
 70: {


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

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

 88:    Not Collective

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

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

 96:    Level: advanced


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

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

115:    Collective on Vec

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

120:    Level: beginner

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

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

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

145:    Collective on Vec

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

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

159:    Level: beginner

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

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

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

182:    Logically Collective on Vec

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

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

190:    Level: advanced

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

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

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


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

222:    Logically Collective on Vec

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

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

230:    Level: advanced

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

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

251:   VecCheckSameSize(w,1,x,2);
252:   VecCheckSameSize(w,1,y,3);
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:   (*w->ops->pointwisemaxabs)(w,x,y);
292:   PetscObjectStateIncrease((PetscObject)w);
293:   return(0);
294: }

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

299:    Logically Collective on Vec

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

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

307:    Level: advanced

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

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

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


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

338:    Collective on Vec

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

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

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

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

353:    Level: beginner

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

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

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

373:    Collective on Vec

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

378:    Level: beginner

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

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

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

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

405:    Collective on Vec

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

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

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

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

423:    Level: intermediate

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

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

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

442:    Collective on Vec

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

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

452:    Level: intermediate

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

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

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

474:    Collective on Vec

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

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

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

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

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

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

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

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

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

525:    Level: beginner


528: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
529:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
530:           PetscRealView(), PetscScalarView(), PetscIntView()
531: @*/
532: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
533: {
534:   PetscErrorCode    ierr;
535:   PetscBool         iascii;
536:   PetscViewerFormat format;
537:   PetscMPIInt       size;

542:   if (!viewer) {
543:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
544:   }
546:   PetscViewerGetFormat(viewer,&format);
547:   MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
548:   if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);

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

552:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
553:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
554:   if (iascii) {
555:     PetscInt rows,bs;

557:     PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
558:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
559:       PetscViewerASCIIPushTab(viewer);
560:       VecGetSize(vec,&rows);
561:       VecGetBlockSize(vec,&bs);
562:       if (bs != 1) {
563:         PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
564:       } else {
565:         PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
566:       }
567:       PetscViewerASCIIPopTab(viewer);
568:     }
569:   }
570:   VecLockReadPush(vec);
571:   if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
572:     (*vec->ops->viewnative)(vec,viewer);
573:   } else {
574:     (*vec->ops->view)(vec,viewer);
575:   }
576:   VecLockReadPop(vec);
577:   PetscLogEventEnd(VEC_View,vec,viewer,0,0);
578:   return(0);
579: }

581: #if defined(PETSC_USE_DEBUG)
582: #include <../src/sys/totalview/tv_data_display.h>
583: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
584: {
585:   const PetscScalar *values;
586:   char              type[32];
587:   PetscErrorCode    ierr;


590:   TV_add_row("Local rows", "int", &v->map->n);
591:   TV_add_row("Global rows", "int", &v->map->N);
592:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
593:   VecGetArrayRead((Vec)v,&values);
594:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
595:   TV_add_row("values",type, values);
596:   VecRestoreArrayRead((Vec)v,&values);
597:   return TV_format_OK;
598: }
599: #endif

601: /*@
602:    VecGetSize - Returns the global number of elements of the vector.

604:    Not Collective

606:    Input Parameter:
607: .  x - the vector

609:    Output Parameters:
610: .  size - the global length of the vector

612:    Level: beginner

614: .seealso: VecGetLocalSize()
615: @*/
616: PetscErrorCode  VecGetSize(Vec x,PetscInt *size)
617: {

624:   (*x->ops->getsize)(x,size);
625:   return(0);
626: }

628: /*@
629:    VecGetLocalSize - Returns the number of elements of the vector stored
630:    in local memory. This routine may be implementation dependent, so use
631:    with care.

633:    Not Collective

635:    Input Parameter:
636: .  x - the vector

638:    Output Parameter:
639: .  size - the length of the local piece of the vector

641:    Level: beginner

643: .seealso: VecGetSize()
644: @*/
645: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
646: {

653:   (*x->ops->getlocalsize)(x,size);
654:   return(0);
655: }

657: /*@C
658:    VecGetOwnershipRange - Returns the range of indices owned by
659:    this processor, assuming that the vectors are laid out with the
660:    first n1 elements on the first processor, next n2 elements on the
661:    second, etc.  For certain parallel layouts this range may not be
662:    well defined.

664:    Not Collective

666:    Input Parameter:
667: .  x - the vector

669:    Output Parameters:
670: +  low - the first local element, pass in NULL if not interested
671: -  high - one more than the last local element, pass in NULL if not interested

673:    Note:
674:    The high argument is one more than the last element stored locally.

676:    Fortran: PETSC_NULL_INTEGER should be used instead of NULL

678:    Level: beginner


681: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
682: @*/
683: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
684: {
690:   if (low)  *low  = x->map->rstart;
691:   if (high) *high = x->map->rend;
692:   return(0);
693: }

695: /*@C
696:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
697:    assuming that the vectors are laid out with the
698:    first n1 elements on the first processor, next n2 elements on the
699:    second, etc.  For certain parallel layouts this range may not be
700:    well defined.

702:    Not Collective

704:    Input Parameter:
705: .  x - the vector

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

710:    Note:
711:    The high argument is one more than the last element stored locally.

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

715:    Level: beginner


718: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
719: @*/
720: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
721: {

727:   PetscLayoutGetRanges(x->map,ranges);
728:   return(0);
729: }

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

734:    Collective on Vec

736:    Input Parameter:
737: +  x - the vector
738: .  op - the option
739: -  flag - turn the option on or off

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

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

760:    Level: intermediate

762: @*/
763: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool flag)
764: {

770:   if (x->ops->setoption) {
771:     (*x->ops->setoption)(x,op,flag);
772:   }
773:   return(0);
774: }

776: /* Default routines for obtaining and releasing; */
777: /* may be used by any implementation */
778: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
779: {
781:   PetscInt       i;

786:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
787:   PetscMalloc1(m,V);
788:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
789:   return(0);
790: }

792: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
793: {
795:   PetscInt       i;

799:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
800:   PetscFree(v);
801:   return(0);
802: }

804: /*@
805:    VecResetArray - Resets a vector to use its default memory. Call this
806:    after the use of VecPlaceArray().

808:    Not Collective

810:    Input Parameters:
811: .  vec - the vector

813:    Level: developer

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

817: @*/
818: PetscErrorCode  VecResetArray(Vec vec)
819: {

825:   if (vec->ops->resetarray) {
826:     (*vec->ops->resetarray)(vec);
827:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
828:   PetscObjectStateIncrease((PetscObject)vec);
829:   return(0);
830: }

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

836:   Collective on PetscViewer

838:   Input Parameters:
839: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
840:            some related function before a call to VecLoad().
841: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
842:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

844:    Level: intermediate

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

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

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

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

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

867:   Notes for advanced users:
868:   Most users should not need to know the details of the binary storage
869:   format, since VecLoad() and VecView() completely hide these details.
870:   But for anyone who's interested, the standard binary vector storage
871:   format is
872: .vb
873:      int    VEC_FILE_CLASSID
874:      int    number of rows
875:      PetscScalar *values of all entries
876: .ve

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

884: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
885: @*/
886: PetscErrorCode  VecLoad(Vec newvec, PetscViewer viewer)
887: {
888:   PetscErrorCode    ierr;
889:   PetscBool         isbinary,ishdf5,isadios,isadios2;
890:   PetscViewerFormat format;

895:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
896:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
897:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
898:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
899:   if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

901:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
902:   if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
903:     VecSetType(newvec, VECSTANDARD);
904:   }
905:   PetscViewerGetFormat(viewer,&format);
906:   if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
907:     (*newvec->ops->loadnative)(newvec,viewer);
908:   } else {
909:     (*newvec->ops->load)(newvec,viewer);
910:   }
911:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
912:   return(0);
913: }


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

919:    Logically Collective on Vec

921:    Input Parameter:
922: .  vec - the vector

924:    Output Parameter:
925: .  vec - the vector reciprocal

927:    Level: intermediate

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

931: @*/
932: PetscErrorCode  VecReciprocal(Vec vec)
933: {

939:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
940:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
941:   (*vec->ops->reciprocal)(vec);
942:   PetscObjectStateIncrease((PetscObject)vec);
943:   return(0);
944: }

946: /*@C
947:     VecSetOperation - Allows user to set a vector operation.

949:    Logically Collective on Vec

951:     Input Parameters:
952: +   vec - the vector
953: .   op - the name of the operation
954: -   f - the function that provides the operation.

956:    Level: advanced

958:     Usage:
959: $      PetscErrorCode userview(Vec,PetscViewer);
960: $      VecCreateMPI(comm,m,M,&x);
961: $      VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);

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

969:     This function is not currently available from Fortran.

971: .seealso: VecCreate(), MatShellSetOperation()
972: @*/
973: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
974: {
977:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
978:     vec->ops->viewnative = vec->ops->view;
979:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
980:     vec->ops->loadnative = vec->ops->load;
981:   }
982:   (((void(**)(void))vec->ops)[(int)op]) = f;
983:   return(0);
984: }


987: /*@
988:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
989:    used during the assembly process to store values that belong to
990:    other processors.

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

994:    Input Parameters:
995: +  vec   - the vector
996: .  size  - the initial size of the stash.
997: -  bsize - the initial size of the block-stash(if used).

999:    Options Database Keys:
1000: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1001: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1003:    Level: intermediate

1005:    Notes:
1006:      The block-stash is used for values set with VecSetValuesBlocked() while
1007:      the stash is used for values set with VecSetValues()

1009:      Run with the option -info and look for output of the form
1010:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1011:      to determine the appropriate value, MM, to use for size and
1012:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1013:      to determine the value, BMM to use for bsize


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

1018: @*/
1019: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1020: {

1025:   VecStashSetInitialSize_Private(&vec->stash,size);
1026:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1027:   return(0);
1028: }

1030: /*@
1031:    VecConjugate - Conjugates a vector.

1033:    Logically Collective on Vec

1035:    Input Parameters:
1036: .  x - the vector

1038:    Level: intermediate

1040: @*/
1041: PetscErrorCode  VecConjugate(Vec x)
1042: {
1043: #if defined(PETSC_USE_COMPLEX)

1049:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1050:   (*x->ops->conjugate)(x);
1051:   /* we need to copy norms here */
1052:   PetscObjectStateIncrease((PetscObject)x);
1053:   return(0);
1054: #else
1055:   return(0);
1056: #endif
1057: }

1059: /*@
1060:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1062:    Logically Collective on Vec

1064:    Input Parameters:
1065: .  x, y  - the vectors

1067:    Output Parameter:
1068: .  w - the result

1070:    Level: advanced

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

1075: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1076: @*/
1077: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1078: {

1090:   VecCheckSameSize(w,1,x,2);
1091:   VecCheckSameSize(w,2,y,3);
1092:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1093:   (*w->ops->pointwisemult)(w,x,y);
1094:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1095:   PetscObjectStateIncrease((PetscObject)w);
1096:   return(0);
1097: }

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

1102:    Logically Collective on Vec

1104:    Input Parameters:
1105: +  x  - the vector
1106: -  rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1107:           it will create one internally.

1109:    Output Parameter:
1110: .  x  - the vector

1112:    Example of Usage:
1113: .vb
1114:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1115:      VecSetRandom(x,rctx);
1116:      PetscRandomDestroy(rctx);
1117: .ve

1119:    Level: intermediate


1122: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1123: @*/
1124: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1125: {
1127:   PetscRandom    randObj = NULL;

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

1135:   if (!rctx) {
1136:     MPI_Comm comm;
1137:     PetscObjectGetComm((PetscObject)x,&comm);
1138:     PetscRandomCreate(comm,&randObj);
1139:     PetscRandomSetFromOptions(randObj);
1140:     rctx = randObj;
1141:   }

1143:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1144:   (*x->ops->setrandom)(x,rctx);
1145:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1147:   PetscRandomDestroy(&randObj);
1148:   PetscObjectStateIncrease((PetscObject)x);
1149:   return(0);
1150: }

1152: /*@
1153:   VecZeroEntries - puts a 0.0 in each element of a vector

1155:   Logically Collective on Vec

1157:   Input Parameter:
1158: . vec - The vector

1160:   Level: beginner

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

1167: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1168: @*/
1169: PetscErrorCode  VecZeroEntries(Vec vec)
1170: {

1174:   VecSet(vec,0);
1175:   return(0);
1176: }

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

1182:   Collective on Vec

1184:   Input Parameter:
1185: . vec - The vector

1187:   Level: intermediate

1189: .seealso: VecSetFromOptions(), VecSetType()
1190: */
1191: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1192: {
1193:   PetscBool      opt;
1194:   VecType        defaultType;
1195:   char           typeName[256];
1196:   PetscMPIInt    size;

1200:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1201:   else {
1202:     MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1203:     if (size > 1) defaultType = VECMPI;
1204:     else defaultType = VECSEQ;
1205:   }

1207:   VecRegisterAll();
1208:   PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1209:   if (opt) {
1210:     VecSetType(vec, typeName);
1211:   } else {
1212:     VecSetType(vec, defaultType);
1213:   }
1214:   return(0);
1215: }

1217: /*@
1218:   VecSetFromOptions - Configures the vector from the options database.

1220:   Collective on Vec

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

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

1229:   Level: beginner


1232: .seealso: VecCreate(), VecSetOptionsPrefix()
1233: @*/
1234: PetscErrorCode  VecSetFromOptions(Vec vec)
1235: {


1241:   PetscObjectOptionsBegin((PetscObject)vec);
1242:   /* Handle vector type options */
1243:   VecSetTypeFromOptions_Private(PetscOptionsObject,vec);

1245:   /* Handle specific vector options */
1246:   if (vec->ops->setfromoptions) {
1247:     (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1248:   }

1250:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1251:   PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1252:   PetscOptionsEnd();
1253:   return(0);
1254: }

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

1259:   Collective on Vec

1261:   Input Parameters:
1262: + v - the vector
1263: . n - the local size (or PETSC_DECIDE to have it set)
1264: - N - the global size (or PETSC_DECIDE)

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

1270:   Level: intermediate

1272: .seealso: VecGetSize(), PetscSplitOwnership()
1273: @*/
1274: PetscErrorCode  VecSetSizes(Vec v, PetscInt n, PetscInt N)
1275: {

1281:   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);
1282:   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);
1283:   v->map->n = n;
1284:   v->map->N = N;
1285:   if (v->ops->create) {
1286:     (*v->ops->create)(v);
1287:     v->ops->create = 0;
1288:   }
1289:   return(0);
1290: }

1292: /*@
1293:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1294:    and VecSetValuesBlockedLocal().

1296:    Logically Collective on Vec

1298:    Input Parameter:
1299: +  v - the vector
1300: -  bs - the blocksize

1302:    Notes:
1303:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1305:    Level: advanced

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

1309: @*/
1310: PetscErrorCode  VecSetBlockSize(Vec v,PetscInt bs)
1311: {

1316:   if (bs < 0 || bs == v->map->bs) return(0);
1318:   PetscLayoutSetBlockSize(v->map,bs);
1319:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1320:   return(0);
1321: }

1323: /*@
1324:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1325:    and VecSetValuesBlockedLocal().

1327:    Not Collective

1329:    Input Parameter:
1330: .  v - the vector

1332:    Output Parameter:
1333: .  bs - the blocksize

1335:    Notes:
1336:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1338:    Level: advanced

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


1343: @*/
1344: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1345: {

1351:   PetscLayoutGetBlockSize(v->map,bs);
1352:   return(0);
1353: }

1355: /*@C
1356:    VecSetOptionsPrefix - Sets the prefix used for searching for all
1357:    Vec options in the database.

1359:    Logically Collective on Vec

1361:    Input Parameter:
1362: +  v - the Vec context
1363: -  prefix - the prefix to prepend to all option names

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

1369:    Level: advanced

1371: .seealso: VecSetFromOptions()
1372: @*/
1373: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1374: {

1379:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1380:   return(0);
1381: }

1383: /*@C
1384:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1385:    Vec options in the database.

1387:    Logically Collective on Vec

1389:    Input Parameters:
1390: +  v - the Vec context
1391: -  prefix - the prefix to prepend to all option names

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

1397:    Level: advanced

1399: .seealso: VecGetOptionsPrefix()
1400: @*/
1401: PetscErrorCode  VecAppendOptionsPrefix(Vec v,const char prefix[])
1402: {

1407:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1408:   return(0);
1409: }

1411: /*@C
1412:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1413:    Vec options in the database.

1415:    Not Collective

1417:    Input Parameter:
1418: .  v - the Vec context

1420:    Output Parameter:
1421: .  prefix - pointer to the prefix string used

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

1427:    Level: advanced

1429: .seealso: VecAppendOptionsPrefix()
1430: @*/
1431: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1432: {

1437:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1438:   return(0);
1439: }

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

1444:    Collective on Vec

1446:    Input Parameters:
1447: .  v - the Vec context

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

1453:    Level: advanced

1455: .seealso: VecCreate(), VecDestroy()
1456: @*/
1457: PetscErrorCode  VecSetUp(Vec v)
1458: {
1459:   PetscMPIInt    size;

1464:   if (!((PetscObject)v)->type_name) {
1465:     MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1466:     if (size == 1) {
1467:       VecSetType(v, VECSEQ);
1468:     } else {
1469:       VecSetType(v, VECMPI);
1470:     }
1471:   }
1472:   return(0);
1473: }

1475: /*
1476:     These currently expose the PetscScalar/PetscReal in updating the
1477:     cached norm. If we push those down into the implementation these
1478:     will become independent of PetscScalar/PetscReal
1479: */

1481: /*@
1482:    VecCopy - Copies a vector. y <- x

1484:    Logically Collective on Vec

1486:    Input Parameter:
1487: .  x - the vector

1489:    Output Parameter:
1490: .  y - the copy

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

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

1501:    Level: beginner

1503: .seealso: VecDuplicate()
1504: @*/
1505: PetscErrorCode  VecCopy(Vec x,Vec y)
1506: {
1507:   PetscBool      flgs[4];
1508:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1510:   PetscInt       i;

1517:   if (x == y) return(0);
1518:   VecCheckSameLocalSize(x,1,y,2);
1519:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1520:   VecSetErrorIfLocked(y,2);

1522: #if !defined(PETSC_USE_MIXED_PRECISION)
1523:   for (i=0; i<4; i++) {
1524:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1525:   }
1526: #endif

1528:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
1529: #if defined(PETSC_USE_MIXED_PRECISION)
1530:   extern PetscErrorCode VecGetArray(Vec,double**);
1531:   extern PetscErrorCode VecRestoreArray(Vec,double**);
1532:   extern PetscErrorCode VecGetArray(Vec,float**);
1533:   extern PetscErrorCode VecRestoreArray(Vec,float**);
1534:   extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1535:   extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1536:   extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1537:   extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1538:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1539:     PetscInt    i,n;
1540:     const float *xx;
1541:     double      *yy;
1542:     VecGetArrayRead(x,&xx);
1543:     VecGetArray(y,&yy);
1544:     VecGetLocalSize(x,&n);
1545:     for (i=0; i<n; i++) yy[i] = xx[i];
1546:     VecRestoreArrayRead(x,&xx);
1547:     VecRestoreArray(y,&yy);
1548:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1549:     PetscInt     i,n;
1550:     float        *yy;
1551:     const double *xx;
1552:     VecGetArrayRead(x,&xx);
1553:     VecGetArray(y,&yy);
1554:     VecGetLocalSize(x,&n);
1555:     for (i=0; i<n; i++) yy[i] = (float) xx[i];
1556:     VecRestoreArrayRead(x,&xx);
1557:     VecRestoreArray(y,&yy);
1558:   } else {
1559:     (*x->ops->copy)(x,y);
1560:   }
1561: #else
1562:   (*x->ops->copy)(x,y);
1563: #endif

1565:   PetscObjectStateIncrease((PetscObject)y);
1566: #if !defined(PETSC_USE_MIXED_PRECISION)
1567:   for (i=0; i<4; i++) {
1568:     if (flgs[i]) {
1569:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1570:     }
1571:   }
1572: #endif

1574:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1575:   return(0);
1576: }

1578: /*@
1579:    VecSwap - Swaps the vectors x and y.

1581:    Logically Collective on Vec

1583:    Input Parameters:
1584: .  x, y  - the vectors

1586:    Level: advanced

1588: @*/
1589: PetscErrorCode  VecSwap(Vec x,Vec y)
1590: {
1591:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1592:   PetscBool      flgxs[4],flgys[4];
1594:   PetscInt       i;

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

1606:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1607:   for (i=0; i<4; i++) {
1608:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1609:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1610:   }
1611:   (*x->ops->swap)(x,y);
1612:   PetscObjectStateIncrease((PetscObject)x);
1613:   PetscObjectStateIncrease((PetscObject)y);
1614:   for (i=0; i<4; i++) {
1615:     if (flgxs[i]) {
1616:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1617:     }
1618:     if (flgys[i]) {
1619:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1620:     }
1621:   }
1622:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1623:   return(0);
1624: }

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

1629:   Collective on VecStash

1631:   Input Parameters:
1632: + obj   - the VecStash object
1633: . bobj - optional other object that provides the prefix
1634: - optionname - option to activate viewing

1636:   Level: intermediate

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

1640: */
1641: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1642: {
1643:   PetscErrorCode    ierr;
1644:   PetscViewer       viewer;
1645:   PetscBool         flg;
1646:   PetscViewerFormat format;
1647:   char              *prefix;

1650:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1651:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1652:   if (flg) {
1653:     PetscViewerPushFormat(viewer,format);
1654:     VecStashView(obj,viewer);
1655:     PetscViewerPopFormat(viewer);
1656:     PetscViewerDestroy(&viewer);
1657:   }
1658:   return(0);
1659: }

1661: /*@
1662:    VecStashView - Prints the entries in the vector stash and block stash.

1664:    Collective on Vec

1666:    Input Parameters:
1667: +  v - the vector
1668: -  viewer - the viewer

1670:    Level: advanced


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

1675: @*/
1676: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1677: {
1679:   PetscMPIInt    rank;
1680:   PetscInt       i,j;
1681:   PetscBool      match;
1682:   VecStash       *s;
1683:   PetscScalar    val;


1690:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1691:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1692:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1693:   MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1694:   s    = &v->bstash;

1696:   /* print block stash */
1697:   PetscViewerASCIIPushSynchronized(viewer);
1698:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1699:   for (i=0; i<s->n; i++) {
1700:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1701:     for (j=0; j<s->bs; j++) {
1702:       val = s->array[i*s->bs+j];
1703: #if defined(PETSC_USE_COMPLEX)
1704:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1705: #else
1706:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1707: #endif
1708:     }
1709:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1710:   }
1711:   PetscViewerFlush(viewer);

1713:   s = &v->stash;

1715:   /* print basic stash */
1716:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1717:   for (i=0; i<s->n; i++) {
1718:     val = s->array[i];
1719: #if defined(PETSC_USE_COMPLEX)
1720:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1721: #else
1722:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1723: #endif
1724:   }
1725:   PetscViewerFlush(viewer);
1726:   PetscViewerASCIIPopSynchronized(viewer);
1727:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1728:   return(0);
1729: }

1731: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1732: {
1733:   PetscInt       i,N,rstart,rend;
1735:   PetscScalar    *xx;
1736:   PetscReal      *xreal;
1737:   PetscBool      iset;

1740:   VecGetOwnershipRange(v,&rstart,&rend);
1741:   VecGetSize(v,&N);
1742:   PetscCalloc1(N,&xreal);
1743:   PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1744:   if (iset) {
1745:     VecGetArray(v,&xx);
1746:     for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1747:     VecRestoreArray(v,&xx);
1748:   }
1749:   PetscFree(xreal);
1750:   if (set) *set = iset;
1751:   return(0);
1752: }

1754: /*@
1755:    VecGetLayout - get PetscLayout describing vector layout

1757:    Not Collective

1759:    Input Arguments:
1760: .  x - the vector

1762:    Output Arguments:
1763: .  map - the layout

1765:    Level: developer

1767: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1768: @*/
1769: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1770: {

1774:   *map = x->map;
1775:   return(0);
1776: }

1778: /*@
1779:    VecSetLayout - set PetscLayout describing vector layout

1781:    Not Collective

1783:    Input Arguments:
1784: +  x - the vector
1785: -  map - the layout

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

1790:    Level: developer

1792: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1793: @*/
1794: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1795: {

1800:   PetscLayoutReference(map,&x->map);
1801:   return(0);
1802: }

1804: PetscErrorCode VecSetInf(Vec xin)
1805: {
1806:   PetscInt       i,n = xin->map->n;
1807:   PetscScalar    *xx;
1808:   PetscScalar    zero=0.0,one=1.0,inf=one/zero;

1812:   VecGetArray(xin,&xx);
1813:   for (i=0; i<n; i++) xx[i] = inf;
1814:   VecRestoreArray(xin,&xx);
1815:   return(0);
1816: }

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

1821:    Input Parameters:
1822: +   v - the vector
1823: -   flg - pin to the CPU if value of PETSC_TRUE

1825: @*/
1826: PetscErrorCode VecPinToCPU(Vec v,PetscBool flg)
1827: {
1829: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)

1832:   if (v->pinnedtocpu == flg) return 0;
1833:   v->pinnedtocpu = flg;
1834:   if (v->ops->pintocpu) {
1835:     (*v->ops->pintocpu)(v,flg);
1836:   }
1837: #endif
1838:   return(0);
1839: }