Actual source code: vector.c

petsc-3.3-p7 2013-05-11
  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;

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

 25:     Not collective

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

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

 36:    Level: advanced

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

 40: @*/
 41: PetscErrorCode  VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
 42: {
 45:   VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
 46:   VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
 47:   return(0);
 48: }

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

 57:    Logically Collective on Vec

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

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

 66:    Level: intermediate

 68:    Concepts: vector^setting values with local numbering

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


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

 91: /*@
 92:    VecSetLocalToGlobalMappingBlock - Sets a local numbering to global numbering used
 93:    by the routine VecSetValuesBlockedLocal() to allow users to insert vector entries
 94:    using a local (per-processor) numbering.

 96:    Logically Collective on Vec

 98:    Input Parameters:
 99: +  x - vector
100: -  mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()

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

105:    Level: intermediate

107:    Concepts: vector^setting values blocked with local numbering

109: .seealso:  VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
110:            VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
111: @*/
112: PetscErrorCode  VecSetLocalToGlobalMappingBlock(Vec x,ISLocalToGlobalMapping mapping)
113: {


120:   PetscLayoutSetISLocalToGlobalMappingBlock(x->map,mapping);
121:   return(0);
122: }

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

129:    Not Collective

131:    Input Parameter:
132: .  X - the vector

134:    Output Parameter:
135: .  mapping - the mapping

137:    Level: advanced

139:    Concepts: vectors^local to global mapping
140:    Concepts: local to global mapping^for vectors

142: .seealso:  VecSetValuesLocal(), VecGetLocalToGlobalMappingBlock()
143: @*/
144: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
145: {
150:   *mapping = X->map->mapping;
151:   return(0);
152: }

156: /*@
157:    VecGetLocalToGlobalMappingBlock - Gets the local-to-global numbering set by VecSetLocalToGlobalMappingBlock()

159:    Not Collective

161:    Input Parameters:
162: .  X - the vector

164:    Output Parameters:
165: .  mapping - the mapping

167:    Level: advanced

169:    Concepts: vectors^local to global mapping blocked
170:    Concepts: local to global mapping^for vectors, blocked

172: .seealso:  VecSetValuesBlockedLocal(), VecGetLocalToGlobalMapping()
173: @*/
174: PetscErrorCode VecGetLocalToGlobalMappingBlock(Vec X,ISLocalToGlobalMapping *mapping)
175: {
180:   *mapping = X->map->bmapping;
181:   return(0);
182: }

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

190:    Collective on Vec

192:    Input Parameter:
193: .  vec - the vector

195:    Level: beginner

197:    Concepts: assembly^vectors

199: .seealso: VecAssemblyEnd(), VecSetValues()
200: @*/
201: PetscErrorCode  VecAssemblyBegin(Vec vec)
202: {
204:   PetscBool      flg = PETSC_FALSE;


210:   PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_stash",&flg,PETSC_NULL);
211:   if (flg) {
212:     PetscViewer viewer;
213:     PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
214:     VecStashView(vec,viewer);
215:   }

217:   PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
218:   if (vec->ops->assemblybegin) {
219:     (*vec->ops->assemblybegin)(vec);
220:   }
221:   PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
222:   PetscObjectStateIncrease((PetscObject)vec);
223:   return(0);
224: }

228: /*
229:   Processes command line options to determine if/how a matrix
230:   is to be viewed. Called by VecAssemblyEnd().

232: .seealso: MatView_Private()

234: */
235: PetscErrorCode  VecView_Private(Vec vec)
236: {
238:   PetscBool      flg = PETSC_FALSE;

241:   PetscObjectOptionsBegin((PetscObject)vec);
242:     PetscOptionsBool("-vec_view_info","Information on vector size","VecView",flg,&flg,PETSC_NULL);
243:     if (flg) {
244:       PetscViewer viewer;

246:       PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
247:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO);
248:       VecView(vec,viewer);
249:       PetscViewerPopFormat(viewer);
250:     }
251:     flg  = PETSC_FALSE;
252:     PetscOptionsBool("-vec_view","Print vector to stdout","VecView",flg,&flg,PETSC_NULL);
253:     if (flg) {
254:       PetscViewer viewer;
255:       PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
256:       VecView(vec,viewer);
257:     }
258:     flg  = PETSC_FALSE;
259:     PetscOptionsBool("-vec_view_matlab","Print vector to stdout in a format MATLAB can read","VecView",flg,&flg,PETSC_NULL);
260:     if (flg) {
261:       PetscViewer viewer;
262:       PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
263:       PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
264:       VecView(vec,viewer);
265:       PetscViewerPopFormat(viewer);
266:     }
267: #if defined(PETSC_HAVE_MATLAB_ENGINE)
268:     flg  = PETSC_FALSE;
269:     PetscOptionsBool("-vec_view_matlab_file","Print vector to matlaboutput.mat format MATLAB can read","VecView",flg,&flg,PETSC_NULL);
270:     if (flg) {
271:       VecView(vec,PETSC_VIEWER_MATLAB_(((PetscObject)vec)->comm));
272:     }
273: #endif
274: #if defined(PETSC_USE_SOCKET_VIEWER)
275:     flg  = PETSC_FALSE;
276:     PetscOptionsBool("-vec_view_socket","Send vector to socket (can be read from matlab)","VecView",flg,&flg,PETSC_NULL);
277:     if (flg) {
278:       VecView(vec,PETSC_VIEWER_SOCKET_(((PetscObject)vec)->comm));
279:       PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)vec)->comm));
280:     }
281: #endif
282:     flg  = PETSC_FALSE;
283:     PetscOptionsBool("-vec_view_binary","Save vector to file in binary format","VecView",flg,&flg,PETSC_NULL);
284:     if (flg) {
285:       VecView(vec,PETSC_VIEWER_BINARY_(((PetscObject)vec)->comm));
286:       PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)vec)->comm));
287:     }
288:   PetscOptionsEnd();
289:   /* These invoke PetscDrawGetDraw which invokes PetscOptionsBegin/End, */
290:   /* hence they should not be inside the above PetscOptionsBegin/End block. */
291:   flg  = PETSC_FALSE;
292:   PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_draw",&flg,PETSC_NULL);
293:   if (flg) {
294:     VecView(vec,PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
295:     PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
296:   }
297:   flg  = PETSC_FALSE;
298:   PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_draw_lg",&flg,PETSC_NULL);
299:   if (flg) {
300:     PetscViewerSetFormat(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm),PETSC_VIEWER_DRAW_LG);
301:     VecView(vec,PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
302:     PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
303:   }
304:   return(0);
305: }

