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[] = "Compute all eigenvalues in an interval of a symmetric-definite problem.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
14 : " -m <m>, where <m> = number of grid subdivisions in y dimension.\n\n";
15 :
16 : #include <slepceps.h>
17 :
18 1 : int main(int argc,char **argv)
19 : {
20 1 : Mat A,B; /* matrices */
21 1 : EPS eps; /* eigenproblem solver context */
22 1 : ST st; /* spectral transformation context */
23 1 : KSP ksp;
24 1 : PC pc;
25 1 : PetscInt N,n=35,m,Istart,Iend,II,nev,i,j,k,*inertias;
26 1 : PetscBool flag,showinertia=PETSC_TRUE;
27 1 : PetscReal int0,int1,*shifts;
28 :
29 1 : PetscFunctionBeginUser;
30 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
31 :
32 1 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-showinertia",&showinertia,NULL));
33 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
34 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
35 1 : if (!flag) m=n;
36 1 : N = n*m;
37 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSymmetric-definite problem with two intervals, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
38 :
39 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
40 : Compute the matrices that define the eigensystem, Ax=kBx
41 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
42 :
43 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
44 1 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
45 1 : PetscCall(MatSetFromOptions(A));
46 :
47 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
48 1 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
49 1 : PetscCall(MatSetFromOptions(B));
50 :
51 1 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
52 1226 : for (II=Istart;II<Iend;II++) {
53 1225 : i = II/n; j = II-i*n;
54 1225 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
55 1225 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
56 1225 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
57 1225 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
58 1225 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
59 1225 : PetscCall(MatSetValue(B,II,II,2.0,INSERT_VALUES));
60 : }
61 1 : if (Istart==0) {
62 1 : PetscCall(MatSetValue(B,0,0,6.0,INSERT_VALUES));
63 1 : PetscCall(MatSetValue(B,0,1,-1.0,INSERT_VALUES));
64 1 : PetscCall(MatSetValue(B,1,0,-1.0,INSERT_VALUES));
65 1 : PetscCall(MatSetValue(B,1,1,1.0,INSERT_VALUES));
66 : }
67 :
68 1 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
69 1 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
70 1 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
71 1 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
72 :
73 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74 : Create the eigensolver and set various options
75 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
76 :
77 1 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
78 1 : PetscCall(EPSSetOperators(eps,A,B));
79 1 : PetscCall(EPSSetProblemType(eps,EPS_GHEP));
80 :
81 : /*
82 : Set first interval and other settings for spectrum slicing
83 : */
84 1 : PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
85 1 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_ALL));
86 1 : PetscCall(EPSSetInterval(eps,1.1,1.3));
87 1 : PetscCall(EPSGetST(eps,&st));
88 1 : PetscCall(STSetType(st,STSINVERT));
89 1 : PetscCall(EPSKrylovSchurGetKSP(eps,&ksp));
90 1 : PetscCall(KSPGetPC(ksp,&pc));
91 1 : PetscCall(KSPSetType(ksp,KSPPREONLY));
92 1 : PetscCall(PCSetType(pc,PCCHOLESKY));
93 1 : PetscCall(EPSSetFromOptions(eps));
94 :
95 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 : Solve for first interval and display info
97 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 :
99 1 : PetscCall(EPSSolve(eps));
100 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
101 1 : PetscCall(EPSGetInterval(eps,&int0,&int1));
102 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1));
103 1 : if (showinertia) {
104 0 : PetscCall(EPSKrylovSchurGetInertias(eps,&k,&shifts,&inertias));
105 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Used %" PetscInt_FMT " shifts (inertia):\n",k));
106 0 : for (i=0;i<k;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
107 0 : PetscCall(PetscFree(shifts));
108 0 : PetscCall(PetscFree(inertias));
109 : }
110 :
111 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
112 : Solve for second interval and display info
113 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
114 1 : PetscCall(EPSSetInterval(eps,1.499,1.6));
115 1 : PetscCall(EPSSolve(eps));
116 1 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
117 1 : PetscCall(EPSGetInterval(eps,&int0,&int1));
118 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1));
119 1 : if (showinertia) {
120 0 : PetscCall(EPSKrylovSchurGetInertias(eps,&k,&shifts,&inertias));
121 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Used %" PetscInt_FMT " shifts (inertia):\n",k));
122 0 : for (i=0;i<k;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
123 0 : PetscCall(PetscFree(shifts));
124 0 : PetscCall(PetscFree(inertias));
125 : }
126 :
127 1 : PetscCall(EPSDestroy(&eps));
128 1 : PetscCall(MatDestroy(&A));
129 1 : PetscCall(MatDestroy(&B));
130 1 : PetscCall(SlepcFinalize());
131 : return 0;
132 : }
133 :
134 : /*TEST
135 :
136 : test:
137 : suffix: 1
138 : args: -showinertia 0 -eps_error_relative
139 : requires: !single
140 :
141 : TEST*/
|
