Actual source code: vecutil.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
7: SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9: */
11: #include <slepc/private/vecimplslepc.h>
13: /*@
14: VecNormalizeComplex - Normalizes a possibly complex vector by the 2-norm.
16: Collective
18: Input Parameters:
19: + xr - the real part of the vector (overwritten on output)
20: . xi - the imaginary part of the vector (not referenced if iscomplex is false)
21: - iscomplex - a flag indicating if the vector is complex
23: Output Parameter:
24: . norm - the vector norm before normalization (can be set to NULL)
26: Level: developer
28: .seealso: `BVNormalize()`
29: @*/
30: PetscErrorCode VecNormalizeComplex(Vec xr,Vec xi,PetscBool iscomplex,PetscReal *norm)
31: {
32: #if !defined(PETSC_USE_COMPLEX)
33: PetscReal normr,normi,alpha;
34: #endif
36: PetscFunctionBegin;
38: #if !defined(PETSC_USE_COMPLEX)
39: if (iscomplex) {
41: PetscCall(VecNormBegin(xr,NORM_2,&normr));
42: PetscCall(VecNormBegin(xi,NORM_2,&normi));
43: PetscCall(VecNormEnd(xr,NORM_2,&normr));
44: PetscCall(VecNormEnd(xi,NORM_2,&normi));
45: alpha = SlepcAbsEigenvalue(normr,normi);
46: if (norm) *norm = alpha;
47: alpha = 1.0 / alpha;
48: PetscCall(VecScale(xr,alpha));
49: PetscCall(VecScale(xi,alpha));
50: } else
51: #endif
52: PetscCall(VecNormalize(xr,norm));
53: PetscFunctionReturn(PETSC_SUCCESS);
54: }
56: static PetscErrorCode VecCheckOrthogonality_Private(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev,PetscBool norm)
57: {
58: PetscInt i,j;
59: PetscScalar *vals;
60: PetscBool isascii;
61: Vec w;
63: PetscFunctionBegin;
64: if (!lev) {
65: if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)*V),&viewer));
67: PetscCheckSameComm(*V,1,viewer,6);
68: PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
69: if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
70: }
72: PetscCall(PetscMalloc1(nv,&vals));
73: if (B) PetscCall(VecDuplicate(V[0],&w));
74: if (lev) *lev = 0.0;
75: for (i=0;i<nw;i++) {
76: if (B) {
77: if (W) PetscCall(MatMult(B,W[i],w));
78: else PetscCall(MatMult(B,V[i],w));
79: } else {
80: if (W) w = W[i];
81: else w = V[i];
82: }
83: PetscCall(VecMDot(w,nv,V,vals));
84: for (j=0;j<nv;j++) {
85: if (lev) {
86: if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
87: else if (norm) {
88: if (PetscRealPart(vals[j])<0.0) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]+PetscRealConstant(1.0))); /* indefinite case */
89: else *lev = PetscMax(*lev,PetscAbsScalar(vals[j]-PetscRealConstant(1.0)));
90: }
91: } else {
92: #if !defined(PETSC_USE_COMPLEX)
93: PetscCall(PetscViewerASCIIPrintf(viewer," %12g ",(double)vals[j]));
94: #else
95: PetscCall(PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j])));
96: #endif
97: }
98: }
99: if (!lev) PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
100: }
101: PetscCall(PetscFree(vals));
102: if (B) PetscCall(VecDestroy(&w));
103: PetscFunctionReturn(PETSC_SUCCESS);
104: }
106: /*@
107: VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
108: of a set of vectors.
110: Collective
112: Input Parameters:
113: + V - a set of vectors
114: . nv - number of V vectors
115: . W - an alternative set of vectors (optional)
116: . nw - number of W vectors
117: . B - Hermitian matrix defining the inner product (optional)
118: - viewer - optional visualization context
120: Output Parameter:
121: . lev - level of orthogonality (optional)
123: Notes:
124: This function computes W'*V and prints the result. It is intended to check
125: the level of bi-orthogonality of the vectors in the two sets. If W is equal
126: to NULL then V is used, thus checking the orthogonality of the V vectors.
128: If matrix B is provided then the check uses the B-inner product, W'*B*V,
129: where B is assumed to be Hermitian.
131: If V, W represent eigenvectors computed by SLEPc, this function will not work
132: correctly if one of the eigenvalues is complex when running with real scalars.
134: If lev is not NULL, it will contain the maximum entry of matrix
135: W'*V - I (in absolute value) omitting the diagonal. Otherwise, the matrix W'*V
136: is printed.
138: Level: developer
140: .seealso: `VecCheckOrthonormality()`
141: @*/
142: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
143: {
144: PetscFunctionBegin;
145: PetscAssertPointer(V,1);
149: if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
150: if (W) {
151: PetscAssertPointer(W,3);
153: PetscCheckSameComm(*V,1,*W,3);
154: }
155: PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
156: PetscFunctionReturn(PETSC_SUCCESS);
157: }
159: /*@
160: VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
161: of a set of vectors.
163: Collective
165: Input Parameters:
166: + V - a set of vectors
167: . nv - number of V vectors
168: . W - an alternative set of vectors (optional)
169: . nw - number of W vectors
170: . B - Hermitian matrix defining the inner product (optional)
171: - viewer - optional visualization context
173: Output Parameter:
174: . lev - level of orthogonality (optional)
176: Notes:
177: This function is equivalent to VecCheckOrthonormality(), but in addition it checks
178: that the diagonal of W'*V (or W'*B*V) is equal to all ones.
180: Level: developer
182: .seealso: `VecCheckOrthogonality()`
183: @*/
184: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
185: {
186: PetscFunctionBegin;
187: PetscAssertPointer(V,1);
191: if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
192: if (W) {
193: PetscAssertPointer(W,3);
195: PetscCheckSameComm(*V,1,*W,3);
196: }
197: PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: /*@
202: VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
203: but without internal array.
205: Collective
207: Input Parameters:
208: . v - a vector to mimic
210: Output Parameter:
211: . newv - location to put new vector
213: Note:
214: This is similar to VecDuplicate(), but the new vector does not have an internal
215: array, so the intended usage is with VecPlaceArray().
217: Level: developer
219: .seealso: `MatCreateVecsEmpty()`
220: @*/
221: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
222: {
223: PetscBool standard,cuda,hip,mpi;
224: PetscInt N,nloc,bs;
226: PetscFunctionBegin;
228: PetscAssertPointer(newv,2);
231: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
232: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
233: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&hip,VECSEQHIP,VECMPIHIP,""));
234: if (standard || cuda || hip) {
235: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
236: PetscCall(VecGetLocalSize(v,&nloc));
237: PetscCall(VecGetSize(v,&N));
238: PetscCall(VecGetBlockSize(v,&bs));
239: if (cuda) {
240: #if defined(PETSC_HAVE_CUDA)
241: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
242: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
243: #endif
244: } else if (hip) {
245: #if defined(PETSC_HAVE_HIP)
246: if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
247: else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
248: #endif
249: } else {
250: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
251: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
252: }
253: } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
254: PetscFunctionReturn(PETSC_SUCCESS);
255: }
257: /*@
258: VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
260: Logically Collective
262: Input Parameters:
263: + v - the vector to be filled with random values
264: . rctx - the random number context (can be NULL)
265: . w1 - first work vector (can be NULL)
266: - w2 - second work vector (can be NULL)
268: Notes:
269: Fills the two work vectors with uniformly distributed random values (VecSetRandom)
270: and then applies the Box-Muller transform to get normally distributed values on v.
272: Level: developer
274: .seealso: `VecSetRandom()`
275: @*/
276: PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
277: {
278: const PetscScalar *x,*y;
279: PetscScalar *z;
280: PetscInt n,i;
281: PetscRandom rand=NULL;
282: Vec v1=NULL,v2=NULL;
284: PetscFunctionBegin;
291: if (!rctx) {
292: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
293: PetscCall(PetscRandomSetFromOptions(rand));
294: rctx = rand;
295: }
296: if (!w1) {
297: PetscCall(VecDuplicate(v,&v1));
298: w1 = v1;
299: }
300: if (!w2) {
301: PetscCall(VecDuplicate(v,&v2));
302: w2 = v2;
303: }
304: PetscCheckSameTypeAndComm(v,1,w1,3);
305: PetscCheckSameTypeAndComm(v,1,w2,4);
307: PetscCall(VecSetRandom(w1,rctx));
308: PetscCall(VecSetRandom(w2,rctx));
309: PetscCall(VecGetLocalSize(v,&n));
310: PetscCall(VecGetArrayWrite(v,&z));
311: PetscCall(VecGetArrayRead(w1,&x));
312: PetscCall(VecGetArrayRead(w2,&y));
313: for (i=0;i<n;i++) {
314: #if defined(PETSC_USE_COMPLEX)
315: z[i] = PetscCMPLX(PetscSqrtReal(-2.0*PetscLogReal(PetscRealPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscRealPart(y[i])),PetscSqrtReal(-2.0*PetscLogReal(PetscImaginaryPart(x[i])))*PetscCosReal(2.0*PETSC_PI*PetscImaginaryPart(y[i])));
316: #else
317: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
318: #endif
319: }
320: PetscCall(VecRestoreArrayWrite(v,&z));
321: PetscCall(VecRestoreArrayRead(w1,&x));
322: PetscCall(VecRestoreArrayRead(w2,&y));
324: PetscCall(VecDestroy(&v1));
325: PetscCall(VecDestroy(&v2));
326: if (!rctx) PetscCall(PetscRandomDestroy(&rand));
327: PetscFunctionReturn(PETSC_SUCCESS);
328: }