309: /*@
310:    VecAssemblyEnd - Completes assembling the vector.  This routine should
311:    be called after VecAssemblyBegin().

313:    Collective on Vec

315:    Input Parameter:
316: .  vec - the vector

318:    Options Database Keys:
319: +  -vec_view - Prints vector in ASCII format
320: .  -vec_view_matlab - Prints vector in ASCII MATLAB format to stdout
321: .  -vec_view_matlab_file - Prints vector in MATLAB format to matlaboutput.mat
322: .  -vec_view_draw - Activates vector viewing using drawing tools
323: .  -display <name> - Sets display name (default is host)
324: .  -draw_pause <sec> - Sets number of seconds to pause after display
325: -  -vec_view_socket - Activates vector viewing using a socket

327:    Level: beginner

329: .seealso: VecAssemblyBegin(), VecSetValues()
330: @*/
331: PetscErrorCode  VecAssemblyEnd(Vec vec)
332: {

337:   PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
339:   if (vec->ops->assemblyend) {
340:     (*vec->ops->assemblyend)(vec);
341:   }
342:   PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
343:   VecView_Private(vec);
344:   return(0);
345: }

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

352:    Logically Collective on Vec

354:    Input Parameters:
355: .  x, y  - the vectors

357:    Output Parameter:
358: .  w - the result

360:    Level: advanced

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

365:    Concepts: vector^pointwise multiply

367: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
368: @*/
369: PetscErrorCode  VecPointwiseMax(Vec w,Vec x,Vec y)
370: {

382:   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");
383:   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");

385:   (*w->ops->pointwisemax)(w,x,y);
386:   PetscObjectStateIncrease((PetscObject)w);
387:   return(0);
388: }


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

396:    Logically Collective on Vec

398:    Input Parameters:
399: .  x, y  - the vectors

401:    Output Parameter:
402: .  w - the result

404:    Level: advanced

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

409:    Concepts: vector^pointwise multiply

411: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
412: @*/
413: PetscErrorCode  VecPointwiseMin(Vec w,Vec x,Vec y)
414: {

426:   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");
427:   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");

429:   (*w->ops->pointwisemin)(w,x,y);
430:   PetscObjectStateIncrease((PetscObject)w);
431:   return(0);
432: }

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

439:    Logically Collective on Vec

441:    Input Parameters:
442: .  x, y  - the vectors

444:    Output Parameter:
445: .  w - the result

447:    Level: advanced

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

451:    Concepts: vector^pointwise multiply

453: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
454: @*/
455: PetscErrorCode  VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
456: {

468:   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");
469:   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");

471:   (*w->ops->pointwisemaxabs)(w,x,y);
472:   PetscObjectStateIncrease((PetscObject)w);
473:   return(0);
474: }

478: /*@
479:    VecPointwiseDivide - Computes the componentwise division w = x/y.

481:    Logically Collective on Vec

483:    Input Parameters:
484: .  x, y  - the vectors

486:    Output Parameter:
487: .  w - the result

489:    Level: advanced

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

493:    Concepts: vector^pointwise divide

495: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
496: @*/
497: PetscErrorCode  VecPointwiseDivide(Vec w,Vec x,Vec y)
498: {

510:   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");
511:   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");

513:   (*w->ops->pointwisedivide)(w,x,y);
514:   PetscObjectStateIncrease((PetscObject)w);
515:   return(0);
516: }


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

