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