Actual source code: vector.c
petsc-3.7.5 2017-01-01
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;
18: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
19: PetscLogEvent VEC_CUDACopyFromGPUSome, VEC_CUDACopyToGPUSome;
21: extern PetscErrorCode VecStashGetInfo_Private(VecStash*,PetscInt*,PetscInt*);
24: /*@
25: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
26: to be communicated to other processors during the VecAssemblyBegin/End() process
28: Not collective
30: Input Parameter:
31: . vec - the vector
33: Output Parameters:
34: + nstash - the size of the stash
35: . reallocs - the number of additional mallocs incurred.
36: . bnstash - the size of the block stash
37: - breallocs - the number of additional mallocs incurred.in the block stash
39: Level: advanced
41: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
43: @*/
44: PetscErrorCode VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
45: {
49: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
50: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
51: return(0);
52: }
56: /*@
57: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
58: by the routine VecSetValuesLocal() to allow users to insert vector entries
59: using a local (per-processor) numbering.
61: Logically Collective on Vec
63: Input Parameters:
64: + x - vector
65: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
67: Notes:
68: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
70: Level: intermediate
72: Concepts: vector^setting values with local numbering
74: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
75: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
76: @*/
77: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
78: {
85: if (x->ops->setlocaltoglobalmapping) {
86: (*x->ops->setlocaltoglobalmapping)(x,mapping);
87: } else {
88: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
89: }
90: return(0);
91: }
95: /*@
96: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
98: Not Collective
100: Input Parameter:
101: . X - the vector
103: Output Parameter:
104: . mapping - the mapping
106: Level: advanced
108: Concepts: vectors^local to global mapping
109: Concepts: local to global mapping^for vectors
111: .seealso: VecSetValuesLocal()
112: @*/
113: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
114: {
119: *mapping = X->map->mapping;
120: return(0);
121: }
125: /*@
126: VecAssemblyBegin - Begins assembling the vector. This routine should
127: be called after completing all calls to VecSetValues().
129: Collective on Vec
131: Input Parameter:
132: . vec - the vector
134: Level: beginner
136: Concepts: assembly^vectors
138: .seealso: VecAssemblyEnd(), VecSetValues()
139: @*/
140: PetscErrorCode VecAssemblyBegin(Vec vec)
141: {
147: VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
148: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
149: if (vec->ops->assemblybegin) {
150: (*vec->ops->assemblybegin)(vec);
151: }
152: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
153: PetscObjectStateIncrease((PetscObject)vec);
154: return(0);
155: }
159: /*@
160: VecAssemblyEnd - Completes assembling the vector. This routine should
161: be called after VecAssemblyBegin().
163: Collective on Vec
165: Input Parameter:
166: . vec - the vector
168: Options Database Keys:
169: + -vec_view - Prints vector in ASCII format
170: . -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
171: . -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
172: . -vec_view draw - Activates vector viewing using drawing tools
173: . -display <name> - Sets display name (default is host)
174: . -draw_pause <sec> - Sets number of seconds to pause after display
175: - -vec_view socket - Activates vector viewing using a socket
177: Level: beginner
179: .seealso: VecAssemblyBegin(), VecSetValues()
180: @*/
181: PetscErrorCode VecAssemblyEnd(Vec vec)
182: {
187: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
189: if (vec->ops->assemblyend) {
190: (*vec->ops->assemblyend)(vec);
191: }
192: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
193: VecViewFromOptions(vec,NULL,"-vec_view");
194: return(0);
195: }
199: /*@
200: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
202: Logically Collective on Vec
204: Input Parameters:
205: . x, y - the vectors
207: Output Parameter:
208: . w - the result
210: Level: advanced
212: Notes: any subset of the x, y, and w may be the same vector.
213: For complex numbers compares only the real part
215: Concepts: vector^pointwise multiply
217: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
218: @*/
219: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
220: {
232: 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");
233: 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");
235: (*w->ops->pointwisemax)(w,x,y);
236: PetscObjectStateIncrease((PetscObject)w);
237: return(0);
238: }
243: /*@
244: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
246: Logically Collective on Vec
248: Input Parameters:
249: . x, y - the vectors
251: Output Parameter:
252: . w - the result
254: Level: advanced
256: Notes: any subset of the x, y, and w may be the same vector.
257: For complex numbers compares only the real part
259: Concepts: vector^pointwise multiply
261: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
262: @*/
263: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
264: {
276: 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");
277: 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");
279: (*w->ops->pointwisemin)(w,x,y);
280: PetscObjectStateIncrease((PetscObject)w);
281: return(0);
282: }
286: /*@
287: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
289: Logically Collective on Vec
291: Input Parameters:
292: . x, y - the vectors
294: Output Parameter:
295: . w - the result
297: Level: advanced
299: Notes: any subset of the x, y, and w may be the same vector.
301: Concepts: vector^pointwise multiply
303: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
304: @*/
305: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
306: {
318: 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");
319: 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");
321: (*w->ops->pointwisemaxabs)(w,x,y);
322: PetscObjectStateIncrease((PetscObject)w);
323: return(0);
324: }
328: /*@
329: VecPointwiseDivide - Computes the componentwise division w = x/y.
331: Logically Collective on Vec
333: Input Parameters:
334: . x, y - the vectors
336: Output Parameter:
337: . w - the result
339: Level: advanced
341: Notes: any subset of the x, y, and w may be the same vector.
343: Concepts: vector^pointwise divide
345: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
346: @*/
347: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
348: {
360: 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");
361: 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");
363: (*w->ops->pointwisedivide)(w,x,y);
364: PetscObjectStateIncrease((PetscObject)w);
365: return(0);
366: }
371: /*@
372: VecDuplicate - Creates a new vector of the same type as an existing vector.
374: Collective on Vec
376: Input Parameters:
377: . v - a vector to mimic
379: Output Parameter:
380: . newv - location to put new vector
382: Notes:
383: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
384: for the new vector. Use VecCopy() to copy a vector.
386: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
387: vectors.
389: Level: beginner
391: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
392: @*/
393: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
394: {
401: (*v->ops->duplicate)(v,newv);
402: PetscObjectStateIncrease((PetscObject)*newv);
403: return(0);
404: }
408: /*@
409: VecDestroy - Destroys a vector.
411: Collective on Vec
413: Input Parameters:
414: . v - the vector
416: Level: beginner
418: .seealso: VecDuplicate(), VecDestroyVecs()
419: @*/
420: PetscErrorCode VecDestroy(Vec *v)
421: {
425: if (!*v) return(0);
427: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
429: PetscObjectSAWsViewOff((PetscObject)*v);
430: /* destroy the internal part */
431: if ((*v)->ops->destroy) {
432: (*(*v)->ops->destroy)(*v);
433: }
434: /* destroy the external/common part */
435: PetscLayoutDestroy(&(*v)->map);
436: PetscHeaderDestroy(v);
437: return(0);
438: }
442: /*@C
443: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
445: Collective on Vec
447: Input Parameters:
448: + m - the number of vectors to obtain
449: - v - a vector to mimic
451: Output Parameter:
452: . V - location to put pointer to array of vectors
454: Notes:
455: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
456: vector.
458: Fortran Note:
459: The Fortran interface is slightly different from that given below, it
460: requires one to pass in V a Vec (integer) array of size at least m.
461: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
463: Level: intermediate
465: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
466: @*/
467: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
468: {
475: (*v->ops->duplicatevecs)(v,m,V);
476: return(0);
477: }
481: /*@C
482: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
484: Collective on Vec
486: Input Parameters:
487: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
488: - m - the number of vectors previously obtained, if zero no vectors are destroyed
490: Fortran Note:
491: The Fortran interface is slightly different from that given below.
492: See the Fortran chapter of the users manual
494: Level: intermediate
496: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
497: @*/
498: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
499: {
504: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
505: if (!m || !*vv) {*vv = NULL; return(0);}
508: (*(**vv)->ops->destroyvecs)(m,*vv);
509: *vv = NULL;
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 - for sequential vectors
527: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
528: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
530: You can change the format the vector is printed using the
531: option PetscViewerPushFormat().
533: The user can open alternative visualization contexts with
534: + PetscViewerASCIIOpen() - Outputs vector to a specified file
535: . PetscViewerBinaryOpen() - Outputs vector in binary to a
536: specified file; corresponding input uses VecLoad()
537: . PetscViewerDrawOpen() - Outputs vector to an X window display
538: - PetscViewerSocketOpen() - Outputs vector to Socket viewer
540: The user can call PetscViewerPushFormat() to specify the output
541: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
542: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
543: + PETSC_VIEWER_DEFAULT - default, prints vector contents
544: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
545: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
546: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
547: format common among all vector types
549: Notes: You can pass any number of vector objects, or other PETSc objects to the same viewer.
551: Notes for binary viewer: If you pass multiply vectors to a binary viewer you can read them back in in the same order
552: $ with VecLoad().
553: $
554: $ If the blocksize of the vector is greater than one then you must provide a unique prefix to
555: $ the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
556: $ vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
557: $ information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
558: $ filename. If you copy the binary file, make sure you copy the associated .info file with it.
560: Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
561: $ for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
562: $ obviously) several times, you must change its name each time before calling the VecView(). The name you use
563: $ here should equal the name that you use in the Vec object that you use with VecLoad().
565: See the manual page for VecLoad() on the exact format the binary viewer stores
566: the values in the file.
568: Level: beginner
570: Concepts: vector^printing
571: Concepts: vector^saving to disk
573: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
574: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
575: PetscRealView(), PetscScalarView(), PetscIntView()
576: @*/
577: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
578: {
579: PetscErrorCode ierr;
580: PetscBool iascii;
581: PetscViewerFormat format;
586: if (!viewer) {
587: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
588: }
591: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
593: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
594: PetscViewerGetFormat(viewer,&format);
595: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
596: if (iascii) {
597: PetscInt rows,bs;
599: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
600: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
601: PetscViewerASCIIPushTab(viewer);
602: VecGetSize(vec,&rows);
603: VecGetBlockSize(vec,&bs);
604: if (bs != 1) {
605: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
606: } else {
607: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
608: }
609: PetscViewerASCIIPopTab(viewer);
610: }
611: }
612: VecLockPush(vec);
613: if (format == PETSC_VIEWER_NATIVE && vec->ops->viewnative) {
614: (*vec->ops->viewnative)(vec,viewer);
615: } else {
616: (*vec->ops->view)(vec,viewer);
617: }
618: VecLockPop(vec);
619: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
620: return(0);
621: }
623: #if defined(PETSC_USE_DEBUG)
624: #include <../src/sys/totalview/tv_data_display.h>
625: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
626: {
627: const PetscScalar *values;
628: char type[32];
629: PetscErrorCode ierr;
632: TV_add_row("Local rows", "int", &v->map->n);
633: TV_add_row("Global rows", "int", &v->map->N);
634: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
635: VecGetArrayRead((Vec)v,&values);
636: PetscSNPrintf(type,32,"double[%d]",v->map->n);
637: TV_add_row("values",type, values);
638: VecRestoreArrayRead((Vec)v,&values);
639: return TV_format_OK;
640: }
641: #endif
645: /*@
646: VecGetSize - Returns the global number of elements of the vector.
648: Not Collective
650: Input Parameter:
651: . x - the vector
653: Output Parameters:
654: . size - the global length of the vector
656: Level: beginner
658: Concepts: vector^local size
660: .seealso: VecGetLocalSize()
661: @*/
662: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
663: {
670: (*x->ops->getsize)(x,size);
671: return(0);
672: }
676: /*@
677: VecGetLocalSize - Returns the number of elements of the vector stored
678: in local memory. This routine may be implementation dependent, so use
679: with care.
681: Not Collective
683: Input Parameter:
684: . x - the vector
686: Output Parameter:
687: . size - the length of the local piece of the vector
689: Level: beginner
691: Concepts: vector^size
693: .seealso: VecGetSize()
694: @*/
695: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
696: {
703: (*x->ops->getlocalsize)(x,size);
704: return(0);
705: }
709: /*@C
710: VecGetOwnershipRange - Returns the range of indices owned by
711: this processor, assuming that the vectors are laid out with the
712: first n1 elements on the first processor, next n2 elements on the
713: second, etc. For certain parallel layouts this range may not be
714: well defined.
716: Not Collective
718: Input Parameter:
719: . x - the vector
721: Output Parameters:
722: + low - the first local element, pass in NULL if not interested
723: - high - one more than the last local element, pass in NULL if not interested
725: Note:
726: The high argument is one more than the last element stored locally.
728: Fortran: NULL_INTEGER should be used instead of NULL
730: Level: beginner
732: Concepts: ownership^of vectors
733: Concepts: vector^ownership of elements
735: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
736: @*/
737: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
738: {
744: if (low) *low = x->map->rstart;
745: if (high) *high = x->map->rend;
746: return(0);
747: }
751: /*@C
752: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
753: assuming that the vectors are laid out with the
754: first n1 elements on the first processor, next n2 elements on the
755: second, etc. For certain parallel layouts this range may not be
756: well defined.
758: Not Collective
760: Input Parameter:
761: . x - the vector
763: Output Parameters:
764: . range - array of length size+1 with the start and end+1 for each process
766: Note:
767: The high argument is one more than the last element stored locally.
769: Fortran: You must PASS in an array of length size+1
771: Level: beginner
773: Concepts: ownership^of vectors
774: Concepts: vector^ownership of elements
776: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
777: @*/
778: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
779: {
785: PetscLayoutGetRanges(x->map,ranges);
786: return(0);
787: }
791: /*@
792: VecSetOption - Sets an option for controling a vector's behavior.
794: Collective on Vec
796: Input Parameter:
797: + x - the vector
798: . op - the option
799: - flag - turn the option on or off
801: Supported Options:
802: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
803: entries destined to be stored on a separate processor. This can be used
804: to eliminate the global reduction in the VecAssemblyXXXX() if you know
805: that you have only used VecSetValues() to set local elements
806: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
807: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
808: ignored.
809: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
810: entries will always be a subset (possibly equal) of the off-process entries set on the
811: first assembly. This reuses the communication pattern, thus avoiding a global reduction.
812: Subsequent assemblies setting off-process values should use the same InsertMode as the
813: first assembly.
815: Developer Note:
816: The InsertMode restriction could be removed by packing the stash messages out of place.
818: Level: intermediate
820: @*/
821: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
822: {
828: if (x->ops->setoption) {
829: (*x->ops->setoption)(x,op,flag);
830: }
831: return(0);
832: }
836: /* Default routines for obtaining and releasing; */
837: /* may be used by any implementation */
838: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
839: {
841: PetscInt i;
846: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
847: PetscMalloc1(m,V);
848: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
849: return(0);
850: }
854: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
855: {
857: PetscInt i;
861: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
862: PetscFree(v);
863: return(0);
864: }
868: /*@
869: VecResetArray - Resets a vector to use its default memory. Call this
870: after the use of VecPlaceArray().
872: Not Collective
874: Input Parameters:
875: . vec - the vector
877: Level: developer
879: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
881: @*/
882: PetscErrorCode VecResetArray(Vec vec)
883: {
889: if (vec->ops->resetarray) {
890: (*vec->ops->resetarray)(vec);
891: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
892: PetscObjectStateIncrease((PetscObject)vec);
893: return(0);
894: }
898: /*@C
899: VecLoad - Loads a vector that has been stored in binary or HDF5 format
900: with VecView().
902: Collective on PetscViewer
904: Input Parameters:
905: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
906: some related function before a call to VecLoad().
907: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
908: HDF5 file viewer, obtained from PetscViewerHDF5Open()
910: Level: intermediate
912: Notes:
913: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
914: before calling this.
916: The input file must contain the full global vector, as
917: written by the routine VecView().
919: If the type or size of newvec is not set before a call to VecLoad, PETSc
920: sets the type and the local and global sizes. If type and/or
921: sizes are already set, then the same are used.
923: If using binary and the blocksize of the vector is greater than one then you must provide a unique prefix to
924: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
925: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
926: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
927: filename. If you copy the binary file, make sure you copy the associated .info file with it.
929: If using HDF5, you must assign the Vec the same name as was used in the Vec
930: that was stored in the file using PetscObjectSetName(). Otherwise you will
931: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
933: Notes for advanced users:
934: Most users should not need to know the details of the binary storage
935: format, since VecLoad() and VecView() completely hide these details.
936: But for anyone who's interested, the standard binary matrix storage
937: format is
938: .vb
939: int VEC_FILE_CLASSID
940: int number of rows
941: PetscScalar *values of all entries
942: .ve
944: In addition, PETSc automatically does the byte swapping for
945: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
946: linux, Windows and the paragon; thus if you write your own binary
947: read/write routines you have to swap the bytes; see PetscBinaryRead()
948: and PetscBinaryWrite() to see how this may be done.
950: Concepts: vector^loading from file
952: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
953: @*/
954: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
955: {
957: PetscBool isbinary,ishdf5;
958: PetscViewerFormat format;
963: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
964: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
965: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
967: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
968: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
969: VecSetType(newvec, VECSTANDARD);
970: }
971: PetscViewerGetFormat(viewer,&format);
972: if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
973: (*newvec->ops->loadnative)(newvec,viewer);
974: } else {
975: (*newvec->ops->load)(newvec,viewer);
976: }
977: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
978: return(0);
979: }
984: /*@
985: VecReciprocal - Replaces each component of a vector by its reciprocal.
987: Logically Collective on Vec
989: Input Parameter:
990: . vec - the vector
992: Output Parameter:
993: . vec - the vector reciprocal
995: Level: intermediate
997: Concepts: vector^reciprocal
999: .seealso: VecLog(), VecExp(), VecSqrtAbs()
1001: @*/
1002: PetscErrorCode VecReciprocal(Vec vec)
1003: {
1009: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1010: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
1011: (*vec->ops->reciprocal)(vec);
1012: PetscObjectStateIncrease((PetscObject)vec);
1013: return(0);
1014: }
1018: /*@C
1019: VecSetOperation - Allows user to set a vector operation.
1021: Logically Collective on Vec
1023: Input Parameters:
1024: + vec - the vector
1025: . op - the name of the operation
1026: - f - the function that provides the operation.
1028: Level: advanced
1030: Usage:
1031: $ PetscErrorCode userview(Vec,PetscViewer);
1032: $ VecCreateMPI(comm,m,M,&x);
1033: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
1035: Notes:
1036: See the file include/petscvec.h for a complete list of matrix
1037: operations, which all have the form VECOP_<OPERATION>, where
1038: <OPERATION> is the name (in all capital letters) of the
1039: user interface routine (e.g., VecView() -> VECOP_VIEW).
1041: This function is not currently available from Fortran.
1043: .keywords: vector, set, operation
1045: .seealso: VecCreate(), MatShellSetOperation()
1046: @*/
1047: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1048: {
1051: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1052: vec->ops->viewnative = vec->ops->view;
1053: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1054: vec->ops->loadnative = vec->ops->load;
1055: }
1056: (((void(**)(void))vec->ops)[(int)op]) = f;
1057: return(0);
1058: }
1063: /*@
1064: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1065: used during the assembly process to store values that belong to
1066: other processors.
1068: Not Collective, different processes can have different size stashes
1070: Input Parameters:
1071: + vec - the vector
1072: . size - the initial size of the stash.
1073: - bsize - the initial size of the block-stash(if used).
1075: Options Database Keys:
1076: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1077: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1079: Level: intermediate
1081: Notes:
1082: The block-stash is used for values set with VecSetValuesBlocked() while
1083: the stash is used for values set with VecSetValues()
1085: Run with the option -info and look for output of the form
1086: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1087: to determine the appropriate value, MM, to use for size and
1088: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1089: to determine the value, BMM to use for bsize
1091: Concepts: vector^stash
1092: Concepts: stash^vector
1094: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1096: @*/
1097: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1098: {
1103: VecStashSetInitialSize_Private(&vec->stash,size);
1104: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1105: return(0);
1106: }
1110: /*@
1111: VecConjugate - Conjugates a vector.
1113: Logically Collective on Vec
1115: Input Parameters:
1116: . x - the vector
1118: Level: intermediate
1120: Concepts: vector^conjugate
1122: @*/
1123: PetscErrorCode VecConjugate(Vec x)
1124: {
1125: #if defined(PETSC_USE_COMPLEX)
1131: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1132: (*x->ops->conjugate)(x);
1133: /* we need to copy norms here */
1134: PetscObjectStateIncrease((PetscObject)x);
1135: return(0);
1136: #else
1137: return(0);
1138: #endif
1139: }
1143: /*@
1144: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1146: Logically Collective on Vec
1148: Input Parameters:
1149: . x, y - the vectors
1151: Output Parameter:
1152: . w - the result
1154: Level: advanced
1156: Notes: any subset of the x, y, and w may be the same vector.
1158: Concepts: vector^pointwise multiply
1160: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1161: @*/
1162: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1163: {
1175: 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");
1177: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1178: (*w->ops->pointwisemult)(w,x,y);
1179: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1180: PetscObjectStateIncrease((PetscObject)w);
1181: return(0);
1182: }
1186: /*@
1187: VecSetRandom - Sets all components of a vector to random numbers.
1189: Logically Collective on Vec
1191: Input Parameters:
1192: + x - the vector
1193: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1194: it will create one internally.
1196: Output Parameter:
1197: . x - the vector
1199: Example of Usage:
1200: .vb
1201: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1202: VecSetRandom(x,rctx);
1203: PetscRandomDestroy(rctx);
1204: .ve
1206: Level: intermediate
1208: Concepts: vector^setting to random
1209: Concepts: random^vector
1211: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1212: @*/
1213: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1214: {
1216: PetscRandom randObj = NULL;
1222: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1224: if (!rctx) {
1225: MPI_Comm comm;
1226: PetscObjectGetComm((PetscObject)x,&comm);
1227: PetscRandomCreate(comm,&randObj);
1228: PetscRandomSetFromOptions(randObj);
1229: rctx = randObj;
1230: }
1232: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1233: (*x->ops->setrandom)(x,rctx);
1234: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1236: PetscRandomDestroy(&randObj);
1237: PetscObjectStateIncrease((PetscObject)x);
1238: return(0);
1239: }
1243: /*@
1244: VecZeroEntries - puts a 0.0 in each element of a vector
1246: Logically Collective on Vec
1248: Input Parameter:
1249: . vec - The vector
1251: Level: beginner
1253: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1254: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1255: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1256: this routine should not exist.
1258: .keywords: Vec, set, options, database
1259: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1260: @*/
1261: PetscErrorCode VecZeroEntries(Vec vec)
1262: {
1266: VecSet(vec,0);
1267: return(0);
1268: }
1272: /*
1273: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1274: processor and a PETSc MPI vector on more than one processor.
1276: Collective on Vec
1278: Input Parameter:
1279: . vec - The vector
1281: Level: intermediate
1283: .keywords: Vec, set, options, database, type
1284: .seealso: VecSetFromOptions(), VecSetType()
1285: */
1286: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1287: {
1288: PetscBool opt;
1289: VecType defaultType;
1290: char typeName[256];
1291: PetscMPIInt size;
1295: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1296: else {
1297: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1298: if (size > 1) defaultType = VECMPI;
1299: else defaultType = VECSEQ;
1300: }
1302: VecRegisterAll();
1303: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1304: if (opt) {
1305: VecSetType(vec, typeName);
1306: } else {
1307: VecSetType(vec, defaultType);
1308: }
1309: return(0);
1310: }
1314: /*@
1315: VecSetFromOptions - Configures the vector from the options database.
1317: Collective on Vec
1319: Input Parameter:
1320: . vec - The vector
1322: Notes: To see all options, run your program with the -help option, or consult the users manual.
1323: Must be called after VecCreate() but before the vector is used.
1325: Level: beginner
1327: Concepts: vectors^setting options
1328: Concepts: vectors^setting type
1330: .keywords: Vec, set, options, database
1331: .seealso: VecCreate(), VecSetOptionsPrefix()
1332: @*/
1333: PetscErrorCode VecSetFromOptions(Vec vec)
1334: {
1340: PetscObjectOptionsBegin((PetscObject)vec);
1341: /* Handle vector type options */
1342: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1344: /* Handle specific vector options */
1345: if (vec->ops->setfromoptions) {
1346: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1347: }
1349: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1350: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1351: PetscOptionsEnd();
1352: return(0);
1353: }
1357: /*@
1358: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1360: Collective on Vec
1362: Input Parameters:
1363: + v - the vector
1364: . n - the local size (or PETSC_DECIDE to have it set)
1365: - N - the global size (or PETSC_DECIDE)
1367: Notes:
1368: n and N cannot be both PETSC_DECIDE
1369: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1371: Level: intermediate
1373: .seealso: VecGetSize(), PetscSplitOwnership()
1374: @*/
1375: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1376: {
1382: 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);
1383: 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);
1384: v->map->n = n;
1385: v->map->N = N;
1386: if (v->ops->create) {
1387: (*v->ops->create)(v);
1388: v->ops->create = 0;
1389: }
1390: return(0);
1391: }
1395: /*@
1396: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1397: and VecSetValuesBlockedLocal().
1399: Logically Collective on Vec
1401: Input Parameter:
1402: + v - the vector
1403: - bs - the blocksize
1405: Notes:
1406: All vectors obtained by VecDuplicate() inherit the same blocksize.
1408: Level: advanced
1410: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1412: Concepts: block size^vectors
1413: @*/
1414: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1415: {
1420: if (bs < 0 || bs == v->map->bs) return(0);
1422: PetscLayoutSetBlockSize(v->map,bs);
1423: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1424: return(0);
1425: }
1429: /*@
1430: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1431: and VecSetValuesBlockedLocal().
1433: Not Collective
1435: Input Parameter:
1436: . v - the vector
1438: Output Parameter:
1439: . bs - the blocksize
1441: Notes:
1442: All vectors obtained by VecDuplicate() inherit the same blocksize.
1444: Level: advanced
1446: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1448: Concepts: vector^block size
1449: Concepts: block^vector
1451: @*/
1452: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1453: {
1459: PetscLayoutGetBlockSize(v->map,bs);
1460: return(0);
1461: }
1465: /*@C
1466: VecSetOptionsPrefix - Sets the prefix used for searching for all
1467: Vec options in the database.
1469: Logically Collective on Vec
1471: Input Parameter:
1472: + v - the Vec context
1473: - prefix - the prefix to prepend to all option names
1475: Notes:
1476: A hyphen (-) must NOT be given at the beginning of the prefix name.
1477: The first character of all runtime options is AUTOMATICALLY the hyphen.
1479: Level: advanced
1481: .keywords: Vec, set, options, prefix, database
1483: .seealso: VecSetFromOptions()
1484: @*/
1485: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1486: {
1491: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1492: return(0);
1493: }
1497: /*@C
1498: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1499: Vec options in the database.
1501: Logically Collective on Vec
1503: Input Parameters:
1504: + v - the Vec context
1505: - prefix - the prefix to prepend to all option names
1507: Notes:
1508: A hyphen (-) must NOT be given at the beginning of the prefix name.
1509: The first character of all runtime options is AUTOMATICALLY the hyphen.
1511: Level: advanced
1513: .keywords: Vec, append, options, prefix, database
1515: .seealso: VecGetOptionsPrefix()
1516: @*/
1517: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1518: {
1523: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1524: return(0);
1525: }
1529: /*@C
1530: VecGetOptionsPrefix - Sets the prefix used for searching for all
1531: Vec options in the database.
1533: Not Collective
1535: Input Parameter:
1536: . v - the Vec context
1538: Output Parameter:
1539: . prefix - pointer to the prefix string used
1541: Notes: On the fortran side, the user should pass in a string 'prefix' of
1542: sufficient length to hold the prefix.
1544: Level: advanced
1546: .keywords: Vec, get, options, prefix, database
1548: .seealso: VecAppendOptionsPrefix()
1549: @*/
1550: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1551: {
1556: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1557: return(0);
1558: }
1562: /*@
1563: VecSetUp - Sets up the internal vector data structures for the later use.
1565: Collective on Vec
1567: Input Parameters:
1568: . v - the Vec context
1570: Notes:
1571: For basic use of the Vec classes the user need not explicitly call
1572: VecSetUp(), since these actions will happen automatically.
1574: Level: advanced
1576: .keywords: Vec, setup
1578: .seealso: VecCreate(), VecDestroy()
1579: @*/
1580: PetscErrorCode VecSetUp(Vec v)
1581: {
1582: PetscMPIInt size;
1587: if (!((PetscObject)v)->type_name) {
1588: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1589: if (size == 1) {
1590: VecSetType(v, VECSEQ);
1591: } else {
1592: VecSetType(v, VECMPI);
1593: }
1594: }
1595: return(0);
1596: }
1598: /*
1599: These currently expose the PetscScalar/PetscReal in updating the
1600: cached norm. If we push those down into the implementation these
1601: will become independent of PetscScalar/PetscReal
1602: */
1606: /*@
1607: VecCopy - Copies a vector. y <- x
1609: Logically Collective on Vec
1611: Input Parameter:
1612: . x - the vector
1614: Output Parameter:
1615: . y - the copy
1617: Notes:
1618: For default parallel PETSc vectors, both x and y must be distributed in
1619: the same manner; local copies are done.
1621: Level: beginner
1623: .seealso: VecDuplicate()
1624: @*/
1625: PetscErrorCode VecCopy(Vec x,Vec y)
1626: {
1627: PetscBool flgs[4];
1628: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1630: PetscInt i;
1637: if (x == y) return(0);
1638: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1639: 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);
1640: VecLocked(y,2);
1642: #if !defined(PETSC_USE_MIXED_PRECISION)
1643: for (i=0; i<4; i++) {
1644: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1645: }
1646: #endif
1648: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1649: #if defined(PETSC_USE_MIXED_PRECISION)
1650: extern PetscErrorCode VecGetArray(Vec,double**);
1651: extern PetscErrorCode VecRestoreArray(Vec,double**);
1652: extern PetscErrorCode VecGetArray(Vec,float**);
1653: extern PetscErrorCode VecRestoreArray(Vec,float**);
1654: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1655: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1656: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1657: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1658: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1659: PetscInt i,n;
1660: const float *xx;
1661: double *yy;
1662: VecGetArrayRead(x,&xx);
1663: VecGetArray(y,&yy);
1664: VecGetLocalSize(x,&n);
1665: for (i=0; i<n; i++) yy[i] = xx[i];
1666: VecRestoreArrayRead(x,&xx);
1667: VecRestoreArray(y,&yy);
1668: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1669: PetscInt i,n;
1670: float *yy;
1671: const double *xx;
1672: VecGetArrayRead(x,&xx);
1673: VecGetArray(y,&yy);
1674: VecGetLocalSize(x,&n);
1675: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1676: VecRestoreArrayRead(x,&xx);
1677: VecRestoreArray(y,&yy);
1678: } else {
1679: (*x->ops->copy)(x,y);
1680: }
1681: #else
1682: (*x->ops->copy)(x,y);
1683: #endif
1685: PetscObjectStateIncrease((PetscObject)y);
1686: #if !defined(PETSC_USE_MIXED_PRECISION)
1687: for (i=0; i<4; i++) {
1688: if (flgs[i]) {
1689: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1690: }
1691: }
1692: #endif
1694: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1695: return(0);
1696: }
1700: /*@
1701: VecSwap - Swaps the vectors x and y.
1703: Logically Collective on Vec
1705: Input Parameters:
1706: . x, y - the vectors
1708: Level: advanced
1710: Concepts: vector^swapping values
1712: @*/
1713: PetscErrorCode VecSwap(Vec x,Vec y)
1714: {
1715: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1716: PetscBool flgxs[4],flgys[4];
1718: PetscInt i;
1726: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1727: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1728: if (x->map->N != y->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1729: if (x->map->n != y->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1731: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1732: for (i=0; i<4; i++) {
1733: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1734: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1735: }
1736: (*x->ops->swap)(x,y);
1737: PetscObjectStateIncrease((PetscObject)x);
1738: PetscObjectStateIncrease((PetscObject)y);
1739: for (i=0; i<4; i++) {
1740: if (flgxs[i]) {
1741: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1742: }
1743: if (flgys[i]) {
1744: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1745: }
1746: }
1747: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1748: return(0);
1749: }
1753: /*
1754: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1756: Collective on VecStash
1758: Input Parameters:
1759: + obj - the VecStash object
1760: . bobj - optional other object that provides the prefix
1761: - optionname - option to activate viewing
1763: Level: intermediate
1765: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1767: */
1768: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1769: {
1770: PetscErrorCode ierr;
1771: PetscViewer viewer;
1772: PetscBool flg;
1773: PetscViewerFormat format;
1774: char *prefix;
1777: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1778: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),prefix,optionname,&viewer,&format,&flg);
1779: if (flg) {
1780: PetscViewerPushFormat(viewer,format);
1781: VecStashView(obj,viewer);
1782: PetscViewerPopFormat(viewer);
1783: PetscViewerDestroy(&viewer);
1784: }
1785: return(0);
1786: }
1790: /*@
1791: VecStashView - Prints the entries in the vector stash and block stash.
1793: Collective on Vec
1795: Input Parameters:
1796: + v - the vector
1797: - viewer - the viewer
1799: Level: advanced
1801: Concepts: vector^stash
1802: Concepts: stash^vector
1804: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1806: @*/
1807: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1808: {
1810: PetscMPIInt rank;
1811: PetscInt i,j;
1812: PetscBool match;
1813: VecStash *s;
1814: PetscScalar val;
1821: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1822: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1823: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1824: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1825: s = &v->bstash;
1827: /* print block stash */
1828: PetscViewerASCIIPushSynchronized(viewer);
1829: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1830: for (i=0; i<s->n; i++) {
1831: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1832: for (j=0; j<s->bs; j++) {
1833: val = s->array[i*s->bs+j];
1834: #if defined(PETSC_USE_COMPLEX)
1835: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1836: #else
1837: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1838: #endif
1839: }
1840: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1841: }
1842: PetscViewerFlush(viewer);
1844: s = &v->stash;
1846: /* print basic stash */
1847: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1848: for (i=0; i<s->n; i++) {
1849: val = s->array[i];
1850: #if defined(PETSC_USE_COMPLEX)
1851: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1852: #else
1853: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1854: #endif
1855: }
1856: PetscViewerFlush(viewer);
1857: PetscViewerASCIIPopSynchronized(viewer);
1858: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1859: return(0);
1860: }
1864: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1865: {
1866: PetscInt i,N,rstart,rend;
1868: PetscScalar *xx;
1869: PetscReal *xreal;
1870: PetscBool iset;
1873: VecGetOwnershipRange(v,&rstart,&rend);
1874: VecGetSize(v,&N);
1875: PetscCalloc1(N,&xreal);
1876: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1877: if (iset) {
1878: VecGetArray(v,&xx);
1879: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1880: VecRestoreArray(v,&xx);
1881: }
1882: PetscFree(xreal);
1883: if (set) *set = iset;
1884: return(0);
1885: }
1889: /*@
1890: VecGetLayout - get PetscLayout describing vector layout
1892: Not Collective
1894: Input Arguments:
1895: . x - the vector
1897: Output Arguments:
1898: . map - the layout
1900: Level: developer
1902: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1903: @*/
1904: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1905: {
1909: *map = x->map;
1910: return(0);
1911: }
1915: /*@
1916: VecSetLayout - set PetscLayout describing vector layout
1918: Not Collective
1920: Input Arguments:
1921: + x - the vector
1922: - map - the layout
1924: Notes:
1925: It is normally only valid to replace the layout with a layout known to be equivalent.
1927: Level: developer
1929: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1930: @*/
1931: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1932: {
1937: PetscLayoutReference(map,&x->map);
1938: return(0);
1939: }
1943: PetscErrorCode VecSetInf(Vec xin)
1944: {
1945: PetscInt i,n = xin->map->n;
1946: PetscScalar *xx;
1947: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1951: VecGetArray(xin,&xx);
1952: for (i=0; i<n; i++) xx[i] = inf;
1953: VecRestoreArray(xin,&xx);
1954: return(0);
1955: }