Actual source code: vector.c
petsc-3.14.3 2021-01-09
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;
19: /*@
20: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
21: to be communicated to other processors during the VecAssemblyBegin/End() process
23: Not collective
25: Input Parameter:
26: . vec - the vector
28: Output Parameters:
29: + nstash - the size of the stash
30: . reallocs - the number of additional mallocs incurred.
31: . bnstash - the size of the block stash
32: - breallocs - the number of additional mallocs incurred.in the block stash
34: Level: advanced
36: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
38: @*/
39: PetscErrorCode VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
40: {
44: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
45: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
46: return(0);
47: }
49: /*@
50: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
51: by the routine VecSetValuesLocal() to allow users to insert vector entries
52: using a local (per-processor) numbering.
54: Logically Collective on Vec
56: Input Parameters:
57: + x - vector
58: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
60: Notes:
61: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
63: Level: intermediate
65: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
66: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
67: @*/
68: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
69: {
76: if (x->ops->setlocaltoglobalmapping) {
77: (*x->ops->setlocaltoglobalmapping)(x,mapping);
78: } else {
79: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
80: }
81: return(0);
82: }
84: /*@
85: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
87: Not Collective
89: Input Parameter:
90: . X - the vector
92: Output Parameter:
93: . mapping - the mapping
95: Level: advanced
98: .seealso: VecSetValuesLocal()
99: @*/
100: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
101: {
106: *mapping = X->map->mapping;
107: return(0);
108: }
110: /*@
111: VecAssemblyBegin - Begins assembling the vector. This routine should
112: be called after completing all calls to VecSetValues().
114: Collective on Vec
116: Input Parameter:
117: . vec - the vector
119: Level: beginner
121: .seealso: VecAssemblyEnd(), VecSetValues()
122: @*/
123: PetscErrorCode VecAssemblyBegin(Vec vec)
124: {
130: VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
131: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
132: if (vec->ops->assemblybegin) {
133: (*vec->ops->assemblybegin)(vec);
134: }
135: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
136: PetscObjectStateIncrease((PetscObject)vec);
137: return(0);
138: }
140: /*@
141: VecAssemblyEnd - Completes assembling the vector. This routine should
142: be called after VecAssemblyBegin().
144: Collective on Vec
146: Input Parameter:
147: . vec - the vector
149: Options Database Keys:
150: + -vec_view - Prints vector in ASCII format
151: . -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
152: . -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
153: . -vec_view draw - Activates vector viewing using drawing tools
154: . -display <name> - Sets display name (default is host)
155: . -draw_pause <sec> - Sets number of seconds to pause after display
156: - -vec_view socket - Activates vector viewing using a socket
158: Level: beginner
160: .seealso: VecAssemblyBegin(), VecSetValues()
161: @*/
162: PetscErrorCode VecAssemblyEnd(Vec vec)
163: {
168: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
170: if (vec->ops->assemblyend) {
171: (*vec->ops->assemblyend)(vec);
172: }
173: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
174: VecViewFromOptions(vec,NULL,"-vec_view");
175: return(0);
176: }
178: /*@
179: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
181: Logically Collective on Vec
183: Input Parameters:
184: . x, y - the vectors
186: Output Parameter:
187: . w - the result
189: Level: advanced
191: Notes:
192: any subset of the x, y, and w may be the same vector.
193: For complex numbers compares only the real part
195: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
196: @*/
197: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
198: {
210: VecCheckSameSize(w,1,x,2);
211: VecCheckSameSize(w,1,y,3);
212: VecSetErrorIfLocked(w,1);
213: (*w->ops->pointwisemax)(w,x,y);
214: PetscObjectStateIncrease((PetscObject)w);
215: return(0);
216: }
218: /*@
219: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
221: Logically Collective on Vec
223: Input Parameters:
224: . x, y - the vectors
226: Output Parameter:
227: . w - the result
229: Level: advanced
231: Notes:
232: any subset of the x, y, and w may be the same vector.
233: For complex numbers compares only the real part
235: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
236: @*/
237: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
238: {
250: VecCheckSameSize(w,1,x,2);
251: VecCheckSameSize(w,1,y,3);
252: VecSetErrorIfLocked(w,1);
253: (*w->ops->pointwisemin)(w,x,y);
254: PetscObjectStateIncrease((PetscObject)w);
255: return(0);
256: }
258: /*@
259: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
261: Logically Collective on Vec
263: Input Parameters:
264: . x, y - the vectors
266: Output Parameter:
267: . w - the result
269: Level: advanced
271: Notes:
272: any subset of the x, y, and w may be the same vector.
274: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
275: @*/
276: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
277: {
289: VecCheckSameSize(w,1,x,2);
290: VecCheckSameSize(w,1,y,3);
291: VecSetErrorIfLocked(w,1);
292: (*w->ops->pointwisemaxabs)(w,x,y);
293: PetscObjectStateIncrease((PetscObject)w);
294: return(0);
295: }
297: /*@
298: VecPointwiseDivide - Computes the componentwise division w = x/y.
300: Logically Collective on Vec
302: Input Parameters:
303: . x, y - the vectors
305: Output Parameter:
306: . w - the result
308: Level: advanced
310: Notes:
311: any subset of the x, y, and w may be the same vector.
313: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
314: @*/
315: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
316: {
328: VecCheckSameSize(w,1,x,2);
329: VecCheckSameSize(w,1,y,3);
330: VecSetErrorIfLocked(w,1);
331: (*w->ops->pointwisedivide)(w,x,y);
332: PetscObjectStateIncrease((PetscObject)w);
333: return(0);
334: }
337: /*@
338: VecDuplicate - Creates a new vector of the same type as an existing vector.
340: Collective on Vec
342: Input Parameters:
343: . v - a vector to mimic
345: Output Parameter:
346: . newv - location to put new vector
348: Notes:
349: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
350: for the new vector. Use VecCopy() to copy a vector.
352: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
353: vectors.
355: Level: beginner
357: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
358: @*/
359: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
360: {
367: (*v->ops->duplicate)(v,newv);
368: PetscObjectStateIncrease((PetscObject)*newv);
369: return(0);
370: }
372: /*@
373: VecDestroy - Destroys a vector.
375: Collective on Vec
377: Input Parameters:
378: . v - the vector
380: Level: beginner
382: .seealso: VecDuplicate(), VecDestroyVecs()
383: @*/
384: PetscErrorCode VecDestroy(Vec *v)
385: {
389: if (!*v) return(0);
391: if (--((PetscObject)(*v))->refct > 0) {*v = NULL; return(0);}
393: PetscObjectSAWsViewOff((PetscObject)*v);
394: /* destroy the internal part */
395: if ((*v)->ops->destroy) {
396: (*(*v)->ops->destroy)(*v);
397: }
398: /* destroy the external/common part */
399: PetscLayoutDestroy(&(*v)->map);
400: PetscHeaderDestroy(v);
401: return(0);
402: }
404: /*@C
405: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
407: Collective on Vec
409: Input Parameters:
410: + m - the number of vectors to obtain
411: - v - a vector to mimic
413: Output Parameter:
414: . V - location to put pointer to array of vectors
416: Notes:
417: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
418: vector.
420: Fortran Note:
421: The Fortran interface is slightly different from that given below, it
422: requires one to pass in V a Vec (integer) array of size at least m.
423: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
425: Level: intermediate
427: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
428: @*/
429: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
430: {
437: (*v->ops->duplicatevecs)(v,m,V);
438: return(0);
439: }
441: /*@C
442: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
444: Collective on Vec
446: Input Parameters:
447: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
448: - m - the number of vectors previously obtained, if zero no vectors are destroyed
450: Fortran Note:
451: The Fortran interface is slightly different from that given below.
452: See the Fortran chapter of the users manual
454: Level: intermediate
456: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
457: @*/
458: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
459: {
464: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
465: if (!m || !*vv) {*vv = NULL; return(0);}
468: (*(**vv)->ops->destroyvecs)(m,*vv);
469: *vv = NULL;
470: return(0);
471: }
473: /*@C
474: VecViewFromOptions - View from Options
476: Collective on Vec
478: Input Parameters:
479: + A - the vector
480: . obj - Optional object
481: - name - command line option
483: Level: intermediate
484: .seealso: Vec, VecView, PetscObjectViewFromOptions(), VecCreate()
485: @*/
486: PetscErrorCode VecViewFromOptions(Vec A,PetscObject obj,const char name[])
487: {
492: PetscObjectViewFromOptions((PetscObject)A,obj,name);
493: return(0);
494: }
496: /*@C
497: VecView - Views a vector object.
499: Collective on Vec
501: Input Parameters:
502: + vec - the vector
503: - viewer - an optional visualization context
505: Notes:
506: The available visualization contexts include
507: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
508: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
509: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
511: You can change the format the vector is printed using the
512: option PetscViewerPushFormat().
514: The user can open alternative viewers with
515: + PetscViewerASCIIOpen() - Outputs vector to a specified file
516: . PetscViewerBinaryOpen() - Outputs vector in binary to a
517: specified file; corresponding input uses VecLoad()
518: . PetscViewerDrawOpen() - Outputs vector to an X window display
519: . PetscViewerSocketOpen() - Outputs vector to Socket viewer
520: - PetscViewerHDF5Open() - Outputs vector to HDF5 file viewer
522: The user can call PetscViewerPushFormat() to specify the output
523: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
524: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
525: + PETSC_VIEWER_DEFAULT - default, prints vector contents
526: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
527: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
528: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
529: format common among all vector types
531: Notes:
532: You can pass any number of vector objects, or other PETSc objects to the same viewer.
534: In the debugger you can do "call VecView(v,0)" to display the vector. (The same holds for any PETSc object viewer).
536: Notes for binary viewer:
537: If you pass multiple vectors to a binary viewer you can read them back in in the same order
538: with VecLoad().
540: If the blocksize of the vector is greater than one then you must provide a unique prefix to
541: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
542: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
543: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
544: filename. If you copy the binary file, make sure you copy the associated .info file with it.
546: See the manual page for VecLoad() on the exact format the binary viewer stores
547: the values in the file.
550: Notes for HDF5 Viewer:
551: The name of the Vec (given with PetscObjectSetName() is the name that is used
552: for the object in the HDF5 file. If you wish to store the same Vec into multiple
553: datasets in the same file (typically with different values), you must change its
554: name each time before calling the VecView(). To load the same vector,
555: the name of the Vec object passed to VecLoad() must be the same.
557: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
558: If the function PetscViewerHDF5SetBaseDimension2()is called then even if the block size is one it will
559: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
560: See also PetscViewerHDF5SetTimestep() which adds an additional complication to reading and writing Vecs
561: with the HDF5 viewer.
563: Level: beginner
566: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
567: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
568: PetscRealView(), PetscScalarView(), PetscIntView(), PetscViewerHDF5SetTimestep()
569: @*/
570: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
571: {
572: PetscErrorCode ierr;
573: PetscBool iascii;
574: PetscViewerFormat format;
575: PetscMPIInt size;
580: if (!viewer) {PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);}
582: PetscViewerGetFormat(viewer,&format);
583: MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
584: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
586: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
588: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
589: if (iascii) {
590: PetscInt rows,bs;
592: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
593: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
594: PetscViewerASCIIPushTab(viewer);
595: VecGetSize(vec,&rows);
596: VecGetBlockSize(vec,&bs);
597: if (bs != 1) {
598: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
599: } else {
600: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
601: }
602: PetscViewerASCIIPopTab(viewer);
603: }
604: }
605: VecLockReadPush(vec);
606: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
607: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
608: (*vec->ops->viewnative)(vec,viewer);
609: } else {
610: (*vec->ops->view)(vec,viewer);
611: }
612: VecLockReadPop(vec);
613: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
614: return(0);
615: }
617: #if defined(PETSC_USE_DEBUG)
618: #include <../src/sys/totalview/tv_data_display.h>
619: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
620: {
621: const PetscScalar *values;
622: char type[32];
623: PetscErrorCode ierr;
626: TV_add_row("Local rows", "int", &v->map->n);
627: TV_add_row("Global rows", "int", &v->map->N);
628: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
629: VecGetArrayRead((Vec)v,&values);
630: PetscSNPrintf(type,32,"double[%d]",v->map->n);
631: TV_add_row("values",type, values);
632: VecRestoreArrayRead((Vec)v,&values);
633: return TV_format_OK;
634: }
635: #endif
637: /*@
638: VecGetSize - Returns the global number of elements of the vector.
640: Not Collective
642: Input Parameter:
643: . x - the vector
645: Output Parameters:
646: . size - the global length of the vector
648: Level: beginner
650: .seealso: VecGetLocalSize()
651: @*/
652: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
653: {
660: (*x->ops->getsize)(x,size);
661: return(0);
662: }
664: /*@
665: VecGetLocalSize - Returns the number of elements of the vector stored
666: in local memory.
668: Not Collective
670: Input Parameter:
671: . x - the vector
673: Output Parameter:
674: . size - the length of the local piece of the vector
676: Level: beginner
678: .seealso: VecGetSize()
679: @*/
680: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
681: {
688: (*x->ops->getlocalsize)(x,size);
689: return(0);
690: }
692: /*@C
693: VecGetOwnershipRange - Returns the range of indices owned by
694: this processor, assuming that the vectors are laid out with the
695: first n1 elements on the first processor, next n2 elements on the
696: second, etc. For certain parallel layouts this range may not be
697: well defined.
699: Not Collective
701: Input Parameter:
702: . x - the vector
704: Output Parameters:
705: + low - the first local element, pass in NULL if not interested
706: - high - one more than the last local element, pass in NULL if not interested
708: Note:
709: The high argument is one more than the last element stored locally.
711: Fortran: PETSC_NULL_INTEGER should be used instead of NULL
713: Level: beginner
716: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
717: @*/
718: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
719: {
725: if (low) *low = x->map->rstart;
726: if (high) *high = x->map->rend;
727: return(0);
728: }
730: /*@C
731: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
732: assuming that the vectors are laid out with the
733: first n1 elements on the first processor, next n2 elements on the
734: second, etc. For certain parallel layouts this range may not be
735: well defined.
737: Not Collective
739: Input Parameter:
740: . x - the vector
742: Output Parameters:
743: . range - array of length size+1 with the start and end+1 for each process
745: Note:
746: The high argument is one more than the last element stored locally.
748: Fortran: You must PASS in an array of length size+1
750: Level: beginner
753: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
754: @*/
755: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
756: {
762: PetscLayoutGetRanges(x->map,ranges);
763: return(0);
764: }
766: /*@
767: VecSetOption - Sets an option for controling a vector's behavior.
769: Collective on Vec
771: Input Parameter:
772: + x - the vector
773: . op - the option
774: - flag - turn the option on or off
776: Supported Options:
777: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
778: entries destined to be stored on a separate processor. This can be used
779: to eliminate the global reduction in the VecAssemblyXXXX() if you know
780: that you have only used VecSetValues() to set local elements
781: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
782: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
783: ignored.
784: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
785: entries will always be a subset (possibly equal) of the off-process entries set on the
786: first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
787: changed this flag afterwards. If this assembly is not such first assembly, then this
788: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
789: a global reduction. Subsequent assemblies setting off-process values should use the same
790: InsertMode as the first assembly.
792: Developer Note:
793: The InsertMode restriction could be removed by packing the stash messages out of place.
795: Level: intermediate
797: @*/
798: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
799: {
805: if (x->ops->setoption) {
806: (*x->ops->setoption)(x,op,flag);
807: }
808: return(0);
809: }
811: /* Default routines for obtaining and releasing; */
812: /* may be used by any implementation */
813: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
814: {
816: PetscInt i;
821: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
822: PetscMalloc1(m,V);
823: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
824: return(0);
825: }
827: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
828: {
830: PetscInt i;
834: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
835: PetscFree(v);
836: return(0);
837: }
839: /*@
840: VecResetArray - Resets a vector to use its default memory. Call this
841: after the use of VecPlaceArray().
843: Not Collective
845: Input Parameters:
846: . vec - the vector
848: Level: developer
850: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
852: @*/
853: PetscErrorCode VecResetArray(Vec vec)
854: {
860: if (vec->ops->resetarray) {
861: (*vec->ops->resetarray)(vec);
862: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
863: PetscObjectStateIncrease((PetscObject)vec);
864: return(0);
865: }
867: /*@C
868: VecLoad - Loads a vector that has been stored in binary or HDF5 format
869: with VecView().
871: Collective on PetscViewer
873: Input Parameters:
874: + vec - the newly loaded vector, this needs to have been created with VecCreate() or
875: some related function before a call to VecLoad().
876: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
877: HDF5 file viewer, obtained from PetscViewerHDF5Open()
879: Level: intermediate
881: Notes:
882: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
883: before calling this.
885: The input file must contain the full global vector, as
886: written by the routine VecView().
888: If the type or size of vec is not set before a call to VecLoad, PETSc
889: sets the type and the local and global sizes. If type and/or
890: sizes are already set, then the same are used.
892: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
893: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
894: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
895: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
896: filename. If you copy the binary file, make sure you copy the associated .info file with it.
898: If using HDF5, you must assign the Vec the same name as was used in the Vec
899: that was stored in the file using PetscObjectSetName(). Otherwise you will
900: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT".
902: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
903: in loading the vector. Hence, for example, using Matlab notation h5create('vector.dat','/Test_Vec',[27 1]);
904: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
905: vectors communicator and the rest of the processes having zero entries
907: Notes for advanced users when using the binary viewer:
908: Most users should not need to know the details of the binary storage
909: format, since VecLoad() and VecView() completely hide these details.
910: But for anyone who's interested, the standard binary vector storage
911: format is
912: .vb
913: PetscInt VEC_FILE_CLASSID
914: PetscInt number of rows
915: PetscScalar *values of all entries
916: .ve
918: In addition, PETSc automatically uses byte swapping to work on all machines; the files
919: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
920: are converted to the small-endian format when they are read in from the file.
921: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
923: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
924: @*/
925: PetscErrorCode VecLoad(Vec vec, PetscViewer viewer)
926: {
927: PetscErrorCode ierr;
928: PetscBool isbinary,ishdf5,isadios,isadios2;
929: PetscViewerFormat format;
935: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
936: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
937: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
938: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
939: if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
941: VecSetErrorIfLocked(vec,1);
942: if (!((PetscObject)vec)->type_name && !vec->ops->create) {
943: VecSetType(vec, VECSTANDARD);
944: }
945: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
946: PetscViewerGetFormat(viewer,&format);
947: if (format == PETSC_VIEWER_NATIVE && vec->ops->loadnative) {
948: (*vec->ops->loadnative)(vec,viewer);
949: } else {
950: (*vec->ops->load)(vec,viewer);
951: }
952: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
953: return(0);
954: }
957: /*@
958: VecReciprocal - Replaces each component of a vector by its reciprocal.
960: Logically Collective on Vec
962: Input Parameter:
963: . vec - the vector
965: Output Parameter:
966: . vec - the vector reciprocal
968: Level: intermediate
970: .seealso: VecLog(), VecExp(), VecSqrtAbs()
972: @*/
973: PetscErrorCode VecReciprocal(Vec vec)
974: {
980: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
981: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
982: VecSetErrorIfLocked(vec,1);
983: (*vec->ops->reciprocal)(vec);
984: PetscObjectStateIncrease((PetscObject)vec);
985: return(0);
986: }
988: /*@C
989: VecSetOperation - Allows user to set a vector operation.
991: Logically Collective on Vec
993: Input Parameters:
994: + vec - the vector
995: . op - the name of the operation
996: - f - the function that provides the operation.
998: Level: advanced
1000: Usage:
1001: $ PetscErrorCode userview(Vec,PetscViewer);
1002: $ VecCreateMPI(comm,m,M,&x);
1003: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
1005: Notes:
1006: See the file include/petscvec.h for a complete list of matrix
1007: operations, which all have the form VECOP_<OPERATION>, where
1008: <OPERATION> is the name (in all capital letters) of the
1009: user interface routine (e.g., VecView() -> VECOP_VIEW).
1011: This function is not currently available from Fortran.
1013: .seealso: VecCreate(), MatShellSetOperation()
1014: @*/
1015: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1016: {
1019: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1020: vec->ops->viewnative = vec->ops->view;
1021: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1022: vec->ops->loadnative = vec->ops->load;
1023: }
1024: (((void(**)(void))vec->ops)[(int)op]) = f;
1025: return(0);
1026: }
1029: /*@
1030: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1031: used during the assembly process to store values that belong to
1032: other processors.
1034: Not Collective, different processes can have different size stashes
1036: Input Parameters:
1037: + vec - the vector
1038: . size - the initial size of the stash.
1039: - bsize - the initial size of the block-stash(if used).
1041: Options Database Keys:
1042: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1043: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1045: Level: intermediate
1047: Notes:
1048: The block-stash is used for values set with VecSetValuesBlocked() while
1049: the stash is used for values set with VecSetValues()
1051: Run with the option -info and look for output of the form
1052: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1053: to determine the appropriate value, MM, to use for size and
1054: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1055: to determine the value, BMM to use for bsize
1058: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1060: @*/
1061: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1062: {
1067: VecStashSetInitialSize_Private(&vec->stash,size);
1068: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1069: return(0);
1070: }
1072: /*@
1073: VecConjugate - Conjugates a vector.
1075: Logically Collective on Vec
1077: Input Parameters:
1078: . x - the vector
1080: Level: intermediate
1082: @*/
1083: PetscErrorCode VecConjugate(Vec x)
1084: {
1085: #if defined(PETSC_USE_COMPLEX)
1091: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1092: VecSetErrorIfLocked(x,1);
1093: (*x->ops->conjugate)(x);
1094: /* we need to copy norms here */
1095: PetscObjectStateIncrease((PetscObject)x);
1096: return(0);
1097: #else
1098: return(0);
1099: #endif
1100: }
1102: /*@
1103: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1105: Logically Collective on Vec
1107: Input Parameters:
1108: . x, y - the vectors
1110: Output Parameter:
1111: . w - the result
1113: Level: advanced
1115: Notes:
1116: any subset of the x, y, and w may be the same vector.
1118: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1119: @*/
1120: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1121: {
1133: VecCheckSameSize(w,1,x,2);
1134: VecCheckSameSize(w,2,y,3);
1135: VecSetErrorIfLocked(w,1);
1136: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1137: (*w->ops->pointwisemult)(w,x,y);
1138: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1139: PetscObjectStateIncrease((PetscObject)w);
1140: return(0);
1141: }
1143: /*@
1144: VecSetRandom - Sets all components of a vector to random numbers.
1146: Logically Collective on Vec
1148: Input Parameters:
1149: + x - the vector
1150: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1151: it will create one internally.
1153: Output Parameter:
1154: . x - the vector
1156: Example of Usage:
1157: .vb
1158: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1159: VecSetRandom(x,rctx);
1160: PetscRandomDestroy(&rctx);
1161: .ve
1163: Level: intermediate
1166: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1167: @*/
1168: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1169: {
1171: PetscRandom randObj = NULL;
1177: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1178: VecSetErrorIfLocked(x,1);
1180: if (!rctx) {
1181: MPI_Comm comm;
1182: PetscObjectGetComm((PetscObject)x,&comm);
1183: PetscRandomCreate(comm,&randObj);
1184: PetscRandomSetFromOptions(randObj);
1185: rctx = randObj;
1186: }
1188: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1189: (*x->ops->setrandom)(x,rctx);
1190: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1192: PetscRandomDestroy(&randObj);
1193: PetscObjectStateIncrease((PetscObject)x);
1194: return(0);
1195: }
1197: /*@
1198: VecZeroEntries - puts a 0.0 in each element of a vector
1200: Logically Collective on Vec
1202: Input Parameter:
1203: . vec - The vector
1205: Level: beginner
1207: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1208: @*/
1209: PetscErrorCode VecZeroEntries(Vec vec)
1210: {
1214: VecSet(vec,0);
1215: return(0);
1216: }
1218: /*
1219: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1220: processor and a PETSc MPI vector on more than one processor.
1222: Collective on Vec
1224: Input Parameter:
1225: . vec - The vector
1227: Level: intermediate
1229: .seealso: VecSetFromOptions(), VecSetType()
1230: */
1231: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1232: {
1233: PetscBool opt;
1234: VecType defaultType;
1235: char typeName[256];
1236: PetscMPIInt size;
1240: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1241: else {
1242: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1243: if (size > 1) defaultType = VECMPI;
1244: else defaultType = VECSEQ;
1245: }
1247: VecRegisterAll();
1248: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1249: if (opt) {
1250: VecSetType(vec, typeName);
1251: } else {
1252: VecSetType(vec, defaultType);
1253: }
1254: return(0);
1255: }
1257: /*@
1258: VecSetFromOptions - Configures the vector from the options database.
1260: Collective on Vec
1262: Input Parameter:
1263: . vec - The vector
1265: Notes:
1266: To see all options, run your program with the -help option, or consult the users manual.
1267: Must be called after VecCreate() but before the vector is used.
1269: Level: beginner
1272: .seealso: VecCreate(), VecSetOptionsPrefix()
1273: @*/
1274: PetscErrorCode VecSetFromOptions(Vec vec)
1275: {
1281: PetscObjectOptionsBegin((PetscObject)vec);
1282: /* Handle vector type options */
1283: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1285: /* Handle specific vector options */
1286: if (vec->ops->setfromoptions) {
1287: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1288: }
1290: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1291: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1292: PetscOptionsEnd();
1293: return(0);
1294: }
1296: /*@
1297: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1299: Collective on Vec
1301: Input Parameters:
1302: + v - the vector
1303: . n - the local size (or PETSC_DECIDE to have it set)
1304: - N - the global size (or PETSC_DECIDE)
1306: Notes:
1307: n and N cannot be both PETSC_DECIDE
1308: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1310: Level: intermediate
1312: .seealso: VecGetSize(), PetscSplitOwnership()
1313: @*/
1314: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1315: {
1321: 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);
1322: 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);
1323: v->map->n = n;
1324: v->map->N = N;
1325: if (v->ops->create) {
1326: (*v->ops->create)(v);
1327: v->ops->create = NULL;
1328: }
1329: return(0);
1330: }
1332: /*@
1333: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1334: and VecSetValuesBlockedLocal().
1336: Logically Collective on Vec
1338: Input Parameter:
1339: + v - the vector
1340: - bs - the blocksize
1342: Notes:
1343: All vectors obtained by VecDuplicate() inherit the same blocksize.
1345: Level: advanced
1347: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1349: @*/
1350: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1351: {
1357: PetscLayoutSetBlockSize(v->map,bs);
1358: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1359: return(0);
1360: }
1362: /*@
1363: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1364: and VecSetValuesBlockedLocal().
1366: Not Collective
1368: Input Parameter:
1369: . v - the vector
1371: Output Parameter:
1372: . bs - the blocksize
1374: Notes:
1375: All vectors obtained by VecDuplicate() inherit the same blocksize.
1377: Level: advanced
1379: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1382: @*/
1383: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1384: {
1390: PetscLayoutGetBlockSize(v->map,bs);
1391: return(0);
1392: }
1394: /*@C
1395: VecSetOptionsPrefix - Sets the prefix used for searching for all
1396: Vec options in the database.
1398: Logically Collective on Vec
1400: Input Parameter:
1401: + v - the Vec context
1402: - prefix - the prefix to prepend to all option names
1404: Notes:
1405: A hyphen (-) must NOT be given at the beginning of the prefix name.
1406: The first character of all runtime options is AUTOMATICALLY the hyphen.
1408: Level: advanced
1410: .seealso: VecSetFromOptions()
1411: @*/
1412: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1413: {
1418: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1419: return(0);
1420: }
1422: /*@C
1423: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1424: Vec options in the database.
1426: Logically Collective on Vec
1428: Input Parameters:
1429: + v - the Vec context
1430: - prefix - the prefix to prepend to all option names
1432: Notes:
1433: A hyphen (-) must NOT be given at the beginning of the prefix name.
1434: The first character of all runtime options is AUTOMATICALLY the hyphen.
1436: Level: advanced
1438: .seealso: VecGetOptionsPrefix()
1439: @*/
1440: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1441: {
1446: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1447: return(0);
1448: }
1450: /*@C
1451: VecGetOptionsPrefix - Sets the prefix used for searching for all
1452: Vec options in the database.
1454: Not Collective
1456: Input Parameter:
1457: . v - the Vec context
1459: Output Parameter:
1460: . prefix - pointer to the prefix string used
1462: Notes:
1463: On the fortran side, the user should pass in a string 'prefix' of
1464: sufficient length to hold the prefix.
1466: Level: advanced
1468: .seealso: VecAppendOptionsPrefix()
1469: @*/
1470: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1471: {
1476: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1477: return(0);
1478: }
1480: /*@
1481: VecSetUp - Sets up the internal vector data structures for the later use.
1483: Collective on Vec
1485: Input Parameters:
1486: . v - the Vec context
1488: Notes:
1489: For basic use of the Vec classes the user need not explicitly call
1490: VecSetUp(), since these actions will happen automatically.
1492: Level: advanced
1494: .seealso: VecCreate(), VecDestroy()
1495: @*/
1496: PetscErrorCode VecSetUp(Vec v)
1497: {
1498: PetscMPIInt size;
1503: if (v->map->n < 0 && v->map->N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1504: if (!((PetscObject)v)->type_name) {
1505: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1506: if (size == 1) {
1507: VecSetType(v, VECSEQ);
1508: } else {
1509: VecSetType(v, VECMPI);
1510: }
1511: }
1512: return(0);
1513: }
1515: /*
1516: These currently expose the PetscScalar/PetscReal in updating the
1517: cached norm. If we push those down into the implementation these
1518: will become independent of PetscScalar/PetscReal
1519: */
1521: /*@
1522: VecCopy - Copies a vector. y <- x
1524: Logically Collective on Vec
1526: Input Parameter:
1527: . x - the vector
1529: Output Parameter:
1530: . y - the copy
1532: Notes:
1533: For default parallel PETSc vectors, both x and y must be distributed in
1534: the same manner; local copies are done.
1536: Developer Notes:
1538: of the vectors to be sequential and one to be parallel so long as both have the same
1539: local sizes. This is used in some internal functions in PETSc.
1541: Level: beginner
1543: .seealso: VecDuplicate()
1544: @*/
1545: PetscErrorCode VecCopy(Vec x,Vec y)
1546: {
1547: PetscBool flgs[4];
1548: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1550: PetscInt i;
1557: if (x == y) return(0);
1558: VecCheckSameLocalSize(x,1,y,2);
1559: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1560: VecSetErrorIfLocked(y,2);
1562: #if !defined(PETSC_USE_MIXED_PRECISION)
1563: for (i=0; i<4; i++) {
1564: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1565: }
1566: #endif
1568: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1569: #if defined(PETSC_USE_MIXED_PRECISION)
1570: extern PetscErrorCode VecGetArray(Vec,double**);
1571: extern PetscErrorCode VecRestoreArray(Vec,double**);
1572: extern PetscErrorCode VecGetArray(Vec,float**);
1573: extern PetscErrorCode VecRestoreArray(Vec,float**);
1574: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1575: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1576: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1577: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1578: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1579: PetscInt i,n;
1580: const float *xx;
1581: double *yy;
1582: VecGetArrayRead(x,&xx);
1583: VecGetArray(y,&yy);
1584: VecGetLocalSize(x,&n);
1585: for (i=0; i<n; i++) yy[i] = xx[i];
1586: VecRestoreArrayRead(x,&xx);
1587: VecRestoreArray(y,&yy);
1588: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1589: PetscInt i,n;
1590: float *yy;
1591: const double *xx;
1592: VecGetArrayRead(x,&xx);
1593: VecGetArray(y,&yy);
1594: VecGetLocalSize(x,&n);
1595: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1596: VecRestoreArrayRead(x,&xx);
1597: VecRestoreArray(y,&yy);
1598: } else {
1599: (*x->ops->copy)(x,y);
1600: }
1601: #else
1602: (*x->ops->copy)(x,y);
1603: #endif
1605: PetscObjectStateIncrease((PetscObject)y);
1606: #if !defined(PETSC_USE_MIXED_PRECISION)
1607: for (i=0; i<4; i++) {
1608: if (flgs[i]) {
1609: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1610: }
1611: }
1612: #endif
1614: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1615: return(0);
1616: }
1618: /*@
1619: VecSwap - Swaps the vectors x and y.
1621: Logically Collective on Vec
1623: Input Parameters:
1624: . x, y - the vectors
1626: Level: advanced
1628: @*/
1629: PetscErrorCode VecSwap(Vec x,Vec y)
1630: {
1631: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1632: PetscBool flgxs[4],flgys[4];
1634: PetscInt i;
1642: VecCheckSameSize(x,1,y,2);
1643: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1644: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1645: VecSetErrorIfLocked(x,1);
1646: VecSetErrorIfLocked(y,2);
1648: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1649: for (i=0; i<4; i++) {
1650: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1651: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1652: }
1653: (*x->ops->swap)(x,y);
1654: PetscObjectStateIncrease((PetscObject)x);
1655: PetscObjectStateIncrease((PetscObject)y);
1656: for (i=0; i<4; i++) {
1657: if (flgxs[i]) {
1658: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1659: }
1660: if (flgys[i]) {
1661: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1662: }
1663: }
1664: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1665: return(0);
1666: }
1668: /*
1669: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1671: Collective on VecStash
1673: Input Parameters:
1674: + obj - the VecStash object
1675: . bobj - optional other object that provides the prefix
1676: - optionname - option to activate viewing
1678: Level: intermediate
1680: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1682: */
1683: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1684: {
1685: PetscErrorCode ierr;
1686: PetscViewer viewer;
1687: PetscBool flg;
1688: PetscViewerFormat format;
1689: char *prefix;
1692: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1693: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1694: if (flg) {
1695: PetscViewerPushFormat(viewer,format);
1696: VecStashView(obj,viewer);
1697: PetscViewerPopFormat(viewer);
1698: PetscViewerDestroy(&viewer);
1699: }
1700: return(0);
1701: }
1703: /*@
1704: VecStashView - Prints the entries in the vector stash and block stash.
1706: Collective on Vec
1708: Input Parameters:
1709: + v - the vector
1710: - viewer - the viewer
1712: Level: advanced
1715: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1717: @*/
1718: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1719: {
1721: PetscMPIInt rank;
1722: PetscInt i,j;
1723: PetscBool match;
1724: VecStash *s;
1725: PetscScalar val;
1732: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1733: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1734: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1735: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1736: s = &v->bstash;
1738: /* print block stash */
1739: PetscViewerASCIIPushSynchronized(viewer);
1740: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1741: for (i=0; i<s->n; i++) {
1742: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1743: for (j=0; j<s->bs; j++) {
1744: val = s->array[i*s->bs+j];
1745: #if defined(PETSC_USE_COMPLEX)
1746: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1747: #else
1748: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1749: #endif
1750: }
1751: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1752: }
1753: PetscViewerFlush(viewer);
1755: s = &v->stash;
1757: /* print basic stash */
1758: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1759: for (i=0; i<s->n; i++) {
1760: val = s->array[i];
1761: #if defined(PETSC_USE_COMPLEX)
1762: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1763: #else
1764: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1765: #endif
1766: }
1767: PetscViewerFlush(viewer);
1768: PetscViewerASCIIPopSynchronized(viewer);
1769: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1770: return(0);
1771: }
1773: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1774: {
1775: PetscInt i,N,rstart,rend;
1777: PetscScalar *xx;
1778: PetscReal *xreal;
1779: PetscBool iset;
1782: VecGetOwnershipRange(v,&rstart,&rend);
1783: VecGetSize(v,&N);
1784: PetscCalloc1(N,&xreal);
1785: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1786: if (iset) {
1787: VecGetArray(v,&xx);
1788: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1789: VecRestoreArray(v,&xx);
1790: }
1791: PetscFree(xreal);
1792: if (set) *set = iset;
1793: return(0);
1794: }
1796: /*@
1797: VecGetLayout - get PetscLayout describing vector layout
1799: Not Collective
1801: Input Arguments:
1802: . x - the vector
1804: Output Arguments:
1805: . map - the layout
1807: Level: developer
1809: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1810: @*/
1811: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1812: {
1816: *map = x->map;
1817: return(0);
1818: }
1820: /*@
1821: VecSetLayout - set PetscLayout describing vector layout
1823: Not Collective
1825: Input Arguments:
1826: + x - the vector
1827: - map - the layout
1829: Notes:
1830: It is normally only valid to replace the layout with a layout known to be equivalent.
1832: Level: developer
1834: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1835: @*/
1836: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1837: {
1842: PetscLayoutReference(map,&x->map);
1843: return(0);
1844: }
1846: PetscErrorCode VecSetInf(Vec xin)
1847: {
1848: PetscInt i,n = xin->map->n;
1849: PetscScalar *xx;
1850: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1854: VecGetArrayWrite(xin,&xx);
1855: for (i=0; i<n; i++) xx[i] = inf;
1856: VecRestoreArrayWrite(xin,&xx);
1857: return(0);
1858: }
1860: /*@
1861: VecBindToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
1863: Input Parameters:
1864: + v - the vector
1865: - flg - bind to the CPU if value of PETSC_TRUE
1867: Level: intermediate
1868: @*/
1869: PetscErrorCode VecBindToCPU(Vec v,PetscBool flg)
1870: {
1871: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1875: if (v->boundtocpu == flg) return(0);
1876: v->boundtocpu = flg;
1877: if (v->ops->bindtocpu) {
1878: (*v->ops->bindtocpu)(v,flg);
1879: }
1880: return(0);
1881: #else
1882: return 0;
1883: #endif
1884: }
1886: /*@C
1887: VecSetPinnedMemoryMin - Set the minimum data size for which pinned memory will be used for host (CPU) allocations.
1889: Logically Collective on Vec
1891: Input Parameters:
1892: + v - the vector
1893: - mbytes - minimum data size in bytes
1895: Options Database Keys:
1897: . -vec_pinned_memory_min <size> - minimum size (in bytes) for an allocation to use pinned memory on host.
1898: Note that this takes a PetscScalar, to accommodate large values;
1899: specifying -1 ensures that pinned memory will never be used.
1901: Level: developer
1903: .seealso: VecGetPinnedMemoryMin()
1904: @*/
1905: PetscErrorCode VecSetPinnedMemoryMin(Vec v,size_t mbytes)
1906: {
1907: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1909: v->minimum_bytes_pinned_memory = mbytes;
1910: return(0);
1911: #else
1912: return 0;
1913: #endif
1914: }
1916: /*@C
1917: VecGetPinnedMemoryMin - Get the minimum data size for which pinned memory will be used for host (CPU) allocations.
1919: Logically Collective on Vec
1921: Input Parameters:
1922: . v - the vector
1924: Output Parameters:
1925: . mbytes - minimum data size in bytes
1927: Level: developer
1929: .seealso: VecSetPinnedMemoryMin()
1930: @*/
1931: PetscErrorCode VecGetPinnedMemoryMin(Vec v,size_t *mbytes)
1932: {
1933: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1935: *mbytes = v->minimum_bytes_pinned_memory;
1936: return(0);
1937: #else
1938: return 0;
1939: #endif
1940: }
1942: /*@
1943: VecGetOffloadMask - Get the offload mask of a Vec.
1945: Not Collective
1947: Input Parameters:
1948: . v - the vector
1950: Output Parameters:
1951: . mask - corresponding PetscOffloadMask enum value.
1953: Level: intermediate
1955: .seealso: VecCreateSeqCUDA(), VecCreateSeqViennaCL(), VecGetArray(), VecGetType()
1956: @*/
1957: PetscErrorCode VecGetOffloadMask(Vec v,PetscOffloadMask* mask)
1958: {
1960: *mask = v->offloadmask;
1961: return(0);
1962: }