524:    Collective on Vec

526:    Input Parameters:
527: .  v - a vector to mimic

529:    Output Parameter:
530: .  newv - location to put new vector

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

536:    Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
537:    vectors.

539:    Level: beginner

541: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
542: @*/
543: PetscErrorCode  VecDuplicate(Vec v,Vec *newv)
544: {

551:   (*v->ops->duplicate)(v,newv);
552:   PetscObjectStateIncrease((PetscObject)*newv);
553:   return(0);
554: }

558: /*@
559:    VecDestroy - Destroys a vector.

561:    Collective on Vec

563:    Input Parameters:
564: .  v  - the vector

566:    Level: beginner

568: .seealso: VecDuplicate(), VecDestroyVecs()
569: @*/
570: PetscErrorCode  VecDestroy(Vec *v)
571: {

575:   if (!*v) return(0);
577:   if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
578:   /* destroy the internal part */
579:   if ((*v)->ops->destroy) {
580:     (*(*v)->ops->destroy)(*v);
581:   }
582:   /* destroy the external/common part */
583:   PetscLayoutDestroy(&(*v)->map);
584:   PetscHeaderDestroy(v);
585:   return(0);
586: }

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

593:    Collective on Vec

595:    Input Parameters:
596: +  m - the number of vectors to obtain
597: -  v - a vector to mimic

599:    Output Parameter:
600: .  V - location to put pointer to array of vectors

602:    Notes:
603:    Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
604:    vector.

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

611:    Level: intermediate

613: .seealso:  VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
614: @*/
615: PetscErrorCode  VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
616: {

623:   (*v->ops->duplicatevecs)(v, m,V);
624:   return(0);
625: }

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

632:    Collective on Vec

634:    Input Parameters:
635: +  vv - pointer to pointer to array of vector pointers
636: -  m - the number of vectors previously obtained

638:    Fortran Note:
639:    The Fortran interface is slightly different from that given below.
640:    See the Fortran chapter of the users manual

642:    Level: intermediate

644: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
645: @*/
646: PetscErrorCode  VecDestroyVecs(PetscInt m,Vec *vv[])
647: {

652:   if (!*vv) return(0);
655:   if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
656:   (*(**vv)->ops->destroyvecs)(m,*vv);
657:   *vv = 0;
658:   return(0);
659: }

661: #undef  __FUNCT__
663: /*@
664:   VecViewFromOptions - This function visualizes the vector based upon user options.

666:   Collective on Vec

668:   Input Parameters:
669: . vec   - The vector
670: . title - The title (currently ignored)

672:   Level: intermediate

674: .keywords: Vec, view, options, database
675: .seealso: VecSetFromOptions(), VecView()
676: @*/
677: PetscErrorCode  VecViewFromOptions(Vec vec, const char *title)
678: {

682:   VecView_Private(vec);
683:   return(0);
684: }

688: /*@C
689:    VecView - Views a vector object.

691:    Collective on Vec

693:    Input Parameters:
694: +  vec - the vector
695: -  viewer - an optional visualization context

697:    Notes:
698:    The available visualization contexts include
699: +     PETSC_VIEWER_STDOUT_SELF - standard output (default)
700: -     PETSC_VIEWER_STDOUT_WORLD - synchronized standard
701:          output where only the first processor opens
702:          the file.  All other processors send their
703:          data to the first processor to print.

705:    You can change the format the vector is printed using the
706:    option PetscViewerSetFormat().

708:    The user can open alternative visualization contexts with
709: +    PetscViewerASCIIOpen() - Outputs vector to a specified file
710: .    PetscViewerBinaryOpen() - Outputs vector in binary to a
711:          specified file; corresponding input uses VecLoad()
712: .    PetscViewerDrawOpen() - Outputs vector to an X window display
713: -    PetscViewerSocketOpen() - Outputs vector to Socket viewer

715:    The user can call PetscViewerSetFormat() to specify the output
716:    format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
717:    PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen).  Available formats include
718: +    PETSC_VIEWER_DEFAULT - default, prints vector contents
719: .    PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
720: .    PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
721: -    PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
722:          format common among all vector types

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

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

732:    Level: beginner

734:    Concepts: vector^printing
735:    Concepts: vector^saving to disk

737: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
738:           PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
739:           PetscRealView(), PetscScalarView(), PetscIntView()
740: @*/
741: PetscErrorCode  VecView(Vec vec,PetscViewer viewer)
742: {
743:   PetscErrorCode    ierr;
744:   PetscBool         iascii;
745:   PetscViewerFormat format;

750:   if (!viewer) {
751:     PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
752:   }
755:   if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");

757:   PetscLogEventBegin(VEC_View,vec,viewer,0,0);
758:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
759:   if (iascii) {
760:     PetscInt rows,bs;

762:     PetscViewerGetFormat(viewer,&format);
763:     if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
764:       PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer,"Vector Object");
765:       PetscViewerASCIIPushTab(viewer);
766:       VecGetSize(vec,&rows);
767:       VecGetBlockSize(vec,&bs);
768:       if (bs != 1) {
769:         PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
770:       } else {
771:         PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
772:       }
773:       PetscViewerASCIIPopTab(viewer);
774:     }
775:   }
776:   (*vec->ops->view)(vec,viewer);
777:   PetscLogEventEnd(VEC_View,vec,viewer,0,0);
778:   return(0);
779: }

