Actual source code: vector.c

petsc-3.6.0 2015-06-09
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>    /*I  "petscvec.h"   I*/

  8: /* Logging support */
  9: PetscClassId  VEC_CLASSID;
 10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_DotBarrier, VEC_Dot, VEC_MDotBarrier, 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_NormBarrier, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
 13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load, VEC_ScatterBarrier;
 14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceBarrier, VEC_ReduceCommunication,VEC_ReduceBegin,VEC_ReduceEnd,VEC_Ops;
 15: PetscLogEvent VEC_DotNormBarrier, VEC_DotNorm, VEC_AXPBYPCZ, VEC_CUSPCopyFromGPU, VEC_CUSPCopyToGPU;
 16: PetscLogEvent VEC_CUSPCopyFromGPUSome, VEC_CUSPCopyToGPUSome;
 17: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;

 19: extern PetscErrorCode VecStashGetInfo_Private(VecStash*,PetscInt*,PetscInt*);
 22: /*@
 23:    VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
 24:        to be communicated to other processors during the VecAssemblyBegin/End() process

 26:     Not collective

 28:    Input Parameter:
 29: .   vec - the vector

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

 37:    Level: advanced

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

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

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

 54: /*@
 55:    VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
 56:    by the routine VecSetValuesLocal() to allow users to insert vector entries
 57:    using a local (per-processor) numbering.

 59:    Logically Collective on Vec

 61:    Input Parameters:
 62: +  x - vector
 63: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

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

 68:    Level: intermediate

 70:    Concepts: vector^setting values with local numbering

 72: seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
 73:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
 74: @*/
 75: PetscErrorCode  VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
 76: {


 83:   if (x->ops->setlocaltoglobalmapping) {
 84:     (*x->ops->setlocaltoglobalmapping)(x,mapping);
 85:   } else {
 86:     PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
 87:   }
 88:   return(0);
 89: }

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

 96:    Not Collective

 98:    Input Parameter:
 99: .  X - the vector

101:    Output Parameter:
102: .  mapping - the mapping

104:    Level: advanced

106:    Concepts: vectors^local to global mapping
107:    Concepts: local to global mapping^for vectors

109: .seealso:  VecSetValuesLocal()
110: @*/
111: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
112: {
117:   *mapping = X->map->mapping;
118:   return(0);
119: }

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

127:    Collective on Vec

129:    Input Parameter:
130: .  vec - the vector

132:    Level: beginner

134:    Concepts: assembly^vectors

136: .seealso: VecAssemblyEnd(), VecSetValues()
137: @*/
138: PetscErrorCode  VecAssemblyBegin(Vec vec)
139: {

145:   VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
146:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
147:   if (vec->ops->assemblybegin) {
148:     (*vec->ops->assemblybegin)(vec);
149:   }
150:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
151:   PetscObjectStateIncrease((PetscObject)vec);
152:   return(0);
153: }

157: /*@
158:    VecAssemblyEnd - Completes assembling the vector.  This routine should
159:    be called after VecAssemblyBegin().

161:    Collective on Vec

163:    Input Parameter:
164: .  vec - the vector

166:    Options Database Keys:
167: +  -vec_view - Prints vector in ASCII format
168: .  -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
169: .  -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
170: .  -vec_view draw - Activates vector viewing using drawing tools
171: .  -display <name> - Sets display name (default is host)
172: .  -draw_pause <sec> - Sets number of seconds to pause after display
173: -  -vec_view socket - Activates vector viewing using a socket

175:    Level: beginner

177: .seealso: VecAssemblyBegin(), VecSetValues()
178: @*/
179: PetscErrorCode  VecAssemblyEnd(Vec vec)
180: {

185:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
187:   if (vec->ops->assemblyend) {
188:     (*vec->ops->assemblyend)(vec);
189:   }
190:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
191:   VecViewFromOptions(vec,NULL,"-vec_view");
192:   return(0);
193: }

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

200:    Logically Collective on Vec

202:    Input Parameters:
203: .  x, y  - the vectors

205:    Output Parameter:
206: .  w - the result

208:    Level: advanced

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

213:    Concepts: vector^pointwise multiply

215: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
216: @*/
217: PetscErrorCode  VecPointwiseMax(Vec w,Vec x,Vec y)
218: {

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

233:   (*w->ops->pointwisemax)(w,x,y);
234:   PetscObjectStateIncrease((PetscObject)w);
235:   return(0);
236: }


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

