Actual source code: vector.c
petsc-3.6.4 2016-04-12
2: /*
3: Provides the interface functions for vector operations that do NOT have PetscScalar/PetscReal in the signature
4: These are the vector functions the user calls.
5: */
6: #include <petsc/private/vecimpl.h> /*I "petscvec.h" I*/
8: /* Logging support */
9: PetscClassId VEC_CLASSID;
10: PetscLogEvent VEC_View, VEC_Max, VEC_Min, VEC_DotBarrier, VEC_Dot, VEC_MDotBarrier, VEC_MDot, VEC_TDot;
11: PetscLogEvent VEC_Norm, VEC_Normalize, VEC_Scale, VEC_Copy, VEC_Set, VEC_AXPY, VEC_AYPX, VEC_WAXPY;
12: PetscLogEvent VEC_MTDot, VEC_NormBarrier, VEC_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load, VEC_ScatterBarrier;
14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceBarrier, VEC_ReduceCommunication,VEC_ReduceBegin,VEC_ReduceEnd,VEC_Ops;
15: PetscLogEvent VEC_DotNormBarrier, VEC_DotNorm, VEC_AXPBYPCZ, VEC_CUSPCopyFromGPU, VEC_CUSPCopyToGPU;
16: PetscLogEvent VEC_CUSPCopyFromGPUSome, VEC_CUSPCopyToGPUSome;
17: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
19: extern PetscErrorCode VecStashGetInfo_Private(VecStash*,PetscInt*,PetscInt*);
22: /*@
23: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
24: to be communicated to other processors during the VecAssemblyBegin/End() process
26: Not collective
28: Input Parameter:
29: . vec - the vector
31: Output Parameters:
32: + nstash - the size of the stash
33: . reallocs - the number of additional mallocs incurred.
34: . bnstash - the size of the block stash
35: - breallocs - the number of additional mallocs incurred.in the block stash
37: Level: advanced
39: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
41: @*/
42: PetscErrorCode VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
43: {
47: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
48: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
49: return(0);
50: }
54: /*@
55: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
56: by the routine VecSetValuesLocal() to allow users to insert vector entries
57: using a local (per-processor) numbering.
59: Logically Collective on Vec
61: Input Parameters:
62: + x - vector
63: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
65: Notes:
66: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
68: Level: intermediate
70: Concepts: vector^setting values with local numbering
72: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
73: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
74: @*/
75: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
76: {
83: if (x->ops->setlocaltoglobalmapping) {
84: (*x->ops->setlocaltoglobalmapping)(x,mapping);
85: } else {
86: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
87: }
88: return(0);
89: }
93: /*@
94: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
96: Not Collective
98: Input Parameter:
99: . X - the vector
101: Output Parameter:
102: . mapping - the mapping
104: Level: advanced
106: Concepts: vectors^local to global mapping
107: Concepts: local to global mapping^for vectors
109: .seealso: VecSetValuesLocal()
110: @*/
111: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
112: {
117: *mapping = X->map->mapping;
118: return(0);
119: }
123: /*@
124: VecAssemblyBegin - Begins assembling the vector. This routine should
125: be called after completing all calls to VecSetValues().
127: Collective on Vec
129: Input Parameter:
130: . vec - the vector
132: Level: beginner
134: Concepts: assembly^vectors
136: .seealso: VecAssemblyEnd(), VecSetValues()
137: @*/
138: PetscErrorCode VecAssemblyBegin(Vec vec)
139: {
145: VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
146: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
147: if (vec->ops->assemblybegin) {
148: (*vec->ops->assemblybegin)(vec);
149: }
150: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
151: PetscObjectStateIncrease((PetscObject)vec);
152: return(0);
153: }
157: /*@
158: VecAssemblyEnd - Completes assembling the vector. This routine should
159: be called after VecAssemblyBegin().
161: Collective on Vec
163: Input Parameter:
164: . vec - the vector
166: Options Database Keys:
167: + -vec_view - Prints vector in ASCII format
168: . -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
169: . -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
170: . -vec_view draw - Activates vector viewing using drawing tools
171: . -display <name> - Sets display name (default is host)
172: . -draw_pause <sec> - Sets number of seconds to pause after display
173: - -vec_view socket - Activates vector viewing using a socket
175: Level: beginner
177: .seealso: VecAssemblyBegin(), VecSetValues()
178: @*/
179: PetscErrorCode VecAssemblyEnd(Vec vec)
180: {
185: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
187: if (vec->ops->assemblyend) {
188: (*vec->ops->assemblyend)(vec);
189: }
190: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
191: VecViewFromOptions(vec,NULL,"-vec_view");
192: return(0);
193: }
197: /*@
198: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
200: Logically Collective on Vec
202: Input Parameters:
203: . x, y - the vectors
205: Output Parameter:
206: . w - the result
208: Level: advanced
210: Notes: any subset of the x, y, and w may be the same vector.
211: For complex numbers compares only the real part
213: Concepts: vector^pointwise multiply
215: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
216: @*/
217: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
218: {
230: 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");
231: 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");
233: (*w->ops->pointwisemax)(w,x,y);
234: PetscObjectStateIncrease((PetscObject)w);
235: return(0);
236: }
241: /*@
242: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
244: Logically Collective on Vec
246: Input Parameters:
247: . x, y - the vectors
249: Output Parameter:
250: . w - the result
252: Level: advanced
254: Notes: any subset of the x, y, and w may be the same vector.
255: For complex numbers compares only the real part
257: Concepts: vector^pointwise multiply
259: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
260: @*/
261: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
262: {
274: 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");
275: 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");
277: (*w->ops->pointwisemin)(w,x,y);
278: PetscObjectStateIncrease((PetscObject)w);
279: return(0);
280: }
284: /*@
285: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
287: Logically Collective on Vec
289: Input Parameters:
290: . x, y - the vectors
292: Output Parameter:
293: . w - the result
295: Level: advanced
297: Notes: any subset of the x, y, and w may be the same vector.
299: Concepts: vector^pointwise multiply
301: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
302: @*/
303: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
304: {
316: 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");
317: 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");
319: (*w->ops->pointwisemaxabs)(w,x,y);
320: PetscObjectStateIncrease((PetscObject)w);
321: return(0);
322: }
326: /*@
327: VecPointwiseDivide - Computes the componentwise division w = x/y.
329: Logically Collective on Vec
331: Input Parameters:
332: . x, y - the vectors
334: Output Parameter:
335: . w - the result
337: Level: advanced
339: Notes: any subset of the x, y, and w may be the same vector.
341: Concepts: vector^pointwise divide
343: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
344: @*/
345: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
346: {
358: 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");
359: 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");
361: (*w->ops->pointwisedivide)(w,x,y);
362: PetscObjectStateIncrease((PetscObject)w);
363: return(0);
364: }
369: /*@
370: VecDuplicate - Creates a new vector of the same type as an existing vector.
372: Collective on Vec
374: Input Parameters:
375: . v - a vector to mimic
377: Output Parameter:
378: . newv - location to put new vector
380: Notes:
381: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
382: for the new vector. Use VecCopy() to copy a vector.
384: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
385: vectors.
387: Level: beginner
389: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
390: @*/
391: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
392: {
399: (*v->ops->duplicate)(v,newv);
400: PetscObjectStateIncrease((PetscObject)*newv);
401: return(0);
402: }
406: /*@
407: VecDestroy - Destroys a vector.
409: Collective on Vec
411: Input Parameters:
412: . v - the vector
414: Level: beginner
416: .seealso: VecDuplicate(), VecDestroyVecs()
417: @*/
418: PetscErrorCode VecDestroy(Vec *v)
419: {
423: if (!*v) return(0);
425: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
427: PetscObjectSAWsViewOff((PetscObject)*v);
428: /* destroy the internal part */
429: if ((*v)->ops->destroy) {
430: (*(*v)->ops->destroy)(*v);
431: }
432: /* destroy the external/common part */
433: PetscLayoutDestroy(&(*v)->map);
434: PetscHeaderDestroy(v);
435: return(0);
436: }
440: /*@C
441: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
443: Collective on Vec
445: Input Parameters:
446: + m - the number of vectors to obtain
447: - v - a vector to mimic
449: Output Parameter:
450: . V - location to put pointer to array of vectors
452: Notes:
453: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
454: vector.
456: Fortran Note:
457: The Fortran interface is slightly different from that given below, it
458: requires one to pass in V a Vec (integer) array of size at least m.
459: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
461: Level: intermediate
463: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
464: @*/
465: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
466: {
473: (*v->ops->duplicatevecs)(v, m,V);
474: return(0);
475: }
479: /*@C
480: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
482: Collective on Vec
484: Input Parameters:
485: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
486: - m - the number of vectors previously obtained, if zero no vectors are destroyed
488: Fortran Note:
489: The Fortran interface is slightly different from that given below.
490: See the Fortran chapter of the users manual
492: Level: intermediate
494: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
495: @*/
496: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
497: {
502: if (!*vv) return(0);
505: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
506: if (!m) return(0);
507: if (!*vv) return(0);
508: (*(**vv)->ops->destroyvecs)(m,*vv);
509: *vv = 0;
510: return(0);
511: }
515: /*@C
516: VecView - Views a vector object.
518: Collective on Vec
520: Input Parameters:
521: + vec - the vector
522: - viewer - an optional visualization context
524: Notes:
525: The available visualization contexts include
526: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
527: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
528: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
530: You can change the format the vector is printed using the
531: option PetscViewerSetFormat().
533: The user can open alternative visualization contexts with
534: + PetscViewerASCIIOpen() - Outputs vector to a specified file
535: . PetscViewerBinaryOpen() - Outputs vector in binary to a
536: specified file; corresponding input uses VecLoad()
537: . PetscViewerDrawOpen() - Outputs vector to an X window display
538: - PetscViewerSocketOpen() - Outputs vector to Socket viewer
540: The user can call PetscViewerSetFormat() to specify the output
541: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
542: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
543: + PETSC_VIEWER_DEFAULT - default, prints vector contents
544: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
545: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
546: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
547: format common among all vector types
549: Notes: You can pass any number of vector objects, or other PETSc objects to the same viewer.
551: Notes for binary viewer: If you pass multiply vectors to a binary viewer you can read them back in in the same order
552: $ with VecLoad().
553: $
554: $ If the blocksize of the vector is greater than one then you must provide a unique prefix to
555: $ the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
556: $ vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
557: $ information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
558: $ filename. If you copy the binary file, make sure you copy the associated .info file with it.
560: Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
561: $ for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
562: $ obviously) several times, you must change its name each time before calling the VecView(). The name you use
563: $ here should equal the name that you use in the Vec object that you use with VecLoad().
565: See the manual page for VecLoad() on the exact format the binary viewer stores
566: the values in the file.
568: Level: beginner
570: Concepts: vector^printing
571: Concepts: vector^saving to disk
573: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
574: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
575: PetscRealView(), PetscScalarView(), PetscIntView()
576: @*/
577: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
578: {
579: PetscErrorCode ierr;
580: PetscBool iascii;
581: PetscViewerFormat format;
586: if (!viewer) {
587: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
588: }
591: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
593: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
594: PetscViewerGetFormat(viewer,&format);
595: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
596: if (iascii) {
597: PetscInt rows,bs;
599: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
600: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
601: PetscViewerASCIIPushTab(viewer);
602: VecGetSize(vec,&rows);
603: VecGetBlockSize(vec,&bs);
604: if (bs != 1) {
605: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
606: } else {
607: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
608: }
609: PetscViewerASCIIPopTab(viewer);
610: }
611: }
612: VecLockPush(vec);
613: if (format == PETSC_VIEWER_NATIVE && vec->ops->viewnative) {
614: (*vec->ops->viewnative)(vec,viewer);
615: } else {
616: (*vec->ops->view)(vec,viewer);
617: }
618: VecLockPop(vec);
619: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
620: return(0);
621: }
623: #if defined(PETSC_USE_DEBUG)
624: #include <../src/sys/totalview/tv_data_display.h>
625: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
626: {
627: const PetscScalar *values;
628: char type[32];
629: PetscErrorCode ierr;
632: TV_add_row("Local rows", "int", &v->map->n);
633: TV_add_row("Global rows", "int", &v->map->N);
634: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
635: VecGetArrayRead((Vec)v,&values);
636: PetscSNPrintf(type,32,"double[%d]",v->map->n);
637: TV_add_row("values",type, values);
638: VecRestoreArrayRead((Vec)v,&values);
639: return TV_format_OK;
640: }
641: #endif
645: /*@
646: VecGetSize - Returns the global number of elements of the vector.
648: Not Collective
650: Input Parameter:
651: . x - the vector
653: Output Parameters:
654: . size - the global length of the vector
656: Level: beginner
658: Concepts: vector^local size
660: .seealso: VecGetLocalSize()
661: @*/
662: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
663: {
670: (*x->ops->getsize)(x,size);
671: return(0);
672: }
676: /*@
677: VecGetLocalSize - Returns the number of elements of the vector stored
678: in local memory. This routine may be implementation dependent, so use
679: with care.
681: Not Collective
683: Input Parameter:
684: . x - the vector
686: Output Parameter:
687: . size - the length of the local piece of the vector
689: Level: beginner
691: Concepts: vector^size
693: .seealso: VecGetSize()
694: @*/
695: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
696: {
703: (*x->ops->getlocalsize)(x,size);
704: return(0);
705: }
709: /*@C
710: VecGetOwnershipRange - Returns the range of indices owned by
711: this processor, assuming that the vectors are laid out with the
712: first n1 elements on the first processor, next n2 elements on the
713: second, etc. For certain parallel layouts this range may not be
714: well defined.
716: Not Collective
718: Input Parameter:
719: . x - the vector
721: Output Parameters:
722: + low - the first local element, pass in NULL if not interested
723: - high - one more than the last local element, pass in NULL if not interested
725: Note:
726: The high argument is one more than the last element stored locally.
728: Fortran: NULL_INTEGER should be used instead of NULL
730: Level: beginner
732: Concepts: ownership^of vectors
733: Concepts: vector^ownership of elements
735: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
736: @*/
737: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
738: {
744: if (low) *low = x->map->rstart;
745: if (high) *high = x->map->rend;
746: return(0);
747: }
751: /*@C
752: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
753: assuming that the vectors are laid out with the
754: first n1 elements on the first processor, next n2 elements on the
755: second, etc. For certain parallel layouts this range may not be
756: well defined.
758: Not Collective
760: Input Parameter:
761: . x - the vector
763: Output Parameters:
764: . range - array of length size+1 with the start and end+1 for each process
766: Note:
767: The high argument is one more than the last element stored locally.
769: Fortran: You must PASS in an array of length size+1
771: Level: beginner
773: Concepts: ownership^of vectors
774: Concepts: vector^ownership of elements
776: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
777: @*/
778: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
779: {
785: PetscLayoutGetRanges(x->map,ranges);
786: return(0);
787: }
791: /*@
792: VecSetOption - Sets an option for controling a vector's behavior.
794: Collective on Vec
796: Input Parameter:
797: + x - the vector
798: . op - the option
799: - flag - turn the option on or off
801: Supported Options:
802: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
803: entries destined to be stored on a separate processor. This can be used
804: to eliminate the global reduction in the VecAssemblyXXXX() if you know
805: that you have only used VecSetValues() to set local elements
806: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
807: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
808: ignored.
810: Level: intermediate
812: @*/
813: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
814: {
820: if (x->ops->setoption) {
821: (*x->ops->setoption)(x,op,flag);
822: }
823: return(0);
824: }
828: /* Default routines for obtaining and releasing; */
829: /* may be used by any implementation */
830: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
831: {
833: PetscInt i;
838: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
839: PetscMalloc1(m,V);
840: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
841: return(0);
842: }
846: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
847: {
849: PetscInt i;
853: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
854: PetscFree(v);
855: return(0);
856: }
860: /*@
861: VecResetArray - Resets a vector to use its default memory. Call this
862: after the use of VecPlaceArray().
864: Not Collective
866: Input Parameters:
867: . vec - the vector
869: Level: developer
871: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
873: @*/
874: PetscErrorCode VecResetArray(Vec vec)
875: {
881: if (vec->ops->resetarray) {
882: (*vec->ops->resetarray)(vec);
883: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
884: PetscObjectStateIncrease((PetscObject)vec);
885: return(0);
886: }
890: /*@C
891: VecLoad - Loads a vector that has been stored in binary or HDF5 format
892: with VecView().
894: Collective on PetscViewer
896: Input Parameters:
897: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
898: some related function before a call to VecLoad().
899: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
900: HDF5 file viewer, obtained from PetscViewerHDF5Open()
902: Level: intermediate
904: Notes:
905: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
906: before calling this.
908: The input file must contain the full global vector, as
909: written by the routine VecView().
911: If the type or size of newvec is not set before a call to VecLoad, PETSc
912: sets the type and the local and global sizes. If type and/or
913: sizes are already set, then the same are used.
915: If using binary and the blocksize of the vector is greater than one then you must provide a unique prefix to
916: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
917: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
918: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
919: filename. If you copy the binary file, make sure you copy the associated .info file with it.
921: If using HDF5, you must assign the Vec the same name as was used in the Vec
922: that was stored in the file using PetscObjectSetName(). Otherwise you will
923: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
925: Notes for advanced users:
926: Most users should not need to know the details of the binary storage
927: format, since VecLoad() and VecView() completely hide these details.
928: But for anyone who's interested, the standard binary matrix storage
929: format is
930: .vb
931: int VEC_FILE_CLASSID
932: int number of rows
933: PetscScalar *values of all entries
934: .ve
936: In addition, PETSc automatically does the byte swapping for
937: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
938: linux, Windows and the paragon; thus if you write your own binary
939: read/write routines you have to swap the bytes; see PetscBinaryRead()
940: and PetscBinaryWrite() to see how this may be done.
942: Concepts: vector^loading from file
944: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
945: @*/
946: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
947: {
949: PetscBool isbinary,ishdf5;
950: PetscViewerFormat format;
955: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
956: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
957: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
959: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
960: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
961: VecSetType(newvec, VECSTANDARD);
962: }
963: PetscViewerGetFormat(viewer,&format);
964: if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
965: (*newvec->ops->loadnative)(newvec,viewer);
966: } else {
967: (*newvec->ops->load)(newvec,viewer);
968: }
969: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
970: return(0);
971: }
976: /*@
977: VecReciprocal - Replaces each component of a vector by its reciprocal.
979: Logically Collective on Vec
981: Input Parameter:
982: . vec - the vector
984: Output Parameter:
985: . vec - the vector reciprocal
987: Level: intermediate
989: Concepts: vector^reciprocal
991: .seealso: VecLog(), VecExp(), VecSqrtAbs()
993: @*/
994: PetscErrorCode VecReciprocal(Vec vec)
995: {
1001: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1002: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
1003: (*vec->ops->reciprocal)(vec);
1004: PetscObjectStateIncrease((PetscObject)vec);
1005: return(0);
1006: }
1010: /*@C
1011: VecSetOperation - Allows user to set a vector operation.
1013: Logically Collective on Vec
1015: Input Parameters:
1016: + vec - the vector
1017: . op - the name of the operation
1018: - f - the function that provides the operation.
1020: Level: advanced
1022: Usage:
1023: $ PetscErrorCode userview(Vec,PetscViewer);
1024: $ VecCreateMPI(comm,m,M,&x);
1025: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
1027: Notes:
1028: See the file include/petscvec.h for a complete list of matrix
1029: operations, which all have the form VECOP_<OPERATION>, where
1030: <OPERATION> is the name (in all capital letters) of the
1031: user interface routine (e.g., VecView() -> VECOP_VIEW).
1033: This function is not currently available from Fortran.
1035: .keywords: vector, set, operation
1037: .seealso: VecCreate(), MatShellSetOperation()
1038: @*/
1039: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1040: {
1043: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1044: vec->ops->viewnative = vec->ops->view;
1045: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1046: vec->ops->loadnative = vec->ops->load;
1047: }
1048: (((void(**)(void))vec->ops)[(int)op]) = f;
1049: return(0);
1050: }
1055: /*@
1056: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1057: used during the assembly process to store values that belong to
1058: other processors.
1060: Not Collective, different processes can have different size stashes
1062: Input Parameters:
1063: + vec - the vector
1064: . size - the initial size of the stash.
1065: - bsize - the initial size of the block-stash(if used).
1067: Options Database Keys:
1068: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1069: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1071: Level: intermediate
1073: Notes:
1074: The block-stash is used for values set with VecSetValuesBlocked() while
1075: the stash is used for values set with VecSetValues()
1077: Run with the option -info and look for output of the form
1078: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1079: to determine the appropriate value, MM, to use for size and
1080: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1081: to determine the value, BMM to use for bsize
1083: Concepts: vector^stash
1084: Concepts: stash^vector
1086: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1088: @*/
1089: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1090: {
1095: VecStashSetInitialSize_Private(&vec->stash,size);
1096: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1097: return(0);
1098: }
1102: /*@
1103: VecConjugate - Conjugates a vector.
1105: Logically Collective on Vec
1107: Input Parameters:
1108: . x - the vector
1110: Level: intermediate
1112: Concepts: vector^conjugate
1114: @*/
1115: PetscErrorCode VecConjugate(Vec x)
1116: {
1117: #if defined(PETSC_USE_COMPLEX)
1123: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1124: (*x->ops->conjugate)(x);
1125: /* we need to copy norms here */
1126: PetscObjectStateIncrease((PetscObject)x);
1127: return(0);
1128: #else
1129: return(0);
1130: #endif
1131: }
1135: /*@
1136: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1138: Logically Collective on Vec
1140: Input Parameters:
1141: . x, y - the vectors
1143: Output Parameter:
1144: . w - the result
1146: Level: advanced
1148: Notes: any subset of the x, y, and w may be the same vector.
1150: Concepts: vector^pointwise multiply
1152: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1153: @*/
1154: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1155: {
1167: 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");
1169: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1170: (*w->ops->pointwisemult)(w,x,y);
1171: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1172: PetscObjectStateIncrease((PetscObject)w);
1173: return(0);
1174: }
1178: /*@
1179: VecSetRandom - Sets all components of a vector to random numbers.
1181: Logically Collective on Vec
1183: Input Parameters:
1184: + x - the vector
1185: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1186: it will create one internally.
1188: Output Parameter:
1189: . x - the vector
1191: Example of Usage:
1192: .vb
1193: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1194: VecSetRandom(x,rctx);
1195: PetscRandomDestroy(rctx);
1196: .ve
1198: Level: intermediate
1200: Concepts: vector^setting to random
1201: Concepts: random^vector
1203: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1204: @*/
1205: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1206: {
1208: PetscRandom randObj = NULL;
1214: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1216: if (!rctx) {
1217: MPI_Comm comm;
1218: PetscObjectGetComm((PetscObject)x,&comm);
1219: PetscRandomCreate(comm,&randObj);
1220: PetscRandomSetFromOptions(randObj);
1221: rctx = randObj;
1222: }
1224: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1225: (*x->ops->setrandom)(x,rctx);
1226: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1228: PetscRandomDestroy(&randObj);
1229: PetscObjectStateIncrease((PetscObject)x);
1230: return(0);
1231: }
1235: /*@
1236: VecZeroEntries - puts a 0.0 in each element of a vector
1238: Logically Collective on Vec
1240: Input Parameter:
1241: . vec - The vector
1243: Level: beginner
1245: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1246: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1247: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1248: this routine should not exist.
1250: .keywords: Vec, set, options, database
1251: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1252: @*/
1253: PetscErrorCode VecZeroEntries(Vec vec)
1254: {
1258: VecSet(vec,0);
1259: return(0);
1260: }
1264: /*
1265: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1266: processor and a PETSc MPI vector on more than one processor.
1268: Collective on Vec
1270: Input Parameter:
1271: . vec - The vector
1273: Level: intermediate
1275: .keywords: Vec, set, options, database, type
1276: .seealso: VecSetFromOptions(), VecSetType()
1277: */
1278: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptions *PetscOptionsObject,Vec vec)
1279: {
1280: PetscBool opt;
1281: VecType defaultType;
1282: char typeName[256];
1283: PetscMPIInt size;
1287: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1288: else {
1289: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1290: if (size > 1) defaultType = VECMPI;
1291: else defaultType = VECSEQ;
1292: }
1294: VecRegisterAll();
1295: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1296: if (opt) {
1297: VecSetType(vec, typeName);
1298: } else {
1299: VecSetType(vec, defaultType);
1300: }
1301: return(0);
1302: }
1306: /*@
1307: VecSetFromOptions - Configures the vector from the options database.
1309: Collective on Vec
1311: Input Parameter:
1312: . vec - The vector
1314: Notes: To see all options, run your program with the -help option, or consult the users manual.
1315: Must be called after VecCreate() but before the vector is used.
1317: Level: beginner
1319: Concepts: vectors^setting options
1320: Concepts: vectors^setting type
1322: .keywords: Vec, set, options, database
1323: .seealso: VecCreate(), VecSetOptionsPrefix()
1324: @*/
1325: PetscErrorCode VecSetFromOptions(Vec vec)
1326: {
1332: PetscObjectOptionsBegin((PetscObject)vec);
1333: /* Handle vector type options */
1334: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1336: /* Handle specific vector options */
1337: if (vec->ops->setfromoptions) {
1338: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1339: }
1341: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1342: PetscObjectProcessOptionsHandlers((PetscObject)vec);
1343: PetscOptionsEnd();
1344: return(0);
1345: }
1349: /*@
1350: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1352: Collective on Vec
1354: Input Parameters:
1355: + v - the vector
1356: . n - the local size (or PETSC_DECIDE to have it set)
1357: - N - the global size (or PETSC_DECIDE)
1359: Notes:
1360: n and N cannot be both PETSC_DECIDE
1361: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1363: Level: intermediate
1365: .seealso: VecGetSize(), PetscSplitOwnership()
1366: @*/
1367: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1368: {
1374: 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);
1375: 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);
1376: v->map->n = n;
1377: v->map->N = N;
1378: if (v->ops->create) {
1379: (*v->ops->create)(v);
1380: v->ops->create = 0;
1381: }
1382: return(0);
1383: }
1387: /*@
1388: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1389: and VecSetValuesBlockedLocal().
1391: Logically Collective on Vec
1393: Input Parameter:
1394: + v - the vector
1395: - bs - the blocksize
1397: Notes:
1398: All vectors obtained by VecDuplicate() inherit the same blocksize.
1400: Level: advanced
1402: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1404: Concepts: block size^vectors
1405: @*/
1406: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1407: {
1412: if (bs < 0 || bs == v->map->bs) return(0);
1414: PetscLayoutSetBlockSize(v->map,bs);
1415: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1416: return(0);
1417: }
1421: /*@
1422: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1423: and VecSetValuesBlockedLocal().
1425: Not Collective
1427: Input Parameter:
1428: . v - the vector
1430: Output Parameter:
1431: . bs - the blocksize
1433: Notes:
1434: All vectors obtained by VecDuplicate() inherit the same blocksize.
1436: Level: advanced
1438: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1440: Concepts: vector^block size
1441: Concepts: block^vector
1443: @*/
1444: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1445: {
1451: PetscLayoutGetBlockSize(v->map,bs);
1452: return(0);
1453: }
1457: /*@C
1458: VecSetOptionsPrefix - Sets the prefix used for searching for all
1459: Vec options in the database.
1461: Logically Collective on Vec
1463: Input Parameter:
1464: + v - the Vec context
1465: - prefix - the prefix to prepend to all option names
1467: Notes:
1468: A hyphen (-) must NOT be given at the beginning of the prefix name.
1469: The first character of all runtime options is AUTOMATICALLY the hyphen.
1471: Level: advanced
1473: .keywords: Vec, set, options, prefix, database
1475: .seealso: VecSetFromOptions()
1476: @*/
1477: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1478: {
1483: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1484: return(0);
1485: }
1489: /*@C
1490: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1491: Vec options in the database.
1493: Logically Collective on Vec
1495: Input Parameters:
1496: + v - the Vec context
1497: - prefix - the prefix to prepend to all option names
1499: Notes:
1500: A hyphen (-) must NOT be given at the beginning of the prefix name.
1501: The first character of all runtime options is AUTOMATICALLY the hyphen.
1503: Level: advanced
1505: .keywords: Vec, append, options, prefix, database
1507: .seealso: VecGetOptionsPrefix()
1508: @*/
1509: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1510: {
1515: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1516: return(0);
1517: }
1521: /*@C
1522: VecGetOptionsPrefix - Sets the prefix used for searching for all
1523: Vec options in the database.
1525: Not Collective
1527: Input Parameter:
1528: . v - the Vec context
1530: Output Parameter:
1531: . prefix - pointer to the prefix string used
1533: Notes: On the fortran side, the user should pass in a string 'prefix' of
1534: sufficient length to hold the prefix.
1536: Level: advanced
1538: .keywords: Vec, get, options, prefix, database
1540: .seealso: VecAppendOptionsPrefix()
1541: @*/
1542: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1543: {
1548: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1549: return(0);
1550: }
1554: /*@
1555: VecSetUp - Sets up the internal vector data structures for the later use.
1557: Collective on Vec
1559: Input Parameters:
1560: . v - the Vec context
1562: Notes:
1563: For basic use of the Vec classes the user need not explicitly call
1564: VecSetUp(), since these actions will happen automatically.
1566: Level: advanced
1568: .keywords: Vec, setup
1570: .seealso: VecCreate(), VecDestroy()
1571: @*/
1572: PetscErrorCode VecSetUp(Vec v)
1573: {
1574: PetscMPIInt size;
1579: if (!((PetscObject)v)->type_name) {
1580: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1581: if (size == 1) {
1582: VecSetType(v, VECSEQ);
1583: } else {
1584: VecSetType(v, VECMPI);
1585: }
1586: }
1587: return(0);
1588: }
1590: /*
1591: These currently expose the PetscScalar/PetscReal in updating the
1592: cached norm. If we push those down into the implementation these
1593: will become independent of PetscScalar/PetscReal
1594: */
1598: /*@
1599: VecCopy - Copies a vector. y <- x
1601: Logically Collective on Vec
1603: Input Parameter:
1604: . x - the vector
1606: Output Parameter:
1607: . y - the copy
1609: Notes:
1610: For default parallel PETSc vectors, both x and y must be distributed in
1611: the same manner; local copies are done.
1613: Level: beginner
1615: .seealso: VecDuplicate()
1616: @*/
1617: PetscErrorCode VecCopy(Vec x,Vec y)
1618: {
1619: PetscBool flgs[4];
1620: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1622: PetscInt i;
1629: if (x == y) return(0);
1630: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1631: 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);
1632: VecLocked(y,2);
1634: #if !defined(PETSC_USE_MIXED_PRECISION)
1635: for (i=0; i<4; i++) {
1636: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1637: }
1638: #endif
1640: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1641: #if defined(PETSC_USE_MIXED_PRECISION)
1642: extern PetscErrorCode VecGetArray(Vec,double**);
1643: extern PetscErrorCode VecRestoreArray(Vec,double**);
1644: extern PetscErrorCode VecGetArray(Vec,float**);
1645: extern PetscErrorCode VecRestoreArray(Vec,float**);
1646: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1647: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1648: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1649: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1650: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1651: PetscInt i,n;
1652: const float *xx;
1653: double *yy;
1654: VecGetArrayRead(x,&xx);
1655: VecGetArray(y,&yy);
1656: VecGetLocalSize(x,&n);
1657: for (i=0; i<n; i++) yy[i] = xx[i];
1658: VecRestoreArrayRead(x,&xx);
1659: VecRestoreArray(y,&yy);
1660: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1661: PetscInt i,n;
1662: float *yy;
1663: const double *xx;
1664: VecGetArrayRead(x,&xx);
1665: VecGetArray(y,&yy);
1666: VecGetLocalSize(x,&n);
1667: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1668: VecRestoreArrayRead(x,&xx);
1669: VecRestoreArray(y,&yy);
1670: } else {
1671: (*x->ops->copy)(x,y);
1672: }
1673: #else
1674: (*x->ops->copy)(x,y);
1675: #endif
1677: PetscObjectStateIncrease((PetscObject)y);
1678: #if !defined(PETSC_USE_MIXED_PRECISION)
1679: for (i=0; i<4; i++) {
1680: if (flgs[i]) {
1681: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1682: }
1683: }
1684: #endif
1686: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1687: return(0);
1688: }
1692: /*@
1693: VecSwap - Swaps the vectors x and y.
1695: Logically Collective on Vec
1697: Input Parameters:
1698: . x, y - the vectors
1700: Level: advanced
1702: Concepts: vector^swapping values
1704: @*/
1705: PetscErrorCode VecSwap(Vec x,Vec y)
1706: {
1707: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1708: PetscBool flgxs[4],flgys[4];
1710: PetscInt i;
1718: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1719: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1720: if (x->map->N != y->map->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector global lengths");
1721: if (x->map->n != y->map->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Incompatible vector local lengths");
1723: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1724: for (i=0; i<4; i++) {
1725: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1726: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1727: }
1728: (*x->ops->swap)(x,y);
1729: PetscObjectStateIncrease((PetscObject)x);
1730: PetscObjectStateIncrease((PetscObject)y);
1731: for (i=0; i<4; i++) {
1732: if (flgxs[i]) {
1733: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1734: }
1735: if (flgys[i]) {
1736: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1737: }
1738: }
1739: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1740: return(0);
1741: }
1745: /*
1746: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1748: Collective on VecStash
1750: Input Parameters:
1751: + obj - the VecStash object
1752: . bobj - optional other object that provides the prefix
1753: - optionname - option to activate viewing
1755: Level: intermediate
1757: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1759: */
1760: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1761: {
1762: PetscErrorCode ierr;
1763: PetscViewer viewer;
1764: PetscBool flg;
1765: PetscViewerFormat format;
1766: char *prefix;
1769: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1770: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),prefix,optionname,&viewer,&format,&flg);
1771: if (flg) {
1772: PetscViewerPushFormat(viewer,format);
1773: VecStashView(obj,viewer);
1774: PetscViewerPopFormat(viewer);
1775: PetscViewerDestroy(&viewer);
1776: }
1777: return(0);
1778: }
1782: /*@
1783: VecStashView - Prints the entries in the vector stash and block stash.
1785: Collective on Vec
1787: Input Parameters:
1788: + v - the vector
1789: - viewer - the viewer
1791: Level: advanced
1793: Concepts: vector^stash
1794: Concepts: stash^vector
1796: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1798: @*/
1799: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1800: {
1802: PetscMPIInt rank;
1803: PetscInt i,j;
1804: PetscBool match;
1805: VecStash *s;
1806: PetscScalar val;
1813: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1814: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1815: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1816: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1817: s = &v->bstash;
1819: /* print block stash */
1820: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
1821: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1822: for (i=0; i<s->n; i++) {
1823: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1824: for (j=0; j<s->bs; j++) {
1825: val = s->array[i*s->bs+j];
1826: #if defined(PETSC_USE_COMPLEX)
1827: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1828: #else
1829: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1830: #endif
1831: }
1832: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1833: }
1834: PetscViewerFlush(viewer);
1836: s = &v->stash;
1838: /* print basic stash */
1839: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1840: for (i=0; i<s->n; i++) {
1841: val = s->array[i];
1842: #if defined(PETSC_USE_COMPLEX)
1843: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1844: #else
1845: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1846: #endif
1847: }
1848: PetscViewerFlush(viewer);
1849: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
1851: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1852: return(0);
1853: }
1857: PetscErrorCode PetscOptionsGetVec(const char prefix[],const char key[],Vec v,PetscBool *set)
1858: {
1859: PetscInt i,N,rstart,rend;
1861: PetscScalar *xx;
1862: PetscReal *xreal;
1863: PetscBool iset;
1866: VecGetOwnershipRange(v,&rstart,&rend);
1867: VecGetSize(v,&N);
1868: PetscCalloc1(N,&xreal);
1869: PetscOptionsGetRealArray(prefix,key,xreal,&N,&iset);
1870: if (iset) {
1871: VecGetArray(v,&xx);
1872: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1873: VecRestoreArray(v,&xx);
1874: }
1875: PetscFree(xreal);
1876: if (set) *set = iset;
1877: return(0);
1878: }
1882: /*@
1883: VecGetLayout - get PetscLayout describing vector layout
1885: Not Collective
1887: Input Arguments:
1888: . x - the vector
1890: Output Arguments:
1891: . map - the layout
1893: Level: developer
1895: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1896: @*/
1897: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1898: {
1902: *map = x->map;
1903: return(0);
1904: }
1908: /*@
1909: VecSetLayout - set PetscLayout describing vector layout
1911: Not Collective
1913: Input Arguments:
1914: + x - the vector
1915: - map - the layout
1917: Notes:
1918: It is normally only valid to replace the layout with a layout known to be equivalent.
1920: Level: developer
1922: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1923: @*/
1924: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1925: {
1930: PetscLayoutReference(map,&x->map);
1931: return(0);
1932: }
1936: PetscErrorCode VecSetInf(Vec xin)
1937: {
1938: PetscInt i,n = xin->map->n;
1939: PetscScalar *xx;
1940: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1944: VecGetArray(xin,&xx);
1945: for (i=0; i<n; i++) xx[i] = inf;
1946: VecRestoreArray(xin,&xx);
1947: return(0);
1948: }