781: #if defined(PETSC_USE_DEBUG)
782: #include <../src/sys/totalview/tv_data_display.h>
783: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
784: {
785:   const PetscScalar *values;
786:   char              type[32];
787:   PetscErrorCode    ierr;


790:   TV_add_row("Local rows", "int", &v->map->n);
791:   TV_add_row("Global rows", "int", &v->map->N);
792:   TV_add_row("Typename", TV_ascii_string_type , ((PetscObject)v)->type_name);
793:   VecGetArrayRead((Vec)v,&values);
794:   PetscSNPrintf(type,32,"double[%d]",v->map->n);
795:   TV_add_row("values",type, values);
796:   VecRestoreArrayRead((Vec)v,&values);
797:   return TV_format_OK;
798: }
799: #endif

803: /*@
804:    VecGetSize - Returns the global number of elements of the vector.

806:    Not Collective

808:    Input Parameter:
809: .  x - the vector

811:    Output Parameters:
812: .  size - the global length of the vector

814:    Level: beginner

816:    Concepts: vector^local size

818: .seealso: VecGetLocalSize()
819: @*/
820: PetscErrorCode  VecGetSize(Vec x,PetscInt *size)
821: {

828:   (*x->ops->getsize)(x,size);
829:   return(0);
830: }

834: /*@
835:    VecGetLocalSize - Returns the number of elements of the vector stored
836:    in local memory. This routine may be implementation dependent, so use
837:    with care.

839:    Not Collective

841:    Input Parameter:
842: .  x - the vector

844:    Output Parameter:
845: .  size - the length of the local piece of the vector

847:    Level: beginner

849:    Concepts: vector^size

851: .seealso: VecGetSize()
852: @*/
853: PetscErrorCode  VecGetLocalSize(Vec x,PetscInt *size)
854: {

861:   (*x->ops->getlocalsize)(x,size);
862:   return(0);
863: }

867: /*@C
868:    VecGetOwnershipRange - Returns the range of indices owned by
869:    this processor, assuming that the vectors are laid out with the
870:    first n1 elements on the first processor, next n2 elements on the
871:    second, etc.  For certain parallel layouts this range may not be
872:    well defined.

874:    Not Collective

876:    Input Parameter:
877: .  x - the vector

879:    Output Parameters:
880: +  low - the first local element, pass in PETSC_NULL if not interested
881: -  high - one more than the last local element, pass in PETSC_NULL if not interested

883:    Note:
884:    The high argument is one more than the last element stored locally.

886:    Fortran: PETSC_NULL_INTEGER should be used instead of PETSC_NULL

888:    Level: beginner

890:    Concepts: ownership^of vectors
891:    Concepts: vector^ownership of elements

893: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
894: @*/
895: PetscErrorCode  VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
896: {
902:   if (low)  *low  = x->map->rstart;
903:   if (high) *high = x->map->rend;
904:   return(0);
905: }

909: /*@C
910:    VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
911:    assuming that the vectors are laid out with the
912:    first n1 elements on the first processor, next n2 elements on the
913:    second, etc.  For certain parallel layouts this range may not be
914:    well defined.

916:    Not Collective

918:    Input Parameter:
919: .  x - the vector

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

924:    Note:
925:    The high argument is one more than the last element stored locally.

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

929:    Level: beginner

931:    Concepts: ownership^of vectors
932:    Concepts: vector^ownership of elements

934: .seealso:   MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
935: @*/
936: PetscErrorCode  VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
937: {

943:   PetscLayoutGetRanges(x->map,ranges);
944:   return(0);
945: }

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

952:    Collective on Vec

954:    Input Parameter:
955: +  x - the vector
956: .  op - the option
957: -  flag - turn the option on or off

959:    Supported Options:
960: +     VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
961:           entries destined to be stored on a separate processor. This can be used
962:           to eliminate the global reduction in the VecAssemblyXXXX() if you know
963:           that you have only used VecSetValues() to set local elements
964: .     VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
965:           in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
966:           ignored.

968:    Level: intermediate

970: @*/
971: PetscErrorCode  VecSetOption(Vec x,VecOption op,PetscBool  flag)
972: {

978:   if (x->ops->setoption) {
979:     (*x->ops->setoption)(x,op,flag);
980:   }
981:   return(0);
982: }

986: /* Default routines for obtaining and releasing; */
987: /* may be used by any implementation */
988: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
989: {
991:   PetscInt       i;

996:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
997:   PetscMalloc(m*sizeof(Vec*),V);
998:   for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
999:   return(0);
1000: }

