Line data Source code
1 : /*
2 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3 : SLEPc - Scalable Library for Eigenvalue Problem Computations
4 : Copyright (c) 2002-, Universitat Politecnica de Valencia, Spain
5 :
6 : This file is part of SLEPc.
7 : SLEPc is distributed under a 2-clause BSD license (see LICENSE).
8 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
9 : */
10 :
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";
16 :
17 : #include <slepceps.h>
18 :
19 2 : int main(int argc,char **argv)
20 : {
21 2 : Mat A,B; /* matrices */
22 2 : EPS eps; /* eigenproblem solver context */
23 2 : ST st; /* spectral transformation context */
24 2 : KSP ksp;
25 2 : PC pc;
26 2 : EPSType type;
27 2 : PetscInt N,n=10,m,Istart,Iend,II,nev,i,j,*inertias,ns;
28 2 : PetscReal inta,intb,*shifts;
29 2 : PetscBool flag,show=PETSC_FALSE,terse;
30 : #if defined(PETSC_HAVE_MUMPS) && !defined(PETSC_USE_COMPLEX)
31 : Mat F;
32 : #endif
33 :
34 2 : PetscFunctionBeginUser;
35 2 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
36 :
37 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
38 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
39 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
40 2 : if (!flag) m=n;
41 2 : N = n*m;
42 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on GHEP, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
43 :
44 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45 : Compute the matrices that define the eigensystem, Ax=kBx
46 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47 :
48 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
49 2 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
50 2 : PetscCall(MatSetFromOptions(A));
51 :
52 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
53 2 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
54 2 : PetscCall(MatSetFromOptions(B));
55 :
56 2 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
57 202 : for (II=Istart;II<Iend;II++) {
58 200 : i = II/n; j = II-i*n;
59 200 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
60 200 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
61 200 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
62 200 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
63 200 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
64 200 : PetscCall(MatSetValue(B,II,II,4.0,INSERT_VALUES));
65 : }
66 :
67 2 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
68 2 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
69 2 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
70 2 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
71 :
72 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
73 : Create the eigensolver and set various options
74 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
75 :
76 : /*
77 : Create eigensolver context
78 : */
79 2 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
80 :
81 : /*
82 : Set operators and set problem type
83 : */
84 2 : PetscCall(EPSSetOperators(eps,A,B));
85 2 : PetscCall(EPSSetProblemType(eps,EPS_GHEP));
86 :
87 : /*
88 : Set interval for spectrum slicing
89 : */
90 2 : inta = 0.1;
91 2 : intb = 0.2;
92 2 : PetscCall(EPSSetInterval(eps,inta,intb));
93 2 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_ALL));
94 :
95 : /*
96 : Spectrum slicing requires Krylov-Schur
97 : */
98 2 : PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
99 :
100 : /*
101 : Set shift-and-invert with Cholesky; select MUMPS if available
102 : */
103 :
104 2 : PetscCall(EPSGetST(eps,&st));
105 2 : PetscCall(STSetType(st,STSINVERT));
106 2 : PetscCall(EPSKrylovSchurGetKSP(eps,&ksp));
107 2 : PetscCall(KSPSetType(ksp,KSPPREONLY));
108 2 : PetscCall(KSPGetPC(ksp,&pc));
109 2 : PetscCall(PCSetType(pc,PCCHOLESKY));
110 :
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)
125 :
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
134 :
135 : /*
136 : Set solver parameters at runtime
137 : */
138 2 : PetscCall(EPSSetFromOptions(eps));
139 :
140 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141 : Solve the eigensystem
142 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143 2 : PetscCall(EPSSetUp(eps));
144 2 : if (show) {
145 0 : PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
146 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
147 0 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
148 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
149 0 : PetscCall(PetscFree(shifts));
150 0 : PetscCall(PetscFree(inertias));
151 : }
152 2 : PetscCall(EPSSolve(eps));
153 2 : if (show) {
154 0 : PetscCall(EPSKrylovSchurGetInertias(eps,&ns,&shifts,&inertias));
155 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
156 0 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
157 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
158 0 : PetscCall(PetscFree(shifts));
159 0 : PetscCall(PetscFree(inertias));
160 : }
161 :
162 : /*
163 : Show eigenvalues in interval and print solution
164 : */
165 2 : PetscCall(EPSGetType(eps,&type));
166 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
167 2 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
168 2 : PetscCall(EPSGetInterval(eps,&inta,&intb));
169 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
170 :
171 : /*
172 : Show detailed info unless -terse option is given by user
173 : */
174 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
175 2 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
176 : else {
177 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
178 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
179 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
180 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
181 : }
182 :
183 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184 : Clean up
185 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
186 2 : PetscCall(EPSDestroy(&eps));
187 2 : PetscCall(MatDestroy(&A));
188 2 : PetscCall(MatDestroy(&B));
189 2 : PetscCall(SlepcFinalize());
190 : return 0;
191 : }
192 :
193 : /*TEST
194 :
195 : testset:
196 : args: -terse
197 : test:
198 : requires: !mumps
199 : test:
200 : requires: mumps !complex
201 :
202 : TEST*/
|