244:    Logically Collective on Vec

246:    Input Parameters:
247: .  x, y  - the vectors

249:    Output Parameter:
250: .  w - the result

252:    Level: advanced

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

257:    Concepts: vector^pointwise multiply

259: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
260: @*/
261: PetscErrorCode  VecPointwiseMin(Vec w,Vec x,Vec y)
262: {

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

277:   (*w->ops->pointwisemin)(w,x,y);
278:   PetscObjectStateIncrease((PetscObject)w);
279:   return(0);
280: }

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

287:    Logically Collective on Vec

289:    Input Parameters:
290: .  x, y  - the vectors

292:    Output Parameter:
293: .  w - the result

295:    Level: advanced

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

299:    Concepts: vector^pointwise multiply

301: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
302: @*/
303: PetscErrorCode  VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
304: {

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

319:   (*w->ops->pointwisemaxabs)(w,x,y);
320:   PetscObjectStateIncrease((PetscObject)w);
321:   return(0);
322: }

326: /*@
327:    VecPointwiseDivide - Computes the componentwise division w = x/y.

329:    Logically Collective on Vec

331:    Input Parameters:
332: .  x, y  - the vectors

334:    Output Parameter:
335: .  w - the result

337:    Level: advanced

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

341:    Concepts: vector^pointwise divide

343: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
344: @*/
345: PetscErrorCode  VecPointwiseDivide(Vec w,Vec x,Vec y)
346: {

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

361:   (*w->ops->pointwisedivide)(w,x,y);
362:   PetscObjectStateIncrease((PetscObject)w);
363:   return(0);
364: }


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

372:    Collective on Vec

374:    Input Parameters:
375: .  v - a vector to mimic

377:    Output Parameter:
378: .  newv - location to put new vector

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

384:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
385:    vectors.

387:    Level: beginner

389: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
390: @*/
391: PetscErrorCode  VecDuplicate(Vec v,Vec *newv)
392: {

399:   (*v->ops->duplicate)(v,newv);
400:   PetscObjectStateIncrease((PetscObject)*newv);
401:   return(0);
402: }

406: /*@
407:    VecDestroy - Destroys a vector.

409:    Collective on Vec

411:    Input Parameters:
412: .  v  - the vector

414:    Level: beginner

416: .seealso: VecDuplicate(), VecDestroyVecs()
417: @*/
418: PetscErrorCode  VecDestroy(Vec *v)
419: {

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

427:   PetscObjectSAWsViewOff((PetscObject)*v);
428:   /* destroy the internal part */
429:   if ((*v)->ops->destroy) {
430:     (*(*v)->ops->destroy)(*v);
431:   }
432:   /* destroy the external/common part */
433:   PetscLayoutDestroy(&(*v)->map);
434:   PetscHeaderDestroy(v);
435:   return(0);
436: }

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

443:    Collective on Vec

445:    Input Parameters:
446: +  m - the number of vectors to obtain
447: -  v - a vector to mimic

449:    Output Parameter:
450: .  V - location to put pointer to array of vectors

452:    Notes:
453:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
454:    vector.

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

461:    Level: intermediate

463: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
464: @*/
465: PetscErrorCode  VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
466: {

473:   (*v->ops->duplicatevecs)(v, m,V);
474:   return(0);
475: }

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

482:    Collective on Vec

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

488:    Fortran Note:
489:    The Fortran interface is slightly different from that given below.
490:    See the Fortran chapter of the users manual

492:    Level: intermediate

494: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
495: @*/
496: PetscErrorCode  VecDestroyVecs(PetscInt m,Vec *vv[])
497: {

502:   if (!*vv) return(0);
505:   if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
506:   if (!m) return(0);
507:   if (!*vv) return(0);
508:   (*(**vv)->ops->destroyvecs)(m,*vv);
509:   *vv  = 0;
510:   return(0);
511: }

