Actual source code: vector.c
1: /*
2: Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
3: These are the vector functions the user calls.
4: */
5: #include <petsc/private/vecimpl.h>
7: /* Logging support */
8: PetscClassId VEC_CLASSID;
9: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_Dot, VEC_MDot, VEC_TDot;
10: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
11: PetscLogEvent VEC_MTDot, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
12: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load;
13: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication,VEC_ReduceBegin,VEC_ReduceEnd,VEC_Ops;
14: PetscLogEvent VEC_DotNorm2, VEC_AXPBYPCZ;
15: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
16: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
17: PetscLogEvent VEC_CUDACopyFromGPUSome, VEC_CUDACopyToGPUSome;
18: PetscLogEvent VEC_HIPCopyFromGPU, VEC_HIPCopyToGPU;
19: PetscLogEvent VEC_HIPCopyFromGPUSome, VEC_HIPCopyToGPUSome;
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: {
46: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
47: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
48: return(0);
49: }
51: /*@
52: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
53: by the routine VecSetValuesLocal() to allow users to insert vector entries
54: using a local (per-processor) numbering.
56: Logically Collective on Vec
58: Input Parameters:
59: + x - vector
60: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
62: Notes:
63: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
65: Level: intermediate
67: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
68: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
69: @*/
70: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
71: {
78: if (x->ops->setlocaltoglobalmapping) {
79: (*x->ops->setlocaltoglobalmapping)(x,mapping);
80: } else {
81: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
82: }
83: return(0);
84: }
86: /*@
87: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
89: Not Collective
91: Input Parameter:
92: . X - the vector
94: Output Parameter:
95: . mapping - the mapping
97: Level: advanced
100: .seealso: VecSetValuesLocal()
101: @*/
102: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
103: {
108: *mapping = X->map->mapping;
109: return(0);
110: }
112: /*@
113: VecAssemblyBegin - Begins assembling the vector. This routine should
114: be called after completing all calls to VecSetValues().
116: Collective on Vec
118: Input Parameter:
119: . vec - the vector
121: Level: beginner
123: .seealso: VecAssemblyEnd(), VecSetValues()
124: @*/
125: PetscErrorCode VecAssemblyBegin(Vec vec)
126: {
132: VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
133: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
134: if (vec->ops->assemblybegin) {
135: (*vec->ops->assemblybegin)(vec);
136: }
137: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
138: PetscObjectStateIncrease((PetscObject)vec);
139: return(0);
140: }
142: /*@
143: VecAssemblyEnd - Completes assembling the vector. This routine should
144: be called after VecAssemblyBegin().
146: Collective on Vec
148: Input Parameter:
149: . vec - the vector
151: Options Database Keys:
152: + -vec_view - Prints vector in ASCII format
153: . -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
154: . -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
155: . -vec_view draw - Activates vector viewing using drawing tools
156: . -display <name> - Sets display name (default is host)
157: . -draw_pause <sec> - Sets number of seconds to pause after display
158: - -vec_view socket - Activates vector viewing using a socket
160: Level: beginner
162: .seealso: VecAssemblyBegin(), VecSetValues()
163: @*/
164: PetscErrorCode VecAssemblyEnd(Vec vec)
165: {
170: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
172: if (vec->ops->assemblyend) {
173: (*vec->ops->assemblyend)(vec);
174: }
175: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
176: VecViewFromOptions(vec,NULL,"-vec_view");
177: return(0);
178: }
180: /*@
181: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
183: Logically Collective on Vec
185: Input Parameters:
186: . x, y - the vectors
188: Output Parameter:
189: . w - the result
191: Level: advanced
193: Notes:
194: any subset of the x, y, and w may be the same vector.
195: For complex numbers compares only the real part
197: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
198: @*/
199: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
200: {
212: VecCheckSameSize(w,1,x,2);
213: VecCheckSameSize(w,1,y,3);
214: VecSetErrorIfLocked(w,1);
215: (*w->ops->pointwisemax)(w,x,y);
216: PetscObjectStateIncrease((PetscObject)w);
217: return(0);
218: }
220: /*@
221: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
223: Logically Collective on Vec
225: Input Parameters:
226: . x, y - the vectors
228: Output Parameter:
229: . w - the result
231: Level: advanced
233: Notes:
234: any subset of the x, y, and w may be the same vector.
235: For complex numbers compares only the real part
237: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
238: @*/
239: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
240: {
252: VecCheckSameSize(w,1,x,2);
253: VecCheckSameSize(w,1,y,3);
254: VecSetErrorIfLocked(w,1);
255: (*w->ops->pointwisemin)(w,x,y);
256: PetscObjectStateIncrease((PetscObject)w);
257: return(0);
258: }
260: /*@
261: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
263: Logically Collective on Vec
265: Input Parameters:
266: . x, y - the vectors
268: Output Parameter:
269: . w - the result
271: Level: advanced
273: Notes:
274: any subset of the x, y, and w may be the same vector.
276: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
277: @*/
278: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
279: {
291: VecCheckSameSize(w,1,x,2);
292: VecCheckSameSize(w,1,y,3);
293: VecSetErrorIfLocked(w,1);
294: (*w->ops->pointwisemaxabs)(w,x,y);
295: PetscObjectStateIncrease((PetscObject)w);
296: return(0);
297: }
299: /*@
300: VecPointwiseDivide - Computes the componentwise division w = x/y.
302: Logically Collective on Vec
304: Input Parameters:
305: . x, y - the vectors
307: Output Parameter:
308: . w - the result
310: Level: advanced
312: Notes:
313: any subset of the x, y, and w may be the same vector.
315: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
316: @*/
317: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
318: {
330: VecCheckSameSize(w,1,x,2);
331: VecCheckSameSize(w,1,y,3);
332: VecSetErrorIfLocked(w,1);
333: (*w->ops->pointwisedivide)(w,x,y);
334: PetscObjectStateIncrease((PetscObject)w);
335: return(0);
336: }
339: /*@
340: VecDuplicate - Creates a new vector of the same type as an existing vector.
342: Collective on Vec
344: Input Parameters:
345: . v - a vector to mimic
347: Output Parameter:
348: . newv - location to put new vector
350: Notes:
351: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
352: for the new vector. Use VecCopy() to copy a vector.
354: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
355: vectors.
357: Level: beginner
359: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
360: @*/
361: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
362: {
369: (*v->ops->duplicate)(v,newv);
370: PetscObjectStateIncrease((PetscObject)*newv);
371: return(0);
372: }
374: /*@C
375: VecDestroy - Destroys a vector.
377: Collective on Vec
379: Input Parameters:
380: . v - the vector
382: Level: beginner
384: .seealso: VecDuplicate(), VecDestroyVecs()
385: @*/
386: PetscErrorCode VecDestroy(Vec *v)
387: {
391: if (!*v) return(0);
393: if (--((PetscObject)(*v))->refct > 0) {*v = NULL; return(0);}
395: PetscObjectSAWsViewOff((PetscObject)*v);
396: /* destroy the internal part */
397: if ((*v)->ops->destroy) {
398: (*(*v)->ops->destroy)(*v);
399: }
400: PetscFree((*v)->defaultrandtype);
401: /* destroy the external/common part */
402: PetscLayoutDestroy(&(*v)->map);
403: PetscHeaderDestroy(v);
404: return(0);
405: }
407: /*@C
408: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
410: Collective on Vec
412: Input Parameters:
413: + m - the number of vectors to obtain
414: - v - a vector to mimic
416: Output Parameter:
417: . V - location to put pointer to array of vectors
419: Notes:
420: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
421: vector.
423: Fortran Note:
424: The Fortran interface is slightly different from that given below, it
425: requires one to pass in V a Vec (integer) array of size at least m.
426: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
428: Level: intermediate
430: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
431: @*/
432: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
433: {
440: (*v->ops->duplicatevecs)(v,m,V);
441: return(0);
442: }
444: /*@C
445: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
447: Collective on Vec
449: Input Parameters:
450: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
451: - m - the number of vectors previously obtained, if zero no vectors are destroyed
453: Fortran Note:
454: The Fortran interface is slightly different from that given below.
455: See the Fortran chapter of the users manual
457: Level: intermediate
459: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
460: @*/
461: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
462: {
467: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
468: if (!m || !*vv) {*vv = NULL; return(0);}
471: (*(**vv)->ops->destroyvecs)(m,*vv);
472: *vv = NULL;
473: return(0);
474: }
476: /*@C
477: VecViewFromOptions - View from Options
479: Collective on Vec
481: Input Parameters:
482: + A - the vector
483: . obj - Optional object
484: - name - command line option
486: Level: intermediate
487: .seealso: Vec, VecView, PetscObjectViewFromOptions(), VecCreate()
488: @*/
489: PetscErrorCode VecViewFromOptions(Vec A,PetscObject obj,const char name[])
490: {
495: PetscObjectViewFromOptions((PetscObject)A,obj,name);
496: return(0);
497: }
499: /*@C
500: VecView - Views a vector object.
502: Collective on Vec
504: Input Parameters:
505: + vec - the vector
506: - viewer - an optional visualization context
508: Notes:
509: The available visualization contexts include
510: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
511: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
512: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
514: You can change the format the vector is printed using the
515: option PetscViewerPushFormat().
517: The user can open alternative viewers with
518: + PetscViewerASCIIOpen() - Outputs vector to a specified file
519: . PetscViewerBinaryOpen() - Outputs vector in binary to a
520: specified file; corresponding input uses VecLoad()
521: . PetscViewerDrawOpen() - Outputs vector to an X window display
522: . PetscViewerSocketOpen() - Outputs vector to Socket viewer
523: - PetscViewerHDF5Open() - Outputs vector to HDF5 file viewer
525: The user can call PetscViewerPushFormat() to specify the output
526: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
527: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
528: + PETSC_VIEWER_DEFAULT - default, prints vector contents
529: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
530: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
531: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
532: format common among all vector types
534: Notes:
535: You can pass any number of vector objects, or other PETSc objects to the same viewer.
537: In the debugger you can do "call VecView(v,0)" to display the vector. (The same holds for any PETSc object viewer).
539: Notes for binary viewer:
540: If you pass multiple vectors to a binary viewer you can read them back in in the same order
541: with VecLoad().
543: If the blocksize of the vector is greater than one then you must provide a unique prefix to
544: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
545: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
546: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
547: filename. If you copy the binary file, make sure you copy the associated .info file with it.
549: See the manual page for VecLoad() on the exact format the binary viewer stores
550: the values in the file.
553: Notes for HDF5 Viewer:
554: The name of the Vec (given with PetscObjectSetName() is the name that is used
555: for the object in the HDF5 file. If you wish to store the same Vec into multiple
556: datasets in the same file (typically with different values), you must change its
557: name each time before calling the VecView(). To load the same vector,
558: the name of the Vec object passed to VecLoad() must be the same.
560: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
561: If the function PetscViewerHDF5SetBaseDimension2()is called then even if the block size is one it will
562: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
563: See also PetscViewerHDF5SetTimestep() which adds an additional complication to reading and writing Vecs
564: with the HDF5 viewer.
566: Level: beginner
569: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
570: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
571: PetscRealView(), PetscScalarView(), PetscIntView(), PetscViewerHDF5SetTimestep()
572: @*/
573: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
574: {
575: PetscErrorCode ierr;
576: PetscBool iascii;
577: PetscViewerFormat format;
578: PetscMPIInt size;
583: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);}
585: PetscViewerGetFormat(viewer,&format);
586: MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
587: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
589: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
591: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
592: if (iascii) {
593: PetscInt rows,bs;
595: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
596: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
597: PetscViewerASCIIPushTab(viewer);
598: VecGetSize(vec,&rows);
599: VecGetBlockSize(vec,&bs);
600: if (bs != 1) {
601: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
602: } else {
603: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
604: }
605: PetscViewerASCIIPopTab(viewer);
606: }
607: }
608: VecLockReadPush(vec);
609: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
610: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
611: (*vec->ops->viewnative)(vec,viewer);
612: } else {
613: (*vec->ops->view)(vec,viewer);
614: }
615: VecLockReadPop(vec);
616: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
617: return(0);
618: }
620: #if defined(PETSC_USE_DEBUG)
621: #include <../src/sys/totalview/tv_data_display.h>
622: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
623: {
624: const PetscScalar *values;
625: char type[32];
626: PetscErrorCode ierr;
629: TV_add_row("Local rows", "int", &v->map->n);
630: TV_add_row("Global rows", "int", &v->map->N);
631: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
632: VecGetArrayRead((Vec)v,&values);
633: PetscSNPrintf(type,32,"double[%d]",v->map->n);
634: TV_add_row("values",type, values);
635: VecRestoreArrayRead((Vec)v,&values);
636: return TV_format_OK;
637: }
638: #endif
640: /*@
641: VecGetSize - Returns the global number of elements of the vector.
643: Not Collective
645: Input Parameter:
646: . x - the vector
648: Output Parameters:
649: . size - the global length of the vector
651: Level: beginner
653: .seealso: VecGetLocalSize()
654: @*/
655: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
656: {
663: (*x->ops->getsize)(x,size);
664: return(0);
665: }
667: /*@
668: VecGetLocalSize - Returns the number of elements of the vector stored
669: in local memory.
671: Not Collective
673: Input Parameter:
674: . x - the vector
676: Output Parameter:
677: . size - the length of the local piece of the vector
679: Level: beginner
681: .seealso: VecGetSize()
682: @*/
683: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
684: {
691: (*x->ops->getlocalsize)(x,size);
692: return(0);
693: }
695: /*@C
696: VecGetOwnershipRange - Returns the range of indices owned by
697: this processor, assuming that the vectors are laid out with the
698: first n1 elements on the first processor, next n2 elements on the
699: second, etc. For certain parallel layouts this range may not be
700: well defined.
702: Not Collective
704: Input Parameter:
705: . x - the vector
707: Output Parameters:
708: + low - the first local element, pass in NULL if not interested
709: - high - one more than the last local element, pass in NULL if not interested
711: Note:
712: The high argument is one more than the last element stored locally.
714: Fortran: PETSC_NULL_INTEGER should be used instead of NULL
716: Level: beginner
718: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
719: @*/
720: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
721: {
727: if (low) *low = x->map->rstart;
728: if (high) *high = x->map->rend;
729: return(0);
730: }
732: /*@C
733: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
734: assuming that the vectors are laid out with the
735: first n1 elements on the first processor, next n2 elements on the
736: second, etc. For certain parallel layouts this range may not be
737: well defined.
739: Not Collective
741: Input Parameter:
742: . x - the vector
744: Output Parameters:
745: . range - array of length size+1 with the start and end+1 for each process
747: Note:
748: The high argument is one more than the last element stored locally.
750: Fortran: You must PASS in an array of length size+1
752: If the ranges are used after all vectors that share the ranges has been destroyed then the program will crash accessing ranges[].
754: Level: beginner
756: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
757: @*/
758: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
759: {
765: PetscLayoutGetRanges(x->map,ranges);
766: return(0);
767: }
769: /*@
770: VecSetOption - Sets an option for controling a vector's behavior.
772: Collective on Vec
774: Input Parameter:
775: + x - the vector
776: . op - the option
777: - flag - turn the option on or off
779: Supported Options:
780: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
781: entries destined to be stored on a separate processor. This can be used
782: to eliminate the global reduction in the VecAssemblyXXXX() if you know
783: that you have only used VecSetValues() to set local elements
784: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
785: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
786: ignored.
787: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
788: entries will always be a subset (possibly equal) of the off-process entries set on the
789: first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
790: changed this flag afterwards. If this assembly is not such first assembly, then this
791: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
792: a global reduction. Subsequent assemblies setting off-process values should use the same
793: InsertMode as the first assembly.
795: Developer Note:
796: The InsertMode restriction could be removed by packing the stash messages out of place.
798: Level: intermediate
800: @*/
801: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
802: {
808: if (x->ops->setoption) {
809: (*x->ops->setoption)(x,op,flag);
810: }
811: return(0);
812: }
814: /* Default routines for obtaining and releasing; */
815: /* may be used by any implementation */
816: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
817: {
819: PetscInt i;
824: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
825: PetscMalloc1(m,V);
826: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
827: return(0);
828: }
830: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
831: {
833: PetscInt i;
837: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
838: PetscFree(v);
839: return(0);
840: }
842: /*@
843: VecResetArray - Resets a vector to use its default memory. Call this
844: after the use of VecPlaceArray().
846: Not Collective
848: Input Parameters:
849: . vec - the vector
851: Level: developer
853: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
855: @*/
856: PetscErrorCode VecResetArray(Vec vec)
857: {
863: if (vec->ops->resetarray) {
864: (*vec->ops->resetarray)(vec);
865: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
866: PetscObjectStateIncrease((PetscObject)vec);
867: return(0);
868: }
870: /*@C
871: VecLoad - Loads a vector that has been stored in binary or HDF5 format
872: with VecView().
874: Collective on PetscViewer
876: Input Parameters:
877: + vec - the newly loaded vector, this needs to have been created with VecCreate() or
878: some related function before a call to VecLoad().
879: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
880: HDF5 file viewer, obtained from PetscViewerHDF5Open()
882: Level: intermediate
884: Notes:
885: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
886: before calling this.
888: The input file must contain the full global vector, as
889: written by the routine VecView().
891: If the type or size of vec is not set before a call to VecLoad, PETSc
892: sets the type and the local and global sizes. If type and/or
893: sizes are already set, then the same are used.
895: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
896: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
897: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
898: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
899: filename. If you copy the binary file, make sure you copy the associated .info file with it.
901: If using HDF5, you must assign the Vec the same name as was used in the Vec
902: that was stored in the file using PetscObjectSetName(). Otherwise you will
903: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT".
905: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
906: in loading the vector. Hence, for example, using Matlab notation h5create('vector.dat','/Test_Vec',[27 1]);
907: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
908: vectors communicator and the rest of the processes having zero entries
910: Notes for advanced users when using the binary viewer:
911: Most users should not need to know the details of the binary storage
912: format, since VecLoad() and VecView() completely hide these details.
913: But for anyone who's interested, the standard binary vector storage
914: format is
915: .vb
916: PetscInt VEC_FILE_CLASSID
917: PetscInt number of rows
918: PetscScalar *values of all entries
919: .ve
921: In addition, PETSc automatically uses byte swapping to work on all machines; the files
922: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
923: are converted to the small-endian format when they are read in from the file.
924: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
926: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
927: @*/
928: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
929: {
930: PetscErrorCode ierr;
931: PetscBool isbinary,ishdf5,isadios,isadios2;
932: PetscViewerFormat format;
938: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
939: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
940: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
941: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
942: if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
944: VecSetErrorIfLocked(vec,1);
945: if (!((PetscObject)vec)->type_name && !vec->ops->create) {
946: VecSetType(vec, VECSTANDARD);
947: }
948: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
949: PetscViewerGetFormat(viewer,&format);
950: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
951: (*vec->ops->loadnative)(vec,viewer);
952: } else {
953: (*vec->ops->load)(vec,viewer);
954: }
955: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
956: return(0);
957: }
960: /*@
961: VecReciprocal - Replaces each component of a vector by its reciprocal.
963: Logically Collective on Vec
965: Input Parameter:
966: . vec - the vector
968: Output Parameter:
969: . vec - the vector reciprocal
971: Level: intermediate
973: .seealso: VecLog(), VecExp(), VecSqrtAbs()
975: @*/
976: PetscErrorCode VecReciprocal(Vec vec)
977: {
983: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
984: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
985: VecSetErrorIfLocked(vec,1);
986: (*vec->ops->reciprocal)(vec);
987: PetscObjectStateIncrease((PetscObject)vec);
988: return(0);
989: }
991: /*@C
992: VecSetOperation - Allows user to set a vector operation.
994: Logically Collective on Vec
996: Input Parameters:
997: + vec - the vector
998: . op - the name of the operation
999: - f - the function that provides the operation.
1001: Level: advanced
1003: Usage:
1004: $ PetscErrorCode userview(Vec,PetscViewer);
1005: $ VecCreateMPI(comm,m,M,&x);
1006: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
1008: Notes:
1009: See the file include/petscvec.h for a complete list of matrix
1010: operations, which all have the form VECOP_<OPERATION>, where
1011: <OPERATION> is the name (in all capital letters) of the
1012: user interface routine (e.g., VecView() -> VECOP_VIEW).
1014: This function is not currently available from Fortran.
1016: .seealso: VecCreate(), MatShellSetOperation()
1017: @*/
1018: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1019: {
1022: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1023: vec->ops->viewnative = vec->ops->view;
1024: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1025: vec->ops->loadnative = vec->ops->load;
1026: }
1027: (((void(**)(void))vec->ops)[(int)op]) = f;
1028: return(0);
1029: }
1032: /*@
1033: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1034: used during the assembly process to store values that belong to
1035: other processors.
1037: Not Collective, different processes can have different size stashes
1039: Input Parameters:
1040: + vec - the vector
1041: . size - the initial size of the stash.
1042: - bsize - the initial size of the block-stash(if used).
1044: Options Database Keys:
1045: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1046: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1048: Level: intermediate
1050: Notes:
1051: The block-stash is used for values set with VecSetValuesBlocked() while
1052: the stash is used for values set with VecSetValues()
1054: Run with the option -info and look for output of the form
1055: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1056: to determine the appropriate value, MM, to use for size and
1057: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1058: to determine the value, BMM to use for bsize
1061: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1063: @*/
1064: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1065: {
1070: VecStashSetInitialSize_Private(&vec->stash,size);
1071: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1072: return(0);
1073: }
1075: /*@
1076: VecConjugate - Conjugates a vector.
1078: Logically Collective on Vec
1080: Input Parameters:
1081: . x - the vector
1083: Level: intermediate
1085: @*/
1086: PetscErrorCode VecConjugate(Vec x)
1087: {
1088: #if defined(PETSC_USE_COMPLEX)
1094: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1095: VecSetErrorIfLocked(x,1);
1096: (*x->ops->conjugate)(x);
1097: /* we need to copy norms here */
1098: PetscObjectStateIncrease((PetscObject)x);
1099: return(0);
1100: #else
1101: return(0);
1102: #endif
1103: }
1105: /*@
1106: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1108: Logically Collective on Vec
1110: Input Parameters:
1111: . x, y - the vectors
1113: Output Parameter:
1114: . w - the result
1116: Level: advanced
1118: Notes:
1119: any subset of the x, y, and w may be the same vector.
1121: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1122: @*/
1123: PetscErrorCode VecPointwiseMult(Vec w,Vec x,Vec y)
1124: {
1136: VecCheckSameSize(w,1,x,2);
1137: VecCheckSameSize(w,2,y,3);
1138: VecSetErrorIfLocked(w,1);
1139: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1140: (*w->ops->pointwisemult)(w,x,y);
1141: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1142: PetscObjectStateIncrease((PetscObject)w);
1143: return(0);
1144: }
1146: /*@
1147: VecSetRandom - Sets all components of a vector to random numbers.
1149: Logically Collective on Vec
1151: Input Parameters:
1152: + x - the vector
1153: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1154: it will create one internally.
1156: Output Parameter:
1157: . x - the vector
1159: Example of Usage:
1160: .vb
1161: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1162: VecSetRandom(x,rctx);
1163: PetscRandomDestroy(&rctx);
1164: .ve
1166: Level: intermediate
1168: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1169: @*/
1170: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1171: {
1173: PetscRandom randObj = NULL;
1179: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1180: VecSetErrorIfLocked(x,1);
1182: if (!rctx) {
1183: PetscRandomCreate(PetscObjectComm((PetscObject)x),&randObj);
1184: PetscRandomSetType(randObj,x->defaultrandtype);
1185: PetscRandomSetFromOptions(randObj);
1186: rctx = randObj;
1187: }
1189: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1190: (*x->ops->setrandom)(x,rctx);
1191: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1193: PetscRandomDestroy(&randObj);
1194: PetscObjectStateIncrease((PetscObject)x);
1195: return(0);
1196: }
1198: /*@
1199: VecZeroEntries - puts a 0.0 in each element of a vector
1201: Logically Collective on Vec
1203: Input Parameter:
1204: . vec - The vector
1206: Level: beginner
1208: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1209: @*/
1210: PetscErrorCode VecZeroEntries(Vec vec)
1211: {
1215: VecSet(vec,0);
1216: return(0);
1217: }
1219: /*
1220: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1221: processor and a PETSc MPI vector on more than one processor.
1223: Collective on Vec
1225: Input Parameter:
1226: . vec - The vector
1228: Level: intermediate
1230: .seealso: VecSetFromOptions(), VecSetType()
1231: */
1232: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1233: {
1234: PetscBool opt;
1235: VecType defaultType;
1236: char typeName[256];
1237: PetscMPIInt size;
1241: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1242: else {
1243: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1244: if (size > 1) defaultType = VECMPI;
1245: else defaultType = VECSEQ;
1246: }
1248: VecRegisterAll();
1249: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1250: if (opt) {
1251: VecSetType(vec, typeName);
1252: } else {
1253: VecSetType(vec, defaultType);
1254: }
1255: return(0);
1256: }
1258: /*@
1259: VecSetFromOptions - Configures the vector from the options database.
1261: Collective on Vec
1263: Input Parameter:
1264: . vec - The vector
1266: Notes:
1267: To see all options, run your program with the -help option, or consult the users manual.
1268: Must be called after VecCreate() but before the vector is used.
1270: Level: beginner
1273: .seealso: VecCreate(), VecSetOptionsPrefix()
1274: @*/
1275: PetscErrorCode VecSetFromOptions(Vec vec)
1276: {
1282: PetscObjectOptionsBegin((PetscObject)vec);
1283: /* Handle vector type options */
1284: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1286: /* Handle specific vector options */
1287: if (vec->ops->setfromoptions) {
1288: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1289: }
1291: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1292: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1293: PetscOptionsEnd();
1294: return(0);
1295: }
1297: /*@
1298: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1300: Collective on Vec
1302: Input Parameters:
1303: + v - the vector
1304: . n - the local size (or PETSC_DECIDE to have it set)
1305: - N - the global size (or PETSC_DECIDE)
1307: Notes:
1308: n and N cannot be both PETSC_DECIDE
1309: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1311: Level: intermediate
1313: .seealso: VecGetSize(), PetscSplitOwnership()
1314: @*/
1315: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1316: {
1322: 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);
1323: 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);
1324: v->map->n = n;
1325: v->map->N = N;
1326: if (v->ops->create) {
1327: (*v->ops->create)(v);
1328: v->ops->create = NULL;
1329: }
1330: return(0);
1331: }
1333: /*@
1334: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1335: and VecSetValuesBlockedLocal().
1337: Logically Collective on Vec
1339: Input Parameter:
1340: + v - the vector
1341: - bs - the blocksize
1343: Notes:
1344: All vectors obtained by VecDuplicate() inherit the same blocksize.
1346: Level: advanced
1348: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1350: @*/
1351: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1352: {
1358: PetscLayoutSetBlockSize(v->map,bs);
1359: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1360: return(0);
1361: }
1363: /*@
1364: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1365: and VecSetValuesBlockedLocal().
1367: Not Collective
1369: Input Parameter:
1370: . v - the vector
1372: Output Parameter:
1373: . bs - the blocksize
1375: Notes:
1376: All vectors obtained by VecDuplicate() inherit the same blocksize.
1378: Level: advanced
1380: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1383: @*/
1384: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1385: {
1391: PetscLayoutGetBlockSize(v->map,bs);
1392: return(0);
1393: }
1395: /*@C
1396: VecSetOptionsPrefix - Sets the prefix used for searching for all
1397: Vec options in the database.
1399: Logically Collective on Vec
1401: Input Parameter:
1402: + v - the Vec context
1403: - prefix - the prefix to prepend to all option names
1405: Notes:
1406: A hyphen (-) must NOT be given at the beginning of the prefix name.
1407: The first character of all runtime options is AUTOMATICALLY the hyphen.
1409: Level: advanced
1411: .seealso: VecSetFromOptions()
1412: @*/
1413: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1414: {
1419: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1420: return(0);
1421: }
1423: /*@C
1424: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1425: Vec options in the database.
1427: Logically Collective on Vec
1429: Input Parameters:
1430: + v - the Vec context
1431: - prefix - the prefix to prepend to all option names
1433: Notes:
1434: A hyphen (-) must NOT be given at the beginning of the prefix name.
1435: The first character of all runtime options is AUTOMATICALLY the hyphen.
1437: Level: advanced
1439: .seealso: VecGetOptionsPrefix()
1440: @*/
1441: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1442: {
1447: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1448: return(0);
1449: }
1451: /*@C
1452: VecGetOptionsPrefix - Sets the prefix used for searching for all
1453: Vec options in the database.
1455: Not Collective
1457: Input Parameter:
1458: . v - the Vec context
1460: Output Parameter:
1461: . prefix - pointer to the prefix string used
1463: Notes:
1464: On the fortran side, the user should pass in a string 'prefix' of
1465: sufficient length to hold the prefix.
1467: Level: advanced
1469: .seealso: VecAppendOptionsPrefix()
1470: @*/
1471: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1472: {
1477: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1478: return(0);
1479: }
1481: /*@
1482: VecSetUp - Sets up the internal vector data structures for the later use.
1484: Collective on Vec
1486: Input Parameters:
1487: . v - the Vec context
1489: Notes:
1490: For basic use of the Vec classes the user need not explicitly call
1491: VecSetUp(), since these actions will happen automatically.
1493: Level: advanced
1495: .seealso: VecCreate(), VecDestroy()
1496: @*/
1497: PetscErrorCode VecSetUp(Vec v)
1498: {
1499: PetscMPIInt size;
1504: if (v->map->n < 0 && v->map->N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1505: if (!((PetscObject)v)->type_name) {
1506: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1507: if (size == 1) {
1508: VecSetType(v, VECSEQ);
1509: } else {
1510: VecSetType(v, VECMPI);
1511: }
1512: }
1513: return(0);
1514: }
1516: /*
1517: These currently expose the PetscScalar/PetscReal in updating the
1518: cached norm. If we push those down into the implementation these
1519: will become independent of PetscScalar/PetscReal
1520: */
1522: /*@
1523: VecCopy - Copies a vector. y <- x
1525: Logically Collective on Vec
1527: Input Parameter:
1528: . x - the vector
1530: Output Parameter:
1531: . y - the copy
1533: Notes:
1534: For default parallel PETSc vectors, both x and y must be distributed in
1535: the same manner; local copies are done.
1537: Developer Notes:
1539: of the vectors to be sequential and one to be parallel so long as both have the same
1540: local sizes. This is used in some internal functions in PETSc.
1542: Level: beginner
1544: .seealso: VecDuplicate()
1545: @*/
1546: PetscErrorCode VecCopy(Vec x,Vec y)
1547: {
1548: PetscBool flgs[4];
1549: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1551: PetscInt i;
1558: if (x == y) return(0);
1559: VecCheckSameLocalSize(x,1,y,2);
1560: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1561: VecSetErrorIfLocked(y,2);
1563: #if !defined(PETSC_USE_MIXED_PRECISION)
1564: for (i=0; i<4; i++) {
1565: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1566: }
1567: #endif
1569: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1570: #if defined(PETSC_USE_MIXED_PRECISION)
1571: extern PetscErrorCode VecGetArray(Vec,double**);
1572: extern PetscErrorCode VecRestoreArray(Vec,double**);
1573: extern PetscErrorCode VecGetArray(Vec,float**);
1574: extern PetscErrorCode VecRestoreArray(Vec,float**);
1575: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1576: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1577: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1578: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1579: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1580: PetscInt i,n;
1581: const float *xx;
1582: double *yy;
1583: VecGetArrayRead(x,&xx);
1584: VecGetArray(y,&yy);
1585: VecGetLocalSize(x,&n);
1586: for (i=0; i<n; i++) yy[i] = xx[i];
1587: VecRestoreArrayRead(x,&xx);
1588: VecRestoreArray(y,&yy);
1589: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1590: PetscInt i,n;
1591: float *yy;
1592: const double *xx;
1593: VecGetArrayRead(x,&xx);
1594: VecGetArray(y,&yy);
1595: VecGetLocalSize(x,&n);
1596: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1597: VecRestoreArrayRead(x,&xx);
1598: VecRestoreArray(y,&yy);
1599: } else {
1600: (*x->ops->copy)(x,y);
1601: }
1602: #else
1603: (*x->ops->copy)(x,y);
1604: #endif
1606: PetscObjectStateIncrease((PetscObject)y);
1607: #if !defined(PETSC_USE_MIXED_PRECISION)
1608: for (i=0; i<4; i++) {
1609: if (flgs[i]) {
1610: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1611: }
1612: }
1613: #endif
1615: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1616: return(0);
1617: }
1619: /*@
1620: VecSwap - Swaps the vectors x and y.
1622: Logically Collective on Vec
1624: Input Parameters:
1625: . x, y - the vectors
1627: Level: advanced
1629: @*/
1630: PetscErrorCode VecSwap(Vec x,Vec y)
1631: {
1632: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1633: PetscBool flgxs[4],flgys[4];
1635: PetscInt i;
1643: VecCheckSameSize(x,1,y,2);
1644: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1645: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1646: VecSetErrorIfLocked(x,1);
1647: VecSetErrorIfLocked(y,2);
1649: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1650: for (i=0; i<4; i++) {
1651: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1652: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1653: }
1654: (*x->ops->swap)(x,y);
1655: PetscObjectStateIncrease((PetscObject)x);
1656: PetscObjectStateIncrease((PetscObject)y);
1657: for (i=0; i<4; i++) {
1658: if (flgxs[i]) {
1659: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1660: }
1661: if (flgys[i]) {
1662: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1663: }
1664: }
1665: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1666: return(0);
1667: }
1669: /*
1670: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1672: Collective on VecStash
1674: Input Parameters:
1675: + obj - the VecStash object
1676: . bobj - optional other object that provides the prefix
1677: - optionname - option to activate viewing
1679: Level: intermediate
1681: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1683: */
1684: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1685: {
1686: PetscErrorCode ierr;
1687: PetscViewer viewer;
1688: PetscBool flg;
1689: PetscViewerFormat format;
1690: char *prefix;
1693: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1694: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1695: if (flg) {
1696: PetscViewerPushFormat(viewer,format);
1697: VecStashView(obj,viewer);
1698: PetscViewerPopFormat(viewer);
1699: PetscViewerDestroy(&viewer);
1700: }
1701: return(0);
1702: }
1704: /*@
1705: VecStashView - Prints the entries in the vector stash and block stash.
1707: Collective on Vec
1709: Input Parameters:
1710: + v - the vector
1711: - viewer - the viewer
1713: Level: advanced
1716: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1718: @*/
1719: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1720: {
1722: PetscMPIInt rank;
1723: PetscInt i,j;
1724: PetscBool match;
1725: VecStash *s;
1726: PetscScalar val;
1733: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1734: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1735: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1736: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1737: s = &v->bstash;
1739: /* print block stash */
1740: PetscViewerASCIIPushSynchronized(viewer);
1741: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1742: for (i=0; i<s->n; i++) {
1743: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1744: for (j=0; j<s->bs; j++) {
1745: val = s->array[i*s->bs+j];
1746: #if defined(PETSC_USE_COMPLEX)
1747: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1748: #else
1749: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1750: #endif
1751: }
1752: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1753: }
1754: PetscViewerFlush(viewer);
1756: s = &v->stash;
1758: /* print basic stash */
1759: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1760: for (i=0; i<s->n; i++) {
1761: val = s->array[i];
1762: #if defined(PETSC_USE_COMPLEX)
1763: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1764: #else
1765: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1766: #endif
1767: }
1768: PetscViewerFlush(viewer);
1769: PetscViewerASCIIPopSynchronized(viewer);
1770: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1771: return(0);
1772: }
1774: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1775: {
1776: PetscInt i,N,rstart,rend;
1778: PetscScalar *xx;
1779: PetscReal *xreal;
1780: PetscBool iset;
1783: VecGetOwnershipRange(v,&rstart,&rend);
1784: VecGetSize(v,&N);
1785: PetscCalloc1(N,&xreal);
1786: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1787: if (iset) {
1788: VecGetArray(v,&xx);
1789: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1790: VecRestoreArray(v,&xx);
1791: }
1792: PetscFree(xreal);
1793: if (set) *set = iset;
1794: return(0);
1795: }
1797: /*@
1798: VecGetLayout - get PetscLayout describing vector layout
1800: Not Collective
1802: Input Arguments:
1803: . x - the vector
1805: Output Arguments:
1806: . map - the layout
1808: Level: developer
1810: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1811: @*/
1812: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1813: {
1817: *map = x->map;
1818: return(0);
1819: }
1821: /*@
1822: VecSetLayout - set PetscLayout describing vector layout
1824: Not Collective
1826: Input Arguments:
1827: + x - the vector
1828: - map - the layout
1830: Notes:
1831: It is normally only valid to replace the layout with a layout known to be equivalent.
1833: Level: developer
1835: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1836: @*/
1837: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1838: {
1843: PetscLayoutReference(map,&x->map);
1844: return(0);
1845: }
1847: PetscErrorCode VecSetInf(Vec xin)
1848: {
1849: PetscInt i,n = xin->map->n;
1850: PetscScalar *xx;
1851: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1855: if (xin->ops->set) { /* can be called by a subset of processes, do not use collective routines */
1856: (*xin->ops->set)(xin,inf);
1857: } else {
1858: VecGetArrayWrite(xin,&xx);
1859: for (i=0; i<n; i++) xx[i] = inf;
1860: VecRestoreArrayWrite(xin,&xx);
1861: }
1862: return(0);
1863: }
1865: /*@
1866: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
1868: Input Parameters:
1869: + v - the vector
1870: - flg - bind to the CPU if value of PETSC_TRUE
1872: Level: intermediate
1873: @*/
1874: PetscErrorCode VecBindToCPU(Vec v,PetscBool flg)
1875: {
1876: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1880: if (v->boundtocpu == flg) return(0);
1881: v->boundtocpu = flg;
1882: if (v->ops->bindtocpu) {
1883: (*v->ops->bindtocpu)(v,flg);
1884: }
1885: return(0);
1886: #else
1887: return 0;
1888: #endif
1889: }
1891: /*@C
1892: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
1894: Logically Collective on Vec
1896: Input Parameters:
1897: + v - the vector
1898: - mbytes - minimum data size in bytes
1900: Options Database Keys:
1902: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.
1903: Note that this takes a PetscScalar, to accommodate large values;
1904: specifying -1 ensures that pinned memory will never be used.
1906: Level: developer
1908: .seealso: VecGetPinnedMemoryMin()
1909: @*/
1910: PetscErrorCode VecSetPinnedMemoryMin(Vec v,size_t mbytes)
1911: {
1912: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1914: v->minimum_bytes_pinned_memory = mbytes;
1915: return(0);
1916: #else
1917: return 0;
1918: #endif
1919: }
1921: /*@C
1922: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
1924: Logically Collective on Vec
1926: Input Parameters:
1927: . v - the vector
1929: Output Parameters:
1930: . mbytes - minimum data size in bytes
1932: Level: developer
1934: .seealso: VecSetPinnedMemoryMin()
1935: @*/
1936: PetscErrorCode VecGetPinnedMemoryMin(Vec v,size_t *mbytes)
1937: {
1938: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA) || defined(PETSC_HAVE_HIP)
1940: *mbytes = v->minimum_bytes_pinned_memory;
1941: return(0);
1942: #else
1943: return 0;
1944: #endif
1945: }
1947: /*@
1948: VecGetOffloadMask - Get the offload mask of a Vec.
1950: Not Collective
1952: Input Parameters:
1953: . v - the vector
1955: Output Parameters:
1956: . mask - corresponding PetscOffloadMask enum value.
1958: Level: intermediate
1960: .seealso: VecCreateSeqCUDA(), VecCreateSeqViennaCL(), VecGetArray(), VecGetType()
1961: @*/
1962: PetscErrorCode VecGetOffloadMask(Vec v,PetscOffloadMask* mask)
1963: {
1965: *mask = v->offloadmask;
1966: return(0);
1967: }
1970: #if !defined(PETSC_HAVE_VIENNACL)
1971: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLContext(Vec v,PETSC_UINTPTR_T* ctx)
1972: {
1973: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to get a Vec's cl_context");
1974: }
1976: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLQueue(Vec v,PETSC_UINTPTR_T* queue)
1977: {
1978: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to get a Vec's cl_command_queue");
1979: }
1981: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMem(Vec v,PETSC_UINTPTR_T* queue)
1982: {
1983: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to get a Vec's cl_mem");
1984: }
1986: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemRead(Vec v,PETSC_UINTPTR_T* queue)
1987: {
1988: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to get a Vec's cl_mem");
1989: }
1991: PETSC_EXTERN PetscErrorCode VecViennaCLGetCLMemWrite(Vec v,PETSC_UINTPTR_T* queue)
1992: {
1993: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to get a Vec's cl_mem");
1994: }
1996: PETSC_EXTERN PetscErrorCode VecViennaCLRestoreCLMemWrite(Vec v)
1997: {
1998: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_LIB,"PETSc must be configured with --with-opencl to restore a Vec's cl_mem");
1999: }
2000: #endif