Actual source code: ex25.c
slepc-3.20.2 2024-03-15
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));
51: PetscCall(MatSetUp(A));
53: PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
54: PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
55: PetscCall(MatSetFromOptions(B));
56: PetscCall(MatSetUp(B));
58: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
59: for (II=Istart;II<Iend;II++) {
60: i = II/n; j = II-i*n;
61: if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
62: if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
63: if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
64: if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
65: PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
66: PetscCall(MatSetValue(B,II,II,4.0,INSERT_VALUES));
67: }
69: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
70: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
71: PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
72: PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create the eigensolver and set various options
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: /*
79: Create eigensolver context
80: */
81: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
83: /*
84: Set operators and set problem type
85: */
86: PetscCall(EPSSetOperators(eps,A,B));
87: PetscCall(EPSSetProblemType(eps,EPS_GHEP));
89: /*
90: Set interval for spectrum slicing
91: */
92: inta = 0.1;
93: intb = 0.2;
94: PetscCall(EPSSetInterval(eps,inta,intb));
95: PetscCall(EPSSetWhichEigenpairs(eps,EPS_ALL));
97: /*
98: Spectrum slicing requires Krylov-Schur
99: */
100: PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
102: /*
103: Set shift-and-invert with Cholesky; select MUMPS if available
104: */
106: PetscCall(EPSGetST(eps,&st));
107: PetscCall(STSetType(st,STSINVERT));
108: PetscCall(EPSKrylovSchurGetKSP(eps,&ksp));
109: PetscCall(KSPSetType(ksp,KSPPREONLY));
110: PetscCall(KSPGetPC(ksp,&pc));
111: PetscCall(PCSetType(pc,PCCHOLESKY));
113: /*
114: Use MUMPS if available.
115: Note that in complex scalars we cannot use MUMPS for spectrum slicing,
116: because MatGetInertia() is not available in that case.
117: */
118: #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
119: PetscCall(EPSKrylovSchurSetDetectZeros(eps,PETSC_TRUE)); /* enforce zero detection */
120: PetscCall(PCFactorSetMatSolverType(pc,MATSOLVERMUMPS));
121: PetscCall(PCFactorSetUpMatSolverType(pc));
122: /*
123: Set several MUMPS options, the corresponding command-line options are:
124: '-st_mat_mumps_icntl_13 1': turn off ScaLAPACK for matrix inertia
125: '-st_mat_mumps_icntl_24 1': detect null pivots in factorization (for the case that a shift is equal to an eigenvalue)
126: '-st_mat_mumps_cntl_3 <tol>': a tolerance used for null pivot detection (must be larger than machine epsilon)
128: Note: depending on the interval, it may be necessary also to increase the workspace:
129: '-st_mat_mumps_icntl_14 <percentage>': increase workspace with a percentage (50, 100 or more)
130: */
131: PetscCall(PCFactorGetMatrix(pc,&F));
132: PetscCall(MatMumpsSetIcntl(F,13,1));
133: PetscCall(MatMumpsSetIcntl(F,24,1));
134: PetscCall(MatMumpsSetCntl(F,3,1e-12));
135: #endif
137: /*
138: Set solver parameters at runtime
139: */
140: PetscCall(EPSSetFromOptions(eps));
142: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143: Solve the eigensystem
144: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145: PetscCall(EPSSetUp(eps));
146: if (show) {
147: PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
148: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
149: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
150: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
151: PetscCall(PetscFree(shifts));
152: PetscCall(PetscFree(inertias));
153: }
154: PetscCall(EPSSolve(eps));
155: if (show) {
156: PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
157: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
158: for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
159: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
160: PetscCall(PetscFree(shifts));
161: PetscCall(PetscFree(inertias));
162: }
164: /*
165: Show eigenvalues in interval and print solution
166: */
167: PetscCall(EPSGetType(eps,&type));
168: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
169: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
170: PetscCall(EPSGetInterval(eps,&inta,&intb));
171: PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
173: /*
174: Show detailed info unless -terse option is given by user
175: */
176: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
177: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
178: else {
179: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
180: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
181: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
182: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
183: }
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Clean up
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
188: PetscCall(EPSDestroy(&eps));
189: PetscCall(MatDestroy(&A));
190: PetscCall(MatDestroy(&B));
191: PetscCall(SlepcFinalize());
192: return 0;
193: }
195: /*TEST
197: testset:
198: args: -terse
199: test:
200: requires: !mumps
201: test:
202: requires: mumps !complex
204: TEST*/