Actual source code: ex40.c

slepc-3.20.1 2023-11-27
Report Typos and Errors
  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: 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";

 17: #include <slepcpep.h>

 19: /*
 20:   This example is based on spring.c, for fixed values mu=1,tau=10,kappa=5

 22:   The transformations are based on the method proposed in [Niendorf and Voss, LAA 2010].
 23: */

 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]);

 30: int main(int argc,char **argv)
 31: {
 32:   Mat            M,C,K,*Op,A[3],At[3],B[3]; /* problem matrices */
 33:   PEP            pep;        /* polynomial eigenproblem solver context */
 34:   ST             st;         /* spectral transformation context */
 35:   KSP            ksp;
 36:   PC             pc;
 37:   PEPProblemType type;
 38:   PetscBool      terse,transform=PETSC_FALSE,nohyp=PETSC_FALSE;
 39:   PetscInt       n=100,Istart,Iend,i,def=0,hyp;
 40:   PetscReal      muu=1,tau=10,kappa=5,inta,intb;
 41:   PetscReal      alpha,beta,xi,mu,at[2]={0.0,0.0},c=.857,s;
 42:   PetscScalar    target,targett,ats[2];

 44:   PetscFunctionBeginUser;
 45:   PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));

 47:   PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
 48:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nPEP example that checks definite property, n=%" PetscInt_FMT "\n\n",n));

 50:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 51:      Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
 52:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 54:   /* K is a tridiagonal */
 55:   PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
 56:   PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
 57:   PetscCall(MatSetFromOptions(K));
 58:   PetscCall(MatSetUp(K));

 60:   PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
 61:   for (i=Istart;i<Iend;i++) {
 62:     if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
 63:     PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
 64:     if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
 65:   }

 67:   PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
 68:   PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));

 70:   /* C is a tridiagonal */
 71:   PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
 72:   PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
 73:   PetscCall(MatSetFromOptions(C));
 74:   PetscCall(MatSetUp(C));

 76:   PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
 77:   for (i=Istart;i<Iend;i++) {
 78:     if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
 79:     PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
 80:     if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
 81:   }

 83:   PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
 84:   PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));

 86:   /* M is a diagonal matrix */
 87:   PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
 88:   PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
 89:   PetscCall(MatSetFromOptions(M));
 90:   PetscCall(MatSetUp(M));
 91:   PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
 92:   for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,muu,INSERT_VALUES));
 93:   PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
 94:   PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));

 96:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-nonhyperbolic",&nohyp,NULL));
 97:   A[0] = K; A[1] = C; A[2] = M;
 98:   if (nohyp) {
 99:     s = c*.6;
100:     PetscCall(TransformMatricesMoebius(A,UNKNOWN_NONZERO_PATTERN,c,s,-s,c,At));
101:     for (i=0;i<3;i++) PetscCall(MatDestroy(&A[i]));
102:     Op = At;
103:   } else Op = A;

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106:                 Create the eigensolver and solve the problem
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

109:   /*
110:      Create eigensolver context
111:   */
112:   PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
113:   PetscCall(PEPSetProblemType(pep,PEP_HERMITIAN));
114:   PetscCall(PEPSetType(pep,PEPSTOAR));
115:   /*
116:      Set operators and set problem type
117:   */
118:   PetscCall(PEPSetOperators(pep,3,Op));

120:   /*
121:      Set shift-and-invert with Cholesky; select MUMPS if available
122:   */
123:   PetscCall(PEPGetST(pep,&st));
124:   PetscCall(STGetKSP(st,&ksp));
125:   PetscCall(KSPSetType(ksp,KSPPREONLY));
126:   PetscCall(KSPGetPC(ksp,&pc));
127:   PetscCall(PCSetType(pc,PCCHOLESKY));

