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