515: /*@C
516:    VecView - Views a vector object.

518:    Collective on Vec

520:    Input Parameters:
521: +  vec - the vector
522: -  viewer - an optional visualization context

524:    Notes:
525:    The available visualization contexts include
526: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
527: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
528:          output where only the first processor opens
529:          the file.  All other processors send their
530:          data to the first processor to print.

532:    You can change the format the vector is printed using the
533:    option PetscViewerSetFormat().

535:    The user can open alternative visualization contexts with
536: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
537: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
538:          specified file; corresponding input uses VecLoad()
539: .    PetscViewerDrawOpen() - Outputs vector to an X window display
540: -    PetscViewerSocketOpen() - Outputs vector to Socket viewer

542:    The user can call PetscViewerSetFormat() to specify the output
543:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
544:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
545: +    PETSC_VIEWER_DEFAULT - default, prints vector contents
546: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
547: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
548: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
549:          format common among all vector types

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

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

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

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

570:    Level: beginner

572:    Concepts: vector^printing
573:    Concepts: vector^saving to disk

575: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
576:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
577:           PetscRealView(), PetscScalarView(), PetscIntView()
578: @*/
579: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
580: {
581:   PetscErrorCode    ierr;
582:   PetscBool         iascii;
583:   PetscViewerFormat format;

588:   if (!viewer) {
589:     PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
590:   }
593:   if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");

595:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
596:   PetscViewerGetFormat(viewer,&format);
597:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
598:   if (iascii) {
599:     PetscInt rows,bs;

601:     PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
602:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
603:       PetscViewerASCIIPushTab(viewer);
604:       VecGetSize(vec,&rows);
605:       VecGetBlockSize(vec,&bs);
606:       if (bs != 1) {
607:         PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
608:       } else {
609:         PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
610:       }
611:       PetscViewerASCIIPopTab(viewer);
612:     }
613:   }
614:   VecLockPush(vec);
615:   if (format == PETSC_VIEWER_NATIVE && vec->ops->viewnative) {
616:     (*vec->ops->viewnative)(vec,viewer);
617:   } else {
618:     (*vec->ops->view)(vec,viewer);
619:   }
620:   VecLockPop(vec);
621:   PetscLogEventEnd(VEC_View,vec,viewer,0,0);
622:   return(0);
623: }

625: #if defined(PETSC_USE_DEBUG)
626: #include <../src/sys/totalview/tv_data_display.h>
627: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
628: {
629:   const PetscScalar *values;
630:   char              type[32];
631:   PetscErrorCode    ierr;


634:   TV_add_row("Local rows", "int", &v->map->n);
635:   TV_add_row("Global rows", "int", &v->map->N);
636:   TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
637:   VecGetArrayRead((Vec)v,&values);
638:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
639:   TV_add_row("values",type, values);
640:   VecRestoreArrayRead((Vec)v,&values);
641:   return TV_format_OK;
642: }
643: #endif

647: /*@
648:    VecGetSize - Returns the global number of elements of the vector.

650:    Not Collective

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

655:    Output Parameters:
656: .  size - the global length of the vector

658:    Level: beginner

660:    Concepts: vector^local size

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

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

678: /*@
679:    VecGetLocalSize - Returns the number of elements of the vector stored
680:    in local memory. This routine may be implementation dependent, so use
681:    with care.

683:    Not Collective

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

688:    Output Parameter:
689: .  size - the length of the local piece of the vector

691:    Level: beginner

693:    Concepts: vector^size

695: .seealso: VecGetSize()
696: @*/
697: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
698: {

705:   (*x->ops->getlocalsize)(x,size);
706:   return(0);
707: }

711: /*@C
712:    VecGetOwnershipRange - Returns the range of indices owned by
713:    this processor, assuming that the vectors are laid out with the
714:    first n1 elements on the first processor, next n2 elements on the
715:    second, etc.  For certain parallel layouts this range may not be
716:    well defined.

718:    Not Collective

720:    Input Parameter:
721: .  x - the vector

723:    Output Parameters:
724: +  low - the first local element, pass in NULL if not interested
725: -  high - one more than the last local element, pass in NULL if not interested

727:    Note:
728:    The high argument is one more than the last element stored locally.

730:    Fortran: NULL_INTEGER should be used instead of NULL

732:    Level: beginner

734:    Concepts: ownership^of vectors
735:    Concepts: vector^ownership of elements

737: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
738: @*/
739: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
740: {
746:   if (low)  *low  = x->map->rstart;
747:   if (high) *high = x->map->rend;
748:   return(0);
749: }

753: /*@C
754:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
755:    assuming that the vectors are laid out with the
756:    first n1 elements on the first processor, next n2 elements on the
757:    second, etc.  For certain parallel layouts this range may not be
758:    well defined.

760:    Not Collective

762:    Input Parameter:
763: .  x - the vector

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

768:    Note:
769:    The high argument is one more than the last element stored locally.

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

773:    Level: beginner

775:    Concepts: ownership^of vectors
776:    Concepts: vector^ownership of elements

778: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
779: @*/
780: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
781: {

787:   PetscLayoutGetRanges(x->map,ranges);
788:   return(0);
789: }

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

