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 !PetscDefined(USE_COMPLEX)
33: PetscReal normr,normi,alpha;
34: #endif
36: PetscFunctionBegin;
38: #if !PetscDefined(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 !PetscDefined(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^*BV$,
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^*BV$) 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 Parameter:
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: `VecDuplicate()`, `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 (mpi) PetscCall(VecCreateMPIWithArrayAndMemType(PetscObjectComm((PetscObject)v),cuda?PETSC_MEMTYPE_CUDA:(hip?PETSC_MEMTYPE_HIP:PETSC_MEMTYPE_HOST),bs,nloc,N,NULL,newv));
240: else PetscCall(VecCreateSeqWithArrayAndMemType(PetscObjectComm((PetscObject)v),cuda?PETSC_MEMTYPE_CUDA:(hip?PETSC_MEMTYPE_HIP:PETSC_MEMTYPE_HOST),bs,N,NULL,newv));
241: } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: /*@
246: VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
248: Logically Collective
250: Input Parameters:
251: + v - the vector to be filled with random values
252: . rctx - the random number context (can be `NULL`)
253: . w1 - first work vector (can be `NULL`)
254: - w2 - second work vector (can be `NULL`)
256: Notes:
257: Fills the two work vectors with uniformly distributed random values (`VecSetRandom()`)
258: and then applies the Box-Muller transform to get normally distributed values on `v`.
260: Level: developer
262: .seealso: `VecSetRandom()`
263: @*/
264: PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
265: {
266: const PetscScalar *x,*y;
267: PetscScalar *z;
268: PetscInt n,i;
269: PetscRandom rand=NULL;
270: Vec v1=NULL,v2=NULL;
272: PetscFunctionBegin;
279: if (!rctx) {
280: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
281: PetscCall(PetscRandomSetFromOptions(rand));
282: rctx = rand;
283: }
284: if (!w1) {
285: PetscCall(VecDuplicate(v,&v1));
286: w1 = v1;
287: }
288: if (!w2) {
289: PetscCall(VecDuplicate(v,&v2));
290: w2 = v2;
291: }
292: PetscCheckSameTypeAndComm(v,1,w1,3);
293: PetscCheckSameTypeAndComm(v,1,w2,4);
295: PetscCall(VecSetRandom(w1,rctx));
296: PetscCall(VecSetRandom(w2,rctx));
297: PetscCall(VecGetLocalSize(v,&n));
298: PetscCall(VecGetArrayWrite(v,&z));
299: PetscCall(VecGetArrayRead(w1,&x));
300: PetscCall(VecGetArrayRead(w2,&y));
301: for (i=0;i<n;i++) {
302: #if PetscDefined(USE_COMPLEX)
303: 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])));
304: #else
305: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
306: #endif
307: }
308: PetscCall(VecRestoreArrayWrite(v,&z));
309: PetscCall(VecRestoreArrayRead(w1,&x));
310: PetscCall(VecRestoreArrayRead(w2,&y));
312: PetscCall(VecDestroy(&v1));
313: PetscCall(VecDestroy(&v2));
314: PetscCall(PetscRandomDestroy(&rand));
315: PetscFunctionReturn(PETSC_SUCCESS);
316: }