Actual source code: ex38.c
slepc-3.16.1 2021-11-17
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2021, 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[] = "Spectrum slicing on quadratic symmetric eigenproblem.\n\n"
12: "The command line options are:\n"
13: " -n <n> ... dimension of the matrices.\n\n";
15: #include <slepcpep.h>
17: int main(int argc,char **argv)
18: {
19: Mat M,C,K,A[3]; /* problem matrices */
20: PEP pep; /* polynomial eigenproblem solver context */
21: ST st; /* spectral transformation context */
22: KSP ksp;
23: PC pc;
24: PEPType type;
25: PetscBool show=PETSC_FALSE,terse;
26: PetscInt n=100,Istart,Iend,nev,i,*inertias,ns;
27: PetscReal mu=1,tau=10,kappa=5,inta,intb,*shifts;
30: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
32: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
33: PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL);
34: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
35: PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%D\n\n",n);
37: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38: Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
39: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
41: /* K is a tridiagonal */
42: MatCreate(PETSC_COMM_WORLD,&K);
43: MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n);
44: MatSetFromOptions(K);
45: MatSetUp(K);
47: MatGetOwnershipRange(K,&Istart,&Iend);
48: for (i=Istart;i<Iend;i++) {
49: if (i>0) {
50: MatSetValue(K,i,i-1,-kappa,INSERT_VALUES);
51: }
52: MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES);
53: if (i<n-1) {
54: MatSetValue(K,i,i+1,-kappa,INSERT_VALUES);
55: }
56: }
58: MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY);
59: MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY);
61: /* C is a tridiagonal */
62: MatCreate(PETSC_COMM_WORLD,&C);
63: MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n);
64: MatSetFromOptions(C);
65: MatSetUp(C);
67: MatGetOwnershipRange(C,&Istart,&Iend);
68: for (i=Istart;i<Iend;i++) {
69: if (i>0) {
70: MatSetValue(C,i,i-1,-tau,INSERT_VALUES);
71: }
72: MatSetValue(C,i,i,tau*3.0,INSERT_VALUES);
73: if (i<n-1) {
74: MatSetValue(C,i,i+1,-tau,INSERT_VALUES);
75: }
76: }
78: MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);
79: MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);
81: /* M is a diagonal matrix */
82: MatCreate(PETSC_COMM_WORLD,&M);
83: MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n);
84: MatSetFromOptions(M);
85: MatSetUp(M);
86: MatGetOwnershipRange(M,&Istart,&Iend);
87: for (i=Istart;i<Iend;i++) {
88: MatSetValue(M,i,i,mu,INSERT_VALUES);
89: }
90: MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY);
91: MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY);
93: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
94: Create the eigensolver and solve the problem
95: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
97: /*
98: Create eigensolver context
99: */
100: PEPCreate(PETSC_COMM_WORLD,&pep);
102: /*
103: Set operators and set problem type
104: */
105: A[0] = K; A[1] = C; A[2] = M;
106: PEPSetOperators(pep,3,A);
107: PEPSetProblemType(pep,PEP_HYPERBOLIC);
109: /*
110: Set interval for spectrum slicing
111: */
112: inta = -11.3;
113: intb = -9.5;
114: PEPSetInterval(pep,inta,intb);
115: PEPSetWhichEigenpairs(pep,PEP_ALL);
117: /*
118: Spectrum slicing requires STOAR
119: */
120: PEPSetType(pep,PEPSTOAR);
122: /*
123: Set shift-and-invert with Cholesky; select MUMPS if available
124: */
125: PEPGetST(pep,&st);
126: STSetType(st,STSINVERT);
128: STGetKSP(st,&ksp);
129: KSPSetType(ksp,KSPPREONLY);
130: KSPGetPC(ksp,&pc);
131: PCSetType(pc,PCCHOLESKY);
133: /*
134: Use MUMPS if available.
135: Note that in complex scalars we cannot use MUMPS for spectrum slicing,
136: because MatGetInertia() is not available in that case.
137: */
138: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
139: PEPSTOARSetDetectZeros(pep,PETSC_TRUE); /* enforce zero detection */
140: PCFactorSetMatSolverType(pc,MATSOLVERMUMPS);
141: /*
142: Add several MUMPS options (see ex43.c for a better way of setting them in program):
143: '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
144: '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
145: '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)
147: Note: depending on the interval, it may be necessary also to increase the workspace:
148: '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
149: */
150: PetscOptionsInsertString(NULL,"-st_mat_mumps_icntl_13 1 -st_mat_mumps_icntl_24 1 -st_mat_mumps_cntl_3 1e-12");
151: #endif
153: /*
154: Set solver parameters at runtime
155: */
156: PEPSetFromOptions(pep);
158: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159: Solve the eigensystem
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161: PEPSetUp(pep);
162: if (show) {
163: PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
164: PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n");
165: for (i=0;i<ns;i++) { PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %D \n",(double)shifts[i],inertias[i]); }
166: PetscPrintf(PETSC_COMM_WORLD,"\n");
167: PetscFree(shifts);
168: PetscFree(inertias);
169: }
170: PEPSolve(pep);
171: if (show && !terse) {
172: PEPSTOARGetInertias(pep,&ns,&shifts,&inertias);
173: PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n");
174: for (i=0;i<ns;i++) { PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %D \n",(double)shifts[i],inertias[i]); }
175: PetscPrintf(PETSC_COMM_WORLD,"\n");
176: PetscFree(shifts);
177: PetscFree(inertias);
178: }
180: /*
181: Show eigenvalues in interval and print solution
182: */
183: PEPGetType(pep,&type);
184: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
185: PEPGetDimensions(pep,&nev,NULL,NULL);
186: PEPGetInterval(pep,&inta,&intb);
187: PetscPrintf(PETSC_COMM_WORLD," %D eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb);
189: /*
190: Show detailed info unless -terse option is given by user
191: */
192: if (terse) {
193: PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL);
194: } else {
195: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
196: PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD);
197: PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD);
198: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
199: }
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Clean up
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: PEPDestroy(&pep);
205: MatDestroy(&M);
206: MatDestroy(&C);
207: MatDestroy(&K);
208: SlepcFinalize();
209: return ierr;
210: }
212: /*TEST
214: testset:
215: requires: !single
216: args: -show_inertias -terse
217: output_file: output/ex38_1.out
218: test:
219: suffix: 1
220: requires: !complex
221: test:
222: suffix: 1_complex
223: requires: complex !mumps
225: test:
226: suffix: 2
227: args: -pep_interval 0,2
229: TEST*/