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 : static char help[] = "Checking the definite property in quadratic symmetric eigenproblem.\n\n"
12 : "The command line options are:\n"
13 : " -n <n> ... dimension of the matrices.\n"
14 : " -transform... whether to transform to a hyperbolic problem or not.\n"
15 : " -nonhyperbolic... to test with a modified (definite) problem that is not hyperbolic.\n\n";
16 :
17 : #include <slepcpep.h>
18 :
19 : /*
20 : This example is based on spring.c, for fixed values mu=1,tau=10,kappa=5
21 :
22 : The transformations are based on the method proposed in [Niendorf and Voss, LAA 2010].
23 : */
24 :
25 : PetscErrorCode QEPDefiniteTransformGetMatrices(PEP,PetscBool,PetscReal,PetscReal,Mat[3]);
26 : PetscErrorCode QEPDefiniteTransformMap(PetscBool,PetscReal,PetscReal,PetscInt,PetscScalar*,PetscBool);
27 : PetscErrorCode QEPDefiniteCheckError(Mat*,PEP,PetscBool,PetscReal,PetscReal);
28 : PetscErrorCode TransformMatricesMoebius(Mat[3],MatStructure,PetscReal,PetscReal,PetscReal,PetscReal,Mat[3]);
29 :
30 2 : int main(int argc,char **argv)
31 : {
32 2 : Mat M,C,K,*Op,A[3],At[3],B[3]; /* problem matrices */
33 2 : PEP pep; /* polynomial eigenproblem solver context */
34 2 : ST st; /* spectral transformation context */
35 2 : KSP ksp;
36 2 : PC pc;
37 2 : PEPProblemType type;
38 2 : PetscBool terse,transform=PETSC_FALSE,nohyp=PETSC_FALSE;
39 2 : PetscInt n=100,Istart,Iend,i,def=0,hyp;
40 2 : PetscReal muu=1,tau=10,kappa=5,inta,intb;
41 2 : PetscReal alpha,beta,xi,mu,at[2]={0.0,0.0},c=.857,s;
42 2 : PetscScalar target,targett,ats[2];
43 :
44 2 : PetscFunctionBeginUser;
45 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
46 :
47 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
48 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPEP example that checks definite property, n=%" PetscInt_FMT "\n\n",n));
49 :
50 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
52 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
53 :
54 : /* K is a tridiagonal */
55 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
56 2 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
57 2 : PetscCall(MatSetFromOptions(K));
58 :
59 2 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
60 202 : for (i=Istart;i<Iend;i++) {
61 200 : if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
62 200 : PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
63 200 : if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
64 : }
65 :
66 2 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
67 2 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
68 :
69 : /* C is a tridiagonal */
70 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
71 2 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
72 2 : PetscCall(MatSetFromOptions(C));
73 :
74 2 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
75 202 : for (i=Istart;i<Iend;i++) {
76 200 : if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
77 200 : PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
78 200 : if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
79 : }
80 :
81 2 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
82 2 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
83 :
84 : /* M is a diagonal matrix */
85 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
86 2 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
87 2 : PetscCall(MatSetFromOptions(M));
88 2 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
89 202 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,muu,INSERT_VALUES));
90 2 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
91 2 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
92 :
93 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-nonhyperbolic",&nohyp,NULL));
94 2 : A[0] = K; A[1] = C; A[2] = M;
95 2 : if (nohyp) {
96 2 : s = c*.6;
97 2 : PetscCall(TransformMatricesMoebius(A,UNKNOWN_NONZERO_PATTERN,c,s,-s,c,At));
98 8 : for (i=0;i<3;i++) PetscCall(MatDestroy(&A[i]));
99 : Op = At;
100 : } else Op = A;
101 :
102 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103 : Create the eigensolver and solve the problem
104 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105 :
106 : /*
107 : Create eigensolver context
108 : */
109 2 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
110 2 : PetscCall(PEPSetProblemType(pep,PEP_HERMITIAN));
111 2 : PetscCall(PEPSetType(pep,PEPSTOAR));
112 : /*
113 : Set operators and set problem type
114 : */
115 2 : PetscCall(PEPSetOperators(pep,3,Op));
116 :
117 : /*
118 : Set shift-and-invert with Cholesky; select MUMPS if available
119 : */
120 2 : PetscCall(PEPGetST(pep,&st));
121 2 : PetscCall(STGetKSP(st,&ksp));
122 2 : PetscCall(KSPSetType(ksp,KSPPREONLY));
123 2 : PetscCall(KSPGetPC(ksp,&pc));
124 2 : PetscCall(PCSetType(pc,PCCHOLESKY));
125 :
126 : /*
127 : Use MUMPS if available.
128 : Note that in complex scalars we cannot use MUMPS for spectrum slicing,
129 : because MatGetInertia() is not available in that case.
130 : */
131 : #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
132 : PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
133 : /*
134 : Add several MUMPS options (see ex43.c for a better way of setting them in program):
135 : '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
136 : */
137 : PetscCall(PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12"));
138 : #endif
139 :
140 : /*
141 : Set solver parameters at runtime
142 : */
143 2 : PetscCall(PEPSetFromOptions(pep));
144 :
145 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-transform",&transform,NULL));
146 2 : if (transform) {
147 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148 : Check if the problem is definite
149 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150 1 : PetscCall(PEPCheckDefiniteQEP(pep,&xi,&mu,&def,&hyp));
151 1 : switch (def) {
152 1 : case 1:
153 1 : if (hyp==1) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Hyperbolic Problem xi=%g\n",(double)xi));
154 1 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Definite Problem xi=%g mu=%g\n",(double)xi,(double)mu));
155 : break;
156 0 : case -1:
157 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Not Definite Problem\n"));
158 : break;
159 0 : default:
160 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cannot determine definiteness\n"));
161 : break;
162 : }
163 :
164 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
165 : Transform the QEP to have a definite inner product in the linearization
166 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167 1 : if (def==1) {
168 1 : PetscCall(QEPDefiniteTransformGetMatrices(pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,B));
169 1 : PetscCall(PEPSetOperators(pep,3,B));
170 1 : PetscCall(PEPGetTarget(pep,&target));
171 1 : targett = target;
172 1 : PetscCall(QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,1,&targett,PETSC_FALSE));
173 1 : PetscCall(PEPSetTarget(pep,targett));
174 1 : PetscCall(PEPGetProblemType(pep,&type));
175 1 : PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
176 1 : PetscCall(PEPSTOARGetLinearization(pep,&alpha,&beta));
177 1 : PetscCall(PEPSTOARSetLinearization(pep,1.0,0.0));
178 1 : PetscCall(PEPGetInterval(pep,&inta,&intb));
179 1 : if (inta!=intb) {
180 0 : ats[0] = inta; ats[1] = intb;
181 0 : PetscCall(QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,2,ats,PETSC_FALSE));
182 0 : at[0] = PetscRealPart(ats[0]); at[1] = PetscRealPart(ats[1]);
183 0 : if (at[0]<at[1]) PetscCall(PEPSetInterval(pep,at[0],at[1]));
184 0 : else PetscCall(PEPSetInterval(pep,PETSC_MIN_REAL,at[1]));
185 : }
186 : }
187 : }
188 :
189 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190 : Solve the eigensystem
191 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
192 2 : PetscCall(PEPSolve(pep));
193 :
194 : /* show detailed info unless -terse option is given by user */
195 2 : if (def!=1) {
196 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
197 1 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
198 : else {
199 1 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
200 1 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
201 1 : PetscCall(PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
202 1 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
203 : }
204 : } else {
205 : /* Map the solution */
206 1 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
207 1 : PetscCall(QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu));
208 4 : for (i=0;i<3;i++) PetscCall(MatDestroy(B+i));
209 : }
210 2 : if (at[0]>at[1]) {
211 0 : PetscCall(PEPSetInterval(pep,at[0],PETSC_MAX_REAL));
212 0 : PetscCall(PEPSolve(pep));
213 0 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
214 : /* Map the solution */
215 0 : PetscCall(QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu));
216 : }
217 2 : if (def==1) {
218 1 : PetscCall(PEPSetTarget(pep,target));
219 1 : PetscCall(PEPSetProblemType(pep,type));
220 1 : PetscCall(PEPSTOARSetLinearization(pep,alpha,beta));
221 1 : if (inta!=intb) PetscCall(PEPSetInterval(pep,inta,intb));
222 : }
223 :
224 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 : Clean up
226 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227 2 : PetscCall(PEPDestroy(&pep));
228 8 : for (i=0;i<3;i++) PetscCall(MatDestroy(Op+i));
229 2 : PetscCall(SlepcFinalize());
230 : return 0;
231 : }
232 :
233 : /* ------------------------------------------------------------------- */
234 : /*
235 : QEPDefiniteTransformMap_Initial - map a scalar value with a certain Moebius transform
236 :
237 : a theta + b
238 : lambda = --------------
239 : c theta + d
240 :
241 : Input:
242 : xi,mu: real values such that Q(xi)<0 and Q(mu)>0
243 : hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
244 : Input/Output:
245 : val (array of length n)
246 : if backtransform=true returns lambda from theta, else returns theta from lambda
247 : */
248 13 : static PetscErrorCode QEPDefiniteTransformMap_Initial(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
249 : {
250 13 : PetscInt i;
251 13 : PetscReal a,b,c,d,s;
252 :
253 13 : PetscFunctionBegin;
254 13 : if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
255 13 : else { a = mu; b = mu*xi-1; c = 1.0; d = xi+mu; }
256 13 : if (!backtransform) { s = a; a = -d; d = -s; }
257 26 : for (i=0;i<n;i++) {
258 13 : if (PetscRealPart(val[i]) >= PETSC_MAX_REAL || PetscRealPart(val[i]) <= PETSC_MIN_REAL) val[i] = a/c;
259 13 : else if (val[i] == -d/c) val[i] = PETSC_MAX_REAL;
260 13 : else val[i] = (a*val[i]+b)/(c*val[i]+d);
261 : }
262 13 : PetscFunctionReturn(PETSC_SUCCESS);
263 : }
264 :
265 : /* ------------------------------------------------------------------- */
266 : /*
267 : QEPDefiniteTransformMap - perform the mapping if the problem is hyperbolic, otherwise
268 : modify the value of xi in advance
269 : */
270 6 : PetscErrorCode QEPDefiniteTransformMap(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
271 : {
272 6 : PetscReal xit;
273 6 : PetscScalar alpha;
274 :
275 6 : PetscFunctionBegin;
276 6 : xit = xi;
277 6 : if (!hyperbolic) {
278 6 : alpha = xi;
279 6 : PetscCall(QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&alpha,PETSC_FALSE));
280 6 : xit = PetscRealPart(alpha);
281 : }
282 6 : PetscCall(QEPDefiniteTransformMap_Initial(hyperbolic,xit,mu,n,val,backtransform));
283 6 : PetscFunctionReturn(PETSC_SUCCESS);
284 : }
285 :
286 : /* ------------------------------------------------------------------- */
287 : /*
288 : TransformMatricesMoebius - transform the coefficient matrices of a QEP
289 :
290 : Input:
291 : A: coefficient matrices of the original QEP
292 : a,b,c,d: parameters of the Moebius transform
293 : str: structure flag for MatAXPY operations
294 : Output:
295 : B: transformed matrices
296 : */
297 3 : PetscErrorCode TransformMatricesMoebius(Mat A[3],MatStructure str,PetscReal a,PetscReal b,PetscReal c,PetscReal d,Mat B[3])
298 : {
299 3 : PetscInt i,k;
300 3 : PetscReal cf[9];
301 :
302 3 : PetscFunctionBegin;
303 12 : for (i=0;i<3;i++) PetscCall(MatDuplicate(A[2],MAT_COPY_VALUES,&B[i]));
304 : /* Ct = b*b*A+b*d*B+d*d*C */
305 3 : cf[0] = d*d; cf[1] = b*d; cf[2] = b*b;
306 : /* Bt = 2*a*b*A+(b*c+a*d)*B+2*c*d*C*/
307 3 : cf[3] = 2*c*d; cf[4] = b*c+a*d; cf[5] = 2*a*b;
308 : /* At = a*a*A+a*c*B+c*c*C */
309 3 : cf[6] = c*c; cf[7] = a*c; cf[8] = a*a;
310 12 : for (k=0;k<3;k++) {
311 9 : PetscCall(MatScale(B[k],cf[k*3+2]));
312 27 : for (i=0;i<2;i++) PetscCall(MatAXPY(B[k],cf[3*k+i],A[i],str));
313 : }
314 3 : PetscFunctionReturn(PETSC_SUCCESS);
315 : }
316 :
317 : /* ------------------------------------------------------------------- */
318 : /*
319 : QEPDefiniteTransformGetMatrices - given a PEP of degree 2, transform the three
320 : matrices with TransformMatricesMoebius
321 :
322 : Input:
323 : pep: polynomial eigenproblem to be transformed, with Q(.) being the quadratic polynomial
324 : xi,mu: real values such that Q(xi)<0 and Q(mu)>0
325 : hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
326 : Output:
327 : T: coefficient matrices of the transformed polynomial
328 : */
329 1 : PetscErrorCode QEPDefiniteTransformGetMatrices(PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu,Mat T[3])
330 : {
331 1 : MatStructure str;
332 1 : ST st;
333 1 : PetscInt i;
334 1 : PetscReal a,b,c,d;
335 1 : PetscScalar xit;
336 1 : Mat A[3];
337 :
338 1 : PetscFunctionBegin;
339 4 : for (i=2;i>=0;i--) PetscCall(PEPGetOperators(pep,i,&A[i]));
340 1 : if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
341 : else {
342 1 : xit = xi;
343 1 : PetscCall(QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&xit,PETSC_FALSE));
344 1 : a = mu; b = mu*PetscRealPart(xit)-1.0; c = 1.0; d = PetscRealPart(xit)+mu;
345 : }
346 1 : PetscCall(PEPGetST(pep,&st));
347 1 : PetscCall(STGetMatStructure(st,&str));
348 1 : PetscCall(TransformMatricesMoebius(A,str,a,b,c,d,T));
349 1 : PetscFunctionReturn(PETSC_SUCCESS);
350 : }
351 :
352 : /* ------------------------------------------------------------------- */
353 : /*
354 : Auxiliary function to compute the residual norm of an eigenpair of a QEP defined
355 : by coefficient matrices A
356 : */
357 5 : static PetscErrorCode PEPResidualNorm(Mat *A,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
358 : {
359 5 : PetscInt i,nmat=3;
360 5 : PetscScalar vals[3];
361 5 : Vec u,w;
362 : #if !defined(PETSC_USE_COMPLEX)
363 : Vec ui,wi;
364 : PetscReal ni;
365 : PetscBool imag;
366 : PetscScalar ivals[3];
367 : #endif
368 :
369 5 : PetscFunctionBegin;
370 5 : u = z[0]; w = z[1];
371 5 : PetscCall(VecSet(u,0.0));
372 : #if !defined(PETSC_USE_COMPLEX)
373 : ui = z[2]; wi = z[3];
374 : #endif
375 5 : vals[0] = 1.0;
376 5 : vals[1] = kr;
377 5 : vals[2] = kr*kr-ki*ki;
378 : #if !defined(PETSC_USE_COMPLEX)
379 : ivals[0] = 0.0;
380 : ivals[1] = ki;
381 : ivals[2] = 2.0*kr*ki;
382 : if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
383 : imag = PETSC_FALSE;
384 : else {
385 : imag = PETSC_TRUE;
386 : PetscCall(VecSet(ui,0.0));
387 : }
388 : #endif
389 20 : for (i=0;i<nmat;i++) {
390 15 : if (vals[i]!=0.0) {
391 15 : PetscCall(MatMult(A[i],xr,w));
392 15 : PetscCall(VecAXPY(u,vals[i],w));
393 : }
394 : #if !defined(PETSC_USE_COMPLEX)
395 : if (imag) {
396 : if (ivals[i]!=0 || vals[i]!=0) {
397 : PetscCall(MatMult(A[i],xi,wi));
398 : if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
399 : }
400 : if (ivals[i]!=0) {
401 : PetscCall(VecAXPY(u,-ivals[i],wi));
402 : PetscCall(VecAXPY(ui,ivals[i],w));
403 : }
404 : if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
405 : }
406 : #endif
407 : }
408 5 : PetscCall(VecNorm(u,NORM_2,norm));
409 : #if !defined(PETSC_USE_COMPLEX)
410 : if (imag) {
411 : PetscCall(VecNorm(ui,NORM_2,&ni));
412 : *norm = SlepcAbsEigenvalue(*norm,ni);
413 : }
414 : #endif
415 5 : PetscFunctionReturn(PETSC_SUCCESS);
416 : }
417 :
418 : /* ------------------------------------------------------------------- */
419 : /*
420 : QEPDefiniteCheckError - check and print the residual norm of a transformed PEP
421 :
422 : Input:
423 : A: coefficient matrices of the original problem
424 : pep: solver containing the computed solution of the transformed problem
425 : xi,mu,hyperbolic: parameters used in transformation
426 : */
427 1 : PetscErrorCode QEPDefiniteCheckError(Mat *A,PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu)
428 : {
429 1 : PetscScalar er,ei;
430 1 : PetscReal re,im,error;
431 1 : Vec vr,vi,w[4];
432 1 : PetscInt i,nconv;
433 1 : BV bv;
434 1 : char ex[30],sep[]=" ---------------------- --------------------\n";
435 :
436 1 : PetscFunctionBegin;
437 1 : PetscCall(PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||"));
438 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s k %s\n%s",sep,ex,sep));
439 1 : PetscCall(PEPGetConverged(pep,&nconv));
440 1 : PetscCall(PEPGetBV(pep,&bv));
441 1 : PetscCall(BVCreateVec(bv,w));
442 1 : PetscCall(VecDuplicate(w[0],&vr));
443 1 : PetscCall(VecDuplicate(w[0],&vi));
444 4 : for (i=1;i<4;i++) PetscCall(VecDuplicate(w[0],w+i));
445 6 : for (i=0;i<nconv;i++) {
446 5 : PetscCall(PEPGetEigenpair(pep,i,&er,&ei,vr,vi));
447 5 : PetscCall(QEPDefiniteTransformMap(hyperbolic,xi,mu,1,&er,PETSC_TRUE));
448 5 : PetscCall(PEPResidualNorm(A,er,0.0,vr,vi,w,&error));
449 5 : error /= SlepcAbsEigenvalue(er,0.0);
450 : #if defined(PETSC_USE_COMPLEX)
451 5 : re = PetscRealPart(er);
452 5 : im = PetscImaginaryPart(ei);
453 : #else
454 : re = er;
455 : im = ei;
456 : #endif
457 5 : if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 9f%+9fi %12g\n",(double)re,(double)im,(double)error));
458 5 : else PetscCall(PetscPrintf(PETSC_COMM_WORLD," % 12f %12g\n",(double)re,(double)error));
459 : }
460 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s",sep));
461 5 : for (i=0;i<4;i++) PetscCall(VecDestroy(w+i));
462 1 : PetscCall(VecDestroy(&vi));
463 1 : PetscCall(VecDestroy(&vr));
464 1 : PetscFunctionReturn(PETSC_SUCCESS);
465 : }
466 :
467 : /*TEST
468 :
469 : testset:
470 : requires: !single
471 : args: -pep_nev 3 -nonhyperbolic -pep_target 2
472 : output_file: output/ex40_1.out
473 : filter: grep -v "Definite" | sed -e "s/iterations [0-9]\([0-9]*\)/iterations xx/g" | sed -e "s/[0-9]\.[0-9]*e[+-]\([0-9]*\)/removed/g"
474 : test:
475 : suffix: 1
476 : requires: !complex
477 : test:
478 : suffix: 1_complex
479 : requires: complex !mumps
480 : test:
481 : suffix: 1_transform
482 : requires: !complex
483 : args: -transform
484 : test:
485 : suffix: 1_transform_complex
486 : requires: complex !mumps
487 : args: -transform
488 :
489 : TEST*/
|
