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 : This file contains some simple default routines for common operations
12 : */
13 :
14 : #include <slepc/private/epsimpl.h> /*I "slepceps.h" I*/
15 : #include <slepcvec.h>
16 :
17 633 : PetscErrorCode EPSBackTransform_Default(EPS eps)
18 : {
19 633 : PetscFunctionBegin;
20 633 : PetscCall(STBackTransform(eps->st,eps->nconv,eps->eigr,eps->eigi));
21 633 : PetscFunctionReturn(PETSC_SUCCESS);
22 : }
23 :
24 : /*
25 : EPSComputeVectors_Hermitian - Copies the Lanczos vectors as eigenvectors
26 : using purification for generalized eigenproblems.
27 : */
28 285 : PetscErrorCode EPSComputeVectors_Hermitian(EPS eps)
29 : {
30 285 : PetscBool iscayley,indef;
31 285 : Mat B,C;
32 :
33 285 : PetscFunctionBegin;
34 285 : if (eps->purify) {
35 72 : PetscCall(EPS_Purify(eps,eps->nconv));
36 72 : PetscCall(BVNormalize(eps->V,NULL));
37 : } else {
38 : /* In the case of Cayley transform, eigenvectors need to be B-normalized */
39 213 : PetscCall(PetscObjectTypeCompare((PetscObject)eps->st,STCAYLEY,&iscayley));
40 213 : if (iscayley && eps->isgeneralized) {
41 1 : PetscCall(STGetMatrix(eps->st,1,&B));
42 1 : PetscCall(BVGetMatrix(eps->V,&C,&indef));
43 1 : PetscCheck(!indef,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_WRONGSTATE,"The inner product should not be indefinite");
44 1 : PetscCall(BVSetMatrix(eps->V,B,PETSC_FALSE));
45 1 : PetscCall(BVNormalize(eps->V,NULL));
46 1 : PetscCall(BVSetMatrix(eps->V,C,PETSC_FALSE)); /* restore original matrix */
47 : }
48 : }
49 285 : PetscFunctionReturn(PETSC_SUCCESS);
50 : }
51 :
52 : /*
53 : EPSComputeVectors_Indefinite - similar to the Schur version but
54 : for indefinite problems
55 : */
56 16 : PetscErrorCode EPSComputeVectors_Indefinite(EPS eps)
57 : {
58 16 : PetscInt n;
59 16 : Mat X;
60 :
61 16 : PetscFunctionBegin;
62 16 : PetscCall(DSGetDimensions(eps->ds,&n,NULL,NULL,NULL));
63 16 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
64 16 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&X));
65 16 : PetscCall(BVMultInPlace(eps->V,X,0,n));
66 16 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&X));
67 :
68 : /* purification */
69 16 : if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
70 :
71 : /* normalization */
72 16 : PetscCall(BVNormalize(eps->V,eps->eigi));
73 16 : PetscFunctionReturn(PETSC_SUCCESS);
74 : }
75 :
76 : /*
77 : EPSComputeVectors_Twosided - Adjust left eigenvectors in generalized problems: y = B^-* y.
78 : */
79 18 : PetscErrorCode EPSComputeVectors_Twosided(EPS eps)
80 : {
81 18 : PetscInt i;
82 18 : Vec w,y;
83 :
84 18 : PetscFunctionBegin;
85 18 : if (!eps->twosided || !eps->isgeneralized) PetscFunctionReturn(PETSC_SUCCESS);
86 3 : PetscCall(EPSSetWorkVecs(eps,1));
87 3 : w = eps->work[0];
88 10 : for (i=0;i<eps->nconv;i++) {
89 7 : PetscCall(BVCopyVec(eps->W,i,w));
90 7 : PetscCall(BVGetColumn(eps->W,i,&y));
91 7 : PetscCall(STMatSolveHermitianTranspose(eps->st,w,y));
92 7 : PetscCall(BVRestoreColumn(eps->W,i,&y));
93 : }
94 3 : PetscFunctionReturn(PETSC_SUCCESS);
95 : }
96 :
97 : /*
98 : EPSComputeVectors_Schur - Compute eigenvectors from the vectors
99 : provided by the eigensolver. This version is intended for solvers
100 : that provide Schur vectors. Given the partial Schur decomposition
101 : OP*V=V*T, the following steps are performed:
102 : 1) compute eigenvectors of T: T*Z=Z*D
103 : 2) compute eigenvectors of OP: X=V*Z
104 : */
105 302 : PetscErrorCode EPSComputeVectors_Schur(EPS eps)
106 : {
107 302 : PetscInt i;
108 302 : Mat Z;
109 302 : Vec z;
110 :
111 302 : PetscFunctionBegin;
112 302 : if (eps->ishermitian) {
113 49 : if (eps->isgeneralized && !eps->ispositive) PetscCall(EPSComputeVectors_Indefinite(eps));
114 49 : else PetscCall(EPSComputeVectors_Hermitian(eps));
115 49 : PetscFunctionReturn(PETSC_SUCCESS);
116 : }
117 :
118 : /* right eigenvectors */
119 253 : PetscCall(DSVectors(eps->ds,DS_MAT_X,NULL,NULL));
120 :
121 : /* V = V * Z */
122 253 : PetscCall(DSGetMat(eps->ds,DS_MAT_X,&Z));
123 253 : PetscCall(BVMultInPlace(eps->V,Z,0,eps->nconv));
124 253 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_X,&Z));
125 :
126 : /* Purify eigenvectors */
127 253 : if (eps->purify) PetscCall(EPS_Purify(eps,eps->nconv));
128 :
129 : /* Fix eigenvectors if balancing was used */
130 253 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
131 74 : for (i=0;i<eps->nconv;i++) {
132 62 : PetscCall(BVGetColumn(eps->V,i,&z));
133 62 : PetscCall(VecPointwiseDivide(z,z,eps->D));
134 62 : PetscCall(BVRestoreColumn(eps->V,i,&z));
135 : }
136 : }
137 :
138 : /* normalize eigenvectors (when using purification or balancing) */
139 253 : if (eps->purify || (eps->balance!=EPS_BALANCE_NONE && eps->D)) PetscCall(BVNormalize(eps->V,eps->eigi));
140 :
141 : /* left eigenvectors */
142 253 : if (eps->twosided) {
143 13 : PetscCall(DSVectors(eps->ds,DS_MAT_Y,NULL,NULL));
144 : /* W = W * Z */
145 13 : PetscCall(DSGetMat(eps->ds,DS_MAT_Y,&Z));
146 13 : PetscCall(BVMultInPlace(eps->W,Z,0,eps->nconv));
147 13 : PetscCall(DSRestoreMat(eps->ds,DS_MAT_Y,&Z));
148 : /* Fix left eigenvectors if balancing was used */
149 13 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
150 10 : for (i=0;i<eps->nconv;i++) {
151 8 : PetscCall(BVGetColumn(eps->W,i,&z));
152 8 : PetscCall(VecPointwiseMult(z,z,eps->D));
153 8 : PetscCall(BVRestoreColumn(eps->W,i,&z));
154 : }
155 : }
156 13 : PetscCall(EPSComputeVectors_Twosided(eps));
157 : /* normalize */
158 13 : PetscCall(BVNormalize(eps->W,eps->eigi));
159 : #if !defined(PETSC_USE_COMPLEX)
160 60 : for (i=0;i<eps->nconv-1;i++) {
161 47 : if (eps->eigi[i] != 0.0) {
162 2 : if (eps->eigi[i] > 0.0) PetscCall(BVScaleColumn(eps->W,i+1,-1.0));
163 2 : i++;
164 : }
165 : }
166 : #endif
167 : }
168 253 : PetscFunctionReturn(PETSC_SUCCESS);
169 : }
170 :
171 : /*@
172 : EPSSetWorkVecs - Sets a number of work vectors into an EPS object.
173 :
174 : Collective
175 :
176 : Input Parameters:
177 : + eps - eigensolver context
178 : - nw - number of work vectors to allocate
179 :
180 : Developer Notes:
181 : This is SLEPC_EXTERN because it may be required by user plugin EPS
182 : implementations.
183 :
184 : Level: developer
185 :
186 : .seealso: EPSSetUp()
187 : @*/
188 4236 : PetscErrorCode EPSSetWorkVecs(EPS eps,PetscInt nw)
189 : {
190 4236 : Vec t;
191 :
192 4236 : PetscFunctionBegin;
193 4236 : PetscValidHeaderSpecific(eps,EPS_CLASSID,1);
194 12708 : PetscValidLogicalCollectiveInt(eps,nw,2);
195 4236 : PetscCheck(nw>0,PetscObjectComm((PetscObject)eps),PETSC_ERR_ARG_OUTOFRANGE,"nw must be > 0: nw = %" PetscInt_FMT,nw);
196 4236 : if (eps->nwork < nw) {
197 705 : PetscCall(VecDestroyVecs(eps->nwork,&eps->work));
198 705 : eps->nwork = nw;
199 705 : PetscCall(BVGetColumn(eps->V,0,&t));
200 705 : PetscCall(VecDuplicateVecs(t,nw,&eps->work));
201 705 : PetscCall(BVRestoreColumn(eps->V,0,&t));
202 : }
203 4236 : PetscFunctionReturn(PETSC_SUCCESS);
204 : }
205 :
206 : /*
207 : EPSSetWhichEigenpairs_Default - Sets the default value for which,
208 : depending on the ST.
209 : */
210 116 : PetscErrorCode EPSSetWhichEigenpairs_Default(EPS eps)
211 : {
212 116 : PetscBool target;
213 :
214 116 : PetscFunctionBegin;
215 116 : PetscCall(PetscObjectTypeCompareAny((PetscObject)eps->st,&target,STSINVERT,STCAYLEY,""));
216 116 : if (target) eps->which = EPS_TARGET_MAGNITUDE;
217 96 : else eps->which = EPS_LARGEST_MAGNITUDE;
218 116 : PetscFunctionReturn(PETSC_SUCCESS);
219 : }
220 :
221 : /*
222 : EPSConvergedRelative - Checks convergence relative to the eigenvalue.
223 : */
224 47430 : PetscErrorCode EPSConvergedRelative(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
225 : {
226 47430 : PetscReal w;
227 :
228 47430 : PetscFunctionBegin;
229 47430 : w = SlepcAbsEigenvalue(eigr,eigi);
230 47430 : *errest = (w!=0.0)? res/w: PETSC_MAX_REAL;
231 47430 : PetscFunctionReturn(PETSC_SUCCESS);
232 : }
233 :
234 : /*
235 : EPSConvergedAbsolute - Checks convergence absolutely.
236 : */
237 4151 : PetscErrorCode EPSConvergedAbsolute(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
238 : {
239 4151 : PetscFunctionBegin;
240 4151 : *errest = res;
241 4151 : PetscFunctionReturn(PETSC_SUCCESS);
242 : }
243 :
244 : /*
245 : EPSConvergedNorm - Checks convergence relative to the eigenvalue and
246 : the matrix norms.
247 : */
248 5026 : PetscErrorCode EPSConvergedNorm(EPS eps,PetscScalar eigr,PetscScalar eigi,PetscReal res,PetscReal *errest,void *ctx)
249 : {
250 5026 : PetscReal w;
251 :
252 5026 : PetscFunctionBegin;
253 5026 : w = SlepcAbsEigenvalue(eigr,eigi);
254 5026 : *errest = res / (eps->nrma + w*eps->nrmb);
255 5026 : PetscFunctionReturn(PETSC_SUCCESS);
256 : }
257 :
258 : /*@C
259 : EPSStoppingBasic - Default routine to determine whether the outer eigensolver
260 : iteration must be stopped.
261 :
262 : Collective
263 :
264 : Input Parameters:
265 : + eps - eigensolver context obtained from EPSCreate()
266 : . its - current number of iterations
267 : . max_it - maximum number of iterations
268 : . nconv - number of currently converged eigenpairs
269 : . nev - number of requested eigenpairs
270 : - ctx - context (not used here)
271 :
272 : Output Parameter:
273 : . reason - result of the stopping test
274 :
275 : Notes:
276 : A positive value of reason indicates that the iteration has finished successfully
277 : (converged), and a negative value indicates an error condition (diverged). If
278 : the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
279 : (zero).
280 :
281 : EPSStoppingBasic() will stop if all requested eigenvalues are converged, or if
282 : the maximum number of iterations has been reached.
283 :
284 : Use EPSSetStoppingTest() to provide your own test instead of using this one.
285 :
286 : Level: advanced
287 :
288 : .seealso: EPSSetStoppingTest(), EPSStoppingThreshold(), EPSConvergedReason, EPSGetConvergedReason()
289 : @*/
290 46776 : PetscErrorCode EPSStoppingBasic(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
291 : {
292 46776 : PetscFunctionBegin;
293 46776 : *reason = EPS_CONVERGED_ITERATING;
294 46776 : if (nconv >= nev) {
295 736 : PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
296 736 : *reason = EPS_CONVERGED_TOL;
297 46040 : } else if (its >= max_it) {
298 14 : *reason = EPS_DIVERGED_ITS;
299 14 : PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
300 : }
301 46776 : PetscFunctionReturn(PETSC_SUCCESS);
302 : }
303 :
304 : /*@C
305 : EPSStoppingThreshold - Routine to determine whether the outer eigenvalue solver
306 : iteration must be stopped, according to some threshold for the computed values.
307 :
308 : Collective
309 :
310 : Input Parameters:
311 : + eps - eigenvalue solver context obtained from EPSCreate()
312 : . its - current number of iterations
313 : . max_it - maximum number of iterations
314 : . nconv - number of currently converged eigenpairs (ignored here)
315 : . nev - number of requested eigenpairs (ignored here)
316 : - ctx - context containing additional data (EPSStoppingCtx)
317 :
318 : Output Parameter:
319 : . reason - result of the stopping test
320 :
321 : Notes:
322 : A positive value of reason indicates that the iteration has finished successfully
323 : (converged), and a negative value indicates an error condition (diverged). If
324 : the iteration needs to be continued, reason must be set to EPS_CONVERGED_ITERATING
325 : (zero).
326 :
327 : EPSStoppingThreshold() will stop when one of the computed eigenvalues is not
328 : above/below the threshold given at EPSSetThreshold(). If a number of wanted
329 : eigenvalues has been specified via EPSSetDimensions() then it is also taken into
330 : account, and the solver will stop when one of the two conditions (threshold or
331 : number of converged values) is met.
332 :
333 : Use EPSSetStoppingTest() to provide your own test instead of using this one.
334 :
335 : Level: advanced
336 :
337 : .seealso: EPSSetStoppingTest(), EPSStoppingBasic(), EPSSetThreshold(), EPSSetDimensions(), EPSConvergedReason, EPSGetConvergedReason()
338 : @*/
339 813 : PetscErrorCode EPSStoppingThreshold(EPS eps,PetscInt its,PetscInt max_it,PetscInt nconv,PetscInt nev,EPSConvergedReason *reason,void *ctx)
340 : {
341 813 : PetscReal thres,firstev,lastev;
342 813 : PetscBool magnit,rel;
343 813 : EPSWhich which;
344 :
345 813 : PetscFunctionBegin;
346 813 : *reason = EPS_CONVERGED_ITERATING;
347 813 : firstev = ((EPSStoppingCtx)ctx)->firstev;
348 813 : lastev = ((EPSStoppingCtx)ctx)->lastev;
349 813 : thres = ((EPSStoppingCtx)ctx)->thres;
350 813 : rel = ((EPSStoppingCtx)ctx)->threlative;
351 813 : which = ((EPSStoppingCtx)ctx)->which;
352 813 : magnit = (which==EPS_SMALLEST_MAGNITUDE || which==EPS_LARGEST_MAGNITUDE || which==EPS_TARGET_MAGNITUDE)? PETSC_TRUE: PETSC_FALSE;
353 813 : if (nconv && magnit && which==EPS_TARGET_MAGNITUDE && ((rel && ((thres>1.0 && lastev>thres*firstev) || (thres<1.0 && lastev<thres*firstev))) || (!rel && lastev>thres))) {
354 3 : if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
355 2 : else if (thres>1.0) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is above the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
356 1 : else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
357 3 : *reason = EPS_CONVERGED_TOL;
358 810 : } else if (nconv && magnit && ((which==EPS_LARGEST_MAGNITUDE && ((rel && lastev<thres*firstev) || (!rel && lastev<thres))) || (which==EPS_SMALLEST_MAGNITUDE && lastev>thres))) {
359 12 : if (which==EPS_SMALLEST_MAGNITUDE) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is above the threshold %g\n",(double)lastev,(double)thres));
360 5 : else if (!rel) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the eigenvalue magnitude %g is below the threshold %g\n",(double)lastev,(double)thres));
361 5 : else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: the ratio %g/%g is below the threshold %g\n",(double)lastev,(double)firstev,(double)thres));
362 12 : *reason = EPS_CONVERGED_TOL;
363 798 : } else if (nconv && !magnit && ((which==EPS_LARGEST_REAL && lastev<thres) || (which==EPS_SMALLEST_REAL && lastev>thres))) {
364 0 : if (which==EPS_LARGEST_REAL) PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is below the threshold %g\n",(double)lastev,(double)thres));
365 0 : else PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: eigenvalue %g is above the threshold %g\n",(double)lastev,(double)thres));
366 0 : *reason = EPS_CONVERGED_TOL;
367 798 : } else if (nev && nconv >= nev) {
368 0 : PetscCall(PetscInfo(eps,"Linear eigensolver finished successfully: %" PetscInt_FMT " eigenpairs converged at iteration %" PetscInt_FMT "\n",nconv,its));
369 0 : *reason = EPS_CONVERGED_TOL;
370 798 : } else if (its >= max_it) {
371 0 : *reason = EPS_DIVERGED_ITS;
372 0 : PetscCall(PetscInfo(eps,"Linear eigensolver iteration reached maximum number of iterations (%" PetscInt_FMT ")\n",its));
373 : }
374 813 : PetscFunctionReturn(PETSC_SUCCESS);
375 : }
376 :
377 : /*
378 : EPSComputeRitzVector - Computes the current Ritz vector.
379 :
380 : Simple case (complex scalars or real scalars with Zi=NULL):
381 : x = V*Zr (V is a basis of nv vectors, Zr has length nv)
382 :
383 : Split case:
384 : x = V*Zr y = V*Zi (Zr and Zi have length nv)
385 : */
386 385 : PetscErrorCode EPSComputeRitzVector(EPS eps,PetscScalar *Zr,PetscScalar *Zi,BV V,Vec x,Vec y)
387 : {
388 385 : PetscInt l,k;
389 385 : PetscReal norm;
390 : #if !defined(PETSC_USE_COMPLEX)
391 385 : Vec z;
392 : #endif
393 :
394 385 : PetscFunctionBegin;
395 : /* compute eigenvector */
396 385 : PetscCall(BVGetActiveColumns(V,&l,&k));
397 385 : PetscCall(BVSetActiveColumns(V,0,k));
398 385 : PetscCall(BVMultVec(V,1.0,0.0,x,Zr));
399 :
400 : /* purify eigenvector if necessary */
401 385 : if (eps->purify) {
402 30 : PetscCall(STApply(eps->st,x,y));
403 30 : if (eps->ishermitian) PetscCall(BVNormVec(eps->V,y,NORM_2,&norm));
404 18 : else PetscCall(VecNorm(y,NORM_2,&norm));
405 30 : PetscCall(VecScale(y,1.0/norm));
406 30 : PetscCall(VecCopy(y,x));
407 : }
408 : /* fix eigenvector if balancing is used */
409 385 : if (!eps->ishermitian && eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(x,x,eps->D));
410 : #if !defined(PETSC_USE_COMPLEX)
411 : /* compute imaginary part of eigenvector */
412 385 : if (Zi) {
413 25 : PetscCall(BVMultVec(V,1.0,0.0,y,Zi));
414 25 : if (eps->ispositive) {
415 0 : PetscCall(BVCreateVec(V,&z));
416 0 : PetscCall(STApply(eps->st,y,z));
417 0 : PetscCall(VecNorm(z,NORM_2,&norm));
418 0 : PetscCall(VecScale(z,1.0/norm));
419 0 : PetscCall(VecCopy(z,y));
420 0 : PetscCall(VecDestroy(&z));
421 : }
422 25 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) PetscCall(VecPointwiseDivide(y,y,eps->D));
423 : } else
424 : #endif
425 360 : PetscCall(VecSet(y,0.0));
426 :
427 : /* normalize eigenvectors (when using balancing) */
428 385 : if (eps->balance!=EPS_BALANCE_NONE && eps->D) {
429 : #if !defined(PETSC_USE_COMPLEX)
430 12 : if (Zi) PetscCall(VecNormalizeComplex(x,y,PETSC_TRUE,NULL));
431 : else
432 : #endif
433 0 : PetscCall(VecNormalize(x,NULL));
434 : }
435 385 : PetscCall(BVSetActiveColumns(V,l,k));
436 385 : PetscFunctionReturn(PETSC_SUCCESS);
437 : }
438 :
439 : /*
440 : EPSBuildBalance_Krylov - uses a Krylov subspace method to compute the
441 : diagonal matrix to be applied for balancing in non-Hermitian problems.
442 : */
443 13 : PetscErrorCode EPSBuildBalance_Krylov(EPS eps)
444 : {
445 13 : Vec z,p,r;
446 13 : PetscInt i,j;
447 13 : PetscReal norma;
448 13 : PetscScalar *pz,*pD;
449 13 : const PetscScalar *pr,*pp;
450 13 : PetscRandom rand;
451 :
452 13 : PetscFunctionBegin;
453 13 : PetscCall(EPSSetWorkVecs(eps,3));
454 13 : PetscCall(BVGetRandomContext(eps->V,&rand));
455 13 : r = eps->work[0];
456 13 : p = eps->work[1];
457 13 : z = eps->work[2];
458 13 : PetscCall(VecSet(eps->D,1.0));
459 :
460 78 : for (j=0;j<eps->balance_its;j++) {
461 :
462 : /* Build a random vector of +-1's */
463 65 : PetscCall(VecSetRandom(z,rand));
464 65 : PetscCall(VecGetArray(z,&pz));
465 5095 : for (i=0;i<eps->nloc;i++) {
466 5030 : if (PetscRealPart(pz[i])<0.5) pz[i]=-1.0;
467 2394 : else pz[i]=1.0;
468 : }
469 65 : PetscCall(VecRestoreArray(z,&pz));
470 :
471 : /* Compute p=DA(D\z) */
472 65 : PetscCall(VecPointwiseDivide(r,z,eps->D));
473 65 : PetscCall(STApply(eps->st,r,p));
474 65 : PetscCall(VecPointwiseMult(p,p,eps->D));
475 65 : if (eps->balance == EPS_BALANCE_TWOSIDE) {
476 40 : if (j==0) {
477 : /* Estimate the matrix inf-norm */
478 8 : PetscCall(VecAbs(p));
479 8 : PetscCall(VecMax(p,NULL,&norma));
480 : }
481 : /* Compute r=D\(A'Dz) */
482 40 : PetscCall(VecPointwiseMult(z,z,eps->D));
483 40 : PetscCall(STApplyHermitianTranspose(eps->st,z,r));
484 40 : PetscCall(VecPointwiseDivide(r,r,eps->D));
485 : }
486 :
487 : /* Adjust values of D */
488 65 : PetscCall(VecGetArrayRead(r,&pr));
489 65 : PetscCall(VecGetArrayRead(p,&pp));
490 65 : PetscCall(VecGetArray(eps->D,&pD));
491 5095 : for (i=0;i<eps->nloc;i++) {
492 5030 : if (eps->balance == EPS_BALANCE_TWOSIDE) {
493 2900 : if (PetscAbsScalar(pp[i])>eps->balance_cutoff*norma && pr[i]!=0.0)
494 2814 : pD[i] *= PetscSqrtReal(PetscAbsScalar(pr[i]/pp[i]));
495 : } else {
496 2130 : if (pp[i]!=0.0) pD[i] /= PetscAbsScalar(pp[i]);
497 : }
498 : }
499 65 : PetscCall(VecRestoreArrayRead(r,&pr));
500 65 : PetscCall(VecRestoreArrayRead(p,&pp));
501 65 : PetscCall(VecRestoreArray(eps->D,&pD));
502 : }
503 13 : PetscFunctionReturn(PETSC_SUCCESS);
504 : }
|
