Actual source code: vector.c
petsc-3.9.4 2018-09-11
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>
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;
16: PetscLogEvent VEC_ViennaCLCopyFromGPU, VEC_ViennaCLCopyToGPU;
17: PetscLogEvent VEC_CUDACopyFromGPU, VEC_CUDACopyToGPU;
18: PetscLogEvent VEC_CUDACopyFromGPUSome, VEC_CUDACopyToGPUSome;
20: /*@
21: VecStashGetInfo - Gets how many values are currently in the vector stash, i.e. need
22: to be communicated to other processors during the VecAssemblyBegin/End() process
24: Not collective
26: Input Parameter:
27: . vec - the vector
29: Output Parameters:
30: + nstash - the size of the stash
31: . reallocs - the number of additional mallocs incurred.
32: . bnstash - the size of the block stash
33: - breallocs - the number of additional mallocs incurred.in the block stash
35: Level: advanced
37: .seealso: VecAssemblyBegin(), VecAssemblyEnd(), Vec, VecStashSetInitialSize(), VecStashView()
39: @*/
40: PetscErrorCode VecStashGetInfo(Vec vec,PetscInt *nstash,PetscInt *reallocs,PetscInt *bnstash,PetscInt *breallocs)
41: {
45: VecStashGetInfo_Private(&vec->stash,nstash,reallocs);
46: VecStashGetInfo_Private(&vec->bstash,bnstash,breallocs);
47: return(0);
48: }
50: /*@
51: VecSetLocalToGlobalMapping - Sets a local numbering to global numbering used
52: by the routine VecSetValuesLocal() to allow users to insert vector entries
53: using a local (per-processor) numbering.
55: Logically Collective on Vec
57: Input Parameters:
58: + x - vector
59: - mapping - mapping created with ISLocalToGlobalMappingCreate() or ISLocalToGlobalMappingCreateIS()
61: Notes:
62: All vectors obtained with VecDuplicate() from this vector inherit the same mapping.
64: Level: intermediate
66: Concepts: vector^setting values with local numbering
68: seealso: VecAssemblyBegin(), VecAssemblyEnd(), VecSetValues(), VecSetValuesLocal(),
69: VecSetLocalToGlobalMapping(), VecSetValuesBlockedLocal()
70: @*/
71: PetscErrorCode VecSetLocalToGlobalMapping(Vec x,ISLocalToGlobalMapping mapping)
72: {
79: if (x->ops->setlocaltoglobalmapping) {
80: (*x->ops->setlocaltoglobalmapping)(x,mapping);
81: } else {
82: PetscLayoutSetISLocalToGlobalMapping(x->map,mapping);
83: }
84: return(0);
85: }
87: /*@
88: VecGetLocalToGlobalMapping - Gets the local-to-global numbering set by VecSetLocalToGlobalMapping()
90: Not Collective
92: Input Parameter:
93: . X - the vector
95: Output Parameter:
96: . mapping - the mapping
98: Level: advanced
100: Concepts: vectors^local to global mapping
101: Concepts: local to global mapping^for vectors
103: .seealso: VecSetValuesLocal()
104: @*/
105: PetscErrorCode VecGetLocalToGlobalMapping(Vec X,ISLocalToGlobalMapping *mapping)
106: {
111: *mapping = X->map->mapping;
112: return(0);
113: }
115: /*@
116: VecAssemblyBegin - Begins assembling the vector. This routine should
117: be called after completing all calls to VecSetValues().
119: Collective on Vec
121: Input Parameter:
122: . vec - the vector
124: Level: beginner
126: Concepts: assembly^vectors
128: .seealso: VecAssemblyEnd(), VecSetValues()
129: @*/
130: PetscErrorCode VecAssemblyBegin(Vec vec)
131: {
137: VecStashViewFromOptions(vec,NULL,"-vec_view_stash");
138: PetscLogEventBegin(VEC_AssemblyBegin,vec,0,0,0);
139: if (vec->ops->assemblybegin) {
140: (*vec->ops->assemblybegin)(vec);
141: }
142: PetscLogEventEnd(VEC_AssemblyBegin,vec,0,0,0);
143: PetscObjectStateIncrease((PetscObject)vec);
144: return(0);
145: }
147: /*@
148: VecAssemblyEnd - Completes assembling the vector. This routine should
149: be called after VecAssemblyBegin().
151: Collective on Vec
153: Input Parameter:
154: . vec - the vector
156: Options Database Keys:
157: + -vec_view - Prints vector in ASCII format
158: . -vec_view ::ascii_matlab - Prints vector in ASCII MATLAB format to stdout
159: . -vec_view matlab:filename - Prints vector in MATLAB format to matlaboutput.mat
160: . -vec_view draw - Activates vector viewing using drawing tools
161: . -display <name> - Sets display name (default is host)
162: . -draw_pause <sec> - Sets number of seconds to pause after display
163: - -vec_view socket - Activates vector viewing using a socket
165: Level: beginner
167: .seealso: VecAssemblyBegin(), VecSetValues()
168: @*/
169: PetscErrorCode VecAssemblyEnd(Vec vec)
170: {
175: PetscLogEventBegin(VEC_AssemblyEnd,vec,0,0,0);
177: if (vec->ops->assemblyend) {
178: (*vec->ops->assemblyend)(vec);
179: }
180: PetscLogEventEnd(VEC_AssemblyEnd,vec,0,0,0);
181: VecViewFromOptions(vec,NULL,"-vec_view");
182: return(0);
183: }
185: /*@
186: VecPointwiseMax - Computes the componentwise maximum w_i = max(x_i, y_i).
188: Logically Collective on Vec
190: Input Parameters:
191: . x, y - the vectors
193: Output Parameter:
194: . w - the result
196: Level: advanced
198: Notes: any subset of the x, y, and w may be the same vector.
199: For complex numbers compares only the real part
201: Concepts: vector^pointwise multiply
203: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
204: @*/
205: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
206: {
218: VecCheckSameSize(w,1,x,2);
219: VecCheckSameSize(w,1,y,3);
220: (*w->ops->pointwisemax)(w,x,y);
221: PetscObjectStateIncrease((PetscObject)w);
222: return(0);
223: }
226: /*@
227: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
229: Logically Collective on Vec
231: Input Parameters:
232: . x, y - the vectors
234: Output Parameter:
235: . w - the result
237: Level: advanced
239: Notes: any subset of the x, y, and w may be the same vector.
240: For complex numbers compares only the real part
242: Concepts: vector^pointwise multiply
244: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
245: @*/
246: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
247: {
259: VecCheckSameSize(w,1,x,2);
260: VecCheckSameSize(w,1,y,3);
261: (*w->ops->pointwisemin)(w,x,y);
262: PetscObjectStateIncrease((PetscObject)w);
263: return(0);
264: }
266: /*@
267: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
269: Logically Collective on Vec
271: Input Parameters:
272: . x, y - the vectors
274: Output Parameter:
275: . w - the result
277: Level: advanced
279: Notes: any subset of the x, y, and w may be the same vector.
281: Concepts: vector^pointwise multiply
283: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
284: @*/
285: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
286: {
298: VecCheckSameSize(w,1,x,2);
299: VecCheckSameSize(w,1,y,3);
300: (*w->ops->pointwisemaxabs)(w,x,y);
301: PetscObjectStateIncrease((PetscObject)w);
302: return(0);
303: }
305: /*@
306: VecPointwiseDivide - Computes the componentwise division w = x/y.
308: Logically Collective on Vec
310: Input Parameters:
311: . x, y - the vectors
313: Output Parameter:
314: . w - the result
316: Level: advanced
318: Notes: any subset of the x, y, and w may be the same vector.
320: Concepts: vector^pointwise divide
322: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
323: @*/
324: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
325: {
337: VecCheckSameSize(w,1,x,2);
338: VecCheckSameSize(w,1,y,3);
339: (*w->ops->pointwisedivide)(w,x,y);
340: PetscObjectStateIncrease((PetscObject)w);
341: return(0);
342: }
345: /*@
346: VecDuplicate - Creates a new vector of the same type as an existing vector.
348: Collective on Vec
350: Input Parameters:
351: . v - a vector to mimic
353: Output Parameter:
354: . newv - location to put new vector
356: Notes:
357: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
358: for the new vector. Use VecCopy() to copy a vector.
360: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
361: vectors.
363: Level: beginner
365: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
366: @*/
367: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
368: {
375: (*v->ops->duplicate)(v,newv);
376: PetscObjectStateIncrease((PetscObject)*newv);
377: return(0);
378: }
380: /*@
381: VecDestroy - Destroys a vector.
383: Collective on Vec
385: Input Parameters:
386: . v - the vector
388: Level: beginner
390: .seealso: VecDuplicate(), VecDestroyVecs()
391: @*/
392: PetscErrorCode VecDestroy(Vec *v)
393: {
397: if (!*v) return(0);
399: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
401: PetscObjectSAWsViewOff((PetscObject)*v);
402: /* destroy the internal part */
403: if ((*v)->ops->destroy) {
404: (*(*v)->ops->destroy)(*v);
405: }
406: /* destroy the external/common part */
407: PetscLayoutDestroy(&(*v)->map);
408: PetscHeaderDestroy(v);
409: return(0);
410: }
412: /*@C
413: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
415: Collective on Vec
417: Input Parameters:
418: + m - the number of vectors to obtain
419: - v - a vector to mimic
421: Output Parameter:
422: . V - location to put pointer to array of vectors
424: Notes:
425: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
426: vector.
428: Fortran Note:
429: The Fortran interface is slightly different from that given below, it
430: requires one to pass in V a Vec (integer) array of size at least m.
431: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
433: Level: intermediate
435: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
436: @*/
437: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
438: {
445: (*v->ops->duplicatevecs)(v,m,V);
446: return(0);
447: }
449: /*@C
450: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
452: Collective on Vec
454: Input Parameters:
455: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
456: - m - the number of vectors previously obtained, if zero no vectors are destroyed
458: Fortran Note:
459: The Fortran interface is slightly different from that given below.
460: See the Fortran chapter of the users manual
462: Level: intermediate
464: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
465: @*/
466: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
467: {
472: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
473: if (!m || !*vv) {*vv = NULL; return(0);}
476: (*(**vv)->ops->destroyvecs)(m,*vv);
477: *vv = NULL;
478: return(0);
479: }
481: /*@C
482: VecView - Views a vector object.
484: Collective on Vec
486: Input Parameters:
487: + vec - the vector
488: - viewer - an optional visualization context
490: Notes:
491: The available visualization contexts include
492: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
493: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
494: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
496: You can change the format the vector is printed using the
497: option PetscViewerPushFormat().
499: The user can open alternative visualization contexts with
500: + PetscViewerASCIIOpen() - Outputs vector to a specified file
501: . PetscViewerBinaryOpen() - Outputs vector in binary to a
502: specified file; corresponding input uses VecLoad()
503: . PetscViewerDrawOpen() - Outputs vector to an X window display
504: - PetscViewerSocketOpen() - Outputs vector to Socket viewer
506: The user can call PetscViewerPushFormat() to specify the output
507: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
508: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
509: + PETSC_VIEWER_DEFAULT - default, prints vector contents
510: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
511: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
512: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
513: format common among all vector types
515: Notes: You can pass any number of vector objects, or other PETSc objects to the same viewer.
517: Notes for binary viewer: If you pass multiply vectors to a binary viewer you can read them back in in the same order
518: $ with VecLoad().
519: $
520: $ If the blocksize of the vector is greater than one then you must provide a unique prefix to
521: $ the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
522: $ vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
523: $ information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
524: $ filename. If you copy the binary file, make sure you copy the associated .info file with it.
526: Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
527: $ for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
528: $ obviously) several times, you must change its name each time before calling the VecView(). The name you use
529: $ here should equal the name that you use in the Vec object that you use with VecLoad().
531: See the manual page for VecLoad() on the exact format the binary viewer stores
532: the values in the file.
534: Level: beginner
536: Concepts: vector^printing
537: Concepts: vector^saving to disk
539: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
540: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
541: PetscRealView(), PetscScalarView(), PetscIntView()
542: @*/
543: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
544: {
545: PetscErrorCode ierr;
546: PetscBool iascii;
547: PetscViewerFormat format;
548: PetscMPIInt size;
553: if (!viewer) {
554: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
555: }
557: PetscViewerGetFormat(viewer,&format);
558: MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
559: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
561: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
563: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
564: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
565: if (iascii) {
566: PetscInt rows,bs;
568: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
569: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
570: PetscViewerASCIIPushTab(viewer);
571: VecGetSize(vec,&rows);
572: VecGetBlockSize(vec,&bs);
573: if (bs != 1) {
574: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
575: } else {
576: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
577: }
578: PetscViewerASCIIPopTab(viewer);
579: }
580: }
581: VecLockPush(vec);
582: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
583: (*vec->ops->viewnative)(vec,viewer);
584: } else {
585: (*vec->ops->view)(vec,viewer);
586: }
587: VecLockPop(vec);
588: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
589: return(0);
590: }
592: #if defined(PETSC_USE_DEBUG)
593: #include <../src/sys/totalview/tv_data_display.h>
594: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
595: {
596: const PetscScalar *values;
597: char type[32];
598: PetscErrorCode ierr;
601: TV_add_row("Local rows", "int", &v->map->n);
602: TV_add_row("Global rows", "int", &v->map->N);
603: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
604: VecGetArrayRead((Vec)v,&values);
605: PetscSNPrintf(type,32,"double[%d]",v->map->n);
606: TV_add_row("values",type, values);
607: VecRestoreArrayRead((Vec)v,&values);
608: return TV_format_OK;
609: }
610: #endif
612: /*@
613: VecGetSize - Returns the global number of elements of the vector.
615: Not Collective
617: Input Parameter:
618: . x - the vector
620: Output Parameters:
621: . size - the global length of the vector
623: Level: beginner
625: Concepts: vector^local size
627: .seealso: VecGetLocalSize()
628: @*/
629: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
630: {
637: (*x->ops->getsize)(x,size);
638: return(0);
639: }
641: /*@
642: VecGetLocalSize - Returns the number of elements of the vector stored
643: in local memory. This routine may be implementation dependent, so use
644: with care.
646: Not Collective
648: Input Parameter:
649: . x - the vector
651: Output Parameter:
652: . size - the length of the local piece of the vector
654: Level: beginner
656: Concepts: vector^size
658: .seealso: VecGetSize()
659: @*/
660: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
661: {
668: (*x->ops->getlocalsize)(x,size);
669: return(0);
670: }
672: /*@C
673: VecGetOwnershipRange - Returns the range of indices owned by
674: this processor, assuming that the vectors are laid out with the
675: first n1 elements on the first processor, next n2 elements on the
676: second, etc. For certain parallel layouts this range may not be
677: well defined.
679: Not Collective
681: Input Parameter:
682: . x - the vector
684: Output Parameters:
685: + low - the first local element, pass in NULL if not interested
686: - high - one more than the last local element, pass in NULL if not interested
688: Note:
689: The high argument is one more than the last element stored locally.
691: Fortran: NULL_INTEGER should be used instead of NULL
693: Level: beginner
695: Concepts: ownership^of vectors
696: Concepts: vector^ownership of elements
698: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
699: @*/
700: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
701: {
707: if (low) *low = x->map->rstart;
708: if (high) *high = x->map->rend;
709: return(0);
710: }
712: /*@C
713: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
714: assuming that the vectors are laid out with the
715: first n1 elements on the first processor, next n2 elements on the
716: second, etc. For certain parallel layouts this range may not be
717: well defined.
719: Not Collective
721: Input Parameter:
722: . x - the vector
724: Output Parameters:
725: . range - array of length size+1 with the start and end+1 for each process
727: Note:
728: The high argument is one more than the last element stored locally.
730: Fortran: You must PASS in an array of length size+1
732: Level: beginner
734: Concepts: ownership^of vectors
735: Concepts: vector^ownership of elements
737: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
738: @*/
739: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
740: {
746: PetscLayoutGetRanges(x->map,ranges);
747: return(0);
748: }
750: /*@
751: VecSetOption - Sets an option for controling a vector's behavior.
753: Collective on Vec
755: Input Parameter:
756: + x - the vector
757: . op - the option
758: - flag - turn the option on or off
760: Supported Options:
761: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
762: entries destined to be stored on a separate processor. This can be used
763: to eliminate the global reduction in the VecAssemblyXXXX() if you know
764: that you have only used VecSetValues() to set local elements
765: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
766: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
767: ignored.
768: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
769: entries will always be a subset (possibly equal) of the off-process entries set on the
770: first assembly. This reuses the communication pattern, thus avoiding a global reduction.
771: Subsequent assemblies setting off-process values should use the same InsertMode as the
772: first assembly.
774: Developer Note:
775: The InsertMode restriction could be removed by packing the stash messages out of place.
777: Level: intermediate
779: @*/
780: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
781: {
787: if (x->ops->setoption) {
788: (*x->ops->setoption)(x,op,flag);
789: }
790: return(0);
791: }
793: /* Default routines for obtaining and releasing; */
794: /* may be used by any implementation */
795: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
796: {
798: PetscInt i;
803: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
804: PetscMalloc1(m,V);
805: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
806: return(0);
807: }
809: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
810: {
812: PetscInt i;
816: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
817: PetscFree(v);
818: return(0);
819: }
821: /*@
822: VecResetArray - Resets a vector to use its default memory. Call this
823: after the use of VecPlaceArray().
825: Not Collective
827: Input Parameters:
828: . vec - the vector
830: Level: developer
832: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
834: @*/
835: PetscErrorCode VecResetArray(Vec vec)
836: {
842: if (vec->ops->resetarray) {
843: (*vec->ops->resetarray)(vec);
844: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
845: PetscObjectStateIncrease((PetscObject)vec);
846: return(0);
847: }
849: /*@C
850: VecLoad - Loads a vector that has been stored in binary or HDF5 format
851: with VecView().
853: Collective on PetscViewer
855: Input Parameters:
856: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
857: some related function before a call to VecLoad().
858: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
859: HDF5 file viewer, obtained from PetscViewerHDF5Open()
861: Level: intermediate
863: Notes:
864: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
865: before calling this.
867: The input file must contain the full global vector, as
868: written by the routine VecView().
870: If the type or size of newvec is not set before a call to VecLoad, PETSc
871: sets the type and the local and global sizes. If type and/or
872: sizes are already set, then the same are used.
874: If using binary and the blocksize of the vector is greater than one then you must provide a unique prefix to
875: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
876: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
877: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
878: filename. If you copy the binary file, make sure you copy the associated .info file with it.
880: If using HDF5, you must assign the Vec the same name as was used in the Vec
881: that was stored in the file using PetscObjectSetName(). Otherwise you will
882: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
884: Notes for advanced users:
885: Most users should not need to know the details of the binary storage
886: format, since VecLoad() and VecView() completely hide these details.
887: But for anyone who's interested, the standard binary vector storage
888: format is
889: .vb
890: int VEC_FILE_CLASSID
891: int number of rows
892: PetscScalar *values of all entries
893: .ve
895: In addition, PETSc automatically does the byte swapping for
896: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
897: linux, Windows and the paragon; thus if you write your own binary
898: read/write routines you have to swap the bytes; see PetscBinaryRead()
899: and PetscBinaryWrite() to see how this may be done.
901: Concepts: vector^loading from file
903: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
904: @*/
905: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
906: {
908: PetscBool isbinary,ishdf5;
909: PetscViewerFormat format;
914: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
915: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
916: if (!isbinary && !ishdf5) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
918: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
919: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
920: VecSetType(newvec, VECSTANDARD);
921: }
922: PetscViewerGetFormat(viewer,&format);
923: if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
924: (*newvec->ops->loadnative)(newvec,viewer);
925: } else {
926: (*newvec->ops->load)(newvec,viewer);
927: }
928: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
929: return(0);
930: }
933: /*@
934: VecReciprocal - Replaces each component of a vector by its reciprocal.
936: Logically Collective on Vec
938: Input Parameter:
939: . vec - the vector
941: Output Parameter:
942: . vec - the vector reciprocal
944: Level: intermediate
946: Concepts: vector^reciprocal
948: .seealso: VecLog(), VecExp(), VecSqrtAbs()
950: @*/
951: PetscErrorCode VecReciprocal(Vec vec)
952: {
958: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
959: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
960: (*vec->ops->reciprocal)(vec);
961: PetscObjectStateIncrease((PetscObject)vec);
962: return(0);
963: }
965: /*@C
966: VecSetOperation - Allows user to set a vector operation.
968: Logically Collective on Vec
970: Input Parameters:
971: + vec - the vector
972: . op - the name of the operation
973: - f - the function that provides the operation.
975: Level: advanced
977: Usage:
978: $ PetscErrorCode userview(Vec,PetscViewer);
979: $ VecCreateMPI(comm,m,M,&x);
980: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
982: Notes:
983: See the file include/petscvec.h for a complete list of matrix
984: operations, which all have the form VECOP_<OPERATION>, where
985: <OPERATION> is the name (in all capital letters) of the
986: user interface routine (e.g., VecView() -> VECOP_VIEW).
988: This function is not currently available from Fortran.
990: .keywords: vector, set, operation
992: .seealso: VecCreate(), MatShellSetOperation()
993: @*/
994: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
995: {
998: if (op == VECOP_VIEW && !vec->ops->viewnative) {
999: vec->ops->viewnative = vec->ops->view;
1000: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1001: vec->ops->loadnative = vec->ops->load;
1002: }
1003: (((void(**)(void))vec->ops)[(int)op]) = f;
1004: return(0);
1005: }
1008: /*@
1009: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1010: used during the assembly process to store values that belong to
1011: other processors.
1013: Not Collective, different processes can have different size stashes
1015: Input Parameters:
1016: + vec - the vector
1017: . size - the initial size of the stash.
1018: - bsize - the initial size of the block-stash(if used).
1020: Options Database Keys:
1021: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1022: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1024: Level: intermediate
1026: Notes:
1027: The block-stash is used for values set with VecSetValuesBlocked() while
1028: the stash is used for values set with VecSetValues()
1030: Run with the option -info and look for output of the form
1031: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1032: to determine the appropriate value, MM, to use for size and
1033: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1034: to determine the value, BMM to use for bsize
1036: Concepts: vector^stash
1037: Concepts: stash^vector
1039: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1041: @*/
1042: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1043: {
1048: VecStashSetInitialSize_Private(&vec->stash,size);
1049: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1050: return(0);
1051: }
1053: /*@
1054: VecConjugate - Conjugates a vector.
1056: Logically Collective on Vec
1058: Input Parameters:
1059: . x - the vector
1061: Level: intermediate
1063: Concepts: vector^conjugate
1065: @*/
1066: PetscErrorCode VecConjugate(Vec x)
1067: {
1068: #if defined(PETSC_USE_COMPLEX)
1074: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1075: (*x->ops->conjugate)(x);
1076: /* we need to copy norms here */
1077: PetscObjectStateIncrease((PetscObject)x);
1078: return(0);
1079: #else
1080: return(0);
1081: #endif
1082: }
1084: /*@
1085: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1087: Logically Collective on Vec
1089: Input Parameters:
1090: . x, y - the vectors
1092: Output Parameter:
1093: . w - the result
1095: Level: advanced
1097: Notes: any subset of the x, y, and w may be the same vector.
1099: Concepts: vector^pointwise multiply
1101: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1102: @*/
1103: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1104: {
1116: VecCheckSameSize(w,1,x,2);
1117: VecCheckSameSize(w,2,y,3);
1118: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1119: (*w->ops->pointwisemult)(w,x,y);
1120: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1121: PetscObjectStateIncrease((PetscObject)w);
1122: return(0);
1123: }
1125: /*@C
1126: VecSetRandom - Sets all components of a vector to random numbers.
1128: Logically Collective on Vec
1130: Input Parameters:
1131: + x - the vector
1132: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1133: it will create one internally.
1135: Output Parameter:
1136: . x - the vector
1138: Example of Usage:
1139: .vb
1140: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1141: VecSetRandom(x,rctx);
1142: PetscRandomDestroy(rctx);
1143: .ve
1145: Level: intermediate
1147: Concepts: vector^setting to random
1148: Concepts: random^vector
1150: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1151: @*/
1152: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1153: {
1155: PetscRandom randObj = NULL;
1161: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1163: if (!rctx) {
1164: MPI_Comm comm;
1165: PetscObjectGetComm((PetscObject)x,&comm);
1166: PetscRandomCreate(comm,&randObj);
1167: PetscRandomSetFromOptions(randObj);
1168: rctx = randObj;
1169: }
1171: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1172: (*x->ops->setrandom)(x,rctx);
1173: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1175: PetscRandomDestroy(&randObj);
1176: PetscObjectStateIncrease((PetscObject)x);
1177: return(0);
1178: }
1180: /*@
1181: VecZeroEntries - puts a 0.0 in each element of a vector
1183: Logically Collective on Vec
1185: Input Parameter:
1186: . vec - The vector
1188: Level: beginner
1190: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1191: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1192: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1193: this routine should not exist.
1195: .keywords: Vec, set, options, database
1196: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1197: @*/
1198: PetscErrorCode VecZeroEntries(Vec vec)
1199: {
1203: VecSet(vec,0);
1204: return(0);
1205: }
1207: /*
1208: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1209: processor and a PETSc MPI vector on more than one processor.
1211: Collective on Vec
1213: Input Parameter:
1214: . vec - The vector
1216: Level: intermediate
1218: .keywords: Vec, set, options, database, type
1219: .seealso: VecSetFromOptions(), VecSetType()
1220: */
1221: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1222: {
1223: PetscBool opt;
1224: VecType defaultType;
1225: char typeName[256];
1226: PetscMPIInt size;
1230: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1231: else {
1232: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1233: if (size > 1) defaultType = VECMPI;
1234: else defaultType = VECSEQ;
1235: }
1237: VecRegisterAll();
1238: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1239: if (opt) {
1240: VecSetType(vec, typeName);
1241: } else {
1242: VecSetType(vec, defaultType);
1243: }
1244: return(0);
1245: }
1247: /*@
1248: VecSetFromOptions - Configures the vector from the options database.
1250: Collective on Vec
1252: Input Parameter:
1253: . vec - The vector
1255: Notes: To see all options, run your program with the -help option, or consult the users manual.
1256: Must be called after VecCreate() but before the vector is used.
1258: Level: beginner
1260: Concepts: vectors^setting options
1261: Concepts: vectors^setting type
1263: .keywords: Vec, set, options, database
1264: .seealso: VecCreate(), VecSetOptionsPrefix()
1265: @*/
1266: PetscErrorCode VecSetFromOptions(Vec vec)
1267: {
1273: PetscObjectOptionsBegin((PetscObject)vec);
1274: /* Handle vector type options */
1275: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1277: /* Handle specific vector options */
1278: if (vec->ops->setfromoptions) {
1279: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1280: }
1282: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1283: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1284: PetscOptionsEnd();
1285: return(0);
1286: }
1288: /*@
1289: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1291: Collective on Vec
1293: Input Parameters:
1294: + v - the vector
1295: . n - the local size (or PETSC_DECIDE to have it set)
1296: - N - the global size (or PETSC_DECIDE)
1298: Notes:
1299: n and N cannot be both PETSC_DECIDE
1300: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1302: Level: intermediate
1304: .seealso: VecGetSize(), PetscSplitOwnership()
1305: @*/
1306: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1307: {
1313: 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);
1314: 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);
1315: v->map->n = n;
1316: v->map->N = N;
1317: if (v->ops->create) {
1318: (*v->ops->create)(v);
1319: v->ops->create = 0;
1320: }
1321: return(0);
1322: }
1324: /*@
1325: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1326: and VecSetValuesBlockedLocal().
1328: Logically Collective on Vec
1330: Input Parameter:
1331: + v - the vector
1332: - bs - the blocksize
1334: Notes:
1335: All vectors obtained by VecDuplicate() inherit the same blocksize.
1337: Level: advanced
1339: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1341: Concepts: block size^vectors
1342: @*/
1343: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1344: {
1349: if (bs < 0 || bs == v->map->bs) return(0);
1351: PetscLayoutSetBlockSize(v->map,bs);
1352: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1353: return(0);
1354: }
1356: /*@
1357: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1358: and VecSetValuesBlockedLocal().
1360: Not Collective
1362: Input Parameter:
1363: . v - the vector
1365: Output Parameter:
1366: . bs - the blocksize
1368: Notes:
1369: All vectors obtained by VecDuplicate() inherit the same blocksize.
1371: Level: advanced
1373: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1375: Concepts: vector^block size
1376: Concepts: block^vector
1378: @*/
1379: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1380: {
1386: PetscLayoutGetBlockSize(v->map,bs);
1387: return(0);
1388: }
1390: /*@C
1391: VecSetOptionsPrefix - Sets the prefix used for searching for all
1392: Vec options in the database.
1394: Logically Collective on Vec
1396: Input Parameter:
1397: + v - the Vec context
1398: - prefix - the prefix to prepend to all option names
1400: Notes:
1401: A hyphen (-) must NOT be given at the beginning of the prefix name.
1402: The first character of all runtime options is AUTOMATICALLY the hyphen.
1404: Level: advanced
1406: .keywords: Vec, set, options, prefix, database
1408: .seealso: VecSetFromOptions()
1409: @*/
1410: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1411: {
1416: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1417: return(0);
1418: }
1420: /*@C
1421: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1422: Vec options in the database.
1424: Logically Collective on Vec
1426: Input Parameters:
1427: + v - the Vec context
1428: - prefix - the prefix to prepend to all option names
1430: Notes:
1431: A hyphen (-) must NOT be given at the beginning of the prefix name.
1432: The first character of all runtime options is AUTOMATICALLY the hyphen.
1434: Level: advanced
1436: .keywords: Vec, append, options, prefix, database
1438: .seealso: VecGetOptionsPrefix()
1439: @*/
1440: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1441: {
1446: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1447: return(0);
1448: }
1450: /*@C
1451: VecGetOptionsPrefix - Sets the prefix used for searching for all
1452: Vec options in the database.
1454: Not Collective
1456: Input Parameter:
1457: . v - the Vec context
1459: Output Parameter:
1460: . prefix - pointer to the prefix string used
1462: Notes: On the fortran side, the user should pass in a string 'prefix' of
1463: sufficient length to hold the prefix.
1465: Level: advanced
1467: .keywords: Vec, get, options, prefix, database
1469: .seealso: VecAppendOptionsPrefix()
1470: @*/
1471: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1472: {
1477: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1478: return(0);
1479: }
1481: /*@
1482: VecSetUp - Sets up the internal vector data structures for the later use.
1484: Collective on Vec
1486: Input Parameters:
1487: . v - the Vec context
1489: Notes:
1490: For basic use of the Vec classes the user need not explicitly call
1491: VecSetUp(), since these actions will happen automatically.
1493: Level: advanced
1495: .keywords: Vec, setup
1497: .seealso: VecCreate(), VecDestroy()
1498: @*/
1499: PetscErrorCode VecSetUp(Vec v)
1500: {
1501: PetscMPIInt size;
1506: if (!((PetscObject)v)->type_name) {
1507: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1508: if (size == 1) {
1509: VecSetType(v, VECSEQ);
1510: } else {
1511: VecSetType(v, VECMPI);
1512: }
1513: }
1514: return(0);
1515: }
1517: /*
1518: These currently expose the PetscScalar/PetscReal in updating the
1519: cached norm. If we push those down into the implementation these
1520: will become independent of PetscScalar/PetscReal
1521: */
1523: /*@
1524: VecCopy - Copies a vector. y <- x
1526: Logically Collective on Vec
1528: Input Parameter:
1529: . x - the vector
1531: Output Parameter:
1532: . y - the copy
1534: Notes:
1535: For default parallel PETSc vectors, both x and y must be distributed in
1536: the same manner; local copies are done.
1538: Developer Notes:
1540: of the vectors to be sequential and one to be parallel so long as both have the same
1541: local sizes. This is used in some internal functions in PETSc.
1543: Level: beginner
1545: .seealso: VecDuplicate()
1546: @*/
1547: PetscErrorCode VecCopy(Vec x,Vec y)
1548: {
1549: PetscBool flgs[4];
1550: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1552: PetscInt i;
1559: if (x == y) return(0);
1560: VecCheckSameLocalSize(x,1,y,2);
1561: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1562: VecLocked(y,2);
1564: #if !defined(PETSC_USE_MIXED_PRECISION)
1565: for (i=0; i<4; i++) {
1566: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1567: }
1568: #endif
1570: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1571: #if defined(PETSC_USE_MIXED_PRECISION)
1572: extern PetscErrorCode VecGetArray(Vec,double**);
1573: extern PetscErrorCode VecRestoreArray(Vec,double**);
1574: extern PetscErrorCode VecGetArray(Vec,float**);
1575: extern PetscErrorCode VecRestoreArray(Vec,float**);
1576: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1577: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1578: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1579: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1580: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1581: PetscInt i,n;
1582: const float *xx;
1583: double *yy;
1584: VecGetArrayRead(x,&xx);
1585: VecGetArray(y,&yy);
1586: VecGetLocalSize(x,&n);
1587: for (i=0; i<n; i++) yy[i] = xx[i];
1588: VecRestoreArrayRead(x,&xx);
1589: VecRestoreArray(y,&yy);
1590: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1591: PetscInt i,n;
1592: float *yy;
1593: const double *xx;
1594: VecGetArrayRead(x,&xx);
1595: VecGetArray(y,&yy);
1596: VecGetLocalSize(x,&n);
1597: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1598: VecRestoreArrayRead(x,&xx);
1599: VecRestoreArray(y,&yy);
1600: } else {
1601: (*x->ops->copy)(x,y);
1602: }
1603: #else
1604: (*x->ops->copy)(x,y);
1605: #endif
1607: PetscObjectStateIncrease((PetscObject)y);
1608: #if !defined(PETSC_USE_MIXED_PRECISION)
1609: for (i=0; i<4; i++) {
1610: if (flgs[i]) {
1611: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1612: }
1613: }
1614: #endif
1616: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1617: return(0);
1618: }
1620: /*@
1621: VecSwap - Swaps the vectors x and y.
1623: Logically Collective on Vec
1625: Input Parameters:
1626: . x, y - the vectors
1628: Level: advanced
1630: Concepts: vector^swapping values
1632: @*/
1633: PetscErrorCode VecSwap(Vec x,Vec y)
1634: {
1635: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1636: PetscBool flgxs[4],flgys[4];
1638: PetscInt i;
1646: VecCheckSameSize(x,1,y,2);
1647: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1648: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1650: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1651: for (i=0; i<4; i++) {
1652: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1653: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1654: }
1655: (*x->ops->swap)(x,y);
1656: PetscObjectStateIncrease((PetscObject)x);
1657: PetscObjectStateIncrease((PetscObject)y);
1658: for (i=0; i<4; i++) {
1659: if (flgxs[i]) {
1660: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1661: }
1662: if (flgys[i]) {
1663: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1664: }
1665: }
1666: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1667: return(0);
1668: }
1670: /*
1671: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1673: Collective on VecStash
1675: Input Parameters:
1676: + obj - the VecStash object
1677: . bobj - optional other object that provides the prefix
1678: - optionname - option to activate viewing
1680: Level: intermediate
1682: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1684: */
1685: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1686: {
1687: PetscErrorCode ierr;
1688: PetscViewer viewer;
1689: PetscBool flg;
1690: PetscViewerFormat format;
1691: char *prefix;
1694: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1695: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),prefix,optionname,&viewer,&format,&flg);
1696: if (flg) {
1697: PetscViewerPushFormat(viewer,format);
1698: VecStashView(obj,viewer);
1699: PetscViewerPopFormat(viewer);
1700: PetscViewerDestroy(&viewer);
1701: }
1702: return(0);
1703: }
1705: /*@
1706: VecStashView - Prints the entries in the vector stash and block stash.
1708: Collective on Vec
1710: Input Parameters:
1711: + v - the vector
1712: - viewer - the viewer
1714: Level: advanced
1716: Concepts: vector^stash
1717: Concepts: stash^vector
1719: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1721: @*/
1722: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1723: {
1725: PetscMPIInt rank;
1726: PetscInt i,j;
1727: PetscBool match;
1728: VecStash *s;
1729: PetscScalar val;
1736: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1737: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1738: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1739: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1740: s = &v->bstash;
1742: /* print block stash */
1743: PetscViewerASCIIPushSynchronized(viewer);
1744: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1745: for (i=0; i<s->n; i++) {
1746: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1747: for (j=0; j<s->bs; j++) {
1748: val = s->array[i*s->bs+j];
1749: #if defined(PETSC_USE_COMPLEX)
1750: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1751: #else
1752: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1753: #endif
1754: }
1755: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1756: }
1757: PetscViewerFlush(viewer);
1759: s = &v->stash;
1761: /* print basic stash */
1762: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1763: for (i=0; i<s->n; i++) {
1764: val = s->array[i];
1765: #if defined(PETSC_USE_COMPLEX)
1766: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1767: #else
1768: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1769: #endif
1770: }
1771: PetscViewerFlush(viewer);
1772: PetscViewerASCIIPopSynchronized(viewer);
1773: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1774: return(0);
1775: }
1777: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1778: {
1779: PetscInt i,N,rstart,rend;
1781: PetscScalar *xx;
1782: PetscReal *xreal;
1783: PetscBool iset;
1786: VecGetOwnershipRange(v,&rstart,&rend);
1787: VecGetSize(v,&N);
1788: PetscCalloc1(N,&xreal);
1789: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1790: if (iset) {
1791: VecGetArray(v,&xx);
1792: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1793: VecRestoreArray(v,&xx);
1794: }
1795: PetscFree(xreal);
1796: if (set) *set = iset;
1797: return(0);
1798: }
1800: /*@
1801: VecGetLayout - get PetscLayout describing vector layout
1803: Not Collective
1805: Input Arguments:
1806: . x - the vector
1808: Output Arguments:
1809: . map - the layout
1811: Level: developer
1813: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1814: @*/
1815: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1816: {
1820: *map = x->map;
1821: return(0);
1822: }
1824: /*@
1825: VecSetLayout - set PetscLayout describing vector layout
1827: Not Collective
1829: Input Arguments:
1830: + x - the vector
1831: - map - the layout
1833: Notes:
1834: It is normally only valid to replace the layout with a layout known to be equivalent.
1836: Level: developer
1838: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1839: @*/
1840: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1841: {
1846: PetscLayoutReference(map,&x->map);
1847: return(0);
1848: }
1850: PetscErrorCode VecSetInf(Vec xin)
1851: {
1852: PetscInt i,n = xin->map->n;
1853: PetscScalar *xx;
1854: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1858: VecGetArray(xin,&xx);
1859: for (i=0; i<n; i++) xx[i] = inf;
1860: VecRestoreArray(xin,&xx);
1861: return(0);
1862: }