796:    Collective on Vec

798:    Input Parameter:
799: +  x - the vector
800: .  op - the option
801: -  flag - turn the option on or off

803:    Supported Options:
804: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
805:           entries destined to be stored on a separate processor. This can be used
806:           to eliminate the global reduction in the VecAssemblyXXXX() if you know
807:           that you have only used VecSetValues() to set local elements
808: .     VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
809:           in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
810:           ignored.

812:    Level: intermediate

814: @*/
815: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool flag)
816: {

822:   if (x->ops->setoption) {
823:     (*x->ops->setoption)(x,op,flag);
824:   }
825:   return(0);
826: }

830: /* Default routines for obtaining and releasing; */
831: /* may be used by any implementation */
832: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
833: {
835:   PetscInt       i;

840:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
841:   PetscMalloc1(m,V);
842:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
843:   return(0);
844: }

848: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
849: {
851:   PetscInt       i;

855:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
856:   PetscFree(v);
857:   return(0);
858: }

862: /*@
863:    VecResetArray - Resets a vector to use its default memory. Call this
864:    after the use of VecPlaceArray().

866:    Not Collective

868:    Input Parameters:
869: .  vec - the vector

871:    Level: developer

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

875: @*/
876: PetscErrorCode  VecResetArray(Vec vec)
877: {

883:   if (vec->ops->resetarray) {
884:     (*vec->ops->resetarray)(vec);
885:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
886:   PetscObjectStateIncrease((PetscObject)vec);
887:   return(0);
888: }

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

896:   Collective on PetscViewer

898:   Input Parameters:
899: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
900:            some related function before a call to VecLoad().
901: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
902:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

904:    Level: intermediate

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

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

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

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

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

927:   Notes for advanced users:
928:   Most users should not need to know the details of the binary storage
929:   format, since VecLoad() and VecView() completely hide these details.
930:   But for anyone who's interested, the standard binary matrix storage
931:   format is
932: .vb
933:      int    VEC_FILE_CLASSID
934:      int    number of rows
935:      PetscScalar *values of all entries
936: .ve

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

944:   Concepts: vector^loading from file

946: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
947: @*/
948: PetscErrorCode  VecLoad(Vec newvec, PetscViewer viewer)
949: {
951:   PetscBool      isbinary,ishdf5;
952:   PetscViewerFormat format;

957:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
958:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
959:   if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");

961:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
962:   if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
963:     VecSetType(newvec, VECSTANDARD);
964:   }
965:   PetscViewerGetFormat(viewer,&format);
966:   if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
967:     (*newvec->ops->loadnative)(newvec,viewer);
968:   } else {
969:     (*newvec->ops->load)(newvec,viewer);
970:   }
971:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
972:   return(0);
973: }


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

981:    Logically Collective on Vec

983:    Input Parameter:
984: .  vec - the vector

986:    Output Parameter:
987: .  vec - the vector reciprocal

989:    Level: intermediate

991:    Concepts: vector^reciprocal

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

995: @*/
996: PetscErrorCode  VecReciprocal(Vec vec)
997: {

1003:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1004:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
1005:   (*vec->ops->reciprocal)(vec);
1006:   PetscObjectStateIncrease((PetscObject)vec);
1007:   return(0);
1008: }

1012: /*@C
1013:     VecSetOperation - Allows user to set a vector operation.

1015:    Logically Collective on Vec

1017:     Input Parameters:
1018: +   vec - the vector
1019: .   op - the name of the operation
1020: -   f - the function that provides the operation.

1022:    Level: advanced

1024:     Usage:
1025: $      PetscErrorCode userview(Vec,PetscViewer);
1026: $      VecCreateMPI(comm,m,M,&x);
1027: $      VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);

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

1035:     This function is not currently available from Fortran.

1037: .keywords: vector, set, operation

1039: .seealso: VecCreate(), MatShellSetOperation()
1040: @*/
1041: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1042: {
1045:   if (op == VECOP_VIEW && !vec->ops->viewnative) {
1046:     vec->ops->viewnative = vec->ops->view;
1047:   } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1048:     vec->ops->loadnative = vec->ops->load;
1049:   }
1050:   (((void(**)(void))vec->ops)[(int)op]) = f;
1051:   return(0);
1052: }