129:   /*
130:      Use MUMPS if available.
131:      Note that in complex scalars we cannot use MUMPS for spectrum slicing,
132:      because MatGetInertia() is not available in that case.
133:   */
134: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
135:   PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
136:   /*
137:      Add several MUMPS options (see ex43.c for a better way of setting them in program):
138:      '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
139:   */
140:   PetscCall(PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12"));
141: #endif

143:   /*
144:      Set solver parameters at runtime
145:   */
146:   PetscCall(PEPSetFromOptions(pep));

148:   PetscCall(PetscOptionsGetBool(NULL,NULL,"-transform",&transform,NULL));
149:   if (transform) {
150:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:                     Check if the problem is definite
152:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153:     PetscCall(PEPCheckDefiniteQEP(pep,&xi,&mu,&def,&hyp));
154:     switch (def) {
155:       case 1:
156:         if (hyp==1) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Hyperbolic Problem xi=%g\n",(double)xi));
157:         else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Definite Problem xi=%g mu=%g\n",(double)xi,(double)mu));
158:         break;
159:       case -1:
160:         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Not Definite Problem\n"));
161:         break;
162:       default:
163:         PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Cannot determine definiteness\n"));
164:         break;
165:     }

167:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
168:       Transform the QEP to have a definite inner product in the linearization
169:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
170:     if (def==1) {
171:       PetscCall(QEPDefiniteTransformGetMatrices(pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,B));
172:       PetscCall(PEPSetOperators(pep,3,B));
173:       PetscCall(PEPGetTarget(pep,&target));
174:       targett = target;
175:       PetscCall(QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,1,&targett,PETSC_FALSE));
176:       PetscCall(PEPSetTarget(pep,targett));
177:       PetscCall(PEPGetProblemType(pep,&type));
178:       PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
179:       PetscCall(PEPSTOARGetLinearization(pep,&alpha,&beta));
180:       PetscCall(PEPSTOARSetLinearization(pep,1.0,0.0));
181:       PetscCall(PEPGetInterval(pep,&inta,&intb));
182:       if (inta!=intb) {
183:         ats[0] = inta; ats[1] = intb;
184:         PetscCall(QEPDefiniteTransformMap(hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu,2,ats,PETSC_FALSE));
185:         at[0] = PetscRealPart(ats[0]); at[1] = PetscRealPart(ats[1]);
186:         if (at[0]<at[1]) PetscCall(PEPSetInterval(pep,at[0],at[1]));
187:         else PetscCall(PEPSetInterval(pep,PETSC_MIN_REAL,at[1]));
188:       }
189:     }
190:   }

192:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193:                       Solve the eigensystem
194:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195:   PetscCall(PEPSolve(pep));

197:   /* show detailed info unless -terse option is given by user */
198:   if (def!=1) {
199:     PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
200:     if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
201:     else {
202:       PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
203:       PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
204:       PetscCall(PEPErrorView(pep,PEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
205:       PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
206:     }
207:   } else {
208:     /* Map the solution */
209:     PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
210:     PetscCall(QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu));
211:     for (i=0;i<3;i++) PetscCall(MatDestroy(B+i));
212:   }
213:   if (at[0]>at[1]) {
214:     PetscCall(PEPSetInterval(pep,at[0],PETSC_MAX_REAL));
215:     PetscCall(PEPSolve(pep));
216:     PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
217:     /* Map the solution */
218:     PetscCall(QEPDefiniteCheckError(Op,pep,hyp==1?PETSC_TRUE:PETSC_FALSE,xi,mu));
219:   }
220:   if (def==1) {
221:     PetscCall(PEPSetTarget(pep,target));
222:     PetscCall(PEPSetProblemType(pep,type));
223:     PetscCall(PEPSTOARSetLinearization(pep,alpha,beta));
224:     if (inta!=intb) PetscCall(PEPSetInterval(pep,inta,intb));
225:   }

227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228:                     Clean up
229:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230:   PetscCall(PEPDestroy(&pep));
231:   for (i=0;i<3;i++) PetscCall(MatDestroy(Op+i));
232:   PetscCall(SlepcFinalize());
233:   return 0;
234: }

