Actual source code: vecutil.c
slepc-3.21.0 2024-03-30
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(MatMultTranspose(B,W[i],w));
78: else PetscCall(MatMultTranspose(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 - 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.
130: If lev is not NULL, it will contain the maximum entry of matrix
131: W'*V - I (in absolute value) omitting the diagonal. Otherwise, the matrix W'*V
132: is printed.
134: Level: developer
136: .seealso: VecCheckOrthonormality()
137: @*/
138: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
139: {
140: PetscFunctionBegin;
141: PetscAssertPointer(V,1);
145: if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
146: if (W) {
147: PetscAssertPointer(W,3);
149: PetscCheckSameComm(*V,1,*W,3);
150: }
151: PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: /*@
156: VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
157: of a set of vectors.
159: Collective
161: Input Parameters:
162: + V - a set of vectors
163: . nv - number of V vectors
164: . W - an alternative set of vectors (optional)
165: . nw - number of W vectors
166: . B - matrix defining the inner product (optional)
167: - viewer - optional visualization context
169: Output Parameter:
170: . lev - level of orthogonality (optional)
172: Notes:
173: This function is equivalent to VecCheckOrthonormality(), but in addition it checks
174: that the diagonal of W'*V (or W'*B*V) is equal to all ones.
176: Level: developer
178: .seealso: VecCheckOrthogonality()
179: @*/
180: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
181: {
182: PetscFunctionBegin;
183: PetscAssertPointer(V,1);
187: if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
188: if (W) {
189: PetscAssertPointer(W,3);
191: PetscCheckSameComm(*V,1,*W,3);
192: }
193: PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@
198: VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
199: but without internal array.
201: Collective
203: Input Parameters:
204: . v - a vector to mimic
206: Output Parameter:
207: . newv - location to put new vector
209: Note:
210: This is similar to VecDuplicate(), but the new vector does not have an internal
211: array, so the intended usage is with VecPlaceArray().
213: Level: developer
215: .seealso: MatCreateVecsEmpty()
216: @*/
217: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
218: {
219: PetscBool standard,cuda,mpi;
220: PetscInt N,nloc,bs;
222: PetscFunctionBegin;
224: PetscAssertPointer(newv,2);
227: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
228: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
229: if (standard || cuda) {
230: PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,""));
231: PetscCall(VecGetLocalSize(v,&nloc));
232: PetscCall(VecGetSize(v,&N));
233: PetscCall(VecGetBlockSize(v,&bs));
234: if (cuda) {
235: #if defined(PETSC_HAVE_CUDA)
236: if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
237: else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
238: #endif
239: } else {
240: if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
241: else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
242: }
243: } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*@
248: VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
250: Logically Collective
252: Input Parameters:
253: + v - the vector to be filled with random values
254: . rctx - the random number context (can be NULL)
255: . w1 - first work vector (can be NULL)
256: - w2 - second work vector (can be NULL)
258: Notes:
259: Fills the two work vectors with uniformly distributed random values (VecSetRandom)
260: and then applies the Box-Muller transform to get normally distributed values on v.
262: Level: developer
264: .seealso: VecSetRandom()
265: @*/
266: PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
267: {
268: const PetscScalar *x,*y;
269: PetscScalar *z;
270: PetscInt n,i;
271: PetscRandom rand=NULL;
272: Vec v1=NULL,v2=NULL;
274: PetscFunctionBegin;
281: if (!rctx) {
282: PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
283: PetscCall(PetscRandomSetFromOptions(rand));
284: rctx = rand;
285: }
286: if (!w1) {
287: PetscCall(VecDuplicate(v,&v1));
288: w1 = v1;
289: }
290: if (!w2) {
291: PetscCall(VecDuplicate(v,&v2));
292: w2 = v2;
293: }
294: PetscCheckSameTypeAndComm(v,1,w1,3);
295: PetscCheckSameTypeAndComm(v,1,w2,4);
297: PetscCall(VecSetRandom(w1,rctx));
298: PetscCall(VecSetRandom(w2,rctx));
299: PetscCall(VecGetLocalSize(v,&n));
300: PetscCall(VecGetArrayWrite(v,&z));
301: PetscCall(VecGetArrayRead(w1,&x));
302: PetscCall(VecGetArrayRead(w2,&y));
303: for (i=0;i<n;i++) {
304: #if defined(PETSC_USE_COMPLEX)
305: 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])));
306: #else
307: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
308: #endif
309: }
310: PetscCall(VecRestoreArrayWrite(v,&z));
311: PetscCall(VecRestoreArrayRead(w1,&x));
312: PetscCall(VecRestoreArrayRead(w2,&y));
314: PetscCall(VecDestroy(&v1));
315: PetscCall(VecDestroy(&v2));
316: if (!rctx) PetscCall(PetscRandomDestroy(&rand));
317: PetscFunctionReturn(PETSC_SUCCESS);
318: }