1057: /*@
1058:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1059:    used during the assembly process to store values that belong to
1060:    other processors.

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

1064:    Input Parameters:
1065: +  vec   - the vector
1066: .  size  - the initial size of the stash.
1067: -  bsize - the initial size of the block-stash(if used).

1069:    Options Database Keys:
1070: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1071: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1073:    Level: intermediate

1075:    Notes:
1076:      The block-stash is used for values set with VecSetValuesBlocked() while
1077:      the stash is used for values set with VecSetValues()

1079:      Run with the option -info and look for output of the form
1080:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1081:      to determine the appropriate value, MM, to use for size and
1082:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1083:      to determine the value, BMM to use for bsize

1085:    Concepts: vector^stash
1086:    Concepts: stash^vector

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

1090: @*/
1091: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1092: {

1097:   VecStashSetInitialSize_Private(&vec->stash,size);
1098:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1099:   return(0);
1100: }

1104: /*@
1105:    VecConjugate - Conjugates a vector.

1107:    Logically Collective on Vec

1109:    Input Parameters:
1110: .  x - the vector

1112:    Level: intermediate

1114:    Concepts: vector^conjugate

1116: @*/
1117: PetscErrorCode  VecConjugate(Vec x)
1118: {
1119: #if defined(PETSC_USE_COMPLEX)

1125:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1126:   (*x->ops->conjugate)(x);
1127:   /* we need to copy norms here */
1128:   PetscObjectStateIncrease((PetscObject)x);
1129:   return(0);
1130: #else
1131:   return(0);
1132: #endif
1133: }

1137: /*@
1138:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1140:    Logically Collective on Vec

1142:    Input Parameters:
1143: .  x, y  - the vectors

1145:    Output Parameter:
1146: .  w - the result

1148:    Level: advanced

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

1152:    Concepts: vector^pointwise multiply

1154: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1155: @*/
1156: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1157: {

1169:   if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1171:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1172:   (*w->ops->pointwisemult)(w,x,y);
1173:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1174:   PetscObjectStateIncrease((PetscObject)w);
1175:   return(0);
1176: }

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

1183:    Logically Collective on Vec

1185:    Input Parameters:
1186: +  x  - the vector
1187: -  rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1188:           it will create one internally.

1190:    Output Parameter:
1191: .  x  - the vector

1193:    Example of Usage:
1194: .vb
1195:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1196:      VecSetRandom(x,rctx);
1197:      PetscRandomDestroy(rctx);
1198: .ve

1200:    Level: intermediate

1202:    Concepts: vector^setting to random
1203:    Concepts: random^vector

1205: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1206: @*/
1207: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1208: {
1210:   PetscRandom    randObj = NULL;

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

1218:   if (!rctx) {
1219:     MPI_Comm comm;
1220:     PetscObjectGetComm((PetscObject)x,&comm);
1221:     PetscRandomCreate(comm,&randObj);
1222:     PetscRandomSetFromOptions(randObj);
1223:     rctx = randObj;
1224:   }

1226:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1227:   (*x->ops->setrandom)(x,rctx);
1228:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1230:   PetscRandomDestroy(&randObj);
1231:   PetscObjectStateIncrease((PetscObject)x);
1232:   return(0);
1233: }

1237: /*@
1238:   VecZeroEntries - puts a 0.0 in each element of a vector

1240:   Logically Collective on Vec

1242:   Input Parameter:
1243: . vec - The vector

1245:   Level: beginner

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

1252: .keywords: Vec, set, options, database
1253: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1254: @*/
1255: PetscErrorCode  VecZeroEntries(Vec vec)
1256: {

1260:   VecSet(vec,0);
1261:   return(0);
1262: }

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

1270:   Collective on Vec

1272:   Input Parameter:
1273: . vec - The vector

1275:   Level: intermediate

1277: .keywords: Vec, set, options, database, type
1278: .seealso: VecSetFromOptions(), VecSetType()
1279: */
1280: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptions *PetscOptionsObject,Vec vec)
1281: {
1282:   PetscBool      opt;
1283:   VecType        defaultType;
1284:   char           typeName[256];
1285:   PetscMPIInt    size;

1289:   if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1290:   else {
1291:     MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1292:     if (size > 1) defaultType = VECMPI;
1293:     else defaultType = VECSEQ;
1294:   }

1296:   VecRegisterAll();
1297:   PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1298:   if (opt) {
1299:     VecSetType(vec, typeName);
1300:   } else {
1301:     VecSetType(vec, defaultType);
1302:   }
1303:   return(0);
1304: }