236: /* ------------------------------------------------------------------- */
237: /*
238:   QEPDefiniteTransformMap_Initial - map a scalar value with a certain Moebius transform

240:                    a theta + b
241:          lambda = --------------
242:                    c theta + d

244:   Input:
245:     xi,mu: real values such that Q(xi)<0 and Q(mu)>0
246:     hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
247:   Input/Output:
248:     val (array of length n)
249:     if backtransform=true returns lambda from theta, else returns theta from lambda
250: */
251: static PetscErrorCode QEPDefiniteTransformMap_Initial(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
252: {
253:   PetscInt  i;
254:   PetscReal a,b,c,d,s;

256:   PetscFunctionBegin;
257:   if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
258:   else { a = mu; b = mu*xi-1; c = 1.0; d = xi+mu; }
259:   if (!backtransform) { s = a; a = -d; d = -s; }
260:   for (i=0;i<n;i++) {
261:     if (PetscRealPart(val[i]) >= PETSC_MAX_REAL || PetscRealPart(val[i]) <= PETSC_MIN_REAL) val[i] = a/c;
262:     else if (val[i] == -d/c) val[i] = PETSC_MAX_REAL;
263:     else val[i] = (a*val[i]+b)/(c*val[i]+d);
264:   }
265:   PetscFunctionReturn(PETSC_SUCCESS);
266: }

268: /* ------------------------------------------------------------------- */
269: /*
270:   QEPDefiniteTransformMap - perform the mapping if the problem is hyperbolic, otherwise
271:   modify the value of xi in advance
272: */
273: PetscErrorCode QEPDefiniteTransformMap(PetscBool hyperbolic,PetscReal xi,PetscReal mu,PetscInt n,PetscScalar *val,PetscBool backtransform)
274: {
275:   PetscReal      xit;
276:   PetscScalar    alpha;

278:   PetscFunctionBegin;
279:   xit = xi;
280:   if (!hyperbolic) {
281:     alpha = xi;
282:     PetscCall(QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&alpha,PETSC_FALSE));
283:     xit = PetscRealPart(alpha);
284:   }
285:   PetscCall(QEPDefiniteTransformMap_Initial(hyperbolic,xit,mu,n,val,backtransform));
286:   PetscFunctionReturn(PETSC_SUCCESS);
287: }

289: /* ------------------------------------------------------------------- */
290: /*
291:   TransformMatricesMoebius - transform the coefficient matrices of a QEP

293:   Input:
294:     A: coefficient matrices of the original QEP
295:     a,b,c,d: parameters of the Moebius transform
296:     str: structure flag for MatAXPY operations
297:   Output:
298:     B: transformed matrices
299: */
300: PetscErrorCode TransformMatricesMoebius(Mat A[3],MatStructure str,PetscReal a,PetscReal b,PetscReal c,PetscReal d,Mat B[3])
301: {
302:   PetscInt       i,k;
303:   PetscReal      cf[9];

305:   PetscFunctionBegin;
306:   for (i=0;i<3;i++) PetscCall(MatDuplicate(A[2],MAT_COPY_VALUES,&B[i]));
307:   /* Ct = b*b*A+b*d*B+d*d*C */
308:   cf[0] = d*d; cf[1] = b*d; cf[2] = b*b;
309:   /* Bt = 2*a*b*A+(b*c+a*d)*B+2*c*d*C*/
310:   cf[3] = 2*c*d; cf[4] = b*c+a*d; cf[5] = 2*a*b;
311:   /* At = a*a*A+a*c*B+c*c*C */
312:   cf[6] = c*c; cf[7] = a*c; cf[8] = a*a;
313:   for (k=0;k<3;k++) {
314:     PetscCall(MatScale(B[k],cf[k*3+2]));
315:     for (i=0;i<2;i++) PetscCall(MatAXPY(B[k],cf[3*k+i],A[i],str));
316:   }
317:   PetscFunctionReturn(PETSC_SUCCESS);
318: }

