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: PetscFunctionReturn(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) PetscFunctionReturn(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) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]-PetscRealConstant(1.0)));
86: } else {
87: #if !defined(PETSC_USE_COMPLEX)
88: PetscViewerASCIIPrintf(viewer," %12g ",(double)vals[j]);
89: #else
90: PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j]));
91: #endif
92: }
93: }
94: if (!lev) PetscViewerASCIIPrintf(viewer,"\n");
95: }
96: PetscFree(vals);
97: if (B) VecDestroy(&w);
98: PetscFunctionReturn(0);
99: }
101: /*@
102: VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
103: of a set of vectors.
105: Collective on V
107: Input Parameters:
108: + V - a set of vectors
109: . nv - number of V vectors
110: . W - an alternative set of vectors (optional)
111: . nw - number of W vectors
112: . B - matrix defining the inner product (optional)
113: - viewer - optional visualization context
115: Output Parameter:
116: . lev - level of orthogonality (optional)
118: Notes:
119: This function computes W'*V and prints the result. It is intended to check
120: the level of bi-orthogonality of the vectors in the two sets. If W is equal
121: to NULL then V is used, thus checking the orthogonality of the V vectors.
123: If matrix B is provided then the check uses the B-inner product, W'*B*V.
125: If lev is not NULL, it will contain the maximum entry of matrix
126: W'*V - I (in absolute value) omitting the diagonal. Otherwise, the matrix W'*V
127: is printed.
129: Level: developer
131: .seealso: VecCheckOrthonormality()
132: @*/
133: PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)134: {
139: if (nv<=0 || nw<=0) PetscFunctionReturn(0);
140: if (W) {
144: }
145: VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE);
146: PetscFunctionReturn(0);
147: }
149: /*@
150: VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
151: of a set of vectors.
153: Collective on V
155: Input Parameters:
156: + V - a set of vectors
157: . nv - number of V vectors
158: . W - an alternative set of vectors (optional)
159: . nw - number of W vectors
160: . B - matrix defining the inner product (optional)
161: - viewer - optional visualization context
163: Output Parameter:
164: . lev - level of orthogonality (optional)
166: Notes:
167: This function is equivalent to VecCheckOrthonormality(), but in addition it checks
168: that the diagonal of W'*V (or W'*B*V) is equal to all ones.
170: Level: developer
172: .seealso: VecCheckOrthogonality()
173: @*/
174: PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)175: {
180: if (nv<=0 || nw<=0) PetscFunctionReturn(0);
181: if (W) {
185: }
186: VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE);
187: PetscFunctionReturn(0);
188: }
190: /*@
191: VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
192: but without internal array.
194: Collective on v
196: Input Parameters:
197: . v - a vector to mimic
199: Output Parameter:
200: . newv - location to put new vector
202: Note:
203: This is similar to VecDuplicate(), but the new vector does not have an internal
204: array, so the intended usage is with VecPlaceArray().
206: Level: developer
208: .seealso: MatCreateVecsEmpty()
209: @*/
210: PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)211: {
212: PetscBool standard,cuda,mpi;
213: PetscInt N,nloc,bs;
219: PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,"");
220: PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,"");
221: if (standard || cuda) {
222: PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,"");
223: VecGetLocalSize(v,&nloc);
224: VecGetSize(v,&N);
225: VecGetBlockSize(v,&bs);
226: if (cuda) {
227: #if defined(PETSC_HAVE_CUDA)
228: if (mpi) VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
229: else VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
230: #endif
231: } else {
232: if (mpi) VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv);
233: else VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv);
234: }
235: } else VecDuplicate(v,newv); /* standard duplicate, with internal array */
236: PetscFunctionReturn(0);
237: }
239: /*@
240: VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
242: Logically Collective on v
244: Input Parameters:
245: + v - the vector to be filled with random values
246: . rctx - the random number context (can be NULL)
247: . w1 - first work vector (can be NULL)
248: - w2 - second work vector (can be NULL)
250: Output Parameter:
251: . v - the vector
253: Notes:
254: Fills the two work vectors with uniformly distributed random values (VecSetRandom)
255: and then applies the Box-Muller transform to get normally distributed values on v.
257: Level: developer
259: .seealso: VecSetRandom()
260: @*/
261: PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)262: {
263: const PetscScalar *x,*y;
264: PetscScalar *z;
265: PetscInt n,i;
266: PetscRandom rand=NULL;
267: Vec v1=NULL,v2=NULL;
275: if (!rctx) {
276: PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand);
277: PetscRandomSetFromOptions(rand);
278: rctx = rand;
279: }
280: if (!w1) {
281: VecDuplicate(v,&v1);
282: w1 = v1;
283: }
284: if (!w2) {
285: VecDuplicate(v,&v2);
286: w2 = v2;
287: }
291: VecSetRandom(w1,rctx);
292: VecSetRandom(w2,rctx);
293: VecGetLocalSize(v,&n);
294: VecGetArrayWrite(v,&z);
295: VecGetArrayRead(w1,&x);
296: VecGetArrayRead(w2,&y);
297: for (i=0;i<n;i++) {
298: #if defined(PETSC_USE_COMPLEX)
299: 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])));
300: #else
301: z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
302: #endif
303: }
304: VecRestoreArrayWrite(v,&z);
305: VecRestoreArrayRead(w1,&x);
306: VecRestoreArrayRead(w2,&y);
308: VecDestroy(&v1);
309: VecDestroy(&v2);
310: if (!rctx) PetscRandomDestroy(&rand);
311: PetscFunctionReturn(0);
312: }