1308: /*@
1309:   VecSetFromOptions - Configures the vector from the options database.

1311:   Collective on Vec

1313:   Input Parameter:
1314: . vec - The vector

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

1319:   Level: beginner

1321:   Concepts: vectors^setting options
1322:   Concepts: vectors^setting type

1324: .keywords: Vec, set, options, database
1325: .seealso: VecCreate(), VecSetOptionsPrefix()
1326: @*/
1327: PetscErrorCode  VecSetFromOptions(Vec vec)
1328: {


1334:   PetscObjectOptionsBegin((PetscObject)vec);
1335:   /* Handle vector type options */
1336:   VecSetTypeFromOptions_Private(PetscOptionsObject,vec);

1338:   /* Handle specific vector options */
1339:   if (vec->ops->setfromoptions) {
1340:     (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1341:   }

1343:   /* process any options handlers added with PetscObjectAddOptionsHandler() */
1344:   PetscObjectProcessOptionsHandlers((PetscObject)vec);
1345:   PetscOptionsEnd();
1346:   return(0);
1347: }

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

1354:   Collective on Vec

1356:   Input Parameters:
1357: + v - the vector
1358: . n - the local size (or PETSC_DECIDE to have it set)
1359: - N - the global size (or PETSC_DECIDE)

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

1365:   Level: intermediate

1367: .seealso: VecGetSize(), PetscSplitOwnership()
1368: @*/
1369: PetscErrorCode  VecSetSizes(Vec v, PetscInt n, PetscInt N)
1370: {

1376:   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);
1377:   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);
1378:   v->map->n = n;
1379:   v->map->N = N;
1380:   if (v->ops->create) {
1381:     (*v->ops->create)(v);
1382:     v->ops->create = 0;
1383:   }
1384:   return(0);
1385: }

1389: /*@
1390:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1391:    and VecSetValuesBlockedLocal().

1393:    Logically Collective on Vec

1395:    Input Parameter:
1396: +  v - the vector
1397: -  bs - the blocksize

1399:    Notes:
1400:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1402:    Level: advanced

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

1406:   Concepts: block size^vectors
1407: @*/
1408: PetscErrorCode  VecSetBlockSize(Vec v,PetscInt bs)
1409: {

1414:   if (bs < 0 || bs == v->map->bs) return(0);
1416:   PetscLayoutSetBlockSize(v->map,bs);
1417:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1418:   return(0);
1419: }

1423: /*@
1424:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1425:    and VecSetValuesBlockedLocal().

1427:    Not Collective

1429:    Input Parameter:
1430: .  v - the vector

1432:    Output Parameter:
1433: .  bs - the blocksize

1435:    Notes:
1436:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1438:    Level: advanced

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

1442:    Concepts: vector^block size
1443:    Concepts: block^vector

1445: @*/
1446: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1447: {

1453:   PetscLayoutGetBlockSize(v->map,bs);
1454:   return(0);
1455: }

1459: /*@C
1460:    VecSetOptionsPrefix - Sets the prefix used for searching for all
1461:    Vec options in the database.

1463:    Logically Collective on Vec

1465:    Input Parameter:
1466: +  v - the Vec context
1467: -  prefix - the prefix to prepend to all option names

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

1473:    Level: advanced

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

1477: .seealso: VecSetFromOptions()
1478: @*/
1479: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1480: {

1485:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1486:   return(0);
1487: }

1491: /*@C
1492:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1493:    Vec options in the database.

1495:    Logically Collective on Vec

1497:    Input Parameters:
1498: +  v - the Vec context
1499: -  prefix - the prefix to prepend to all option names

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

1505:    Level: advanced

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

1509: .seealso: VecGetOptionsPrefix()
1510: @*/
1511: PetscErrorCode  VecAppendOptionsPrefix(Vec v,const char prefix[])
1512: {

1517:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1518:   return(0);
1519: }

1523: /*@C
1524:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1525:    Vec options in the database.

1527:    Not Collective

1529:    Input Parameter:
1530: .  v - the Vec context

1532:    Output Parameter:
1533: .  prefix - pointer to the prefix string used

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

1538:    Level: advanced

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

1542: .seealso: VecAppendOptionsPrefix()
1543: @*/
1544: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1545: {

1550:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1551:   return(0);
1552: }

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

