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 262 : static PetscErrorCode VecCheckOrthogonality_Private(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev,PetscBool norm)
57 : {
58 262 : PetscInt i,j;
59 262 : PetscScalar *vals;
60 262 : PetscBool isascii;
61 262 : Vec w;
62 :
63 262 : PetscFunctionBegin;
64 262 : 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 262 : PetscCall(PetscMalloc1(nv,&vals));
73 262 : if (B) PetscCall(VecDuplicate(V[0],&w));
74 262 : if (lev) *lev = 0.0;
75 3238 : for (i=0;i<nw;i++) {
76 2976 : if (B) {
77 1210 : if (W) PetscCall(MatMultTranspose(B,W[i],w));
78 1204 : else PetscCall(MatMultTranspose(B,V[i],w));
79 : } else {
80 1766 : if (W) w = W[i];
81 1481 : else w = V[i];
82 : }
83 2976 : PetscCall(VecMDot(w,nv,V,vals));
84 187692 : for (j=0;j<nv;j++) {
85 184716 : if (lev) {
86 365971 : if (i!=j) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]));
87 2976 : else if (norm) {
88 2685 : if (PetscRealPart(vals[j])<0.0) *lev = PetscMax(*lev,PetscAbsScalar(vals[j]+PetscRealConstant(1.0))); /* indefinite case */
89 4934 : 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 184716 : PetscCall(PetscViewerASCIIPrintf(viewer," %12g%+12gi ",(double)PetscRealPart(vals[j]),(double)PetscImaginaryPart(vals[j])));
96 : #endif
97 : }
98 : }
99 2976 : if (!lev) PetscCall(PetscViewerASCIIPrintf(viewer,"\n"));
100 : }
101 262 : PetscCall(PetscFree(vals));
102 262 : if (B) PetscCall(VecDestroy(&w));
103 262 : 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 22 : PetscErrorCode VecCheckOrthogonality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
139 : {
140 22 : PetscFunctionBegin;
141 22 : PetscAssertPointer(V,1);
142 22 : PetscValidHeaderSpecific(*V,VEC_CLASSID,1);
143 66 : PetscValidLogicalCollectiveInt(*V,nv,2);
144 66 : PetscValidLogicalCollectiveInt(*V,nw,4);
145 22 : if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
146 22 : if (W) {
147 22 : PetscAssertPointer(W,3);
148 22 : PetscValidHeaderSpecific(*W,VEC_CLASSID,3);
149 22 : PetscCheckSameComm(*V,1,*W,3);
150 : }
151 22 : PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_FALSE));
152 22 : 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 240 : PetscErrorCode VecCheckOrthonormality(Vec V[],PetscInt nv,Vec W[],PetscInt nw,Mat B,PetscViewer viewer,PetscReal *lev)
181 : {
182 240 : PetscFunctionBegin;
183 240 : PetscAssertPointer(V,1);
184 240 : PetscValidHeaderSpecific(*V,VEC_CLASSID,1);
185 720 : PetscValidLogicalCollectiveInt(*V,nv,2);
186 720 : PetscValidLogicalCollectiveInt(*V,nw,4);
187 240 : if (nv<=0 || nw<=0) PetscFunctionReturn(PETSC_SUCCESS);
188 240 : if (W) {
189 0 : PetscAssertPointer(W,3);
190 0 : PetscValidHeaderSpecific(*W,VEC_CLASSID,3);
191 0 : PetscCheckSameComm(*V,1,*W,3);
192 : }
193 240 : PetscCall(VecCheckOrthogonality_Private(V,nv,W,nw,B,viewer,lev,PETSC_TRUE));
194 240 : 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,hip,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 : PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&hip,VECSEQHIP,VECMPIHIP,""));
230 10 : if (standard || cuda || hip) {
231 10 : PetscCall(PetscObjectTypeCompareAny((PetscObject)v,&mpi,VECMPI,VECMPICUDA,VECMPIHIP,""));
232 10 : PetscCall(VecGetLocalSize(v,&nloc));
233 10 : PetscCall(VecGetSize(v,&N));
234 10 : PetscCall(VecGetBlockSize(v,&bs));
235 10 : if (cuda) {
236 : #if defined(PETSC_HAVE_CUDA)
237 : if (mpi) PetscCall(VecCreateMPICUDAWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
238 : else PetscCall(VecCreateSeqCUDAWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
239 : #endif
240 10 : } else if (hip) {
241 : #if defined(PETSC_HAVE_HIP)
242 : if (mpi) PetscCall(VecCreateMPIHIPWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
243 : else PetscCall(VecCreateSeqHIPWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
244 : #endif
245 : } else {
246 10 : if (mpi) PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)v),bs,nloc,N,NULL,newv));
247 10 : else PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)v),bs,N,NULL,newv));
248 : }
249 0 : } else PetscCall(VecDuplicate(v,newv)); /* standard duplicate, with internal array */
250 10 : PetscFunctionReturn(PETSC_SUCCESS);
251 : }
252 :
253 : /*@
254 : VecSetRandomNormal - Sets all components of a vector to normally distributed random values.
255 :
256 : Logically Collective
257 :
258 : Input Parameters:
259 : + v - the vector to be filled with random values
260 : . rctx - the random number context (can be NULL)
261 : . w1 - first work vector (can be NULL)
262 : - w2 - second work vector (can be NULL)
263 :
264 : Notes:
265 : Fills the two work vectors with uniformly distributed random values (VecSetRandom)
266 : and then applies the Box-Muller transform to get normally distributed values on v.
267 :
268 : Level: developer
269 :
270 : .seealso: VecSetRandom()
271 : @*/
272 109 : PetscErrorCode VecSetRandomNormal(Vec v,PetscRandom rctx,Vec w1,Vec w2)
273 : {
274 109 : const PetscScalar *x,*y;
275 109 : PetscScalar *z;
276 109 : PetscInt n,i;
277 109 : PetscRandom rand=NULL;
278 109 : Vec v1=NULL,v2=NULL;
279 :
280 109 : PetscFunctionBegin;
281 109 : PetscValidHeaderSpecific(v,VEC_CLASSID,1);
282 109 : PetscValidType(v,1);
283 109 : if (rctx) PetscValidHeaderSpecific(rctx,PETSC_RANDOM_CLASSID,2);
284 109 : if (w1) PetscValidHeaderSpecific(w1,VEC_CLASSID,3);
285 109 : if (w2) PetscValidHeaderSpecific(w2,VEC_CLASSID,4);
286 :
287 109 : if (!rctx) {
288 0 : PetscCall(PetscRandomCreate(PetscObjectComm((PetscObject)v),&rand));
289 0 : PetscCall(PetscRandomSetFromOptions(rand));
290 0 : rctx = rand;
291 : }
292 109 : if (!w1) {
293 0 : PetscCall(VecDuplicate(v,&v1));
294 0 : w1 = v1;
295 : }
296 109 : if (!w2) {
297 0 : PetscCall(VecDuplicate(v,&v2));
298 0 : w2 = v2;
299 : }
300 109 : PetscCheckSameTypeAndComm(v,1,w1,3);
301 109 : PetscCheckSameTypeAndComm(v,1,w2,4);
302 :
303 109 : PetscCall(VecSetRandom(w1,rctx));
304 109 : PetscCall(VecSetRandom(w2,rctx));
305 109 : PetscCall(VecGetLocalSize(v,&n));
306 109 : PetscCall(VecGetArrayWrite(v,&z));
307 109 : PetscCall(VecGetArrayRead(w1,&x));
308 109 : PetscCall(VecGetArrayRead(w2,&y));
309 6513 : for (i=0;i<n;i++) {
310 : #if defined(PETSC_USE_COMPLEX)
311 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])));
312 : #else
313 : z[i] = PetscSqrtReal(-2.0*PetscLogReal(x[i]))*PetscCosReal(2.0*PETSC_PI*y[i]);
314 : #endif
315 : }
316 109 : PetscCall(VecRestoreArrayWrite(v,&z));
317 109 : PetscCall(VecRestoreArrayRead(w1,&x));
318 109 : PetscCall(VecRestoreArrayRead(w2,&y));
319 :
320 109 : PetscCall(VecDestroy(&v1));
321 109 : PetscCall(VecDestroy(&v2));
322 109 : if (!rctx) PetscCall(PetscRandomDestroy(&rand));
323 109 : PetscFunctionReturn(PETSC_SUCCESS);
324 : }
|
