Actual source code: ex25.c
slepc-3.21.0 2024-03-30
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 generalized symmetric eigenproblem.\n\n"
12: "The problem is similar to ex13.c.\n\n"
13: "The command line options are:\n"
14: " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
15: " -m <m>, where <m> = number of grid subdivisions in y dimension.\n";
17: #include <slepceps.h>
19: int main(int argc,char **argv)
20: {
21: Mat A,B; /* matrices */
22: EPS eps; /* eigenproblem solver context */
23: ST st; /* spectral transformation context */
24: KSP ksp;
25: PC pc;
26: EPSType type;
27: PetscInt N,n=10,m,Istart,Iend,II,nev,i,j,*inertias,ns;
28: PetscReal inta,intb,*shifts;
29: PetscBool flag,show=PETSC_FALSE,terse;
30: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
31: Mat F;
32: #endif
34: PetscFunctionBeginUser;
35: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
37: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
38: PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
39: PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
40: if (!flag) m=n;
41: N = n*m;
42: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on GHEP, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
44: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45: Compute the matrices that define the eigensystem, Ax=kBx
46: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
48: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
49: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
50: PetscCall(MatSetFromOptions(A));
52: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
53: PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
54: PetscCall(MatSetFromOptions(B));
56: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
57: for (II=Istart;II<Iend;II++) {
58: i = II/n; j = II-i*n;
59: if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
60: if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
61: if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
62: if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
63: PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
64: PetscCall(MatSetValue(B,II,II,4.0,INSERT_VALUES));
65: }
67: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
68: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
69: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
72: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73: Create the eigensolver and set various options
74: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76: /*
77: Create eigensolver context
78: */
79: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
81: /*
82: Set operators and set problem type
83: */
84: PetscCall(EPSSetOperators(eps,A,B));
85: PetscCall(EPSSetProblemType(eps,EPS_GHEP));
87: /*
88: Set interval for spectrum slicing
89: */
90: inta = 0.1;
91: intb = 0.2;
92: PetscCall(EPSSetInterval(eps,inta,intb));
93: PetscCall(EPSSetWhichEigenpairs(eps,EPS_ALL));
95: /*
96: Spectrum slicing requires Krylov-Schur
97: */
98: PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
100: /*
101: Set shift-and-invert with Cholesky; select MUMPS if available
102: */
104: PetscCall(EPSGetST(eps,&st));
105: PetscCall(STSetType(st,STSINVERT));
106: PetscCall(EPSKrylovSchurGetKSP(eps,&ksp));
107: PetscCall(KSPSetType(ksp,KSPPREONLY));
108: PetscCall(KSPGetPC(ksp,&pc));
109: PetscCall(PCSetType(pc,PCCHOLESKY));
111: /*
112: Use MUMPS if available.
113: Note that in complex scalars we cannot use MUMPS for spectrum slicing,
114: because MatGetInertia() is not available in that case.
115: */
116: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
117: PetscCall(EPSKrylovSchurSetDetectZeros(eps,PETSC_TRUE)); /* enforce zero detection */
118: PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
119: PetscCall(PCFactorSetUpMatSolverType(pc));
120: /*
121: Set several MUMPS options, the corresponding command-line options are:
122: '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
123: '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
124: '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)
126: Note: depending on the interval, it may be necessary also to increase the workspace:
127: '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
128: */
129: PetscCall(PCFactorGetMatrix(pc,&F));
130: PetscCall(MatMumpsSetIcntl(F,13,1));
131: PetscCall(MatMumpsSetIcntl(F,24,1));
132: PetscCall(MatMumpsSetCntl(F,3,1e-12));
133: #endif
135: /*
136: Set solver parameters at runtime
137: */
138: PetscCall(EPSSetFromOptions(eps));
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Solve the eigensystem
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: PetscCall(EPSSetUp(eps));
144: if (show) {
145: PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
146: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
147: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
148: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
149: PetscCall(PetscFree(shifts));
150: PetscCall(PetscFree(inertias));
151: }
152: PetscCall(EPSSolve(eps));
153: if (show) {
154: PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
155: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
156: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
157: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
158: PetscCall(PetscFree(shifts));
159: PetscCall(PetscFree(inertias));
160: }
162: /*
163: Show eigenvalues in interval and print solution
164: */
165: PetscCall(EPSGetType(eps,&type));
166: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
167: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
168: PetscCall(EPSGetInterval(eps,&inta,&intb));
169: PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
171: /*
172: Show detailed info unless -terse option is given by user
173: */
174: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
175: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
176: else {
177: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
178: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
179: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
180: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
181: }
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Clean up
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186: PetscCall(EPSDestroy(&eps));
187: PetscCall(MatDestroy(&A));
188: PetscCall(MatDestroy(&B));
189: PetscCall(SlepcFinalize());
190: return 0;
191: }
193: /*TEST
195: testset:
196: args: -terse
197: test:
198: requires: !mumps
199: test:
200: requires: mumps !complex
202: TEST*/