1559:    Collective on Vec

1561:    Input Parameters:
1562: .  v - the Vec context

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

1568:    Level: advanced

1570: .keywords: Vec, setup

1572: .seealso: VecCreate(), VecDestroy()
1573: @*/
1574: PetscErrorCode  VecSetUp(Vec v)
1575: {
1576:   PetscMPIInt    size;

1581:   if (!((PetscObject)v)->type_name) {
1582:     MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1583:     if (size == 1) {
1584:       VecSetType(v, VECSEQ);
1585:     } else {
1586:       VecSetType(v, VECMPI);
1587:     }
1588:   }
1589:   return(0);
1590: }

1592: /*
1593:     These currently expose the PetscScalar/PetscReal in updating the
1594:     cached norm. If we push those down into the implementation these
1595:     will become independent of PetscScalar/PetscReal
1596: */

1600: /*@
1601:    VecCopy - Copies a vector. y <- x

1603:    Logically Collective on Vec

1605:    Input Parameter:
1606: .  x - the vector

1608:    Output Parameter:
1609: .  y - the copy

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

1615:    Level: beginner

1617: .seealso: VecDuplicate()
1618: @*/
1619: PetscErrorCode  VecCopy(Vec x,Vec y)
1620: {
1621:   PetscBool      flgs[4];
1622:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1624:   PetscInt       i;

1631:   if (x == y) return(0);
1632:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1633:   if (x->map->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths %d != %d", x->map->n, y->map->n);
1634:   VecLocked(y,2);

1636: #if !defined(PETSC_USE_MIXED_PRECISION)
1637:   for (i=0; i<4; i++) {
1638:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1639:   }
1640: #endif

1642:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
1643: #if defined(PETSC_USE_MIXED_PRECISION)
1644:   extern PetscErrorCode VecGetArray(Vec,double**);
1645:   extern PetscErrorCode VecRestoreArray(Vec,double**);
1646:   extern PetscErrorCode VecGetArray(Vec,float**);
1647:   extern PetscErrorCode VecRestoreArray(Vec,float**);
1648:   extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1649:   extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1650:   extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1651:   extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1652:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1653:     PetscInt    i,n;
1654:     const float *xx;
1655:     double      *yy;
1656:     VecGetArrayRead(x,&xx);
1657:     VecGetArray(y,&yy);
1658:     VecGetLocalSize(x,&n);
1659:     for (i=0; i<n; i++) yy[i] = xx[i];
1660:     VecRestoreArrayRead(x,&xx);
1661:     VecRestoreArray(y,&yy);
1662:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1663:     PetscInt     i,n;
1664:     float        *yy;
1665:     const double *xx;
1666:     VecGetArrayRead(x,&xx);
1667:     VecGetArray(y,&yy);
1668:     VecGetLocalSize(x,&n);
1669:     for (i=0; i<n; i++) yy[i] = (float) xx[i];
1670:     VecRestoreArrayRead(x,&xx);
1671:     VecRestoreArray(y,&yy);
1672:   } else {
1673:     (*x->ops->copy)(x,y);
1674:   }
1675: #else
1676:   (*x->ops->copy)(x,y);
1677: #endif

1679:   PetscObjectStateIncrease((PetscObject)y);
1680: #if !defined(PETSC_USE_MIXED_PRECISION)
1681:   for (i=0; i<4; i++) {
1682:     if (flgs[i]) {
1683:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1684:     }
1685:   }
1686: #endif

1688:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1689:   return(0);
1690: }