320: /* ------------------------------------------------------------------- */
321: /*
322:   QEPDefiniteTransformGetMatrices - given a PEP of degree 2, transform the three
323:   matrices with TransformMatricesMoebius

325:   Input:
326:     pep: polynomial eigenproblem to be transformed, with Q(.) being the quadratic polynomial
327:     xi,mu: real values such that Q(xi)<0 and Q(mu)>0
328:     hyperbolic: if true the problem is assumed hyperbolic (mu is not used)
329:   Output:
330:     T: coefficient matrices of the transformed polynomial
331: */
332: PetscErrorCode QEPDefiniteTransformGetMatrices(PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu,Mat T[3])
333: {
334:   MatStructure   str;
335:   ST             st;
336:   PetscInt       i;
337:   PetscReal      a,b,c,d;
338:   PetscScalar    xit;
339:   Mat            A[3];

341:   PetscFunctionBegin;
342:   for (i=2;i>=0;i--) PetscCall(PEPGetOperators(pep,i,&A[i]));
343:   if (hyperbolic) { a = 1.0; b = xi; c =0.0; d = 1.0; }
344:   else {
345:     xit = xi;
346:     PetscCall(QEPDefiniteTransformMap_Initial(PETSC_FALSE,0.0,mu,1,&xit,PETSC_FALSE));
347:     a = mu; b = mu*PetscRealPart(xit)-1.0; c = 1.0; d = PetscRealPart(xit)+mu;
348:   }
349:   PetscCall(PEPGetST(pep,&st));
350:   PetscCall(STGetMatStructure(st,&str));
351:   PetscCall(TransformMatricesMoebius(A,str,a,b,c,d,T));
352:   PetscFunctionReturn(PETSC_SUCCESS);
353: }

355: /* ------------------------------------------------------------------- */
356: /*
357:   Auxiliary function to compute the residual norm of an eigenpair of a QEP defined
358:   by coefficient matrices A
359: */
360: static PetscErrorCode PEPResidualNorm(Mat *A,PetscScalar kr,PetscScalar ki,Vec xr,Vec xi,Vec *z,PetscReal *norm)
361: {
362:   PetscInt       i,nmat=3;
363:   PetscScalar    vals[3];
364:   Vec            u,w;
365: #if !defined(PETSC_USE_COMPLEX)
366:   Vec            ui,wi;
367:   PetscReal      ni;
368:   PetscBool      imag;
369:   PetscScalar    ivals[3];
370: #endif

372:   PetscFunctionBegin;
373:   u = z[0]; w = z[1];
374:   PetscCall(VecSet(u,0.0));
375: #if !defined(PETSC_USE_COMPLEX)
376:   ui = z[2]; wi = z[3];
377: #endif
378:   vals[0] = 1.0;
379:   vals[1] = kr;
380:   vals[2] = kr*kr-ki*ki;
381: #if !defined(PETSC_USE_COMPLEX)
382:   ivals[0] = 0.0;
383:   ivals[1] = ki;
384:   ivals[2] = 2.0*kr*ki;
385:   if (ki == 0 || PetscAbsScalar(ki) < PetscAbsScalar(kr*PETSC_MACHINE_EPSILON))
386:     imag = PETSC_FALSE;
387:   else {
388:     imag = PETSC_TRUE;
389:     PetscCall(VecSet(ui,0.0));
390:   }
391: #endif
392:   for (i=0;i<nmat;i++) {
393:     if (vals[i]!=0.0) {
394:       PetscCall(MatMult(A[i],xr,w));
395:       PetscCall(VecAXPY(u,vals[i],w));
396:     }
397: #if !defined(PETSC_USE_COMPLEX)
398:     if (imag) {
399:       if (ivals[i]!=0 || vals[i]!=0) {
400:         PetscCall(MatMult(A[i],xi,wi));
401:         if (vals[i]==0) PetscCall(MatMult(A[i],xr,w));
402:       }
403:       if (ivals[i]!=0) {
404:         PetscCall(VecAXPY(u,-ivals[i],wi));
405:         PetscCall(VecAXPY(ui,ivals[i],w));
406:       }
407:       if (vals[i]!=0) PetscCall(VecAXPY(ui,vals[i],wi));
408:     }
409: #endif
410:   }
411:   PetscCall(VecNorm(u,NORM_2,norm));
412: #if !defined(PETSC_USE_COMPLEX)
413:   if (imag) {
414:     PetscCall(VecNorm(ui,NORM_2,&ni));
415:     *norm = SlepcAbsEigenvalue(*norm,ni);
416:   }
417: #endif
418:   PetscFunctionReturn(PETSC_SUCCESS);
419: }

