Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
11 : #include <slepc/private/vecimplslepc.h> /*I "slepcvec.h" I*/
12 :
13 : /*@
14 : VecNormalizeComplex - Normalizes a possibly complex vector by the 2-norm.
15 :
16 : Collective
17 :
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
22 :
23 : Output Parameter:
24 : . norm - the vector norm before normalization (can be set to NULL)
25 :
26 : Level: developer
27 :
28 : .seealso: BVNormalize()
29 : @*/
30 34 : PetscErrorCode VecNormalizeComplex(Vec xr,Vec xi,PetscBool iscomplex,PetscReal *norm)
31 : {
32 : #if !defined(PETSC_USE_COMPLEX)
33 : PetscReal normr,normi,alpha;
34 : #endif
35 :
36 34 : PetscFunctionBegin;
37 34 : PetscValidHeaderSpecific(xr,VEC_CLASSID,1);
38 : #if !defined(PETSC_USE_COMPLEX)
39 : if (iscomplex) {
40 : PetscValidHeaderSpecific(xi,VEC_CLASSID,2);
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 34 : PetscCall(VecNormalize(xr,norm));
53 34 : PetscFunctionReturn(PETSC_SUCCESS);
54 : }
55 :
56 243 : static PetscErrorCode VecCheckOrthogonality_Private(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev,PetscBool norm)
57 : {
58 243 : PetscInt i,j;
59 243 : PetscScalar *vals;
60 243 : PetscBool isascii;
61 243 : Vec w;
62 :
63 243 : PetscFunctionBegin;
64 243 : if (!lev) {
65 0 : if (!viewer) PetscCall(PetscViewerASCIIGetStdout(PetscObjectComm((PetscObject)*V),&viewer));
66 0 : PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,6);
67 0 : PetscCheckSameComm(*V,1,viewer,6);
68 0 : PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&isascii));
69 0 : if (!isascii) PetscFunctionReturn(PETSC_SUCCESS);
70 : }
71 :
72 243 : PetscCall(PetscMalloc1(nv,&vals));
73 243 : if (B) PetscCall(VecDuplicate(V[0],&w));
74 243 : if (lev) *lev = 0.0;
75 2967 : for (i=0;i<nw;i++) {
76 2724 : if (B) {
77 1208 : if (W) PetscCall(MatMultTranspose(B,W[i],w));
78 1202 : else PetscCall(MatMultTranspose(B,V[i],w));
79 : } else {
80 1516 : if (W) w = W[i];
81 1499 : else w = V[i];
82 : }
83 2724 : PetscCall(VecMDot(w,nv,V,vals));
84 181532 : for (j=0;j<nv;j++) {
85 178808 : if (lev) {
86 354542 : if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
87 2724 : else if (norm) {
88 2701 : if (PetscRealPart(vals[j])<0.0) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]+PetscRealConstant(1.0))); /* indefinite case */
89 4968 : 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 178808 : PetscCall(PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j])));
96 : #endif
97 : }
98 : }
99 2724 : if (!lev) PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
100 : }
101 243 : PetscCall(PetscFree(vals));
102 243 : if (B) PetscCall(VecDestroy(&w));
103 243 : PetscFunctionReturn(PETSC_SUCCESS);
104 : }
105 :
106 : /*@
107 : VecCheckOrthogonality - Checks (or prints) the level of (bi-)orthogonality
108 : of a set of vectors.
109 :
110 : Collective
111 :
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
119 :
120 : Output Parameter:
121 : . lev - level of orthogonality (optional)
122 :
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.
127 :
128 : If matrix B is provided then the check uses the B-inner product, W'*B*V.
129 :
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.
133 :
134 : Level: developer
135 :
136 : .seealso: VecCheckOrthonormality()
137 : @*/
138 5 : PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
139 : {
140 5 : PetscFunctionBegin;
141 5 : PetscAssertPointer(V,1);
142 5 : PetscValidHeaderSpecific(*V,VEC_CLASSID,1);
143 20 : PetscValidLogicalCollectiveInt(*V,nv,2);
144 20 : PetscValidLogicalCollectiveInt(*V,nw,4);
145 5 : if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
146 5 : if (W) {
147 5 : PetscAssertPointer(W,3);
148 5 : PetscValidHeaderSpecific(*W,VEC_CLASSID,3);
149 5 : PetscCheckSameComm(*V,1,*W,3);
150 : }
151 5 : PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
152 5 : PetscFunctionReturn(PETSC_SUCCESS);
153 : }
154 :
155 : /*@
156 : VecCheckOrthonormality - Checks (or prints) the level of (bi-)orthonormality
157 : of a set of vectors.
158 :
159 : Collective
160 :
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
168 :
169 : Output Parameter:
170 : . lev - level of orthogonality (optional)
171 :
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.
175 :
176 : Level: developer
177 :
178 : .seealso: VecCheckOrthogonality()
179 : @*/
180 238 : PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
181 : {
182 238 : PetscFunctionBegin;
183 238 : PetscAssertPointer(V,1);
184 238 : PetscValidHeaderSpecific(*V,VEC_CLASSID,1);
185 952 : PetscValidLogicalCollectiveInt(*V,nv,2);
186 952 : PetscValidLogicalCollectiveInt(*V,nw,4);
187 238 : if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
188 238 : if (W) {
189 0 : PetscAssertPointer(W,3);
190 0 : PetscValidHeaderSpecific(*W,VEC_CLASSID,3);
191 0 : PetscCheckSameComm(*V,1,*W,3);
192 : }
193 238 : PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
194 238 : PetscFunctionReturn(PETSC_SUCCESS);
195 : }
196 :
197 : /*@
198 : VecDuplicateEmpty - Creates a new vector of the same type as an existing vector,
199 : but without internal array.
200 :
201 : Collective
202 :
203 : Input Parameters:
204 : . v - a vector to mimic
205 :
206 : Output Parameter:
207 : . newv - location to put new vector
208 :
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().
212 :
213 : Level: developer
214 :
215 : .seealso: MatCreateVecsEmpty()
216 : @*/
217 10 : PetscErrorCode VecDuplicateEmpty(Vec v,Vec *newv)
218 : {
219 10 : PetscBool standard,cuda,mpi;
220 10 : PetscInt N,nloc,bs;
221 :
222 10 : PetscFunctionBegin;
223 10 : PetscValidHeaderSpecific(v,VEC_CLASSID,1);
224 10 : PetscAssertPointer(newv,2);
225 10 : PetscValidType(v,1);
226 :
227 10 : PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&standard,VECSEQ,VECMPI,""));
228 10 : PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&cuda,VECSEQCUDA,VECMPICUDA,""));
229 10 : if (standard || cuda) {
230 10 : PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,""));
231 10 : PetscCall(VecGetLocalSize(v,&nloc));
232 10 : PetscCall(VecGetSize(v,&N));
233 10 : PetscCall(VecGetBlockSize(v,&bs));
234 10 : 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 10 : if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
241 10 : else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
242 : }
243 0 : } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
244 10 : PetscFunctionReturn(PETSC_SUCCESS);
245 : }
246 :
247 : /*@
248 : VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
249 :
250 : Logically Collective
251 :
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)
257 :
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.
261 :
262 : Level: developer
263 :
264 : .seealso: VecSetRandom()
265 : @*/
266 109 : PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
267 : {
268 109 : const PetscScalar *x,*y;
269 109 : PetscScalar *z;
270 109 : PetscInt n,i;
271 109 : PetscRandom rand=NULL;
272 109 : Vec v1=NULL,v2=NULL;
273 :
274 109 : PetscFunctionBegin;
275 109 : PetscValidHeaderSpecific(v,VEC_CLASSID,1);
276 109 : PetscValidType(v,1);
277 109 : if (rctx) PetscValidHeaderSpecific(rctx,PETSC_RANDOM_CLASSID,2);
278 109 : if (w1) PetscValidHeaderSpecific(w1,VEC_CLASSID,3);
279 109 : if (w2) PetscValidHeaderSpecific(w2,VEC_CLASSID,4);
280 :
281 109 : if (!rctx) {
282 0 : PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
283 0 : PetscCall(PetscRandomSetFromOptions(rand));
284 0 : rctx = rand;
285 : }
286 109 : if (!w1) {
287 0 : PetscCall(VecDuplicate(v,&v1));
288 0 : w1 = v1;
289 : }
290 109 : if (!w2) {
291 0 : PetscCall(VecDuplicate(v,&v2));
292 0 : w2 = v2;
293 : }
294 109 : PetscCheckSameTypeAndComm(v,1,w1,3);
295 109 : PetscCheckSameTypeAndComm(v,1,w2,4);
296 :
297 109 : PetscCall(VecSetRandom(w1,rctx));
298 109 : PetscCall(VecSetRandom(w2,rctx));
299 109 : PetscCall(VecGetLocalSize(v,&n));
300 109 : PetscCall(VecGetArrayWrite(v,&z));
301 109 : PetscCall(VecGetArrayRead(w1,&x));
302 109 : PetscCall(VecGetArrayRead(w2,&y));
303 6513 : for (i=0;i<n;i++) {
304 : #if defined(PETSC_USE_COMPLEX)
305 6404 : 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 109 : PetscCall(VecRestoreArrayWrite(v,&z));
311 109 : PetscCall(VecRestoreArrayRead(w1,&x));
312 109 : PetscCall(VecRestoreArrayRead(w2,&y));
313 :
314 109 : PetscCall(VecDestroy(&v1));
315 109 : PetscCall(VecDestroy(&v2));
316 109 : if (!rctx) PetscCall(PetscRandomDestroy(&rand));
317 109 : PetscFunctionReturn(PETSC_SUCCESS);
318 : }
|