1694: /*@
1695:    VecSwap - Swaps the vectors x and y.

1697:    Logically Collective on Vec

1699:    Input Parameters:
1700: .  x, y  - the vectors

1702:    Level: advanced

1704:    Concepts: vector^swapping values

1706: @*/
1707: PetscErrorCode  VecSwap(Vec x,Vec y)
1708: {
1709:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1710:   PetscBool      flgxs[4],flgys[4];
1712:   PetscInt       i;

1720:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1721:   if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1722:   if (x->map->N != y->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1723:   if (x->map->n != y->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");

1725:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1726:   for (i=0; i<4; i++) {
1727:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1728:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1729:   }
1730:   (*x->ops->swap)(x,y);
1731:   PetscObjectStateIncrease((PetscObject)x);
1732:   PetscObjectStateIncrease((PetscObject)y);
1733:   for (i=0; i<4; i++) {
1734:     if (flgxs[i]) {
1735:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1736:     }
1737:     if (flgys[i]) {
1738:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1739:     }
1740:   }
1741:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1742:   return(0);
1743: }

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

1750:   Collective on VecStash

1752:   Input Parameters:
1753: + obj   - the VecStash object
1754: . bobj - optional other object that provides the prefix
1755: - optionname - option to activate viewing

1757:   Level: intermediate

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

1761: */
1762: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1763: {
1764:   PetscErrorCode    ierr;
1765:   PetscViewer       viewer;
1766:   PetscBool         flg;
1767:   PetscViewerFormat format;
1768:   char              *prefix;

1771:   prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1772:   PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),prefix,optionname,&viewer,&format,&flg);
1773:   if (flg) {
1774:     PetscViewerPushFormat(viewer,format);
1775:     VecStashView(obj,viewer);
1776:     PetscViewerPopFormat(viewer);
1777:     PetscViewerDestroy(&viewer);
1778:   }
1779:   return(0);
1780: }

1784: /*@
1785:    VecStashView - Prints the entries in the vector stash and block stash.

1787:    Collective on Vec

1789:    Input Parameters:
1790: +  v - the vector
1791: -  viewer - the viewer

1793:    Level: advanced

1795:    Concepts: vector^stash
1796:    Concepts: stash^vector

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

1800: @*/
1801: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1802: {
1804:   PetscMPIInt    rank;
1805:   PetscInt       i,j;
1806:   PetscBool      match;
1807:   VecStash       *s;
1808:   PetscScalar    val;


1815:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1816:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1817:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1818:   MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1819:   s    = &v->bstash;

1821:   /* print block stash */
1822:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1823:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1824:   for (i=0; i<s->n; i++) {
1825:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1826:     for (j=0; j<s->bs; j++) {
1827:       val = s->array[i*s->bs+j];
1828: #if defined(PETSC_USE_COMPLEX)
1829:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1830: #else
1831:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1832: #endif
1833:     }
1834:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1835:   }
1836:   PetscViewerFlush(viewer);

1838:   s = &v->stash;

1840:   /* print basic stash */
1841:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1842:   for (i=0; i<s->n; i++) {
1843:     val = s->array[i];
1844: #if defined(PETSC_USE_COMPLEX)
1845:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1846: #else
1847:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1848: #endif
1849:   }
1850:   PetscViewerFlush(viewer);
1851:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);

1853:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1854:   return(0);
1855: }

1859: PetscErrorCode PetscOptionsGetVec(const char prefix[],const char key[],Vec v,PetscBool *set)
1860: {
1861:   PetscInt       i,N,rstart,rend;
1863:   PetscScalar    *xx;
1864:   PetscReal      *xreal;
1865:   PetscBool      iset;

1868:   VecGetOwnershipRange(v,&rstart,&rend);
1869:   VecGetSize(v,&N);
1870:   PetscCalloc1(N,&xreal);
1871:   PetscOptionsGetRealArray(prefix,key,xreal,&N,&iset);
1872:   if (iset) {
1873:     VecGetArray(v,&xx);
1874:     for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1875:     VecRestoreArray(v,&xx);
1876:   }
1877:   PetscFree(xreal);
1878:   if (set) *set = iset;
1879:   return(0);
1880: }

1884: /*@
1885:    VecGetLayout - get PetscLayout describing vector layout

1887:    Not Collective

1889:    Input Arguments:
1890: .  x - the vector

1892:    Output Arguments:
1893: .  map - the layout

1895:    Level: developer

1897: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1898: @*/
1899: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1900: {

1904:   *map = x->map;
1905:   return(0);
1906: }

1910: /*@
1911:    VecSetLayout - set PetscLayout describing vector layout

1913:    Not Collective

1915:    Input Arguments:
1916: +  x - the vector
1917: -  map - the layout

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

1922:    Level: developer

1924: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1925: @*/
1926: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1927: {

1932:   PetscLayoutReference(map,&x->map);
1933:   return(0);
1934: }

1938: PetscErrorCode VecSetInf(Vec xin)
1939: {
1940:   PetscInt       i,n = xin->map->n;
1941:   PetscScalar    *xx;
1942:   PetscScalar    zero=0.0,one=1.0,inf=one/zero;

1946:   VecGetArray(xin,&xx);
1947:   for (i=0; i<n; i++) xx[i] = inf;
1948:   VecRestoreArray(xin,&xx);
1949:   return(0);
1950: }