421: /* ------------------------------------------------------------------- */
422: /*
423:   QEPDefiniteCheckError - check and print the residual norm of a transformed PEP

425:   Input:
426:     A: coefficient matrices of the original problem
427:     pep: solver containing the computed solution of the transformed problem
428:     xi,mu,hyperbolic: parameters used in transformation
429: */
430: PetscErrorCode QEPDefiniteCheckError(Mat *A,PEP pep,PetscBool hyperbolic,PetscReal xi,PetscReal mu)
431: {
432:   PetscScalar    er,ei;
433:   PetscReal      re,im,error;
434:   Vec            vr,vi,w[4];
435:   PetscInt       i,nconv;
436:   BV             bv;
437:   char           ex[30],sep[]=" ---------------------- --------------------\n";

439:   PetscFunctionBegin;
440:   PetscCall(PetscSNPrintf(ex,sizeof(ex),"||P(k)x||/||kx||"));
441:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s            k             %s\n%s",sep,ex,sep));
442:   PetscCall(PEPGetConverged(pep,&nconv));
443:   PetscCall(PEPGetBV(pep,&bv));
444:   PetscCall(BVCreateVec(bv,w));
445:   PetscCall(VecDuplicate(w[0],&vr));
446:   PetscCall(VecDuplicate(w[0],&vi));
447:   for (i=1;i<4;i++) PetscCall(VecDuplicate(w[0],w+i));
448:   for (i=0;i<nconv;i++) {
449:     PetscCall(PEPGetEigenpair(pep,i,&er,&ei,vr,vi));
450:     PetscCall(QEPDefiniteTransformMap(hyperbolic,xi,mu,1,&er,PETSC_TRUE));
451:     PetscCall(PEPResidualNorm(A,er,0.0,vr,vi,w,&error));
452:     error /= SlepcAbsEigenvalue(er,0.0);
453: #if defined(PETSC_USE_COMPLEX)
454:     re = PetscRealPart(er);
455:     im = PetscImaginaryPart(ei);
456: #else
457:     re = er;
458:     im = ei;
459: #endif
460:     if (im!=0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"  % 9f%+9fi      %12g\n",(double)re,(double)im,(double)error));
461:     else PetscCall(PetscPrintf(PETSC_COMM_WORLD,"    % 12f           %12g\n",(double)re,(double)error));
462:   }
463:   PetscCall(PetscPrintf(PETSC_COMM_WORLD,"%s",sep));
464:   for (i=0;i<4;i++) PetscCall(VecDestroy(w+i));
465:   PetscCall(VecDestroy(&vi));
466:   PetscCall(VecDestroy(&vr));
467:   PetscFunctionReturn(PETSC_SUCCESS);
468: }

470: /*TEST

472:    testset:
473:       requires: !single
474:       args: -pep_nev 3 -nonhyperbolic -pep_target 2
475:       output_file: output/ex40_1.out
476:       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"
477:       test:
478:          suffix: 1
479:          requires: !complex
480:       test:
481:          suffix: 1_complex
482:          requires: complex !mumps
483:       test:
484:          suffix: 1_transform
485:          requires: !complex
486:          args: -transform
487:       test:
488:          suffix: 1_transform_complex
489:          requires: complex !mumps
490:          args: -transform

492: TEST*/