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 quadratic symmetric eigenproblem.\n\n"
12 : "The command line options are:\n"
13 : " -n <n> ... dimension of the matrices.\n\n";
14 :
15 : #include <slepcpep.h>
16 :
17 2 : int main(int argc,char **argv)
18 : {
19 2 : Mat M,C,K,A[3]; /* problem matrices */
20 2 : PEP pep; /* polynomial eigenproblem solver context */
21 2 : ST st; /* spectral transformation context */
22 2 : KSP ksp;
23 2 : PC pc;
24 2 : PEPType type;
25 2 : PetscBool show=PETSC_FALSE,terse;
26 2 : PetscInt n=100,Istart,Iend,nev,i,*inertias,ns;
27 2 : PetscReal mu=1,tau=10,kappa=5,inta,intb,*shifts;
28 :
29 2 : PetscFunctionBeginUser;
30 2 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
31 :
32 2 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
33 2 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-show_inertias",&show,NULL));
34 2 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
35 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\n\n",n));
36 :
37 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
38 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
39 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
40 :
41 : /* K is a tridiagonal */
42 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
43 2 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
44 2 : PetscCall(MatSetFromOptions(K));
45 :
46 2 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
47 202 : for (i=Istart;i<Iend;i++) {
48 200 : if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
49 200 : PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
50 200 : if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
51 : }
52 :
53 2 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
54 2 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
55 :
56 : /* C is a tridiagonal */
57 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
58 2 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
59 2 : PetscCall(MatSetFromOptions(C));
60 :
61 2 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
62 202 : for (i=Istart;i<Iend;i++) {
63 200 : if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
64 200 : PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
65 200 : if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
66 : }
67 :
68 2 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
69 2 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
70 :
71 : /* M is a diagonal matrix */
72 2 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
73 2 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
74 2 : PetscCall(MatSetFromOptions(M));
75 2 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
76 202 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,mu,INSERT_VALUES));
77 2 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
78 2 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
79 :
80 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
81 : Create the eigensolver and solve the problem
82 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
83 :
84 : /*
85 : Create eigensolver context
86 : */
87 2 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
88 :
89 : /*
90 : Set operators and set problem type
91 : */
92 2 : A[0] = K; A[1] = C; A[2] = M;
93 2 : PetscCall(PEPSetOperators(pep,3,A));
94 2 : PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
95 :
96 : /*
97 : Set interval for spectrum slicing
98 : */
99 2 : inta = -11.3;
100 2 : intb = -9.5;
101 2 : PetscCall(PEPSetInterval(pep,inta,intb));
102 2 : PetscCall(PEPSetWhichEigenpairs(pep,PEP_ALL));
103 :
104 : /*
105 : Spectrum slicing requires STOAR
106 : */
107 2 : PetscCall(PEPSetType(pep,PEPSTOAR));
108 :
109 : /*
110 : Set shift-and-invert with Cholesky; select MUMPS if available
111 : */
112 2 : PetscCall(PEPGetST(pep,&st));
113 2 : PetscCall(STSetType(st,STSINVERT));
114 :
115 2 : PetscCall(STGetKSP(st,&ksp));
116 2 : PetscCall(KSPSetType(ksp,KSPPREONLY));
117 2 : PetscCall(KSPGetPC(ksp,&pc));
118 2 : PetscCall(PCSetType(pc,PCCHOLESKY));
119 :
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)
133 :
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
139 :
140 : /*
141 : Set solver parameters at runtime
142 : */
143 2 : PetscCall(PEPSetFromOptions(pep));
144 :
145 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146 : Solve the eigensystem
147 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148 2 : PetscCall(PEPSetUp(pep));
149 2 : if (show) {
150 1 : PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
151 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Subintervals (after setup):\n"));
152 3 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
153 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
154 1 : PetscCall(PetscFree(shifts));
155 1 : PetscCall(PetscFree(inertias));
156 : }
157 2 : PetscCall(PEPSolve(pep));
158 2 : if (show && !terse) {
159 0 : PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
160 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"All shifts (after solve):\n"));
161 0 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Shift %g Inertia %" PetscInt_FMT " \n",(double)shifts[i],inertias[i]));
162 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
163 0 : PetscCall(PetscFree(shifts));
164 0 : PetscCall(PetscFree(inertias));
165 : }
166 :
167 : /*
168 : Show eigenvalues in interval and print solution
169 : */
170 2 : PetscCall(PEPGetType(pep,&type));
171 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
172 2 : PetscCall(PEPGetDimensions(pep,&nev,NULL,NULL));
173 2 : PetscCall(PEPGetInterval(pep,&inta,&intb));
174 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," %" PetscInt_FMT " eigenvalues found in [%g, %g]\n",nev,(double)inta,(double)intb));
175 :
176 : /*
177 : Show detailed info unless -terse option is given by user
178 : */
179 2 : if (terse) PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
180 : else {
181 1 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
182 1 : PetscCall(PEPConvergedReasonView(pep,PETSC_VIEWER_STDOUT_WORLD));
183 1 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,PETSC_VIEWER_STDOUT_WORLD));
184 1 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
185 : }
186 :
187 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188 : Clean up
189 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
190 2 : PetscCall(PEPDestroy(&pep));
191 2 : PetscCall(MatDestroy(&M));
192 2 : PetscCall(MatDestroy(&C));
193 2 : PetscCall(MatDestroy(&K));
194 2 : PetscCall(SlepcFinalize());
195 : return 0;
196 : }
197 :
198 : /*TEST
199 :
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
210 :
211 : test:
212 : suffix: 2
213 : args: -pep_interval 0,2
214 :
215 : TEST*/
|