1004: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
1005: {
1007:   PetscInt       i;

1011:   if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
1012:   for (i=0; i<m; i++) {VecDestroy(&v[i]);}
1013:   PetscFree(v);
1014:   return(0);
1015: }

1019: /*@
1020:    VecResetArray - Resets a vector to use its default memory. Call this
1021:    after the use of VecPlaceArray().

1023:    Not Collective

1025:    Input Parameters:
1026: .  vec - the vector

1028:    Level: developer

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

1032: @*/
1033: PetscErrorCode  VecResetArray(Vec vec)
1034: {

1040:   if (vec->ops->resetarray) {
1041:     (*vec->ops->resetarray)(vec);
1042:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
1043:   PetscObjectStateIncrease((PetscObject)vec);
1044:   return(0);
1045: }

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

1053:   Collective on PetscViewer 

1055:   Input Parameters:
1056: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
1057:            some related function before a call to VecLoad(). 
1058: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
1059:            HDF5 file viewer, obtained from PetscViewerHDF5Open()

1061:    Level: intermediate

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

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

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

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

1078:   Notes for advanced users:
1079:   Most users should not need to know the details of the binary storage
1080:   format, since VecLoad() and VecView() completely hide these details.
1081:   But for anyone who's interested, the standard binary matrix storage
1082:   format is
1083: .vb
1084:      int    VEC_FILE_CLASSID
1085:      int    number of rows
1086:      PetscScalar *values of all entries
1087: .ve

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

1095:   Concepts: vector^loading from file

1097: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad() 
1098: @*/
1099: PetscErrorCode  VecLoad(Vec newvec, PetscViewer viewer)
1100: {


1107:   PetscLogEventBegin(VEC_Load,viewer,0,0,0);
1108:   if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
1109:     VecSetType(newvec, VECSTANDARD);
1110:   }
1111:   (*newvec->ops->load)(newvec,viewer);
1112:   PetscLogEventEnd(VEC_Load,viewer,0,0,0);
1113:   return(0);
1114: }


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

1122:    Logically Collective on Vec

1124:    Input Parameter:
1125: .  vec - the vector

1127:    Output Parameter:
1128: .  vec - the vector reciprocal

1130:    Level: intermediate

1132:    Concepts: vector^reciprocal

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

1136: @*/
1137: PetscErrorCode  VecReciprocal(Vec vec)
1138: {

1144:   if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1145:   if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
1146:   (*vec->ops->reciprocal)(vec);
1147:   PetscObjectStateIncrease((PetscObject)vec);
1148:   return(0);
1149: }

1153: PetscErrorCode  VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1154: {
1157:   (((void(**)(void))vec->ops)[(int)op]) = f;
1158:   return(0);
1159: }


1164: /*@
1165:    VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1166:    used during the assembly process to store values that belong to
1167:    other processors.

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

1171:    Input Parameters:
1172: +  vec   - the vector
1173: .  size  - the initial size of the stash.
1174: -  bsize - the initial size of the block-stash(if used).

1176:    Options Database Keys:
1177: +   -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1178: -   -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>

1180:    Level: intermediate

1182:    Notes:
1183:      The block-stash is used for values set with VecSetValuesBlocked() while
1184:      the stash is used for values set with VecSetValues()

1186:      Run with the option -info and look for output of the form
1187:      VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1188:      to determine the appropriate value, MM, to use for size and
1189:      VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1190:      to determine the value, BMM to use for bsize

1192:    Concepts: vector^stash
1193:    Concepts: stash^vector

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

1197: @*/
1198: PetscErrorCode  VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1199: {

1204:   VecStashSetInitialSize_Private(&vec->stash,size);
1205:   VecStashSetInitialSize_Private(&vec->bstash,bsize);
1206:   return(0);
1207: }

1211: /*@
1212:    VecConjugate - Conjugates a vector.

1214:    Logically Collective on Vec

1216:    Input Parameters:
1217: .  x - the vector

1219:    Level: intermediate

1221:    Concepts: vector^conjugate

1223: @*/
1224: PetscErrorCode  VecConjugate(Vec x)
1225: {
1226: #ifdef PETSC_USE_COMPLEX

1232:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1233:   (*x->ops->conjugate)(x);
1234:   /* we need to copy norms here */
1235:   PetscObjectStateIncrease((PetscObject)x);
1236:   return(0);
1237: #else
1238:   return(0);
1239: #endif
1240: }

1244: /*@
1245:    VecPointwiseMult - Computes the componentwise multiplication w = x*y.

1247:    Logically Collective on Vec

1249:    Input Parameters:
1250: .  x, y  - the vectors

1252:    Output Parameter:
1253: .  w - the result

1255:    Level: advanced

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

1259:    Concepts: vector^pointwise multiply

1261: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1262: @*/
1263: PetscErrorCode  VecPointwiseMult(Vec w, Vec x,Vec y)
1264: {

1276:   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");

1278:   PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1279:   (*w->ops->pointwisemult)(w,x,y);
1280:   PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1281:   PetscObjectStateIncrease((PetscObject)w);
1282:   return(0);
1283: }

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

