Actual source code: vector.c
2: /*
3: Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
4: These are the vector functions the user calls.
5: */
6: #include <private/vecimpl.h> /*I "petscvec.h" I*/
8: /* Logging support */
9: PetscClassId VEC_CLASSID;
10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_DotBarrier, VEC_Dot, VEC_MDotBarrier, VEC_MDot, VEC_TDot;
11: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
12: PetscLogEvent VEC_MTDot, VEC_NormBarrier, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load, VEC_ScatterBarrier;
14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceBarrier, VEC_ReduceCommunication,VEC_Ops;
15: PetscLogEvent VEC_DotNormBarrier, VEC_DotNorm, VEC_AXPBYPCZ, VEC_CUSPCopyFromGPU, VEC_CUSPCopyToGPU;
16: PetscLogEvent VEC_CUSPCopyFromGPUSome, VEC_CUSPCopyToGPUSome;
18: extern PetscErrorCode VecStashGetInfo_Private(VecStash*,PetscInt*,PetscInt*);
21: /*@
22: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
23: to be communicated to other processors during the VecAssemblyBegin/End() process
25: Not collective
27: Input Parameter:
28: . vec - the vector
30: Output Parameters:
31: + nstash - the size of the stash
32: . reallocs - the number of additional mallocs incurred.
33: . bnstash - the size of the block stash
34: - breallocs - the number of additional mallocs incurred.in the block stash
36: Level: advanced
38: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
40: @*/
41: PetscErrorCode VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
42: {
45: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
46: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
47: return(0);
48: }
52: /*@
53: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
54: by the routine VecSetValuesLocal() to allow users to insert vector entries
55: using a local (per-processor) numbering.
57: Logically Collective on Vec
59: Input Parameters:
60: + x - vector
61: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
63: Notes:
64: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
66: Level: intermediate
68: Concepts: vector^setting values with local numbering
70: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
71: VecSetLocalToGlobalMappingBlock(), VecSetValuesBlockedLocal()
72: @*/
73: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
74: {
81: if (x->ops->setlocaltoglobalmapping) {
82: (*x->ops->setlocaltoglobalmapping)(x,mapping);
83: } else {
84: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
85: }
86: return(0);
87: }
91: /*@
92: VecSetLocalToGlobalMappingBlock - Sets a local numbering to global numbering used
93: by the routine VecSetValuesBlockedLocal() to allow users to insert vector entries
94: using a local (per-processor) numbering.
96: Logically Collective on Vec
98: Input Parameters:
99: + x - vector
100: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
102: Notes:
103: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
105: Level: intermediate
107: Concepts: vector^setting values blocked with local numbering
109: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
110: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
111: @*/
112: PetscErrorCode VecSetLocalToGlobalMappingBlock(Vec x,ISLocalToGlobalMapping mapping)
113: {
120: PetscLayoutSetISLocalToGlobalMappingBlock(x->map,mapping);
121: return(0);
122: }
126: /*@
127: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
129: Not Collective
131: Input Parameter:
132: . X - the vector
134: Output Parameter:
135: . mapping - the mapping
137: Level: advanced
139: Concepts: vectors^local to global mapping
140: Concepts: local to global mapping^for vectors
142: .seealso: VecSetValuesLocal(), VecGetLocalToGlobalMappingBlock()
143: @*/
144: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
145: {
150: *mapping = X->map->mapping;
151: return(0);
152: }
156: /*@
157: VecGetLocalToGlobalMappingBlock - Gets the local-to-global numbering set by VecSetLocalToGlobalMappingBlock()
159: Not Collective
161: Input Parameters:
162: . X - the vector
164: Output Parameters:
165: . mapping - the mapping
167: Level: advanced
169: Concepts: vectors^local to global mapping blocked
170: Concepts: local to global mapping^for vectors, blocked
172: .seealso: VecSetValuesBlockedLocal(), VecGetLocalToGlobalMapping()
173: @*/
174: PetscErrorCode VecGetLocalToGlobalMappingBlock(Vec X,ISLocalToGlobalMapping *mapping)
175: {
180: *mapping = X->map->bmapping;
181: return(0);
182: }
186: /*@
187: VecAssemblyBegin - Begins assembling the vector. This routine should
188: be called after completing all calls to VecSetValues().
190: Collective on Vec
192: Input Parameter:
193: . vec - the vector
195: Level: beginner
197: Concepts: assembly^vectors
199: .seealso: VecAssemblyEnd(), VecSetValues()
200: @*/
201: PetscErrorCode VecAssemblyBegin(Vec vec)
202: {
204: PetscBool flg = PETSC_FALSE;
210: PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_stash",&flg,PETSC_NULL);
211: if (flg) {
212: PetscViewer viewer;
213: PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
214: VecStashView(vec,viewer);
215: }
217: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
218: if (vec->ops->assemblybegin) {
219: (*vec->ops->assemblybegin)(vec);
220: }
221: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
222: PetscObjectStateIncrease((PetscObject)vec);
223: return(0);
224: }
228: /*
229: Processes command line options to determine if/how a matrix
230: is to be viewed. Called by VecAssemblyEnd().
232: .seealso: MatView_Private()
234: */
235: PetscErrorCode VecView_Private(Vec vec)
236: {
238: PetscBool flg = PETSC_FALSE;
241: PetscObjectOptionsBegin((PetscObject)vec);
242: PetscOptionsBool("-vec_view","Print vector to stdout","VecView",flg,&flg,PETSC_NULL);
243: if (flg) {
244: PetscViewer viewer;
245: PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
246: VecView(vec,viewer);
247: }
248: flg = PETSC_FALSE;
249: PetscOptionsBool("-vec_view_matlab","Print vector to stdout in a format MATLAB can read","VecView",flg,&flg,PETSC_NULL);
250: if (flg) {
251: PetscViewer viewer;
252: PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
253: PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_MATLAB);
254: VecView(vec,viewer);
255: PetscViewerPopFormat(viewer);
256: }
257: #if defined(PETSC_HAVE_MATLAB_ENGINE)
258: flg = PETSC_FALSE;
259: PetscOptionsBool("-vec_view_matlab_file","Print vector to matlaboutput.mat format MATLAB can read","VecView",flg,&flg,PETSC_NULL);
260: if (flg) {
261: VecView(vec,PETSC_VIEWER_MATLAB_(((PetscObject)vec)->comm));
262: }
263: #endif
264: #if defined(PETSC_USE_SOCKET_VIEWER)
265: flg = PETSC_FALSE;
266: PetscOptionsBool("-vec_view_socket","Send vector to socket (can be read from matlab)","VecView",flg,&flg,PETSC_NULL);
267: if (flg) {
268: VecView(vec,PETSC_VIEWER_SOCKET_(((PetscObject)vec)->comm));
269: PetscViewerFlush(PETSC_VIEWER_SOCKET_(((PetscObject)vec)->comm));
270: }
271: #endif
272: flg = PETSC_FALSE;
273: PetscOptionsBool("-vec_view_binary","Save vector to file in binary format","VecView",flg,&flg,PETSC_NULL);
274: if (flg) {
275: VecView(vec,PETSC_VIEWER_BINARY_(((PetscObject)vec)->comm));
276: PetscViewerFlush(PETSC_VIEWER_BINARY_(((PetscObject)vec)->comm));
277: }
278: PetscOptionsEnd();
279: /* These invoke PetscDrawGetDraw which invokes PetscOptionsBegin/End, */
280: /* hence they should not be inside the above PetscOptionsBegin/End block. */
281: flg = PETSC_FALSE;
282: PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_draw",&flg,PETSC_NULL);
283: if (flg) {
284: VecView(vec,PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
285: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
286: }
287: flg = PETSC_FALSE;
288: PetscOptionsGetBool(((PetscObject)vec)->prefix,"-vec_view_draw_lg",&flg,PETSC_NULL);
289: if (flg) {
290: PetscViewerSetFormat(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm),PETSC_VIEWER_DRAW_LG);
291: VecView(vec,PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
292: PetscViewerFlush(PETSC_VIEWER_DRAW_(((PetscObject)vec)->comm));
293: }
294: return(0);
295: }
299: /*@
300: VecAssemblyEnd - Completes assembling the vector. This routine should
301: be called after VecAssemblyBegin().
303: Collective on Vec
305: Input Parameter:
306: . vec - the vector
308: Options Database Keys:
309: + -vec_view - Prints vector in ASCII format
310: . -vec_view_matlab - Prints vector in ASCII MATLAB format to stdout
311: . -vec_view_matlab_file - Prints vector in MATLAB format to matlaboutput.mat
312: . -vec_view_draw - Activates vector viewing using drawing tools
313: . -display <name> - Sets display name (default is host)
314: . -draw_pause <sec> - Sets number of seconds to pause after display
315: - -vec_view_socket - Activates vector viewing using a socket
317: Level: beginner
319: .seealso: VecAssemblyBegin(), VecSetValues()
320: @*/
321: PetscErrorCode VecAssemblyEnd(Vec vec)
322: {
327: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
329: if (vec->ops->assemblyend) {
330: (*vec->ops->assemblyend)(vec);
331: }
332: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
333: VecView_Private(vec);
334: return(0);
335: }
339: /*@
340: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
342: Logically Collective on Vec
344: Input Parameters:
345: . x, y - the vectors
347: Output Parameter:
348: . w - the result
350: Level: advanced
352: Notes: any subset of the x, y, and w may be the same vector.
353: For complex numbers compares only the real part
355: Concepts: vector^pointwise multiply
357: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
358: @*/
359: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
360: {
372: if (x->map->N != y->map->N || x->map->N != w->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
373: if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
375: (*w->ops->pointwisemax)(w,x,y);
376: PetscObjectStateIncrease((PetscObject)w);
377: return(0);
378: }
383: /*@
384: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
386: Logically Collective on Vec
388: Input Parameters:
389: . x, y - the vectors
391: Output Parameter:
392: . w - the result
394: Level: advanced
396: Notes: any subset of the x, y, and w may be the same vector.
397: For complex numbers compares only the real part
399: Concepts: vector^pointwise multiply
401: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
402: @*/
403: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
404: {
416: if (x->map->N != y->map->N || x->map->N != w->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
417: if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
419: (*w->ops->pointwisemin)(w,x,y);
420: PetscObjectStateIncrease((PetscObject)w);
421: return(0);
422: }
426: /*@
427: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
429: Logically Collective on Vec
431: Input Parameters:
432: . x, y - the vectors
434: Output Parameter:
435: . w - the result
437: Level: advanced
439: Notes: any subset of the x, y, and w may be the same vector.
441: Concepts: vector^pointwise multiply
443: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
444: @*/
445: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
446: {
458: if (x->map->N != y->map->N || x->map->N != w->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
459: if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
461: (*w->ops->pointwisemaxabs)(w,x,y);
462: PetscObjectStateIncrease((PetscObject)w);
463: return(0);
464: }
468: /*@
469: VecPointwiseDivide - Computes the componentwise division w = x/y.
471: Logically Collective on Vec
473: Input Parameters:
474: . x, y - the vectors
476: Output Parameter:
477: . w - the result
479: Level: advanced
481: Notes: any subset of the x, y, and w may be the same vector.
483: Concepts: vector^pointwise divide
485: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
486: @*/
487: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
488: {
500: if (x->map->N != y->map->N || x->map->N != w->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
501: if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
503: (*w->ops->pointwisedivide)(w,x,y);
504: PetscObjectStateIncrease((PetscObject)w);
505: return(0);
506: }
511: /*@
512: VecDuplicate - Creates a new vector of the same type as an existing vector.
514: Collective on Vec
516: Input Parameters:
517: . v - a vector to mimic
519: Output Parameter:
520: . newv - location to put new vector
522: Notes:
523: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
524: for the new vector. Use VecCopy() to copy a vector.
526: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
527: vectors.
529: Level: beginner
531: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
532: @*/
533: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
534: {
541: (*v->ops->duplicate)(v,newv);
542: PetscObjectStateIncrease((PetscObject)*newv);
543: return(0);
544: }
548: /*@
549: VecDestroy - Destroys a vector.
551: Collective on Vec
553: Input Parameters:
554: . v - the vector
556: Level: beginner
558: .seealso: VecDuplicate(), VecDestroyVecs()
559: @*/
560: PetscErrorCode VecDestroy(Vec *v)
561: {
565: if (!*v) return(0);
567: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
568: /* destroy the internal part */
569: if ((*v)->ops->destroy) {
570: (*(*v)->ops->destroy)(*v);
571: }
572: /* destroy the external/common part */
573: PetscLayoutDestroy(&(*v)->map);
574: PetscHeaderDestroy(v);
575: return(0);
576: }
580: /*@C
581: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
583: Collective on Vec
585: Input Parameters:
586: + m - the number of vectors to obtain
587: - v - a vector to mimic
589: Output Parameter:
590: . V - location to put pointer to array of vectors
592: Notes:
593: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
594: vector.
596: Fortran Note:
597: The Fortran interface is slightly different from that given below, it
598: requires one to pass in V a Vec (integer) array of size at least m.
599: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
601: Level: intermediate
603: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
604: @*/
605: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
606: {
613: (*v->ops->duplicatevecs)(v, m,V);
614: return(0);
615: }
619: /*@C
620: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
622: Collective on Vec
624: Input Parameters:
625: + vv - pointer to pointer to array of vector pointers
626: - m - the number of vectors previously obtained
628: Fortran Note:
629: The Fortran interface is slightly different from that given below.
630: See the Fortran chapter of the users manual
632: Level: intermediate
634: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
635: @*/
636: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
637: {
642: if (!*vv) return(0);
645: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
646: (*(**vv)->ops->destroyvecs)(m,*vv);
647: *vv = 0;
648: return(0);
649: }
651: #undef __FUNCT__
653: /*@
654: VecViewFromOptions - This function visualizes the vector based upon user options.
656: Collective on Vec
658: Input Parameters:
659: . vec - The vector
660: . title - The title (currently ignored)
662: Level: intermediate
664: .keywords: Vec, view, options, database
665: .seealso: VecSetFromOptions(), VecView()
666: @*/
667: PetscErrorCode VecViewFromOptions(Vec vec, const char *title)
668: {
672: VecView_Private(vec);
673: return(0);
674: }
678: /*@C
679: VecView - Views a vector object.
681: Collective on Vec
683: Input Parameters:
684: + vec - the vector
685: - viewer - an optional visualization context
687: Notes:
688: The available visualization contexts include
689: + PETSC_VIEWER_STDOUT_SELF - standard output (default)
690: - PETSC_VIEWER_STDOUT_WORLD - synchronized standard
691: output where only the first processor opens
692: the file. All other processors send their
693: data to the first processor to print.
695: You can change the format the vector is printed using the
696: option PetscViewerSetFormat().
698: The user can open alternative visualization contexts with
699: + PetscViewerASCIIOpen() - Outputs vector to a specified file
700: . PetscViewerBinaryOpen() - Outputs vector in binary to a
701: specified file; corresponding input uses VecLoad()
702: . PetscViewerDrawOpen() - Outputs vector to an X window display
703: - PetscViewerSocketOpen() - Outputs vector to Socket viewer
705: The user can call PetscViewerSetFormat() to specify the output
706: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
707: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
708: + PETSC_VIEWER_DEFAULT - default, prints vector contents
709: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
710: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
711: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
712: format common among all vector types
714: Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
715: for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
716: obviously) several times, you must change its name each time before calling the VecView(). The name you use
717: here should equal the name that you use in the Vec object that you use with VecLoad().
719: See the manual page for VecLoad() on the exact format the binary viewer stores
720: the values in the file.
722: Level: beginner
724: Concepts: vector^printing
725: Concepts: vector^saving to disk
727: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
728: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
729: PetscRealView(), PetscScalarView(), PetscIntView()
730: @*/
731: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
732: {
733: PetscErrorCode ierr;
738: if (!viewer) {
739: PetscViewerASCIIGetStdout(((PetscObject)vec)->comm,&viewer);
740: }
743: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
745: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
746: (*vec->ops->view)(vec,viewer);
747: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
748: return(0);
749: }
751: #if defined(PETSC_USE_DEBUG)
752: #include <../src/sys/totalview/tv_data_display.h>
753: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
754: {
755: const PetscScalar *values;
756: char type[32];
757: PetscErrorCode ierr;
760: TV_add_row("Local rows", "int", &v->map->n);
761: TV_add_row("Global rows", "int", &v->map->N);
762: TV_add_row("Typename", TV_ascii_string_type , ((PetscObject)v)->type_name);
763: VecGetArrayRead((Vec)v,&values);
764: PetscSNPrintf(type,32,"double[%d]",v->map->n);
765: TV_add_row("values",type, values);
766: VecRestoreArrayRead((Vec)v,&values);
767: return TV_format_OK;
768: }
769: #endif
773: /*@
774: VecGetSize - Returns the global number of elements of the vector.
776: Not Collective
778: Input Parameter:
779: . x - the vector
781: Output Parameters:
782: . size - the global length of the vector
784: Level: beginner
786: Concepts: vector^local size
788: .seealso: VecGetLocalSize()
789: @*/
790: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
791: {
798: (*x->ops->getsize)(x,size);
799: return(0);
800: }
804: /*@
805: VecGetLocalSize - Returns the number of elements of the vector stored
806: in local memory. This routine may be implementation dependent, so use
807: with care.
809: Not Collective
811: Input Parameter:
812: . x - the vector
814: Output Parameter:
815: . size - the length of the local piece of the vector
817: Level: beginner
819: Concepts: vector^size
821: .seealso: VecGetSize()
822: @*/
823: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
824: {
831: (*x->ops->getlocalsize)(x,size);
832: return(0);
833: }
837: /*@C
838: VecGetOwnershipRange - Returns the range of indices owned by
839: this processor, assuming that the vectors are laid out with the
840: first n1 elements on the first processor, next n2 elements on the
841: second, etc. For certain parallel layouts this range may not be
842: well defined.
844: Not Collective
846: Input Parameter:
847: . x - the vector
849: Output Parameters:
850: + low - the first local element, pass in PETSC_NULL if not interested
851: - high - one more than the last local element, pass in PETSC_NULL if not interested
853: Note:
854: The high argument is one more than the last element stored locally.
856: Fortran: PETSC_NULL_INTEGER should be used instead of PETSC_NULL
858: Level: beginner
860: Concepts: ownership^of vectors
861: Concepts: vector^ownership of elements
863: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
864: @*/
865: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
866: {
872: if (low) *low = x->map->rstart;
873: if (high) *high = x->map->rend;
874: return(0);
875: }
879: /*@C
880: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
881: assuming that the vectors are laid out with the
882: first n1 elements on the first processor, next n2 elements on the
883: second, etc. For certain parallel layouts this range may not be
884: well defined.
886: Not Collective
888: Input Parameter:
889: . x - the vector
891: Output Parameters:
892: . range - array of length size+1 with the start and end+1 for each process
894: Note:
895: The high argument is one more than the last element stored locally.
897: Fortran: You must PASS in an array of length size+1
899: Level: beginner
901: Concepts: ownership^of vectors
902: Concepts: vector^ownership of elements
904: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
905: @*/
906: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
907: {
913: PetscLayoutGetRanges(x->map,ranges);
914: return(0);
915: }
919: /*@
920: VecSetOption - Sets an option for controling a vector's behavior.
922: Collective on Vec
924: Input Parameter:
925: + x - the vector
926: . op - the option
927: - flag - turn the option on or off
929: Supported Options:
930: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
931: entries destined to be stored on a separate processor. This can be used
932: to eliminate the global reduction in the VecAssemblyXXXX() if you know
933: that you have only used VecSetValues() to set local elements
934: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
935: in ix in calls to VecSetValues or VecGetValues. These rows are simply
936: ignored.
938: Level: intermediate
940: @*/
941: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
942: {
948: if (x->ops->setoption) {
949: (*x->ops->setoption)(x,op,flag);
950: }
951: return(0);
952: }
956: /* Default routines for obtaining and releasing; */
957: /* may be used by any implementation */
958: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
959: {
961: PetscInt i;
966: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
967: PetscMalloc(m*sizeof(Vec*),V);
968: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
969: return(0);
970: }
974: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
975: {
977: PetscInt i;
981: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
982: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
983: PetscFree(v);
984: return(0);
985: }
989: /*@
990: VecResetArray - Resets a vector to use its default memory. Call this
991: after the use of VecPlaceArray().
993: Not Collective
995: Input Parameters:
996: . vec - the vector
998: Level: developer
1000: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
1002: @*/
1003: PetscErrorCode VecResetArray(Vec vec)
1004: {
1010: if (vec->ops->resetarray) {
1011: (*vec->ops->resetarray)(vec);
1012: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
1013: PetscObjectStateIncrease((PetscObject)vec);
1014: return(0);
1015: }
1019: /*@C
1020: VecLoad - Loads a vector that has been stored in binary or HDF5 format
1021: with VecView().
1023: Collective on PetscViewer
1025: Input Parameters:
1026: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
1027: some related function before a call to VecLoad().
1028: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
1029: HDF5 file viewer, obtained from PetscViewerHDF5Open()
1031: Level: intermediate
1033: Notes:
1034: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
1035: before calling this.
1037: The input file must contain the full global vector, as
1038: written by the routine VecView().
1040: If the type or size of newvec is not set before a call to VecLoad, PETSc
1041: sets the type and the local and global sizes.If type and/or
1042: sizes are already set, then the same are used.
1044: IF using HDF5, you must assign the Vec the same name as was used in the Vec
1045: that was stored in the file using PetscObjectSetName(). Otherwise you will
1046: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
1048: Notes for advanced users:
1049: Most users should not need to know the details of the binary storage
1050: format, since VecLoad() and VecView() completely hide these details.
1051: But for anyone who's interested, the standard binary matrix storage
1052: format is
1053: .vb
1054: int VEC_FILE_CLASSID
1055: int number of rows
1056: PetscScalar *values of all entries
1057: .ve
1059: In addition, PETSc automatically does the byte swapping for
1060: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
1061: linux, Windows and the paragon; thus if you write your own binary
1062: read/write routines you have to swap the bytes; see PetscBinaryRead()
1063: and PetscBinaryWrite() to see how this may be done.
1065: Concepts: vector^loading from file
1067: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
1068: @*/
1069: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
1070: {
1077: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
1078: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
1079: VecSetType(newvec, VECSTANDARD);
1080: }
1081: (*newvec->ops->load)(newvec,viewer);
1082: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
1083: return(0);
1084: }
1089: /*@
1090: VecReciprocal - Replaces each component of a vector by its reciprocal.
1092: Logically Collective on Vec
1094: Input Parameter:
1095: . vec - the vector
1097: Output Parameter:
1098: . vec - the vector reciprocal
1100: Level: intermediate
1102: Concepts: vector^reciprocal
1104: .seealso: VecLog(), VecExp(), VecSqrtAbs()
1106: @*/
1107: PetscErrorCode VecReciprocal(Vec vec)
1108: {
1114: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1115: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
1116: (*vec->ops->reciprocal)(vec);
1117: PetscObjectStateIncrease((PetscObject)vec);
1118: return(0);
1119: }
1123: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1124: {
1127: (((void(**)(void))vec->ops)[(int)op]) = f;
1128: return(0);
1129: }
1134: /*@
1135: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1136: used during the assembly process to store values that belong to
1137: other processors.
1139: Not Collective, different processes can have different size stashes
1141: Input Parameters:
1142: + vec - the vector
1143: . size - the initial size of the stash.
1144: - bsize - the initial size of the block-stash(if used).
1146: Options Database Keys:
1147: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1148: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1150: Level: intermediate
1152: Notes:
1153: The block-stash is used for values set with VecSetValuesBlocked() while
1154: the stash is used for values set with VecSetValues()
1156: Run with the option -info and look for output of the form
1157: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1158: to determine the appropriate value, MM, to use for size and
1159: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1160: to determine the value, BMM to use for bsize
1162: Concepts: vector^stash
1163: Concepts: stash^vector
1165: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1167: @*/
1168: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1169: {
1174: VecStashSetInitialSize_Private(&vec->stash,size);
1175: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1176: return(0);
1177: }
1181: /*@
1182: VecConjugate - Conjugates a vector.
1184: Logically Collective on Vec
1186: Input Parameters:
1187: . x - the vector
1189: Level: intermediate
1191: Concepts: vector^conjugate
1193: @*/
1194: PetscErrorCode VecConjugate(Vec x)
1195: {
1196: #ifdef PETSC_USE_COMPLEX
1202: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1203: (*x->ops->conjugate)(x);
1204: /* we need to copy norms here */
1205: PetscObjectStateIncrease((PetscObject)x);
1206: return(0);
1207: #else
1208: return(0);
1209: #endif
1210: }
1214: /*@
1215: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1217: Logically Collective on Vec
1219: Input Parameters:
1220: . x, y - the vectors
1222: Output Parameter:
1223: . w - the result
1225: Level: advanced
1227: Notes: any subset of the x, y, and w may be the same vector.
1229: Concepts: vector^pointwise multiply
1231: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1232: @*/
1233: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1234: {
1246: if (x->map->n != y->map->n || x->map->n != w->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1248: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1249: (*w->ops->pointwisemult)(w,x,y);
1250: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1251: PetscObjectStateIncrease((PetscObject)w);
1252: return(0);
1253: }
1257: /*@
1258: VecSetRandom - Sets all components of a vector to random numbers.
1260: Logically Collective on Vec
1262: Input Parameters:
1263: + x - the vector
1264: - rctx - the random number context, formed by PetscRandomCreate(), or PETSC_NULL and
1265: it will create one internally.
1267: Output Parameter:
1268: . x - the vector
1270: Example of Usage:
1271: .vb
1272: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1273: VecSetRandom(x,rctx);
1274: PetscRandomDestroy(rctx);
1275: .ve
1277: Level: intermediate
1279: Concepts: vector^setting to random
1280: Concepts: random^vector
1282: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1283: @*/
1284: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1285: {
1287: PetscRandom randObj = PETSC_NULL;
1293: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1295: if (!rctx) {
1296: MPI_Comm comm;
1297: PetscObjectGetComm((PetscObject)x,&comm);
1298: PetscRandomCreate(comm,&randObj);
1299: PetscRandomSetFromOptions(randObj);
1300: rctx = randObj;
1301: }
1303: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1304: (*x->ops->setrandom)(x,rctx);
1305: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1307: PetscRandomDestroy(&randObj);
1308: PetscObjectStateIncrease((PetscObject)x);
1309: return(0);
1310: }
1314: /*@
1315: VecZeroEntries - puts a 0.0 in each element of a vector
1317: Logically Collective on Vec
1319: Input Parameter:
1320: . vec - The vector
1322: Level: beginner
1324: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1325: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1326: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1327: this routine should not exist.
1329: .keywords: Vec, set, options, database
1330: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1331: @*/
1332: PetscErrorCode VecZeroEntries(Vec vec)
1333: {
1336: VecSet(vec,0);
1337: return(0);
1338: }
1342: /*
1343: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1344: processor and a PETSc MPI vector on more than one processor.
1346: Collective on Vec
1348: Input Parameter:
1349: . vec - The vector
1351: Level: intermediate
1353: .keywords: Vec, set, options, database, type
1354: .seealso: VecSetFromOptions(), VecSetType()
1355: */
1356: static PetscErrorCode VecSetTypeFromOptions_Private(Vec vec)
1357: {
1358: PetscBool opt;
1359: const VecType defaultType;
1360: char typeName[256];
1361: PetscMPIInt size;
1365: if (((PetscObject)vec)->type_name) {
1366: defaultType = ((PetscObject)vec)->type_name;
1367: } else {
1368: MPI_Comm_size(((PetscObject)vec)->comm, &size);
1369: if (size > 1) {
1370: defaultType = VECMPI;
1371: } else {
1372: defaultType = VECSEQ;
1373: }
1374: }
1376: if (!VecRegisterAllCalled) {VecRegisterAll(PETSC_NULL);}
1377: PetscOptionsList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1378: if (opt) {
1379: VecSetType(vec, typeName);
1380: } else {
1381: VecSetType(vec, defaultType);
1382: }
1383: return(0);
1384: }
1388: /*@
1389: VecSetFromOptions - Configures the vector from the options database.
1391: Collective on Vec
1393: Input Parameter:
1394: . vec - The vector
1396: Notes: To see all options, run your program with the -help option, or consult the users manual.
1397: Must be called after VecCreate() but before the vector is used.
1399: Level: beginner
1401: Concepts: vectors^setting options
1402: Concepts: vectors^setting type
1404: .keywords: Vec, set, options, database
1405: .seealso: VecCreate(), VecSetOptionsPrefix()
1406: @*/
1407: PetscErrorCode VecSetFromOptions(Vec vec)
1408: {
1414: PetscObjectOptionsBegin((PetscObject)vec);
1415: /* Handle vector type options */
1416: VecSetTypeFromOptions_Private(vec);
1418: /* Handle specific vector options */
1419: if (vec->ops->setfromoptions) {
1420: (*vec->ops->setfromoptions)(vec);
1421: }
1423: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1424: PetscObjectProcessOptionsHandlers((PetscObject)vec);
1425: PetscOptionsEnd();
1427: VecViewFromOptions(vec, ((PetscObject)vec)->name);
1428: return(0);
1429: }
1433: /*@
1434: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1436: Collective on Vec
1438: Input Parameters:
1439: + v - the vector
1440: . n - the local size (or PETSC_DECIDE to have it set)
1441: - N - the global size (or PETSC_DECIDE)
1443: Notes:
1444: n and N cannot be both PETSC_DECIDE
1445: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1447: Level: intermediate
1449: .seealso: VecGetSize(), PetscSplitOwnership()
1450: @*/
1451: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1452: {
1457: 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);
1458: 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);
1459: v->map->n = n;
1460: v->map->N = N;
1461: if (v->ops->create) {
1462: (*v->ops->create)(v);
1463: v->ops->create = 0;
1464: }
1465: return(0);
1466: }
1470: /*@
1471: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1472: and VecSetValuesBlockedLocal().
1474: Logically Collective on Vec
1476: Input Parameter:
1477: + v - the vector
1478: - bs - the blocksize
1480: Notes:
1481: All vectors obtained by VecDuplicate() inherit the same blocksize.
1483: Level: advanced
1485: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlock(), VecGetBlockSize()
1487: Concepts: block size^vectors
1488: @*/
1489: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1490: {
1493: if (bs <= 0) bs = 1;
1494: if (bs == v->map->bs) return(0);
1495: if (v->map->N != -1 && v->map->N % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Vector length not divisible by blocksize %D %D",v->map->N,bs);
1496: if (v->map->n != -1 && v->map->n % bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Local vector length not divisible by blocksize %D %D\n\
1497: Try setting blocksize before setting the vector type",v->map->n,bs);
1500: v->map->bs = bs;
1501: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1502: return(0);
1503: }
1507: /*@
1508: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1509: and VecSetValuesBlockedLocal().
1511: Not Collective
1513: Input Parameter:
1514: . v - the vector
1516: Output Parameter:
1517: . bs - the blocksize
1519: Notes:
1520: All vectors obtained by VecDuplicate() inherit the same blocksize.
1522: Level: advanced
1524: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMappingBlock(), VecSetBlockSize()
1526: Concepts: vector^block size
1527: Concepts: block^vector
1529: @*/
1530: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1531: {
1535: *bs = v->map->bs;
1536: return(0);
1537: }
1541: /*@C
1542: VecSetOptionsPrefix - Sets the prefix used for searching for all
1543: Vec options in the database.
1545: Logically Collective on Vec
1547: Input Parameter:
1548: + v - the Vec context
1549: - prefix - the prefix to prepend to all option names
1551: Notes:
1552: A hyphen (-) must NOT be given at the beginning of the prefix name.
1553: The first character of all runtime options is AUTOMATICALLY the hyphen.
1555: Level: advanced
1557: .keywords: Vec, set, options, prefix, database
1559: .seealso: VecSetFromOptions()
1560: @*/
1561: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1562: {
1567: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1568: return(0);
1569: }
1573: /*@C
1574: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1575: Vec options in the database.
1577: Logically Collective on Vec
1579: Input Parameters:
1580: + v - the Vec context
1581: - prefix - the prefix to prepend to all option names
1583: Notes:
1584: A hyphen (-) must NOT be given at the beginning of the prefix name.
1585: The first character of all runtime options is AUTOMATICALLY the hyphen.
1587: Level: advanced
1589: .keywords: Vec, append, options, prefix, database
1591: .seealso: VecGetOptionsPrefix()
1592: @*/
1593: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1594: {
1599: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1600: return(0);
1601: }
1605: /*@C
1606: VecGetOptionsPrefix - Sets the prefix used for searching for all
1607: Vec options in the database.
1609: Not Collective
1611: Input Parameter:
1612: . v - the Vec context
1614: Output Parameter:
1615: . prefix - pointer to the prefix string used
1617: Notes: On the fortran side, the user should pass in a string 'prefix' of
1618: sufficient length to hold the prefix.
1620: Level: advanced
1622: .keywords: Vec, get, options, prefix, database
1624: .seealso: VecAppendOptionsPrefix()
1625: @*/
1626: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1627: {
1632: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1633: return(0);
1634: }
1638: /*@
1639: VecSetUp - Sets up the internal vector data structures for the later use.
1641: Collective on Vec
1643: Input Parameters:
1644: . v - the Vec context
1646: Notes:
1647: For basic use of the Vec classes the user need not explicitly call
1648: VecSetUp(), since these actions will happen automatically.
1650: Level: advanced
1652: .keywords: Vec, setup
1654: .seealso: VecCreate(), VecDestroy()
1655: @*/
1656: PetscErrorCode VecSetUp(Vec v)
1657: {
1658: PetscMPIInt size;
1663: if (!((PetscObject)v)->type_name) {
1664: MPI_Comm_size(((PetscObject)v)->comm, &size);
1665: if (size == 1) {
1666: VecSetType(v, VECSEQ);
1667: } else {
1668: VecSetType(v, VECMPI);
1669: }
1670: }
1671: return(0);
1672: }
1674: /*
1675: These currently expose the PetscScalar/PetscReal in updating the
1676: cached norm. If we push those down into the implementation these
1677: will become independent of PetscScalar/PetscReal
1678: */
1682: /*@
1683: VecCopy - Copies a vector. y <- x
1685: Logically Collective on Vec
1687: Input Parameter:
1688: . x - the vector
1690: Output Parameter:
1691: . y - the copy
1693: Notes:
1694: For default parallel PETSc vectors, both x and y must be distributed in
1695: the same manner; local copies are done.
1697: Level: beginner
1699: .seealso: VecDuplicate()
1700: @*/
1701: PetscErrorCode VecCopy(Vec x,Vec y)
1702: {
1703: PetscBool flgs[4];
1704: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1706: PetscInt i;
1713: if (x == y) return(0);
1714: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1715: if (x->map->n != y->map->n) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths %d != %d", x->map->n, y->map->n);
1717: #if !defined(PETSC_USE_MIXED_PRECISION)
1718: for (i=0; i<4; i++) {
1719: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1720: }
1721: #endif
1723: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1724: #if defined(PETSC_USE_MIXED_PRECISION)
1725: extern PetscErrorCode VecGetArray(Vec,double **);
1726: extern PetscErrorCode VecRestoreArray(Vec,double **);
1727: extern PetscErrorCode VecGetArray(Vec,float **);
1728: extern PetscErrorCode VecRestoreArray(Vec,float **);
1729: extern PetscErrorCode VecGetArrayRead(Vec,const double **);
1730: extern PetscErrorCode VecRestoreArrayRead(Vec,const double **);
1731: extern PetscErrorCode VecGetArrayRead(Vec,const float **);
1732: extern PetscErrorCode VecRestoreArrayRead(Vec,const float **);
1733: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1734: PetscInt i,n;
1735: const float *xx;
1736: double *yy;
1737: VecGetArrayRead(x,&xx);
1738: VecGetArray(y,&yy);
1739: VecGetLocalSize(x,&n);
1740: for (i=0; i<n; i++) {
1741: yy[i] = xx[i];
1742: }
1743: VecRestoreArrayRead(x,&xx);
1744: VecRestoreArray(y,&yy);
1745: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1746: PetscInt i,n;
1747: float *yy;
1748: const double *xx;
1749: VecGetArrayRead(x,&xx);
1750: VecGetArray(y,&yy);
1751: VecGetLocalSize(x,&n);
1752: for (i=0; i<n; i++) {
1753: yy[i] = (float) xx[i];
1754: }
1755: VecRestoreArrayRead(x,&xx);
1756: VecRestoreArray(y,&yy);
1757: } else {
1758: (*x->ops->copy)(x,y);
1759: }
1760: #else
1761: (*x->ops->copy)(x,y);
1762: #endif
1764: PetscObjectStateIncrease((PetscObject)y);
1765: #if !defined(PETSC_USE_MIXED_PRECISION)
1766: for (i=0; i<4; i++) {
1767: if (flgs[i]) {
1768: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1769: }
1770: }
1771: #endif
1773: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1774: return(0);
1775: }
1779: /*@
1780: VecSwap - Swaps the vectors x and y.
1782: Logically Collective on Vec
1784: Input Parameters:
1785: . x, y - the vectors
1787: Level: advanced
1789: Concepts: vector^swapping values
1791: @*/
1792: PetscErrorCode VecSwap(Vec x,Vec y)
1793: {
1794: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1795: PetscBool flgxs[4],flgys[4];
1797: PetscInt i;
1805: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1806: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1807: if (x->map->N != y->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1808: if (x->map->n != y->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1810: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1811: for (i=0; i<4; i++) {
1812: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1813: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1814: }
1815: (*x->ops->swap)(x,y);
1816: PetscObjectStateIncrease((PetscObject)x);
1817: PetscObjectStateIncrease((PetscObject)y);
1818: for (i=0; i<4; i++) {
1819: if (flgxs[i]) {
1820: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1821: }
1822: if (flgys[i]) {
1823: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1824: }
1825: }
1826: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1827: return(0);
1828: }
1832: /*@
1833: VecStashView - Prints the entries in the vector stash and block stash.
1835: Collective on Vec
1837: Input Parameters:
1838: + v - the vector
1839: - viewer - the viewer
1841: Level: advanced
1843: Concepts: vector^stash
1844: Concepts: stash^vector
1846: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1848: @*/
1849: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1850: {
1852: PetscMPIInt rank;
1853: PetscInt i,j;
1854: PetscBool match;
1855: VecStash *s;
1856: PetscScalar val;
1863: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1864: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1865: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1866: MPI_Comm_rank(((PetscObject)v)->comm,&rank);
1867: s = &v->bstash;
1869: /* print block stash */
1870: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1871: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1872: for (i=0; i<s->n; i++) {
1873: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1874: for (j=0; j<s->bs; j++) {
1875: val = s->array[i*s->bs+j];
1876: #if defined(PETSC_USE_COMPLEX)
1877: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1878: #else
1879: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1880: #endif
1881: }
1882: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1883: }
1884: PetscViewerFlush(viewer);
1886: s = &v->stash;
1888: /* print basic stash */
1889: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1890: for (i=0; i<s->n; i++) {
1891: val = s->array[i];
1892: #if defined(PETSC_USE_COMPLEX)
1893: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1894: #else
1895: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1896: #endif
1897: }
1898: PetscViewerFlush(viewer);
1899: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1901: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1902: return(0);
1903: }