Actual source code: vecutil.c
slepc-3.18.2 2023-01-26
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 on xr
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
37: #if !defined(PETSC_USE_COMPLEX)
38: if (iscomplex) {
40: VecNormBegin(xr,NORM_2,&normr);
41: VecNormBegin(xi,NORM_2,&normi);
42: VecNormEnd(xr,NORM_2,&normr);
43: VecNormEnd(xi,NORM_2,&normi);
44: alpha = SlepcAbsEigenvalue(normr,normi);
45: if (norm) *norm = alpha;
46: alpha = 1.0 / alpha;
47: VecScale(xr,alpha);
48: VecScale(xi,alpha);
49: } else
50: #endif
51: VecNormalize(xr,norm);
52: return 0;
53: }
55: static PetscErrorCode VecCheckOrthogonality_Private(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev,PetscBool norm)
56: {
57: PetscInt i,j;
58: PetscScalar *vals;
59: PetscBool isascii;
60: Vec w;
62: if (!lev) {
63: if (!viewer) PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)*V),&viewer);
66: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii);
67: if (!isascii) return 0;
68: }
70: PetscMalloc1(nv,&vals);
71: if (B) VecDuplicate(V[0],&w);
72: if (lev) *lev = 0.0;
73: for (i=0;i<nw;i++) {
74: if (B) {
75: if (W) MatMultTranspose(B,W[i],w);
76: else MatMultTranspose(B,V[i],w);
77: } else {
78: if (W) w = W[i];
79: else w = V[i];
80: }
81: VecMDot(w,nv,V,vals);
82: for (j=0;j<nv;j++) {
83: if (lev) {
84: if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
85: else if (norm) {
86: if (PetscRealPart(vals[j])<0.0) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]+PetscRealConstant(1.0))); /* indefinite case */
87: else *lev = PetscMax(*lev,PetscAbsScalar(vals[j]-PetscRealConstant(1.0)));
88: }
89: } else {
90: #if !defined(PETSC_USE_COMPLEX)
91: PetscViewerASCIIPrintf(viewer," %12g ",(double)vals[j]);
92: #else
93: PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j]));
94: #endif
95: }
96: }
97: if (!lev) PetscViewerASCIIPrintf(viewer,"\n");
98: }
99: PetscFree(vals);
100: if (B) VecDestroy(&w);
101: return 0;
102: }
104: /*@
105: VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
106: of a set of vectors.
108: Collective on V
110: Input Parameters:
111: + V - a set of vectors
112: . nv - number of V vectors
113: . W - an alternative set of vectors (optional)
114: . nw - number of W vectors
115: . B - matrix defining the inner product (optional)
116: - viewer - optional visualization context
118: Output Parameter:
119: . lev - level of orthogonality (optional)
121: Notes:
122: This function computes W'*V and prints the result. It is intended to check
123: the level of bi-orthogonality of the vectors in the two sets. If W is equal
124: to NULL then V is used, thus checking the orthogonality of the V vectors.
126: If matrix B is provided then the check uses the B-inner product, W'*B*V.
128: If lev is not NULL, it will contain the maximum entry of matrix
129: W'*V - I (in absolute value) omitting the diagonal. Otherwise, the matrix W'*V
130: is printed.
132: Level: developer
134: .seealso: VecCheckOrthonormality()
135: @*/
136: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
137: {
142: if (nv<=0 || nw<=0) return 0;
143: if (W) {
147: }
148: VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE);
149: return 0;
150: }
152: /*@
153: VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
154: of a set of vectors.
156: Collective on V
158: Input Parameters:
159: + V - a set of vectors
160: . nv - number of V vectors
161: . W - an alternative set of vectors (optional)
162: . nw - number of W vectors
163: . B - matrix defining the inner product (optional)
164: - viewer - optional visualization context
166: Output Parameter:
167: . lev - level of orthogonality (optional)
169: Notes:
170: This function is equivalent to VecCheckOrthonormality(), but in addition it checks
171: that the diagonal of W'*V (or W'*B*V) is equal to all ones.
173: Level: developer
175: .seealso: VecCheckOrthogonality()
176: @*/
177: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
178: {
183: if (nv<=0 || nw<=0) return 0;
184: if (W) {
188: }
189: VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE);
190: return 0;
191: }
193: /*@
194: VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
195: but without internal array.
197: Collective on v
199: Input Parameters:
200: . v - a vector to mimic
202: Output Parameter:
203: . newv - location to put new vector
205: Note:
206: This is similar to VecDuplicate(), but the new vector does not have an internal
207: array, so the intended usage is with VecPlaceArray().
209: Level: developer
211: .seealso: MatCreateVecsEmpty()
212: @*/
213: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
214: {
215: PetscBool standard,cuda,mpi;
216: PetscInt N,nloc,bs;
222: PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,"");
223: PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,"");
224: if (standard || cuda) {
225: PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,"");
226: VecGetLocalSize(v,&nloc);
227: VecGetSize(v,&N);
228: VecGetBlockSize(v,&bs);
229: if (cuda) {
230: #if defined(PETSC_HAVE_CUDA)
231: if (mpi) VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
232: else VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
233: #endif
234: } else {
235: if (mpi) VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
236: else VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
237: }
238: } else VecDuplicate(v,newv); /* standard duplicate, with internal array */
239: return 0;
240: }
242: /*@
243: VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
245: Logically Collective on v
247: Input Parameters:
248: + v - the vector to be filled with random values
249: . rctx - the random number context (can be NULL)
250: . w1 - first work vector (can be NULL)
251: - w2 - second work vector (can be NULL)
253: Output Parameter:
254: . v - the vector
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;
278: if (!rctx) {
279: PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand);
280: PetscRandomSetFromOptions(rand);
281: rctx = rand;
282: }
283: if (!w1) {
284: VecDuplicate(v,&v1);
285: w1 = v1;
286: }
287: if (!w2) {
288: VecDuplicate(v,&v2);
289: w2 = v2;
290: }
294: VecSetRandom(w1,rctx);
295: VecSetRandom(w2,rctx);
296: VecGetLocalSize(v,&n);
297: VecGetArrayWrite(v,&z);
298: VecGetArrayRead(w1,&x);
299: VecGetArrayRead(w2,&y);
300: for (i=0;i<n;i++) {
301: #if defined(PETSC_USE_COMPLEX)
302: 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])));
303: #else
304: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
305: #endif
306: }
307: VecRestoreArrayWrite(v,&z);
308: VecRestoreArrayRead(w1,&x);
309: VecRestoreArrayRead(w2,&y);
311: VecDestroy(&v1);
312: VecDestroy(&v2);
313: if (!rctx) PetscRandomDestroy(&rand);
314: return 0;
315: }