Actual source code: vector.c
petsc-3.11.4 2019-09-28
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_Dot, 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_MAXPY, VEC_Swap, VEC_AssemblyBegin, VEC_ScatterBegin, VEC_ScatterEnd;
13: PetscLogEvent VEC_AssemblyEnd, VEC_PointwiseMult, VEC_SetValues, VEC_Load;
14: PetscLogEvent VEC_SetRandom, VEC_ReduceArithmetic, VEC_ReduceCommunication,VEC_ReduceBegin,VEC_ReduceEnd,VEC_Ops;
15: PetscLogEvent VEC_DotNorm2, 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:
199: any subset of the x, y, and w may be the same vector.
200: For complex numbers compares only the real part
202: Concepts: vector^pointwise multiply
204: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
205: @*/
206: PetscErrorCode VecPointwiseMax(Vec w,Vec x,Vec y)
207: {
219: VecCheckSameSize(w,1,x,2);
220: VecCheckSameSize(w,1,y,3);
221: (*w->ops->pointwisemax)(w,x,y);
222: PetscObjectStateIncrease((PetscObject)w);
223: return(0);
224: }
227: /*@
228: VecPointwiseMin - Computes the componentwise minimum w_i = min(x_i, y_i).
230: Logically Collective on Vec
232: Input Parameters:
233: . x, y - the vectors
235: Output Parameter:
236: . w - the result
238: Level: advanced
240: Notes:
241: any subset of the x, y, and w may be the same vector.
242: For complex numbers compares only the real part
244: Concepts: vector^pointwise multiply
246: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
247: @*/
248: PetscErrorCode VecPointwiseMin(Vec w,Vec x,Vec y)
249: {
261: VecCheckSameSize(w,1,x,2);
262: VecCheckSameSize(w,1,y,3);
263: (*w->ops->pointwisemin)(w,x,y);
264: PetscObjectStateIncrease((PetscObject)w);
265: return(0);
266: }
268: /*@
269: VecPointwiseMaxAbs - Computes the componentwise maximum of the absolute values w_i = max(abs(x_i), abs(y_i)).
271: Logically Collective on Vec
273: Input Parameters:
274: . x, y - the vectors
276: Output Parameter:
277: . w - the result
279: Level: advanced
281: Notes:
282: any subset of the x, y, and w may be the same vector.
284: Concepts: vector^pointwise multiply
286: .seealso: VecPointwiseDivide(), VecPointwiseMult(), VecPointwiseMin(), VecPointwiseMax(), VecMaxPointwiseDivide()
287: @*/
288: PetscErrorCode VecPointwiseMaxAbs(Vec w,Vec x,Vec y)
289: {
301: VecCheckSameSize(w,1,x,2);
302: VecCheckSameSize(w,1,y,3);
303: (*w->ops->pointwisemaxabs)(w,x,y);
304: PetscObjectStateIncrease((PetscObject)w);
305: return(0);
306: }
308: /*@
309: VecPointwiseDivide - Computes the componentwise division w = x/y.
311: Logically Collective on Vec
313: Input Parameters:
314: . x, y - the vectors
316: Output Parameter:
317: . w - the result
319: Level: advanced
321: Notes:
322: any subset of the x, y, and w may be the same vector.
324: Concepts: vector^pointwise divide
326: .seealso: VecPointwiseMult(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
327: @*/
328: PetscErrorCode VecPointwiseDivide(Vec w,Vec x,Vec y)
329: {
341: VecCheckSameSize(w,1,x,2);
342: VecCheckSameSize(w,1,y,3);
343: (*w->ops->pointwisedivide)(w,x,y);
344: PetscObjectStateIncrease((PetscObject)w);
345: return(0);
346: }
349: /*@
350: VecDuplicate - Creates a new vector of the same type as an existing vector.
352: Collective on Vec
354: Input Parameters:
355: . v - a vector to mimic
357: Output Parameter:
358: . newv - location to put new vector
360: Notes:
361: VecDuplicate() DOES NOT COPY the vector entries, but rather allocates storage
362: for the new vector. Use VecCopy() to copy a vector.
364: Use VecDestroy() to free the space. Use VecDuplicateVecs() to get several
365: vectors.
367: Level: beginner
369: .seealso: VecDestroy(), VecDuplicateVecs(), VecCreate(), VecCopy()
370: @*/
371: PetscErrorCode VecDuplicate(Vec v,Vec *newv)
372: {
379: (*v->ops->duplicate)(v,newv);
380: PetscObjectStateIncrease((PetscObject)*newv);
381: return(0);
382: }
384: /*@
385: VecDestroy - Destroys a vector.
387: Collective on Vec
389: Input Parameters:
390: . v - the vector
392: Level: beginner
394: .seealso: VecDuplicate(), VecDestroyVecs()
395: @*/
396: PetscErrorCode VecDestroy(Vec *v)
397: {
401: if (!*v) return(0);
403: if (--((PetscObject)(*v))->refct > 0) {*v = 0; return(0);}
405: PetscObjectSAWsViewOff((PetscObject)*v);
406: /* destroy the internal part */
407: if ((*v)->ops->destroy) {
408: (*(*v)->ops->destroy)(*v);
409: }
410: /* destroy the external/common part */
411: PetscLayoutDestroy(&(*v)->map);
412: PetscHeaderDestroy(v);
413: return(0);
414: }
416: /*@C
417: VecDuplicateVecs - Creates several vectors of the same type as an existing vector.
419: Collective on Vec
421: Input Parameters:
422: + m - the number of vectors to obtain
423: - v - a vector to mimic
425: Output Parameter:
426: . V - location to put pointer to array of vectors
428: Notes:
429: Use VecDestroyVecs() to free the space. Use VecDuplicate() to form a single
430: vector.
432: Fortran Note:
433: The Fortran interface is slightly different from that given below, it
434: requires one to pass in V a Vec (integer) array of size at least m.
435: See the Fortran chapter of the users manual and petsc/src/vec/vec/examples for details.
437: Level: intermediate
439: .seealso: VecDestroyVecs(), VecDuplicate(), VecCreate(), VecDuplicateVecsF90()
440: @*/
441: PetscErrorCode VecDuplicateVecs(Vec v,PetscInt m,Vec *V[])
442: {
449: (*v->ops->duplicatevecs)(v,m,V);
450: return(0);
451: }
453: /*@C
454: VecDestroyVecs - Frees a block of vectors obtained with VecDuplicateVecs().
456: Collective on Vec
458: Input Parameters:
459: + vv - pointer to pointer to array of vector pointers, if NULL no vectors are destroyed
460: - m - the number of vectors previously obtained, if zero no vectors are destroyed
462: Fortran Note:
463: The Fortran interface is slightly different from that given below.
464: See the Fortran chapter of the users manual
466: Level: intermediate
468: .seealso: VecDuplicateVecs(), VecDestroyVecsf90()
469: @*/
470: PetscErrorCode VecDestroyVecs(PetscInt m,Vec *vv[])
471: {
476: if (m < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Trying to destroy negative number of vectors %D",m);
477: if (!m || !*vv) {*vv = NULL; return(0);}
480: (*(**vv)->ops->destroyvecs)(m,*vv);
481: *vv = NULL;
482: return(0);
483: }
485: /*@C
486: VecView - Views a vector object.
488: Collective on Vec
490: Input Parameters:
491: + vec - the vector
492: - viewer - an optional visualization context
494: Notes:
495: The available visualization contexts include
496: + PETSC_VIEWER_STDOUT_SELF - for sequential vectors
497: . PETSC_VIEWER_STDOUT_WORLD - for parallel vectors created on PETSC_COMM_WORLD
498: - PETSC_VIEWER_STDOUT_(comm) - for parallel vectors created on MPI communicator comm
500: You can change the format the vector is printed using the
501: option PetscViewerPushFormat().
503: The user can open alternative visualization contexts with
504: + PetscViewerASCIIOpen() - Outputs vector to a specified file
505: . PetscViewerBinaryOpen() - Outputs vector in binary to a
506: specified file; corresponding input uses VecLoad()
507: . PetscViewerDrawOpen() - Outputs vector to an X window display
508: - PetscViewerSocketOpen() - Outputs vector to Socket viewer
510: The user can call PetscViewerPushFormat() to specify the output
511: format of ASCII printed objects (when using PETSC_VIEWER_STDOUT_SELF,
512: PETSC_VIEWER_STDOUT_WORLD and PetscViewerASCIIOpen). Available formats include
513: + PETSC_VIEWER_DEFAULT - default, prints vector contents
514: . PETSC_VIEWER_ASCII_MATLAB - prints vector contents in MATLAB format
515: . PETSC_VIEWER_ASCII_INDEX - prints vector contents, including indices of vector elements
516: - PETSC_VIEWER_ASCII_COMMON - prints vector contents, using a
517: format common among all vector types
519: Notes:
520: You can pass any number of vector objects, or other PETSc objects to the same viewer.
522: Notes for binary viewer: If you pass multiply vectors to a binary viewer you can read them back in in the same order
523: $ with VecLoad().
524: $
525: $ If the blocksize of the vector is greater than one then you must provide a unique prefix to
526: $ the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
527: $ vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
528: $ information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
529: $ filename. If you copy the binary file, make sure you copy the associated .info file with it.
531: Notes for HDF5 Viewer: the name of the Vec (given with PetscObjectSetName() is the name that is used
532: $ for the object in the HDF5 file. If you wish to store the same vector to the HDF5 viewer (with different values,
533: $ obviously) several times, you must change its name each time before calling the VecView(). The name you use
534: $ here should equal the name that you use in the Vec object that you use with VecLoad().
536: See the manual page for VecLoad() on the exact format the binary viewer stores
537: the values in the file.
539: Level: beginner
541: Concepts: vector^printing
542: Concepts: vector^saving to disk
544: .seealso: PetscViewerASCIIOpen(), PetscViewerDrawOpen(), PetscDrawLGCreate(),
545: PetscViewerSocketOpen(), PetscViewerBinaryOpen(), VecLoad(), PetscViewerCreate(),
546: PetscRealView(), PetscScalarView(), PetscIntView()
547: @*/
548: PetscErrorCode VecView(Vec vec,PetscViewer viewer)
549: {
550: PetscErrorCode ierr;
551: PetscBool iascii;
552: PetscViewerFormat format;
553: PetscMPIInt size;
558: if (!viewer) {
559: PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)vec),&viewer);
560: }
562: PetscViewerGetFormat(viewer,&format);
563: MPI_Comm_size(PetscObjectComm((PetscObject)vec),&size);
564: if (size == 1 && format == PETSC_VIEWER_LOAD_BALANCE) return(0);
566: if (vec->stash.n || vec->bstash.n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Must call VecAssemblyBegin/End() before viewing this vector");
568: PetscLogEventBegin(VEC_View,vec,viewer,0,0);
569: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
570: if (iascii) {
571: PetscInt rows,bs;
573: PetscObjectPrintClassNamePrefixType((PetscObject)vec,viewer);
574: if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) {
575: PetscViewerASCIIPushTab(viewer);
576: VecGetSize(vec,&rows);
577: VecGetBlockSize(vec,&bs);
578: if (bs != 1) {
579: PetscViewerASCIIPrintf(viewer,"length=%D, bs=%D\n",rows,bs);
580: } else {
581: PetscViewerASCIIPrintf(viewer,"length=%D\n",rows);
582: }
583: PetscViewerASCIIPopTab(viewer);
584: }
585: }
586: VecLockReadPush(vec);
587: if ((format == PETSC_VIEWER_NATIVE || format == PETSC_VIEWER_LOAD_BALANCE) && vec->ops->viewnative) {
588: (*vec->ops->viewnative)(vec,viewer);
589: } else {
590: (*vec->ops->view)(vec,viewer);
591: }
592: VecLockReadPop(vec);
593: PetscLogEventEnd(VEC_View,vec,viewer,0,0);
594: return(0);
595: }
597: #if defined(PETSC_USE_DEBUG)
598: #include <../src/sys/totalview/tv_data_display.h>
599: PETSC_UNUSED static int TV_display_type(const struct _p_Vec *v)
600: {
601: const PetscScalar *values;
602: char type[32];
603: PetscErrorCode ierr;
606: TV_add_row("Local rows", "int", &v->map->n);
607: TV_add_row("Global rows", "int", &v->map->N);
608: TV_add_row("Typename", TV_ascii_string_type, ((PetscObject)v)->type_name);
609: VecGetArrayRead((Vec)v,&values);
610: PetscSNPrintf(type,32,"double[%d]",v->map->n);
611: TV_add_row("values",type, values);
612: VecRestoreArrayRead((Vec)v,&values);
613: return TV_format_OK;
614: }
615: #endif
617: /*@
618: VecGetSize - Returns the global number of elements of the vector.
620: Not Collective
622: Input Parameter:
623: . x - the vector
625: Output Parameters:
626: . size - the global length of the vector
628: Level: beginner
630: Concepts: vector^local size
632: .seealso: VecGetLocalSize()
633: @*/
634: PetscErrorCode VecGetSize(Vec x,PetscInt *size)
635: {
642: (*x->ops->getsize)(x,size);
643: return(0);
644: }
646: /*@
647: VecGetLocalSize - Returns the number of elements of the vector stored
648: in local memory. This routine may be implementation dependent, so use
649: with care.
651: Not Collective
653: Input Parameter:
654: . x - the vector
656: Output Parameter:
657: . size - the length of the local piece of the vector
659: Level: beginner
661: Concepts: vector^size
663: .seealso: VecGetSize()
664: @*/
665: PetscErrorCode VecGetLocalSize(Vec x,PetscInt *size)
666: {
673: (*x->ops->getlocalsize)(x,size);
674: return(0);
675: }
677: /*@C
678: VecGetOwnershipRange - Returns the range of indices owned by
679: this processor, assuming that the vectors are laid out with the
680: first n1 elements on the first processor, next n2 elements on the
681: second, etc. For certain parallel layouts this range may not be
682: well defined.
684: Not Collective
686: Input Parameter:
687: . x - the vector
689: Output Parameters:
690: + low - the first local element, pass in NULL if not interested
691: - high - one more than the last local element, pass in NULL if not interested
693: Note:
694: The high argument is one more than the last element stored locally.
696: Fortran: PETSC_NULL_INTEGER should be used instead of NULL
698: Level: beginner
700: Concepts: ownership^of vectors
701: Concepts: vector^ownership of elements
703: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRanges()
704: @*/
705: PetscErrorCode VecGetOwnershipRange(Vec x,PetscInt *low,PetscInt *high)
706: {
712: if (low) *low = x->map->rstart;
713: if (high) *high = x->map->rend;
714: return(0);
715: }
717: /*@C
718: VecGetOwnershipRanges - Returns the range of indices owned by EACH processor,
719: assuming that the vectors are laid out with the
720: first n1 elements on the first processor, next n2 elements on the
721: second, etc. For certain parallel layouts this range may not be
722: well defined.
724: Not Collective
726: Input Parameter:
727: . x - the vector
729: Output Parameters:
730: . range - array of length size+1 with the start and end+1 for each process
732: Note:
733: The high argument is one more than the last element stored locally.
735: Fortran: You must PASS in an array of length size+1
737: Level: beginner
739: Concepts: ownership^of vectors
740: Concepts: vector^ownership of elements
742: .seealso: MatGetOwnershipRange(), MatGetOwnershipRanges(), VecGetOwnershipRange()
743: @*/
744: PetscErrorCode VecGetOwnershipRanges(Vec x,const PetscInt *ranges[])
745: {
751: PetscLayoutGetRanges(x->map,ranges);
752: return(0);
753: }
755: /*@
756: VecSetOption - Sets an option for controling a vector's behavior.
758: Collective on Vec
760: Input Parameter:
761: + x - the vector
762: . op - the option
763: - flag - turn the option on or off
765: Supported Options:
766: + VEC_IGNORE_OFF_PROC_ENTRIES, which causes VecSetValues() to ignore
767: entries destined to be stored on a separate processor. This can be used
768: to eliminate the global reduction in the VecAssemblyXXXX() if you know
769: that you have only used VecSetValues() to set local elements
770: . VEC_IGNORE_NEGATIVE_INDICES, which means you can pass negative indices
771: in ix in calls to VecSetValues() or VecGetValues(). These rows are simply
772: ignored.
773: - VEC_SUBSET_OFF_PROC_ENTRIES, which causes VecAssemblyBegin() to assume that the off-process
774: entries will always be a subset (possibly equal) of the off-process entries set on the
775: first assembly which had a true VEC_SUBSET_OFF_PROC_ENTRIES and the vector has not
776: changed this flag afterwards. If this assembly is not such first assembly, then this
777: assembly can reuse the communication pattern setup in that first assembly, thus avoiding
778: a global reduction. Subsequent assemblies setting off-process values should use the same
779: InsertMode as the first assembly.
781: Developer Note:
782: The InsertMode restriction could be removed by packing the stash messages out of place.
784: Level: intermediate
786: @*/
787: PetscErrorCode VecSetOption(Vec x,VecOption op,PetscBool flag)
788: {
794: if (x->ops->setoption) {
795: (*x->ops->setoption)(x,op,flag);
796: }
797: return(0);
798: }
800: /* Default routines for obtaining and releasing; */
801: /* may be used by any implementation */
802: PetscErrorCode VecDuplicateVecs_Default(Vec w,PetscInt m,Vec *V[])
803: {
805: PetscInt i;
810: if (m <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"m must be > 0: m = %D",m);
811: PetscMalloc1(m,V);
812: for (i=0; i<m; i++) {VecDuplicate(w,*V+i);}
813: return(0);
814: }
816: PetscErrorCode VecDestroyVecs_Default(PetscInt m,Vec v[])
817: {
819: PetscInt i;
823: for (i=0; i<m; i++) {VecDestroy(&v[i]);}
824: PetscFree(v);
825: return(0);
826: }
828: /*@
829: VecResetArray - Resets a vector to use its default memory. Call this
830: after the use of VecPlaceArray().
832: Not Collective
834: Input Parameters:
835: . vec - the vector
837: Level: developer
839: .seealso: VecGetArray(), VecRestoreArray(), VecReplaceArray(), VecPlaceArray()
841: @*/
842: PetscErrorCode VecResetArray(Vec vec)
843: {
849: if (vec->ops->resetarray) {
850: (*vec->ops->resetarray)(vec);
851: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Cannot reset array in this type of vector");
852: PetscObjectStateIncrease((PetscObject)vec);
853: return(0);
854: }
856: /*@C
857: VecLoad - Loads a vector that has been stored in binary or HDF5 format
858: with VecView().
860: Collective on PetscViewer
862: Input Parameters:
863: + newvec - the newly loaded vector, this needs to have been created with VecCreate() or
864: some related function before a call to VecLoad().
865: - viewer - binary file viewer, obtained from PetscViewerBinaryOpen() or
866: HDF5 file viewer, obtained from PetscViewerHDF5Open()
868: Level: intermediate
870: Notes:
871: Defaults to the standard Seq or MPI Vec, if you want some other type of Vec call VecSetFromOptions()
872: before calling this.
874: The input file must contain the full global vector, as
875: written by the routine VecView().
877: If the type or size of newvec is not set before a call to VecLoad, PETSc
878: sets the type and the local and global sizes. If type and/or
879: sizes are already set, then the same are used.
881: If using binary and the blocksize of the vector is greater than one then you must provide a unique prefix to
882: the vector with PetscObjectSetOptionsPrefix((PetscObject)vec,"uniqueprefix"); BEFORE calling VecView() on the
883: vector to be stored and then set that same unique prefix on the vector that you pass to VecLoad(). The blocksize
884: information is stored in an ASCII file with the same name as the binary file plus a ".info" appended to the
885: filename. If you copy the binary file, make sure you copy the associated .info file with it.
887: If using HDF5, you must assign the Vec the same name as was used in the Vec
888: that was stored in the file using PetscObjectSetName(). Otherwise you will
889: get the error message: "Cannot H5DOpen2() with Vec name NAMEOFOBJECT"
891: Notes for advanced users:
892: Most users should not need to know the details of the binary storage
893: format, since VecLoad() and VecView() completely hide these details.
894: But for anyone who's interested, the standard binary vector storage
895: format is
896: .vb
897: int VEC_FILE_CLASSID
898: int number of rows
899: PetscScalar *values of all entries
900: .ve
902: In addition, PETSc automatically does the byte swapping for
903: machines that store the bytes reversed, e.g. DEC alpha, freebsd,
904: linux, Windows and the paragon; thus if you write your own binary
905: read/write routines you have to swap the bytes; see PetscBinaryRead()
906: and PetscBinaryWrite() to see how this may be done.
908: Concepts: vector^loading from file
910: .seealso: PetscViewerBinaryOpen(), VecView(), MatLoad(), VecLoad()
911: @*/
912: PetscErrorCode VecLoad(Vec newvec, PetscViewer viewer)
913: {
914: PetscErrorCode ierr;
915: PetscBool isbinary,ishdf5,isadios,isadios2;
916: PetscViewerFormat format;
921: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
922: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5,&ishdf5);
923: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios);
924: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERADIOS,&isadios2);
925: if (!isbinary && !ishdf5 && !isadios && !isadios2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Invalid viewer; open viewer with PetscViewerBinaryOpen()");
927: PetscLogEventBegin(VEC_Load,viewer,0,0,0);
928: if (!((PetscObject)newvec)->type_name && !newvec->ops->create) {
929: VecSetType(newvec, VECSTANDARD);
930: }
931: PetscViewerGetFormat(viewer,&format);
932: if (format == PETSC_VIEWER_NATIVE && newvec->ops->loadnative) {
933: (*newvec->ops->loadnative)(newvec,viewer);
934: } else {
935: (*newvec->ops->load)(newvec,viewer);
936: }
937: PetscLogEventEnd(VEC_Load,viewer,0,0,0);
938: return(0);
939: }
942: /*@
943: VecReciprocal - Replaces each component of a vector by its reciprocal.
945: Logically Collective on Vec
947: Input Parameter:
948: . vec - the vector
950: Output Parameter:
951: . vec - the vector reciprocal
953: Level: intermediate
955: Concepts: vector^reciprocal
957: .seealso: VecLog(), VecExp(), VecSqrtAbs()
959: @*/
960: PetscErrorCode VecReciprocal(Vec vec)
961: {
967: if (vec->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
968: if (!vec->ops->reciprocal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Vector does not support reciprocal operation");
969: (*vec->ops->reciprocal)(vec);
970: PetscObjectStateIncrease((PetscObject)vec);
971: return(0);
972: }
974: /*@C
975: VecSetOperation - Allows user to set a vector operation.
977: Logically Collective on Vec
979: Input Parameters:
980: + vec - the vector
981: . op - the name of the operation
982: - f - the function that provides the operation.
984: Level: advanced
986: Usage:
987: $ PetscErrorCode userview(Vec,PetscViewer);
988: $ VecCreateMPI(comm,m,M,&x);
989: $ VecSetOperation(x,VECOP_VIEW,(void(*)(void))userview);
991: Notes:
992: See the file include/petscvec.h for a complete list of matrix
993: operations, which all have the form VECOP_<OPERATION>, where
994: <OPERATION> is the name (in all capital letters) of the
995: user interface routine (e.g., VecView() -> VECOP_VIEW).
997: This function is not currently available from Fortran.
999: .keywords: vector, set, operation
1001: .seealso: VecCreate(), MatShellSetOperation()
1002: @*/
1003: PetscErrorCode VecSetOperation(Vec vec,VecOperation op, void (*f)(void))
1004: {
1007: if (op == VECOP_VIEW && !vec->ops->viewnative) {
1008: vec->ops->viewnative = vec->ops->view;
1009: } else if (op == VECOP_LOAD && !vec->ops->loadnative) {
1010: vec->ops->loadnative = vec->ops->load;
1011: }
1012: (((void(**)(void))vec->ops)[(int)op]) = f;
1013: return(0);
1014: }
1017: /*@
1018: VecStashSetInitialSize - sets the sizes of the vec-stash, that is
1019: used during the assembly process to store values that belong to
1020: other processors.
1022: Not Collective, different processes can have different size stashes
1024: Input Parameters:
1025: + vec - the vector
1026: . size - the initial size of the stash.
1027: - bsize - the initial size of the block-stash(if used).
1029: Options Database Keys:
1030: + -vecstash_initial_size <size> or <size0,size1,...sizep-1>
1031: - -vecstash_block_initial_size <bsize> or <bsize0,bsize1,...bsizep-1>
1033: Level: intermediate
1035: Notes:
1036: The block-stash is used for values set with VecSetValuesBlocked() while
1037: the stash is used for values set with VecSetValues()
1039: Run with the option -info and look for output of the form
1040: VecAssemblyBegin_MPIXXX:Stash has MM entries, uses nn mallocs.
1041: to determine the appropriate value, MM, to use for size and
1042: VecAssemblyBegin_MPIXXX:Block-Stash has BMM entries, uses nn mallocs.
1043: to determine the value, BMM to use for bsize
1045: Concepts: vector^stash
1046: Concepts: stash^vector
1048: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked(), VecStashView()
1050: @*/
1051: PetscErrorCode VecStashSetInitialSize(Vec vec,PetscInt size,PetscInt bsize)
1052: {
1057: VecStashSetInitialSize_Private(&vec->stash,size);
1058: VecStashSetInitialSize_Private(&vec->bstash,bsize);
1059: return(0);
1060: }
1062: /*@
1063: VecConjugate - Conjugates a vector.
1065: Logically Collective on Vec
1067: Input Parameters:
1068: . x - the vector
1070: Level: intermediate
1072: Concepts: vector^conjugate
1074: @*/
1075: PetscErrorCode VecConjugate(Vec x)
1076: {
1077: #if defined(PETSC_USE_COMPLEX)
1083: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1084: (*x->ops->conjugate)(x);
1085: /* we need to copy norms here */
1086: PetscObjectStateIncrease((PetscObject)x);
1087: return(0);
1088: #else
1089: return(0);
1090: #endif
1091: }
1093: /*@
1094: VecPointwiseMult - Computes the componentwise multiplication w = x*y.
1096: Logically Collective on Vec
1098: Input Parameters:
1099: . x, y - the vectors
1101: Output Parameter:
1102: . w - the result
1104: Level: advanced
1106: Notes:
1107: any subset of the x, y, and w may be the same vector.
1109: Concepts: vector^pointwise multiply
1111: .seealso: VecPointwiseDivide(), VecPointwiseMax(), VecPointwiseMin(), VecPointwiseMaxAbs(), VecMaxPointwiseDivide()
1112: @*/
1113: PetscErrorCode VecPointwiseMult(Vec w, Vec x,Vec y)
1114: {
1126: VecCheckSameSize(w,1,x,2);
1127: VecCheckSameSize(w,2,y,3);
1128: PetscLogEventBegin(VEC_PointwiseMult,x,y,w,0);
1129: (*w->ops->pointwisemult)(w,x,y);
1130: PetscLogEventEnd(VEC_PointwiseMult,x,y,w,0);
1131: PetscObjectStateIncrease((PetscObject)w);
1132: return(0);
1133: }
1135: /*@
1136: VecSetRandom - Sets all components of a vector to random numbers.
1138: Logically Collective on Vec
1140: Input Parameters:
1141: + x - the vector
1142: - rctx - the random number context, formed by PetscRandomCreate(), or NULL and
1143: it will create one internally.
1145: Output Parameter:
1146: . x - the vector
1148: Example of Usage:
1149: .vb
1150: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
1151: VecSetRandom(x,rctx);
1152: PetscRandomDestroy(rctx);
1153: .ve
1155: Level: intermediate
1157: Concepts: vector^setting to random
1158: Concepts: random^vector
1160: .seealso: VecSet(), VecSetValues(), PetscRandomCreate(), PetscRandomDestroy()
1161: @*/
1162: PetscErrorCode VecSetRandom(Vec x,PetscRandom rctx)
1163: {
1165: PetscRandom randObj = NULL;
1171: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1173: if (!rctx) {
1174: MPI_Comm comm;
1175: PetscObjectGetComm((PetscObject)x,&comm);
1176: PetscRandomCreate(comm,&randObj);
1177: PetscRandomSetFromOptions(randObj);
1178: rctx = randObj;
1179: }
1181: PetscLogEventBegin(VEC_SetRandom,x,rctx,0,0);
1182: (*x->ops->setrandom)(x,rctx);
1183: PetscLogEventEnd(VEC_SetRandom,x,rctx,0,0);
1185: PetscRandomDestroy(&randObj);
1186: PetscObjectStateIncrease((PetscObject)x);
1187: return(0);
1188: }
1190: /*@
1191: VecZeroEntries - puts a 0.0 in each element of a vector
1193: Logically Collective on Vec
1195: Input Parameter:
1196: . vec - The vector
1198: Level: beginner
1200: Developer Note: This routine does not need to exist since the exact functionality is obtained with
1201: VecSet(vec,0); I guess someone added it to mirror the functionality of MatZeroEntries() but Mat is nothing
1202: like a Vec (one is an operator and one is an element of a vector space, yeah yeah dual blah blah blah) so
1203: this routine should not exist.
1205: .keywords: Vec, set, options, database
1206: .seealso: VecCreate(), VecSetOptionsPrefix(), VecSet(), VecSetValues()
1207: @*/
1208: PetscErrorCode VecZeroEntries(Vec vec)
1209: {
1213: VecSet(vec,0);
1214: return(0);
1215: }
1217: /*
1218: VecSetTypeFromOptions_Private - Sets the type of vector from user options. Defaults to a PETSc sequential vector on one
1219: processor and a PETSc MPI vector on more than one processor.
1221: Collective on Vec
1223: Input Parameter:
1224: . vec - The vector
1226: Level: intermediate
1228: .keywords: Vec, set, options, database, type
1229: .seealso: VecSetFromOptions(), VecSetType()
1230: */
1231: static PetscErrorCode VecSetTypeFromOptions_Private(PetscOptionItems *PetscOptionsObject,Vec vec)
1232: {
1233: PetscBool opt;
1234: VecType defaultType;
1235: char typeName[256];
1236: PetscMPIInt size;
1240: if (((PetscObject)vec)->type_name) defaultType = ((PetscObject)vec)->type_name;
1241: else {
1242: MPI_Comm_size(PetscObjectComm((PetscObject)vec), &size);
1243: if (size > 1) defaultType = VECMPI;
1244: else defaultType = VECSEQ;
1245: }
1247: VecRegisterAll();
1248: PetscOptionsFList("-vec_type","Vector type","VecSetType",VecList,defaultType,typeName,256,&opt);
1249: if (opt) {
1250: VecSetType(vec, typeName);
1251: } else {
1252: VecSetType(vec, defaultType);
1253: }
1254: return(0);
1255: }
1257: /*@
1258: VecSetFromOptions - Configures the vector from the options database.
1260: Collective on Vec
1262: Input Parameter:
1263: . vec - The vector
1265: Notes:
1266: To see all options, run your program with the -help option, or consult the users manual.
1267: Must be called after VecCreate() but before the vector is used.
1269: Level: beginner
1271: Concepts: vectors^setting options
1272: Concepts: vectors^setting type
1274: .keywords: Vec, set, options, database
1275: .seealso: VecCreate(), VecSetOptionsPrefix()
1276: @*/
1277: PetscErrorCode VecSetFromOptions(Vec vec)
1278: {
1284: PetscObjectOptionsBegin((PetscObject)vec);
1285: /* Handle vector type options */
1286: VecSetTypeFromOptions_Private(PetscOptionsObject,vec);
1288: /* Handle specific vector options */
1289: if (vec->ops->setfromoptions) {
1290: (*vec->ops->setfromoptions)(PetscOptionsObject,vec);
1291: }
1293: /* process any options handlers added with PetscObjectAddOptionsHandler() */
1294: PetscObjectProcessOptionsHandlers(PetscOptionsObject,(PetscObject)vec);
1295: PetscOptionsEnd();
1296: return(0);
1297: }
1299: /*@
1300: VecSetSizes - Sets the local and global sizes, and checks to determine compatibility
1302: Collective on Vec
1304: Input Parameters:
1305: + v - the vector
1306: . n - the local size (or PETSC_DECIDE to have it set)
1307: - N - the global size (or PETSC_DECIDE)
1309: Notes:
1310: n and N cannot be both PETSC_DECIDE
1311: If one processor calls this with N of PETSC_DECIDE then all processors must, otherwise the program will hang.
1313: Level: intermediate
1315: .seealso: VecGetSize(), PetscSplitOwnership()
1316: @*/
1317: PetscErrorCode VecSetSizes(Vec v, PetscInt n, PetscInt N)
1318: {
1324: 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);
1325: 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);
1326: v->map->n = n;
1327: v->map->N = N;
1328: if (v->ops->create) {
1329: (*v->ops->create)(v);
1330: v->ops->create = 0;
1331: }
1332: return(0);
1333: }
1335: /*@
1336: VecSetBlockSize - Sets the blocksize for future calls to VecSetValuesBlocked()
1337: and VecSetValuesBlockedLocal().
1339: Logically Collective on Vec
1341: Input Parameter:
1342: + v - the vector
1343: - bs - the blocksize
1345: Notes:
1346: All vectors obtained by VecDuplicate() inherit the same blocksize.
1348: Level: advanced
1350: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecGetBlockSize()
1352: Concepts: block size^vectors
1353: @*/
1354: PetscErrorCode VecSetBlockSize(Vec v,PetscInt bs)
1355: {
1360: if (bs < 0 || bs == v->map->bs) return(0);
1362: PetscLayoutSetBlockSize(v->map,bs);
1363: v->bstash.bs = bs; /* use the same blocksize for the vec's block-stash */
1364: return(0);
1365: }
1367: /*@
1368: VecGetBlockSize - Gets the blocksize for the vector, i.e. what is used for VecSetValuesBlocked()
1369: and VecSetValuesBlockedLocal().
1371: Not Collective
1373: Input Parameter:
1374: . v - the vector
1376: Output Parameter:
1377: . bs - the blocksize
1379: Notes:
1380: All vectors obtained by VecDuplicate() inherit the same blocksize.
1382: Level: advanced
1384: .seealso: VecSetValuesBlocked(), VecSetLocalToGlobalMapping(), VecSetBlockSize()
1386: Concepts: vector^block size
1387: Concepts: block^vector
1389: @*/
1390: PetscErrorCode VecGetBlockSize(Vec v,PetscInt *bs)
1391: {
1397: PetscLayoutGetBlockSize(v->map,bs);
1398: return(0);
1399: }
1401: /*@C
1402: VecSetOptionsPrefix - Sets the prefix used for searching for all
1403: Vec options in the database.
1405: Logically Collective on Vec
1407: Input Parameter:
1408: + v - the Vec context
1409: - prefix - the prefix to prepend to all option names
1411: Notes:
1412: A hyphen (-) must NOT be given at the beginning of the prefix name.
1413: The first character of all runtime options is AUTOMATICALLY the hyphen.
1415: Level: advanced
1417: .keywords: Vec, set, options, prefix, database
1419: .seealso: VecSetFromOptions()
1420: @*/
1421: PetscErrorCode VecSetOptionsPrefix(Vec v,const char prefix[])
1422: {
1427: PetscObjectSetOptionsPrefix((PetscObject)v,prefix);
1428: return(0);
1429: }
1431: /*@C
1432: VecAppendOptionsPrefix - Appends to the prefix used for searching for all
1433: Vec options in the database.
1435: Logically Collective on Vec
1437: Input Parameters:
1438: + v - the Vec context
1439: - prefix - the prefix to prepend to all option names
1441: Notes:
1442: A hyphen (-) must NOT be given at the beginning of the prefix name.
1443: The first character of all runtime options is AUTOMATICALLY the hyphen.
1445: Level: advanced
1447: .keywords: Vec, append, options, prefix, database
1449: .seealso: VecGetOptionsPrefix()
1450: @*/
1451: PetscErrorCode VecAppendOptionsPrefix(Vec v,const char prefix[])
1452: {
1457: PetscObjectAppendOptionsPrefix((PetscObject)v,prefix);
1458: return(0);
1459: }
1461: /*@C
1462: VecGetOptionsPrefix - Sets the prefix used for searching for all
1463: Vec options in the database.
1465: Not Collective
1467: Input Parameter:
1468: . v - the Vec context
1470: Output Parameter:
1471: . prefix - pointer to the prefix string used
1473: Notes:
1474: On the fortran side, the user should pass in a string 'prefix' of
1475: sufficient length to hold the prefix.
1477: Level: advanced
1479: .keywords: Vec, get, options, prefix, database
1481: .seealso: VecAppendOptionsPrefix()
1482: @*/
1483: PetscErrorCode VecGetOptionsPrefix(Vec v,const char *prefix[])
1484: {
1489: PetscObjectGetOptionsPrefix((PetscObject)v,prefix);
1490: return(0);
1491: }
1493: /*@
1494: VecSetUp - Sets up the internal vector data structures for the later use.
1496: Collective on Vec
1498: Input Parameters:
1499: . v - the Vec context
1501: Notes:
1502: For basic use of the Vec classes the user need not explicitly call
1503: VecSetUp(), since these actions will happen automatically.
1505: Level: advanced
1507: .keywords: Vec, setup
1509: .seealso: VecCreate(), VecDestroy()
1510: @*/
1511: PetscErrorCode VecSetUp(Vec v)
1512: {
1513: PetscMPIInt size;
1518: if (!((PetscObject)v)->type_name) {
1519: MPI_Comm_size(PetscObjectComm((PetscObject)v), &size);
1520: if (size == 1) {
1521: VecSetType(v, VECSEQ);
1522: } else {
1523: VecSetType(v, VECMPI);
1524: }
1525: }
1526: return(0);
1527: }
1529: /*
1530: These currently expose the PetscScalar/PetscReal in updating the
1531: cached norm. If we push those down into the implementation these
1532: will become independent of PetscScalar/PetscReal
1533: */
1535: /*@
1536: VecCopy - Copies a vector. y <- x
1538: Logically Collective on Vec
1540: Input Parameter:
1541: . x - the vector
1543: Output Parameter:
1544: . y - the copy
1546: Notes:
1547: For default parallel PETSc vectors, both x and y must be distributed in
1548: the same manner; local copies are done.
1550: Developer Notes:
1552: of the vectors to be sequential and one to be parallel so long as both have the same
1553: local sizes. This is used in some internal functions in PETSc.
1555: Level: beginner
1557: .seealso: VecDuplicate()
1558: @*/
1559: PetscErrorCode VecCopy(Vec x,Vec y)
1560: {
1561: PetscBool flgs[4];
1562: PetscReal norms[4] = {0.0,0.0,0.0,0.0};
1564: PetscInt i;
1571: if (x == y) return(0);
1572: VecCheckSameLocalSize(x,1,y,2);
1573: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1574: VecSetErrorIfLocked(y,2);
1576: #if !defined(PETSC_USE_MIXED_PRECISION)
1577: for (i=0; i<4; i++) {
1578: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],norms[i],flgs[i]);
1579: }
1580: #endif
1582: PetscLogEventBegin(VEC_Copy,x,y,0,0);
1583: #if defined(PETSC_USE_MIXED_PRECISION)
1584: extern PetscErrorCode VecGetArray(Vec,double**);
1585: extern PetscErrorCode VecRestoreArray(Vec,double**);
1586: extern PetscErrorCode VecGetArray(Vec,float**);
1587: extern PetscErrorCode VecRestoreArray(Vec,float**);
1588: extern PetscErrorCode VecGetArrayRead(Vec,const double**);
1589: extern PetscErrorCode VecRestoreArrayRead(Vec,const double**);
1590: extern PetscErrorCode VecGetArrayRead(Vec,const float**);
1591: extern PetscErrorCode VecRestoreArrayRead(Vec,const float**);
1592: if ((((PetscObject)x)->precision == PETSC_PRECISION_SINGLE) && (((PetscObject)y)->precision == PETSC_PRECISION_DOUBLE)) {
1593: PetscInt i,n;
1594: const float *xx;
1595: double *yy;
1596: VecGetArrayRead(x,&xx);
1597: VecGetArray(y,&yy);
1598: VecGetLocalSize(x,&n);
1599: for (i=0; i<n; i++) yy[i] = xx[i];
1600: VecRestoreArrayRead(x,&xx);
1601: VecRestoreArray(y,&yy);
1602: } else if ((((PetscObject)x)->precision == PETSC_PRECISION_DOUBLE) && (((PetscObject)y)->precision == PETSC_PRECISION_SINGLE)) {
1603: PetscInt i,n;
1604: float *yy;
1605: const double *xx;
1606: VecGetArrayRead(x,&xx);
1607: VecGetArray(y,&yy);
1608: VecGetLocalSize(x,&n);
1609: for (i=0; i<n; i++) yy[i] = (float) xx[i];
1610: VecRestoreArrayRead(x,&xx);
1611: VecRestoreArray(y,&yy);
1612: } else {
1613: (*x->ops->copy)(x,y);
1614: }
1615: #else
1616: (*x->ops->copy)(x,y);
1617: #endif
1619: PetscObjectStateIncrease((PetscObject)y);
1620: #if !defined(PETSC_USE_MIXED_PRECISION)
1621: for (i=0; i<4; i++) {
1622: if (flgs[i]) {
1623: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],norms[i]);
1624: }
1625: }
1626: #endif
1628: PetscLogEventEnd(VEC_Copy,x,y,0,0);
1629: return(0);
1630: }
1632: /*@
1633: VecSwap - Swaps the vectors x and y.
1635: Logically Collective on Vec
1637: Input Parameters:
1638: . x, y - the vectors
1640: Level: advanced
1642: Concepts: vector^swapping values
1644: @*/
1645: PetscErrorCode VecSwap(Vec x,Vec y)
1646: {
1647: PetscReal normxs[4]={0.0,0.0,0.0,0.0},normys[4]={0.0,0.0,0.0,0.0};
1648: PetscBool flgxs[4],flgys[4];
1650: PetscInt i;
1658: VecCheckSameSize(x,1,y,2);
1659: if (x->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1660: if (y->stash.insertmode != NOT_SET_VALUES) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled vector");
1662: PetscLogEventBegin(VEC_Swap,x,y,0,0);
1663: for (i=0; i<4; i++) {
1664: PetscObjectComposedDataGetReal((PetscObject)x,NormIds[i],normxs[i],flgxs[i]);
1665: PetscObjectComposedDataGetReal((PetscObject)y,NormIds[i],normys[i],flgys[i]);
1666: }
1667: (*x->ops->swap)(x,y);
1668: PetscObjectStateIncrease((PetscObject)x);
1669: PetscObjectStateIncrease((PetscObject)y);
1670: for (i=0; i<4; i++) {
1671: if (flgxs[i]) {
1672: PetscObjectComposedDataSetReal((PetscObject)y,NormIds[i],normxs[i]);
1673: }
1674: if (flgys[i]) {
1675: PetscObjectComposedDataSetReal((PetscObject)x,NormIds[i],normys[i]);
1676: }
1677: }
1678: PetscLogEventEnd(VEC_Swap,x,y,0,0);
1679: return(0);
1680: }
1682: /*
1683: VecStashViewFromOptions - Processes command line options to determine if/how an VecStash object is to be viewed.
1685: Collective on VecStash
1687: Input Parameters:
1688: + obj - the VecStash object
1689: . bobj - optional other object that provides the prefix
1690: - optionname - option to activate viewing
1692: Level: intermediate
1694: Developer Note: This cannot use PetscObjectViewFromOptions() because it takes a Vec as an argument but does not use VecView
1696: */
1697: PetscErrorCode VecStashViewFromOptions(Vec obj,PetscObject bobj,const char optionname[])
1698: {
1699: PetscErrorCode ierr;
1700: PetscViewer viewer;
1701: PetscBool flg;
1702: PetscViewerFormat format;
1703: char *prefix;
1706: prefix = bobj ? bobj->prefix : ((PetscObject)obj)->prefix;
1707: PetscOptionsGetViewer(PetscObjectComm((PetscObject)obj),((PetscObject)obj)->options,prefix,optionname,&viewer,&format,&flg);
1708: if (flg) {
1709: PetscViewerPushFormat(viewer,format);
1710: VecStashView(obj,viewer);
1711: PetscViewerPopFormat(viewer);
1712: PetscViewerDestroy(&viewer);
1713: }
1714: return(0);
1715: }
1717: /*@
1718: VecStashView - Prints the entries in the vector stash and block stash.
1720: Collective on Vec
1722: Input Parameters:
1723: + v - the vector
1724: - viewer - the viewer
1726: Level: advanced
1728: Concepts: vector^stash
1729: Concepts: stash^vector
1731: .seealso: VecSetBlockSize(), VecSetValues(), VecSetValuesBlocked()
1733: @*/
1734: PetscErrorCode VecStashView(Vec v,PetscViewer viewer)
1735: {
1737: PetscMPIInt rank;
1738: PetscInt i,j;
1739: PetscBool match;
1740: VecStash *s;
1741: PetscScalar val;
1748: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&match);
1749: if (!match) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"Stash viewer only works with ASCII viewer not %s\n",((PetscObject)v)->type_name);
1750: PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
1751: MPI_Comm_rank(PetscObjectComm((PetscObject)v),&rank);
1752: s = &v->bstash;
1754: /* print block stash */
1755: PetscViewerASCIIPushSynchronized(viewer);
1756: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector Block stash size %D block size %D\n",rank,s->n,s->bs);
1757: for (i=0; i<s->n; i++) {
1758: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D ",rank,s->idx[i]);
1759: for (j=0; j<s->bs; j++) {
1760: val = s->array[i*s->bs+j];
1761: #if defined(PETSC_USE_COMPLEX)
1762: PetscViewerASCIISynchronizedPrintf(viewer,"(%18.16e %18.16e) ",PetscRealPart(val),PetscImaginaryPart(val));
1763: #else
1764: PetscViewerASCIISynchronizedPrintf(viewer,"%18.16e ",val);
1765: #endif
1766: }
1767: PetscViewerASCIISynchronizedPrintf(viewer,"\n");
1768: }
1769: PetscViewerFlush(viewer);
1771: s = &v->stash;
1773: /* print basic stash */
1774: PetscViewerASCIISynchronizedPrintf(viewer,"[%d]Vector stash size %D\n",rank,s->n);
1775: for (i=0; i<s->n; i++) {
1776: val = s->array[i];
1777: #if defined(PETSC_USE_COMPLEX)
1778: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D (%18.16e %18.16e) ",rank,s->idx[i],PetscRealPart(val),PetscImaginaryPart(val));
1779: #else
1780: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] Element %D %18.16e\n",rank,s->idx[i],val);
1781: #endif
1782: }
1783: PetscViewerFlush(viewer);
1784: PetscViewerASCIIPopSynchronized(viewer);
1785: PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
1786: return(0);
1787: }
1789: PetscErrorCode PetscOptionsGetVec(PetscOptions options,const char prefix[],const char key[],Vec v,PetscBool *set)
1790: {
1791: PetscInt i,N,rstart,rend;
1793: PetscScalar *xx;
1794: PetscReal *xreal;
1795: PetscBool iset;
1798: VecGetOwnershipRange(v,&rstart,&rend);
1799: VecGetSize(v,&N);
1800: PetscCalloc1(N,&xreal);
1801: PetscOptionsGetRealArray(options,prefix,key,xreal,&N,&iset);
1802: if (iset) {
1803: VecGetArray(v,&xx);
1804: for (i=rstart; i<rend; i++) xx[i-rstart] = xreal[i];
1805: VecRestoreArray(v,&xx);
1806: }
1807: PetscFree(xreal);
1808: if (set) *set = iset;
1809: return(0);
1810: }
1812: /*@
1813: VecGetLayout - get PetscLayout describing vector layout
1815: Not Collective
1817: Input Arguments:
1818: . x - the vector
1820: Output Arguments:
1821: . map - the layout
1823: Level: developer
1825: .seealso: VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1826: @*/
1827: PetscErrorCode VecGetLayout(Vec x,PetscLayout *map)
1828: {
1832: *map = x->map;
1833: return(0);
1834: }
1836: /*@
1837: VecSetLayout - set PetscLayout describing vector layout
1839: Not Collective
1841: Input Arguments:
1842: + x - the vector
1843: - map - the layout
1845: Notes:
1846: It is normally only valid to replace the layout with a layout known to be equivalent.
1848: Level: developer
1850: .seealso: VecGetLayout(), VecGetSizes(), VecGetOwnershipRange(), VecGetOwnershipRanges()
1851: @*/
1852: PetscErrorCode VecSetLayout(Vec x,PetscLayout map)
1853: {
1858: PetscLayoutReference(map,&x->map);
1859: return(0);
1860: }
1862: PetscErrorCode VecSetInf(Vec xin)
1863: {
1864: PetscInt i,n = xin->map->n;
1865: PetscScalar *xx;
1866: PetscScalar zero=0.0,one=1.0,inf=one/zero;
1870: VecGetArray(xin,&xx);
1871: for (i=0; i<n; i++) xx[i] = inf;
1872: VecRestoreArray(xin,&xx);
1873: return(0);
1874: }