Actual source code: vector.c
petsc-3.12.5 2020-03-29
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: (*w->ops->pointwisemax)(w,x,y);
213: PetscObjectStateIncrease((PetscObject)w);
214: return(0);
215: }
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: (*w->ops->pointwisemin)(w,x,y);
253: PetscObjectStateIncrease((PetscObject)w);
254: return(0);
255: }
257: /*@
258: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
260: Logically Collective on Vec
262: Input Parameters:
263: . x, y - the vectors
265: Output Parameter:
266: . w - the result
268: Level: advanced
270: Notes:
271: any subset of the x, y, and w may be the same vector.
273: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
274: @*/
275: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
276: {
288: VecCheckSameSize(w,1,x,2);
289: VecCheckSameSize(w,1,y,3);
290: (*w->ops->pointwisemaxabs)(w,x,y);
291: PetscObjectStateIncrease((PetscObject)w);
292: return(0);
293: }
295: /*@
296: VecPointwiseDivide - Computes the componentwise division w = x/y.
298: Logically Collective on Vec
300: Input Parameters:
301: . x, y - the vectors
303: Output Parameter:
304: . w - the result
306: Level: advanced
308: Notes:
309: any subset of the x, y, and w may be the same vector.
311: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
312: @*/
313: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
314: {
326: VecCheckSameSize(w,1,x,2);
327: VecCheckSameSize(w,1,y,3);
328: (*w->ops->pointwisedivide)(w,x,y);
329: PetscObjectStateIncrease((PetscObject)w);
330: return(0);
331: }
334: /*@
335: VecDuplicate - Creates a new vector of the same type as an existing vector.
337: Collective on Vec
339: Input Parameters:
340: . v - a vector to mimic
342: Output Parameter:
343: . newv - location to put new vector
345: Notes:
346: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
347: for the new vector. Use VecCopy() to copy a vector.
349: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
350: vectors.
352: Level: beginner
354: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
355: @*/
356: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
357: {
364: (*v->ops->duplicate)(v,newv);
365: PetscObjectStateIncrease((PetscObject)*newv);
366: return(0);
367: }
369: /*@
370: VecDestroy - Destroys a vector.
372: Collective on Vec
374: Input Parameters:
375: . v - the vector
377: Level: beginner
379: .seealso: VecDuplicate(), VecDestroyVecs()
380: @*/
381: PetscErrorCode VecDestroy(Vec *v)
382: {
386: if (!*v) return(0);
388: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
390: PetscObjectSAWsViewOff((PetscObject)*v);
391: /* destroy the internal part */
392: if ((*v)->ops->destroy) {
393: (*(*v)->ops->destroy)(*v);
394: }
395: /* destroy the external/common part */
396: PetscLayoutDestroy(&(*v)->map);
397: PetscHeaderDestroy(v);
398: return(0);
399: }
401: /*@C
402: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
404: Collective on Vec
406: Input Parameters:
407: + m - the number of vectors to obtain
408: - v - a vector to mimic
410: Output Parameter:
411: . V - location to put pointer to array of vectors
413: Notes:
414: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
415: vector.
417: Fortran Note:
418: The Fortran interface is slightly different from that given below, it
419: requires one to pass in V a Vec (integer) array of size at least m.
420: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
422: Level: intermediate
424: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
425: @*/
426: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
427: {
434: (*v->ops->duplicatevecs)(v,m,V);
435: return(0);
436: }
438: /*@C
439: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
441: Collective on Vec
443: Input Parameters:
444: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
445: - m - the number of vectors previously obtained, if zero no vectors are destroyed
447: Fortran Note:
448: The Fortran interface is slightly different from that given below.
449: See the Fortran chapter of the users manual
451: Level: intermediate
453: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
454: @*/
455: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
456: {
461: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
462: if (!m || !*vv) {*vv = NULL; return(0);}
465: (*(**vv)->ops->destroyvecs)(m,*vv);
466: *vv = NULL;
467: return(0);
468: }
470: /*@C
471: VecView - Views a vector object.
473: Collective on Vec
475: Input Parameters:
476: + vec - the vector
477: - viewer - an optional visualization context
479: Notes:
480: The available visualization contexts include
481: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
482: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
483: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
485: You can change the format the vector is printed using the
486: option PetscViewerPushFormat().
488: The user can open alternative viewers with
489: + PetscViewerASCIIOpen() - Outputs vector to a specified file
490: . PetscViewerBinaryOpen() - Outputs vector in binary to a
491: specified file; corresponding input uses VecLoad()
492: . PetscViewerDrawOpen() - Outputs vector to an X window display
493: . PetscViewerSocketOpen() - Outputs vector to Socket viewer
494: - PetscViewerHDF5Open() - Outputs vector to HDF5 file viewer
496: The user can call PetscViewerPushFormat() to specify the output
497: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
498: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
499: + PETSC_VIEWER_DEFAULT - default, prints vector contents
500: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
501: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
502: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
503: format common among all vector types
505: Notes:
506: You can pass any number of vector objects, or other PETSc objects to the same viewer.
508: Notes for binary viewer:
509: If you pass multiple vectors to a binary viewer you can read them back in in the same order
510: with VecLoad().
512: If the blocksize of the vector is greater than one then you must provide a unique prefix to
513: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
514: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
515: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
516: filename. If you copy the binary file, make sure you copy the associated .info file with it.
517:
518: See the manual page for VecLoad() on the exact format the binary viewer stores
519: the values in the file.
522: Notes for HDF5 Viewer:
523: The name of the Vec (given with PetscObjectSetName() is the name that is used
524: for the object in the HDF5 file. If you wish to store the same Vec into multiple
525: datasets in the same file (typically with different values), you must change its
526: name each time before calling the VecView(). To load the same vector,
527: the name of the Vec object passed to VecLoad() must be the same.
528:
529: If the block size of the vector is greater than 1 then it is used as the first dimension in the HDF5 array.
530: If the function PetscViewerHDF5SetBaseDimension2()is called then even if the block size is one it will
531: be used as the first dimension in the HDF5 array (that is the HDF5 array will always be two dimensional)
532: See also PetscViewerHDF5SetTimestep() which adds an additional complication to reading and writing Vecs
533: with the HDF5 viewer.
534:
535: Level: beginner
538: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
539: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
540: PetscRealView(), PetscScalarView(), PetscIntView(), PetscViewerHDF5SetTimestep()
541: @*/
542: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
543: {
544: PetscErrorCode ierr;
545: PetscBool iascii;
546: PetscViewerFormat format;
547: PetscMPIInt size;
552: if (!viewer) {
553: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
554: }
556: PetscViewerGetFormat(viewer,&format);
557: MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
558: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
560: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
562: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
563: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
564: if (iascii) {
565: PetscInt rows,bs;
567: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
568: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
569: PetscViewerASCIIPushTab(viewer);
570: VecGetSize(vec,&rows);
571: VecGetBlockSize(vec,&bs);
572: if (bs != 1) {
573: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
574: } else {
575: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
576: }
577: PetscViewerASCIIPopTab(viewer);
578: }
579: }
580: VecLockReadPush(vec);
581: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
582: (*vec->ops->viewnative)(vec,viewer);
583: } else {
584: (*vec->ops->view)(vec,viewer);
585: }
586: VecLockReadPop(vec);
587: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
588: return(0);
589: }
591: #if defined(PETSC_USE_DEBUG)
592: #include <../src/sys/totalview/tv_data_display.h>
593: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
594: {
595: const PetscScalar *values;
596: char type[32];
597: PetscErrorCode ierr;
600: TV_add_row("Local rows", "int", &v->map->n);
601: TV_add_row("Global rows", "int", &v->map->N);
602: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
603: VecGetArrayRead((Vec)v,&values);
604: PetscSNPrintf(type,32,"double[%d]",v->map->n);
605: TV_add_row("values",type, values);
606: VecRestoreArrayRead((Vec)v,&values);
607: return TV_format_OK;
608: }
609: #endif
611: /*@
612: VecGetSize - Returns the global number of elements of the vector.
614: Not Collective
616: Input Parameter:
617: . x - the vector
619: Output Parameters:
620: . size - the global length of the vector
622: Level: beginner
624: .seealso: VecGetLocalSize()
625: @*/
626: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
627: {
634: (*x->ops->getsize)(x,size);
635: return(0);
636: }
638: /*@
639: VecGetLocalSize - Returns the number of elements of the vector stored
640: in local memory. This routine may be implementation dependent, so use
641: with care.
643: Not Collective
645: Input Parameter:
646: . x - the vector
648: Output Parameter:
649: . size - the length of the local piece of the vector
651: Level: beginner
653: .seealso: VecGetSize()
654: @*/
655: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
656: {
663: (*x->ops->getlocalsize)(x,size);
664: return(0);
665: }
667: /*@C
668: VecGetOwnershipRange - Returns the range of indices owned by
669: this processor, assuming that the vectors are laid out with the
670: first n1 elements on the first processor, next n2 elements on the
671: second, etc. For certain parallel layouts this range may not be
672: well defined.
674: Not Collective
676: Input Parameter:
677: . x - the vector
679: Output Parameters:
680: + low - the first local element, pass in NULL if not interested
681: - high - one more than the last local element, pass in NULL if not interested
683: Note:
684: The high argument is one more than the last element stored locally.
686: Fortran: PETSC_NULL_INTEGER should be used instead of NULL
688: Level: beginner
691: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
692: @*/
693: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
694: {
700: if (low) *low = x->map->rstart;
701: if (high) *high = x->map->rend;
702: return(0);
703: }
705: /*@C
706: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
707: assuming that the vectors are laid out with the
708: first n1 elements on the first processor, next n2 elements on the
709: second, etc. For certain parallel layouts this range may not be
710: well defined.
712: Not Collective
714: Input Parameter:
715: . x - the vector
717: Output Parameters:
718: . range - array of length size+1 with the start and end+1 for each process
720: Note:
721: The high argument is one more than the last element stored locally.
723: Fortran: You must PASS in an array of length size+1
725: Level: beginner
728: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
729: @*/
730: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
731: {
737: PetscLayoutGetRanges(x->map,ranges);
738: return(0);
739: }
741: /*@
742: VecSetOption - Sets an option for controling a vector's behavior.
744: Collective on Vec
746: Input Parameter:
747: + x - the vector
748: . op - the option
749: - flag - turn the option on or off
751: Supported Options:
752: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
753: entries destined to be stored on a separate processor. This can be used
754: to eliminate the global reduction in the VecAssemblyXXXX() if you know
755: that you have only used VecSetValues() to set local elements
756: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
757: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
758: ignored.
759: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
760: entries will always be a subset (possibly equal) of the off-process entries set on the
761: first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
762: changed this flag afterwards. If this assembly is not such first assembly, then this
763: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
764: a global reduction. Subsequent assemblies setting off-process values should use the same
765: InsertMode as the first assembly.
767: Developer Note:
768: The InsertMode restriction could be removed by packing the stash messages out of place.
770: Level: intermediate
772: @*/
773: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
774: {
780: if (x->ops->setoption) {
781: (*x->ops->setoption)(x,op,flag);
782: }
783: return(0);
784: }
786: /* Default routines for obtaining and releasing; */
787: /* may be used by any implementation */
788: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
789: {
791: PetscInt i;
796: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
797: PetscMalloc1(m,V);
798: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
799: return(0);
800: }
802: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
803: {
805: PetscInt i;
809: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
810: PetscFree(v);
811: return(0);
812: }
814: /*@
815: VecResetArray - Resets a vector to use its default memory. Call this
816: after the use of VecPlaceArray().
818: Not Collective
820: Input Parameters:
821: . vec - the vector
823: Level: developer
825: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
827: @*/
828: PetscErrorCode VecResetArray(Vec vec)
829: {
835: if (vec->ops->resetarray) {
836: (*vec->ops->resetarray)(vec);
837: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
838: PetscObjectStateIncrease((PetscObject)vec);
839: return(0);
840: }
842: /*@C
843: VecLoad - Loads a vector that has been stored in binary or HDF5 format
844: with VecView().
846: Collective on PetscViewer
848: Input Parameters:
849: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
850: some related function before a call to VecLoad().
851: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
852: HDF5 file viewer, obtained from PetscViewerHDF5Open()
854: Level: intermediate
856: Notes:
857: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
858: before calling this.
860: The input file must contain the full global vector, as
861: written by the routine VecView().
863: If the type or size of newvec is not set before a call to VecLoad, PETSc
864: sets the type and the local and global sizes. If type and/or
865: sizes are already set, then the same are used.
867: If using the binary viewer and the blocksize of the vector is greater than one then you must provide a unique prefix to
868: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
869: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
870: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
871: filename. If you copy the binary file, make sure you copy the associated .info file with it.
873: If using HDF5, you must assign the Vec the same name as was used in the Vec
874: that was stored in the file using PetscObjectSetName(). Otherwise you will
875: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT".
877: If the HDF5 file contains a two dimensional array the first dimension is treated as the block size
878: in loading the vector. Hence, for example, using Matlab notation h5create('vector.dat','/Test_Vec',[27 1]);
879: will load a vector of size 27 and block size 27 thus resulting in all 27 entries being on the first process of
880: vectors communicator and the rest of the processes having zero entries
882: Notes for advanced users when using the binary viewer:
883: Most users should not need to know the details of the binary storage
884: format, since VecLoad() and VecView() completely hide these details.
885: But for anyone who's interested, the standard binary vector storage
886: format is
887: .vb
888: PetscInt VEC_FILE_CLASSID
889: PetscInt number of rows
890: PetscScalar *values of all entries
891: .ve
893: In addition, PETSc automatically uses byte swapping to work on all machines; the files
894: are written ALWAYS using big-endian ordering. On small-endian machines the numbers
895: are converted to the small-endian format when they are read in from the file.
896: See PetscBinaryRead() and PetscBinaryWrite() to see how this may be done.
898: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
899: @*/
900: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
901: {
902: PetscErrorCode ierr;
903: PetscBool isbinary,ishdf5,isadios,isadios2;
904: PetscViewerFormat format;
909: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
910: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
911: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
912: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
913: if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
915: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
916: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
917: VecSetType(newvec, VECSTANDARD);
918: }
919: PetscViewerGetFormat(viewer,&format);
920: if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
921: (*newvec->ops->loadnative)(newvec,viewer);
922: } else {
923: (*newvec->ops->load)(newvec,viewer);
924: }
925: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
926: return(0);
927: }
930: /*@
931: VecReciprocal - Replaces each component of a vector by its reciprocal.
933: Logically Collective on Vec
935: Input Parameter:
936: . vec - the vector
938: Output Parameter:
939: . vec - the vector reciprocal
941: Level: intermediate
943: .seealso: VecLog(), VecExp(), VecSqrtAbs()
945: @*/
946: PetscErrorCode VecReciprocal(Vec vec)
947: {
953: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
954: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
955: (*vec->ops->reciprocal)(vec);
956: PetscObjectStateIncrease((PetscObject)vec);
957: return(0);
958: }
960: /*@C
961: VecSetOperation - Allows user to set a vector operation.
963: Logically Collective on Vec
965: Input Parameters:
966: + vec - the vector
967: . op - the name of the operation
968: - f - the function that provides the operation.
970: Level: advanced
972: Usage:
973: $ PetscErrorCode userview(Vec,PetscViewer);
974: $ VecCreateMPI(comm,m,M,&x);
975: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
977: Notes:
978: See the file include/petscvec.h for a complete list of matrix
979: operations, which all have the form VECOP_<OPERATION>, where
980: <OPERATION> is the name (in all capital letters) of the
981: user interface routine (e.g., VecView() -> VECOP_VIEW).
983: This function is not currently available from Fortran.
985: .seealso: VecCreate(), MatShellSetOperation()
986: @*/
987: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
988: {
991: if (op == VECOP_VIEW && !vec->ops->viewnative) {
992: vec->ops->viewnative = vec->ops->view;
993: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
994: vec->ops->loadnative = vec->ops->load;
995: }
996: (((void(**)(void))vec->ops)[(int)op]) = f;
997: return(0);
998: }
1001: /*@
1002: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1003: used during the assembly process to store values that belong to
1004: other processors.
1006: Not Collective, different processes can have different size stashes
1008: Input Parameters:
1009: + vec - the vector
1010: . size - the initial size of the stash.
1011: - bsize - the initial size of the block-stash(if used).
1013: Options Database Keys:
1014: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1015: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1017: Level: intermediate
1019: Notes:
1020: The block-stash is used for values set with VecSetValuesBlocked() while
1021: the stash is used for values set with VecSetValues()
1023: Run with the option -info and look for output of the form
1024: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1025: to determine the appropriate value, MM, to use for size and
1026: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1027: to determine the value, BMM to use for bsize
1030: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1032: @*/
1033: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1034: {
1039: VecStashSetInitialSize_Private(&vec->stash,size);
1040: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1041: return(0);
1042: }
1044: /*@
1045: VecConjugate - Conjugates a vector.
1047: Logically Collective on Vec
1049: Input Parameters:
1050: . x - the vector
1052: Level: intermediate
1054: @*/
1055: PetscErrorCode VecConjugate(Vec x)
1056: {
1057: #if defined(PETSC_USE_COMPLEX)
1063: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1064: (*x->ops->conjugate)(x);
1065: /* we need to copy norms here */
1066: PetscObjectStateIncrease((PetscObject)x);
1067: return(0);
1068: #else
1069: return(0);
1070: #endif
1071: }
1073: /*@
1074: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1076: Logically Collective on Vec
1078: Input Parameters:
1079: . x, y - the vectors
1081: Output Parameter:
1082: . w - the result
1084: Level: advanced
1086: Notes:
1087: any subset of the x, y, and w may be the same vector.
1089: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1090: @*/
1091: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1092: {
1104: VecCheckSameSize(w,1,x,2);
1105: VecCheckSameSize(w,2,y,3);
1106: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1107: (*w->ops->pointwisemult)(w,x,y);
1108: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1109: PetscObjectStateIncrease((PetscObject)w);
1110: return(0);
1111: }
1113: /*@
1114: VecSetRandom - Sets all components of a vector to random numbers.
1116: Logically Collective on Vec
1118: Input Parameters:
1119: + x - the vector
1120: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1121: it will create one internally.
1123: Output Parameter:
1124: . x - the vector
1126: Example of Usage:
1127: .vb
1128: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1129: VecSetRandom(x,rctx);
1130: PetscRandomDestroy(rctx);
1131: .ve
1133: Level: intermediate
1136: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1137: @*/
1138: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1139: {
1141: PetscRandom randObj = NULL;
1147: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1149: if (!rctx) {
1150: MPI_Comm comm;
1151: PetscObjectGetComm((PetscObject)x,&comm);
1152: PetscRandomCreate(comm,&randObj);
1153: PetscRandomSetFromOptions(randObj);
1154: rctx = randObj;
1155: }
1157: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1158: (*x->ops->setrandom)(x,rctx);
1159: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1161: PetscRandomDestroy(&randObj);
1162: PetscObjectStateIncrease((PetscObject)x);
1163: return(0);
1164: }
1166: /*@
1167: VecZeroEntries - puts a 0.0 in each element of a vector
1169: Logically Collective on Vec
1171: Input Parameter:
1172: . vec - The vector
1174: Level: beginner
1176: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1177: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1178: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1179: this routine should not exist.
1181: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1182: @*/
1183: PetscErrorCode VecZeroEntries(Vec vec)
1184: {
1188: VecSet(vec,0);
1189: return(0);
1190: }
1192: /*
1193: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1194: processor and a PETSc MPI vector on more than one processor.
1196: Collective on Vec
1198: Input Parameter:
1199: . vec - The vector
1201: Level: intermediate
1203: .seealso: VecSetFromOptions(), VecSetType()
1204: */
1205: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1206: {
1207: PetscBool opt;
1208: VecType defaultType;
1209: char typeName[256];
1210: PetscMPIInt size;
1214: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1215: else {
1216: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1217: if (size > 1) defaultType = VECMPI;
1218: else defaultType = VECSEQ;
1219: }
1221: VecRegisterAll();
1222: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1223: if (opt) {
1224: VecSetType(vec, typeName);
1225: } else {
1226: VecSetType(vec, defaultType);
1227: }
1228: return(0);
1229: }
1231: /*@
1232: VecSetFromOptions - Configures the vector from the options database.
1234: Collective on Vec
1236: Input Parameter:
1237: . vec - The vector
1239: Notes:
1240: To see all options, run your program with the -help option, or consult the users manual.
1241: Must be called after VecCreate() but before the vector is used.
1243: Level: beginner
1246: .seealso: VecCreate(), VecSetOptionsPrefix()
1247: @*/
1248: PetscErrorCode VecSetFromOptions(Vec vec)
1249: {
1255: PetscObjectOptionsBegin((PetscObject)vec);
1256: /* Handle vector type options */
1257: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1259: /* Handle specific vector options */
1260: if (vec->ops->setfromoptions) {
1261: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1262: }
1264: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1265: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1266: PetscOptionsEnd();
1267: return(0);
1268: }
1270: /*@
1271: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1273: Collective on Vec
1275: Input Parameters:
1276: + v - the vector
1277: . n - the local size (or PETSC_DECIDE to have it set)
1278: - N - the global size (or PETSC_DECIDE)
1280: Notes:
1281: n and N cannot be both PETSC_DECIDE
1282: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1284: Level: intermediate
1286: .seealso: VecGetSize(), PetscSplitOwnership()
1287: @*/
1288: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1289: {
1295: 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);
1296: 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);
1297: v->map->n = n;
1298: v->map->N = N;
1299: if (v->ops->create) {
1300: (*v->ops->create)(v);
1301: v->ops->create = 0;
1302: }
1303: return(0);
1304: }
1306: /*@
1307: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1308: and VecSetValuesBlockedLocal().
1310: Logically Collective on Vec
1312: Input Parameter:
1313: + v - the vector
1314: - bs - the blocksize
1316: Notes:
1317: All vectors obtained by VecDuplicate() inherit the same blocksize.
1319: Level: advanced
1321: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1323: @*/
1324: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1325: {
1330: if (bs < 0 || bs == v->map->bs) return(0);
1332: PetscLayoutSetBlockSize(v->map,bs);
1333: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1334: return(0);
1335: }
1337: /*@
1338: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1339: and VecSetValuesBlockedLocal().
1341: Not Collective
1343: Input Parameter:
1344: . v - the vector
1346: Output Parameter:
1347: . bs - the blocksize
1349: Notes:
1350: All vectors obtained by VecDuplicate() inherit the same blocksize.
1352: Level: advanced
1354: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1357: @*/
1358: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1359: {
1365: PetscLayoutGetBlockSize(v->map,bs);
1366: return(0);
1367: }
1369: /*@C
1370: VecSetOptionsPrefix - Sets the prefix used for searching for all
1371: Vec options in the database.
1373: Logically Collective on Vec
1375: Input Parameter:
1376: + v - the Vec context
1377: - prefix - the prefix to prepend to all option names
1379: Notes:
1380: A hyphen (-) must NOT be given at the beginning of the prefix name.
1381: The first character of all runtime options is AUTOMATICALLY the hyphen.
1383: Level: advanced
1385: .seealso: VecSetFromOptions()
1386: @*/
1387: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1388: {
1393: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1394: return(0);
1395: }
1397: /*@C
1398: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1399: Vec options in the database.
1401: Logically Collective on Vec
1403: Input Parameters:
1404: + v - the Vec context
1405: - prefix - the prefix to prepend to all option names
1407: Notes:
1408: A hyphen (-) must NOT be given at the beginning of the prefix name.
1409: The first character of all runtime options is AUTOMATICALLY the hyphen.
1411: Level: advanced
1413: .seealso: VecGetOptionsPrefix()
1414: @*/
1415: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1416: {
1421: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1422: return(0);
1423: }
1425: /*@C
1426: VecGetOptionsPrefix - Sets the prefix used for searching for all
1427: Vec options in the database.
1429: Not Collective
1431: Input Parameter:
1432: . v - the Vec context
1434: Output Parameter:
1435: . prefix - pointer to the prefix string used
1437: Notes:
1438: On the fortran side, the user should pass in a string 'prefix' of
1439: sufficient length to hold the prefix.
1441: Level: advanced
1443: .seealso: VecAppendOptionsPrefix()
1444: @*/
1445: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1446: {
1451: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1452: return(0);
1453: }
1455: /*@
1456: VecSetUp - Sets up the internal vector data structures for the later use.
1458: Collective on Vec
1460: Input Parameters:
1461: . v - the Vec context
1463: Notes:
1464: For basic use of the Vec classes the user need not explicitly call
1465: VecSetUp(), since these actions will happen automatically.
1467: Level: advanced
1469: .seealso: VecCreate(), VecDestroy()
1470: @*/
1471: PetscErrorCode VecSetUp(Vec v)
1472: {
1473: PetscMPIInt size;
1478: if (v->map->n < 0 && v->map->N < 0) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Sizes not set");
1479: if (!((PetscObject)v)->type_name) {
1480: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1481: if (size == 1) {
1482: VecSetType(v, VECSEQ);
1483: } else {
1484: VecSetType(v, VECMPI);
1485: }
1486: }
1487: return(0);
1488: }
1490: /*
1491: These currently expose the PetscScalar/PetscReal in updating the
1492: cached norm. If we push those down into the implementation these
1493: will become independent of PetscScalar/PetscReal
1494: */
1496: /*@
1497: VecCopy - Copies a vector. y <- x
1499: Logically Collective on Vec
1501: Input Parameter:
1502: . x - the vector
1504: Output Parameter:
1505: . y - the copy
1507: Notes:
1508: For default parallel PETSc vectors, both x and y must be distributed in
1509: the same manner; local copies are done.
1511: Developer Notes:
1513: of the vectors to be sequential and one to be parallel so long as both have the same
1514: local sizes. This is used in some internal functions in PETSc.
1516: Level: beginner
1518: .seealso: VecDuplicate()
1519: @*/
1520: PetscErrorCode VecCopy(Vec x,Vec y)
1521: {
1522: PetscBool flgs[4];
1523: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1525: PetscInt i;
1532: if (x == y) return(0);
1533: VecCheckSameLocalSize(x,1,y,2);
1534: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1535: VecSetErrorIfLocked(y,2);
1537: #if !defined(PETSC_USE_MIXED_PRECISION)
1538: for (i=0; i<4; i++) {
1539: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1540: }
1541: #endif
1543: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1544: #if defined(PETSC_USE_MIXED_PRECISION)
1545: extern PetscErrorCode VecGetArray(Vec,double**);
1546: extern PetscErrorCode VecRestoreArray(Vec,double**);
1547: extern PetscErrorCode VecGetArray(Vec,float**);
1548: extern PetscErrorCode VecRestoreArray(Vec,float**);
1549: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1550: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1551: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1552: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1553: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1554: PetscInt i,n;
1555: const float *xx;
1556: double *yy;
1557: VecGetArrayRead(x,&xx);
1558: VecGetArray(y,&yy);
1559: VecGetLocalSize(x,&n);
1560: for (i=0; i<n; i++) yy[i] = xx[i];
1561: VecRestoreArrayRead(x,&xx);
1562: VecRestoreArray(y,&yy);
1563: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1564: PetscInt i,n;
1565: float *yy;
1566: const double *xx;
1567: VecGetArrayRead(x,&xx);
1568: VecGetArray(y,&yy);
1569: VecGetLocalSize(x,&n);
1570: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1571: VecRestoreArrayRead(x,&xx);
1572: VecRestoreArray(y,&yy);
1573: } else {
1574: (*x->ops->copy)(x,y);
1575: }
1576: #else
1577: (*x->ops->copy)(x,y);
1578: #endif
1580: PetscObjectStateIncrease((PetscObject)y);
1581: #if !defined(PETSC_USE_MIXED_PRECISION)
1582: for (i=0; i<4; i++) {
1583: if (flgs[i]) {
1584: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1585: }
1586: }
1587: #endif
1589: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1590: return(0);
1591: }
1593: /*@
1594: VecSwap - Swaps the vectors x and y.
1596: Logically Collective on Vec
1598: Input Parameters:
1599: . x, y - the vectors
1601: Level: advanced
1603: @*/
1604: PetscErrorCode VecSwap(Vec x,Vec y)
1605: {
1606: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1607: PetscBool flgxs[4],flgys[4];
1609: PetscInt i;
1617: VecCheckSameSize(x,1,y,2);
1618: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1619: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1621: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1622: for (i=0; i<4; i++) {
1623: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1624: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1625: }
1626: (*x->ops->swap)(x,y);
1627: PetscObjectStateIncrease((PetscObject)x);
1628: PetscObjectStateIncrease((PetscObject)y);
1629: for (i=0; i<4; i++) {
1630: if (flgxs[i]) {
1631: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1632: }
1633: if (flgys[i]) {
1634: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1635: }
1636: }
1637: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1638: return(0);
1639: }
1641: /*
1642: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1644: Collective on VecStash
1646: Input Parameters:
1647: + obj - the VecStash object
1648: . bobj - optional other object that provides the prefix
1649: - optionname - option to activate viewing
1651: Level: intermediate
1653: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1655: */
1656: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1657: {
1658: PetscErrorCode ierr;
1659: PetscViewer viewer;
1660: PetscBool flg;
1661: PetscViewerFormat format;
1662: char *prefix;
1665: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1666: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1667: if (flg) {
1668: PetscViewerPushFormat(viewer,format);
1669: VecStashView(obj,viewer);
1670: PetscViewerPopFormat(viewer);
1671: PetscViewerDestroy(&viewer);
1672: }
1673: return(0);
1674: }
1676: /*@
1677: VecStashView - Prints the entries in the vector stash and block stash.
1679: Collective on Vec
1681: Input Parameters:
1682: + v - the vector
1683: - viewer - the viewer
1685: Level: advanced
1688: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1690: @*/
1691: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1692: {
1694: PetscMPIInt rank;
1695: PetscInt i,j;
1696: PetscBool match;
1697: VecStash *s;
1698: PetscScalar val;
1705: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1706: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1707: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1708: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1709: s = &v->bstash;
1711: /* print block stash */
1712: PetscViewerASCIIPushSynchronized(viewer);
1713: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1714: for (i=0; i<s->n; i++) {
1715: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1716: for (j=0; j<s->bs; j++) {
1717: val = s->array[i*s->bs+j];
1718: #if defined(PETSC_USE_COMPLEX)
1719: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1720: #else
1721: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1722: #endif
1723: }
1724: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1725: }
1726: PetscViewerFlush(viewer);
1728: s = &v->stash;
1730: /* print basic stash */
1731: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1732: for (i=0; i<s->n; i++) {
1733: val = s->array[i];
1734: #if defined(PETSC_USE_COMPLEX)
1735: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1736: #else
1737: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1738: #endif
1739: }
1740: PetscViewerFlush(viewer);
1741: PetscViewerASCIIPopSynchronized(viewer);
1742: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1743: return(0);
1744: }
1746: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1747: {
1748: PetscInt i,N,rstart,rend;
1750: PetscScalar *xx;
1751: PetscReal *xreal;
1752: PetscBool iset;
1755: VecGetOwnershipRange(v,&rstart,&rend);
1756: VecGetSize(v,&N);
1757: PetscCalloc1(N,&xreal);
1758: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1759: if (iset) {
1760: VecGetArray(v,&xx);
1761: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1762: VecRestoreArray(v,&xx);
1763: }
1764: PetscFree(xreal);
1765: if (set) *set = iset;
1766: return(0);
1767: }
1769: /*@
1770: VecGetLayout - get PetscLayout describing vector layout
1772: Not Collective
1774: Input Arguments:
1775: . x - the vector
1777: Output Arguments:
1778: . map - the layout
1780: Level: developer
1782: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1783: @*/
1784: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1785: {
1789: *map = x->map;
1790: return(0);
1791: }
1793: /*@
1794: VecSetLayout - set PetscLayout describing vector layout
1796: Not Collective
1798: Input Arguments:
1799: + x - the vector
1800: - map - the layout
1802: Notes:
1803: It is normally only valid to replace the layout with a layout known to be equivalent.
1805: Level: developer
1807: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1808: @*/
1809: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1810: {
1815: PetscLayoutReference(map,&x->map);
1816: return(0);
1817: }
1819: PetscErrorCode VecSetInf(Vec xin)
1820: {
1821: PetscInt i,n = xin->map->n;
1822: PetscScalar *xx;
1823: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1827: VecGetArray(xin,&xx);
1828: for (i=0; i<n; i++) xx[i] = inf;
1829: VecRestoreArray(xin,&xx);
1830: return(0);
1831: }
1833: /*@
1834: VecPinToCPU - marks a vector to temporarily stay on the CPU and perform computations on the CPU
1836: Input Parameters:
1837: + v - the vector
1838: - flg - pin to the CPU if value of PETSC_TRUE
1840: Level: intermediate
1841: @*/
1842: PetscErrorCode VecPinToCPU(Vec v,PetscBool flg)
1843: {
1844: #if defined(PETSC_HAVE_VIENNACL) || defined(PETSC_HAVE_CUDA)
1848: if (v->pinnedtocpu == flg) return(0);
1849: v->pinnedtocpu = flg;
1850: if (v->ops->pintocpu) {
1851: (*v->ops->pintocpu)(v,flg);
1852: }
1853: return(0);
1854: #else
1855: return 0;
1856: #endif
1857: }