1290:    Logically Collective on Vec

1292:    Input Parameters:
1293: +  x  - the vector
1294: -  rctx - the random number context, formed by PetscRandomCreate(), or PETSC_NULL and
1295:           it will create one internally.

1297:    Output Parameter:
1298: .  x  - the vector

1300:    Example of Usage:
1301: .vb
1302:      PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1303:      VecSetRandom(x,rctx);
1304:      PetscRandomDestroy(rctx);
1305: .ve

1307:    Level: intermediate

1309:    Concepts: vector^setting to random
1310:    Concepts: random^vector

1312: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1313: @*/
1314: PetscErrorCode  VecSetRandom(Vec x,PetscRandom rctx)
1315: {
1317:   PetscRandom    randObj = PETSC_NULL;

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

1325:   if (!rctx) {
1326:     MPI_Comm    comm;
1327:     PetscObjectGetComm((PetscObject)x,&comm);
1328:     PetscRandomCreate(comm,&randObj);
1329:     PetscRandomSetFromOptions(randObj);
1330:     rctx = randObj;
1331:   }

1333:   PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1334:   (*x->ops->setrandom)(x,rctx);
1335:   PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);

1337:   PetscRandomDestroy(&randObj);
1338:   PetscObjectStateIncrease((PetscObject)x);
1339:   return(0);
1340: }

1344: /*@
1345:   VecZeroEntries - puts a 0.0 in each element of a vector

1347:   Logically Collective on Vec

1349:   Input Parameter:
1350: . vec - The vector

1352:   Level: beginner

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

1359: .keywords: Vec, set, options, database
1360: .seealso: VecCreate(),  VecSetOptionsPrefix(), VecSet(), VecSetValues()
1361: @*/
1362: PetscErrorCode  VecZeroEntries(Vec vec)
1363: {
1366:   VecSet(vec,0);
1367:   return(0);
1368: }

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

1376:   Collective on Vec

1378:   Input Parameter:
1379: . vec - The vector

1381:   Level: intermediate

1383: .keywords: Vec, set, options, database, type
1384: .seealso: VecSetFromOptions(), VecSetType()
1385: */
1386: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec)
1387: {
1388:   PetscBool      opt;
1389:   const VecType  defaultType;
1390:   char           typeName[256];
1391:   PetscMPIInt    size;

1395:   if (((PetscObject)vec)->type_name) {
1396:     defaultType = ((PetscObject)vec)->type_name;
1397:   } else {
1398:     MPI_Comm_size(((PetscObject)vec)->comm, &size);
1399:     if (size > 1) {
1400:       defaultType = VECMPI;
1401:     } else {
1402:       defaultType = VECSEQ;
1403:     }
1404:   }

1406:   if (!VecRegisterAllCalled) {VecRegisterAll(PETSC_NULL);}
1407:   PetscOptionsList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1408:   if (opt) {
1409:     VecSetType(vec, typeName);
1410:   } else {
1411:     VecSetType(vec, defaultType);
1412:   }
1413:   return(0);
1414: }

1418: /*@
1419:   VecSetFromOptions - Configures the vector from the options database.

1421:   Collective on Vec

1423:   Input Parameter:
1424: . vec - The vector

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

1429:   Level: beginner

1431:   Concepts: vectors^setting options
1432:   Concepts: vectors^setting type

1434: .keywords: Vec, set, options, database
1435: .seealso: VecCreate(), VecSetOptionsPrefix()
1436: @*/
1437: PetscErrorCode  VecSetFromOptions(Vec vec)
1438: {


1444:   PetscObjectOptionsBegin((PetscObject)vec);
1445:     /* Handle vector type options */
1446:     VecSetTypeFromOptions_Private(vec);

1448:     /* Handle specific vector options */
1449:     if (vec->ops->setfromoptions) {
1450:       (*vec->ops->setfromoptions)(vec);
1451:     }

1453:     /* process any options handlers added with PetscObjectAddOptionsHandler() */
1454:     PetscObjectProcessOptionsHandlers((PetscObject)vec);
1455:   PetscOptionsEnd();

1457:   VecViewFromOptions(vec, ((PetscObject)vec)->name);
1458:   return(0);
1459: }

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

1466:   Collective on Vec

1468:   Input Parameters:
1469: + v - the vector
1470: . n - the local size (or PETSC_DECIDE to have it set)
1471: - N - the global size (or PETSC_DECIDE)

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

1477:   Level: intermediate

1479: .seealso: VecGetSize(), PetscSplitOwnership()
1480: @*/
1481: PetscErrorCode  VecSetSizes(Vec v, PetscInt n, PetscInt N)
1482: {

1487:   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);
1488:   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);
1489:   v->map->n = n;
1490:   v->map->N = N;
1491:   if (v->ops->create) {
1492:     (*v->ops->create)(v);
1493:     v->ops->create = 0;
1494:   }
1495:   return(0);
1496: }

