Actual source code: ex38.c
slepc-main 2024-11-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[] = "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;
29: PetscFunctionBeginUser;
30: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
32: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
33: PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
34: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
35: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\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: PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
43: PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
44: PetscCall(MatSetFromOptions(K));
46: PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
47: for (i=Istart;i<Iend;i++) {
48: if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
49: PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
50: if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
51: }
53: PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
54: PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
56: /* C is a tridiagonal */
57: PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
58: PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
59: PetscCall(MatSetFromOptions(C));
61: PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
62: for (i=Istart;i<Iend;i++) {
63: if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
64: PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
65: if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
66: }
68: PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
69: PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
71: /* M is a diagonal matrix */
72: PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
73: PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
74: PetscCall(MatSetFromOptions(M));
75: PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
76: for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,mu,INSERT_VALUES));
77: PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
78: PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
80: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81: Create the eigensolver and solve the problem
82: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84: /*
85: Create eigensolver context
86: */
87: PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
89: /*
90: Set operators and set problem type
91: */
92: A[0] = K; A[1] = C; A[2] = M;
93: PetscCall(PEPSetOperators(pep,3,A));
94: PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
96: /*
97: Set interval for spectrum slicing
98: */
99: inta = -11.3;
100: intb = -9.5;
101: PetscCall(PEPSetInterval(pep,inta,intb));
102: PetscCall(PEPSetWhichEigenpairs(pep,PEP_ALL));
104: /*
105: Spectrum slicing requires STOAR
106: */
107: PetscCall(PEPSetType(pep,PEPSTOAR));
109: /*
110: Set shift-and-invert with Cholesky; select MUMPS if available
111: */
112: PetscCall(PEPGetST(pep,&st));
113: PetscCall(STSetType(st,STSINVERT));
115: PetscCall(STGetKSP(st,&ksp));
116: PetscCall(KSPSetType(ksp,KSPPREONLY));
117: PetscCall(KSPGetPC(ksp,&pc));
118: PetscCall(PCSetType(pc,PCCHOLESKY));
120: /*
121: Use MUMPS if available.
122: Note that in complex scalars we cannot use MUMPS for spectrum slicing,
123: because MatGetInertia() is not available in that case.
124: */
125: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
126: PetscCall(PEPSTOARSetDetectZeros(pep,PETSC_TRUE)); /* enforce zero detection */
127: PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
128: /*
129: Add several MUMPS options (see ex43.c for a better way of setting them in program):
130: '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
131: '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
132: '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)
134: Note: depending on the interval, it may be necessary also to increase the workspace:
135: '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
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
140: /*
141: Set solver parameters at runtime
142: */
143: PetscCall(PEPSetFromOptions(pep));
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Solve the eigensystem
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: PetscCall(PEPSetUp(pep));
149: if (show) {
150: PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
151: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
152: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
153: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
154: PetscCall(PetscFree(shifts));
155: PetscCall(PetscFree(inertias));
156: }
157: PetscCall(PEPSolve(pep));
158: if (show && !terse) {
159: PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
160: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
161: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
162: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
163: PetscCall(PetscFree(shifts));
164: PetscCall(PetscFree(inertias));
165: }
167: /*
168: Show eigenvalues in interval and print solution
169: */
170: PetscCall(PEPGetType(pep,&type));
171: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
172: PetscCall(PEPGetDimensions(pep,&nev,NULL,NULL));
173: PetscCall(PEPGetInterval(pep,&inta,&intb));
174: PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
176: /*
177: Show detailed info unless -terse option is given by user
178: */
179: if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
180: else {
181: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
182: PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
183: PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
184: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
185: }
187: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: Clean up
189: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190: PetscCall(PEPDestroy(&pep));
191: PetscCall(MatDestroy(&M));
192: PetscCall(MatDestroy(&C));
193: PetscCall(MatDestroy(&K));
194: PetscCall(SlepcFinalize());
195: return 0;
196: }
198: /*TEST
200: testset:
201: requires: !single
202: args: -show_inertias -terse
203: output_file: output/ex38_1.out
204: test:
205: suffix: 1
206: requires: !complex
207: test:
208: suffix: 1_complex
209: requires: complex !mumps
211: test:
212: suffix: 2
213: args: -pep_interval 0,2
215: TEST*/