1500: /*@
1501:    VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1502:    and VecSetValuesBlockedLocal().

1504:    Logically Collective on Vec

1506:    Input Parameter:
1507: +  v - the vector
1508: -  bs - the blocksize

1510:    Notes:
1511:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1513:    Level: advanced

1515: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlock(), VecGetBlockSize()

1517:   Concepts: block size^vectors
1518: @*/
1519: PetscErrorCode  VecSetBlockSize(Vec v,PetscInt bs)
1520: {

1525:   if (bs == v->map->bs) return(0);
1526:   PetscLayoutSetBlockSize(v->map,bs);
1527:   v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1528:   return(0);
1529: }

1533: /*@
1534:    VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1535:    and VecSetValuesBlockedLocal().

1537:    Not Collective

1539:    Input Parameter:
1540: .  v - the vector

1542:    Output Parameter:
1543: .  bs - the blocksize

1545:    Notes:
1546:    All vectors obtained by VecDuplicate() inherit the same blocksize.

1548:    Level: advanced

1550: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlock(), VecSetBlockSize()

1552:    Concepts: vector^block size
1553:    Concepts: block^vector

1555: @*/
1556: PetscErrorCode  VecGetBlockSize(Vec v,PetscInt *bs)
1557: {
1562:   PetscLayoutGetBlockSize(v->map,bs);
1563:   return(0);
1564: }

1568: /*@C
1569:    VecSetOptionsPrefix - Sets the prefix used for searching for all
1570:    Vec options in the database.

1572:    Logically Collective on Vec

1574:    Input Parameter:
1575: +  v - the Vec context
1576: -  prefix - the prefix to prepend to all option names

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

1582:    Level: advanced

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

1586: .seealso: VecSetFromOptions()
1587: @*/
1588: PetscErrorCode  VecSetOptionsPrefix(Vec v,const char prefix[])
1589: {

1594:   PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1595:   return(0);
1596: }

1600: /*@C
1601:    VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1602:    Vec options in the database.

1604:    Logically Collective on Vec

1606:    Input Parameters:
1607: +  v - the Vec context
1608: -  prefix - the prefix to prepend to all option names

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

1614:    Level: advanced

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

1618: .seealso: VecGetOptionsPrefix()
1619: @*/
1620: PetscErrorCode  VecAppendOptionsPrefix(Vec v,const char prefix[])
1621: {

1626:   PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1627:   return(0);
1628: }

1632: /*@C
1633:    VecGetOptionsPrefix - Sets the prefix used for searching for all
1634:    Vec options in the database.

1636:    Not Collective

1638:    Input Parameter:
1639: .  v - the Vec context

1641:    Output Parameter:
1642: .  prefix - pointer to the prefix string used

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

1647:    Level: advanced

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

1651: .seealso: VecAppendOptionsPrefix()
1652: @*/
1653: PetscErrorCode  VecGetOptionsPrefix(Vec v,const char *prefix[])
1654: {

1659:   PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1660:   return(0);
1661: }

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

1668:    Collective on Vec

1670:    Input Parameters:
1671: .  v - the Vec context

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

1677:    Level: advanced

1679: .keywords: Vec, setup

1681: .seealso: VecCreate(), VecDestroy()
1682: @*/
1683: PetscErrorCode  VecSetUp(Vec v)
1684: {
1685:   PetscMPIInt    size;

1690:   if (!((PetscObject)v)->type_name) {
1691:     MPI_Comm_size(((PetscObject)v)->comm, &size);
1692:     if (size == 1) {
1693:       VecSetType(v, VECSEQ);
1694:     } else {
1695:       VecSetType(v, VECMPI);
1696:     }
1697:   }
1698:   return(0);
1699: }

1701: /*
1702:     These currently expose the PetscScalar/PetscReal in updating the
1703:     cached norm. If we push those down into the implementation these
1704:     will become independent of PetscScalar/PetscReal
1705: */

1709: /*@
1710:    VecCopy - Copies a vector. y <- x

1712:    Logically Collective on Vec

1714:    Input Parameter:
1715: .  x - the vector

1717:    Output Parameter:
1718: .  y - the copy

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

1724:    Level: beginner

1726: .seealso: VecDuplicate()
1727: @*/
1728: PetscErrorCode  VecCopy(Vec x,Vec y)
1729: {
1730:   PetscBool      flgs[4];
1731:   PetscReal      norms[4] = {0.0,0.0,0.0,0.0};
1733:   PetscInt       i;

1740:   if (x == y) return(0);
1741:   if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1742:   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);

1744: #if !defined(PETSC_USE_MIXED_PRECISION)
1745:   for (i=0; i<4; i++) {
1746:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1747:   }
1748: #endif

1750:   PetscLogEventBegin(VEC_Copy,x,y,0,0);
1751: #if defined(PETSC_USE_MIXED_PRECISION)
1752:   extern PetscErrorCode VecGetArray(Vec,double **);
1753:   extern PetscErrorCode VecRestoreArray(Vec,double **);
1754:   extern PetscErrorCode VecGetArray(Vec,float **);
1755:   extern PetscErrorCode VecRestoreArray(Vec,float **);
1756:   extern PetscErrorCode VecGetArrayRead(Vec,const double **);
1757:   extern PetscErrorCode VecRestoreArrayRead(Vec,const double **);
1758:   extern PetscErrorCode VecGetArrayRead(Vec,const float **);
1759:   extern PetscErrorCode VecRestoreArrayRead(Vec,const float **);
1760:   if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1761:     PetscInt    i,n;
1762:     const float *xx;
1763:     double      *yy;
1764:     VecGetArrayRead(x,&xx);
1765:     VecGetArray(y,&yy);
1766:     VecGetLocalSize(x,&n);
1767:     for (i=0; i<n; i++) {
1768:       yy[i] = xx[i];
1769:     }
1770:     VecRestoreArrayRead(x,&xx);
1771:     VecRestoreArray(y,&yy);
1772:   } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1773:     PetscInt     i,n;
1774:     float        *yy;
1775:     const double *xx;
1776:     VecGetArrayRead(x,&xx);
1777:     VecGetArray(y,&yy);
1778:     VecGetLocalSize(x,&n);
1779:     for (i=0; i<n; i++) {
1780:       yy[i] = (float) xx[i];
1781:     }
1782:     VecRestoreArrayRead(x,&xx);
1783:     VecRestoreArray(y,&yy);
1784:   } else {
1785:     (*x->ops->copy)(x,y);
1786:   }
1787: #else
1788:   (*x->ops->copy)(x,y);
1789: #endif

1791:   PetscObjectStateIncrease((PetscObject)y);
1792: #if !defined(PETSC_USE_MIXED_PRECISION)
1793:   for (i=0; i<4; i++) {
1794:     if (flgs[i]) {
1795:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1796:     }
1797:   }
1798: #endif

1800:   PetscLogEventEnd(VEC_Copy,x,y,0,0);
1801:   return(0);
1802: }

1806: /*@
1807:    VecSwap - Swaps the vectors x and y.

1809:    Logically Collective on Vec

1811:    Input Parameters:
1812: .  x, y  - the vectors

1814:    Level: advanced

1816:    Concepts: vector^swapping values

1818: @*/
1819: PetscErrorCode  VecSwap(Vec x,Vec y)
1820: {
1821:   PetscReal      normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1822:   PetscBool      flgxs[4],flgys[4];
1824:   PetscInt       i;

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

1837:   PetscLogEventBegin(VEC_Swap,x,y,0,0);
1838:   for (i=0; i<4; i++) {
1839:     PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1840:     PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1841:   }
1842:   (*x->ops->swap)(x,y);
1843:   PetscObjectStateIncrease((PetscObject)x);
1844:   PetscObjectStateIncrease((PetscObject)y);
1845:   for (i=0; i<4; i++) {
1846:     if (flgxs[i]) {
1847:       PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1848:     }
1849:     if (flgys[i]) {
1850:       PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1851:     }
1852:   }
1853:   PetscLogEventEnd(VEC_Swap,x,y,0,0);
1854:   return(0);
1855: }

1859: /*@
1860:    VecStashView - Prints the entries in the vector stash and block stash.

1862:    Collective on Vec

1864:    Input Parameters:
1865: +  v - the vector
1866: -  viewer - the viewer

1868:    Level: advanced

1870:    Concepts: vector^stash
1871:    Concepts: stash^vector

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

1875: @*/
1876: PetscErrorCode  VecStashView(Vec v,PetscViewer viewer)
1877: {
1879:   PetscMPIInt    rank;
1880:   PetscInt       i,j;
1881:   PetscBool      match;
1882:   VecStash       *s;
1883:   PetscScalar    val;


1890:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1891:   if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1892:   PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1893:   MPI_Comm_rank(((PetscObject)v)->comm,&rank);
1894:   s = &v->bstash;

1896:   /* print block stash */
1897:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1898:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1899:   for (i=0; i<s->n; i++) {
1900:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1901:     for (j=0; j<s->bs; j++) {
1902:       val = s->array[i*s->bs+j];
1903: #if defined(PETSC_USE_COMPLEX)
1904:       PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1905: #else
1906:       PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1907: #endif
1908:     }
1909:     PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1910:   }
1911:   PetscViewerFlush(viewer);

1913:   s = &v->stash;

1915:   /* print basic stash */
1916:   PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1917:   for (i=0; i<s->n; i++) {
1918:     val = s->array[i];
1919: #if defined(PETSC_USE_COMPLEX)
1920:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1921: #else
1922:     PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1923: #endif
1924:   }
1925:   PetscViewerFlush(viewer);
1926:   PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);

1928:   PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1929:   return(0);
1930: }