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[] = "Test interface functions of spectrum-slicing STOAR.\n\n"
12 : "This is based on ex38.c. The command line options are:\n"
13 : " -n <n> ... dimension of the matrices.\n\n";
14 :
15 : #include <slepcpep.h>
16 :
17 1 : int main(int argc,char **argv)
18 : {
19 1 : Mat M,C,K,A[3]; /* problem matrices */
20 1 : PEP pep; /* polynomial eigenproblem solver context */
21 1 : ST st; /* spectral transformation context */
22 1 : KSP ksp;
23 1 : PC pc;
24 1 : PetscBool showinertia=PETSC_TRUE,lock,detect,checket;
25 1 : PetscInt n=100,Istart,Iend,i,*inertias,ns,nev,ncv,mpd;
26 1 : PetscReal mu=1.0,tau=10.0,kappa=5.0,int0,int1,*shifts;
27 :
28 1 : PetscFunctionBeginUser;
29 1 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
30 :
31 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
32 1 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-showinertia",&showinertia,NULL));
33 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum slicing on PEP, n=%" PetscInt_FMT "\n\n",n));
34 :
35 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36 : Compute the matrices that define the eigensystem, (k^2*M+k*C+K)x=0
37 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
38 :
39 : /* K is a tridiagonal */
40 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&K));
41 1 : PetscCall(MatSetSizes(K,PETSC_DECIDE,PETSC_DECIDE,n,n));
42 1 : PetscCall(MatSetFromOptions(K));
43 :
44 1 : PetscCall(MatGetOwnershipRange(K,&Istart,&Iend));
45 101 : for (i=Istart;i<Iend;i++) {
46 100 : if (i>0) PetscCall(MatSetValue(K,i,i-1,-kappa,INSERT_VALUES));
47 100 : PetscCall(MatSetValue(K,i,i,kappa*3.0,INSERT_VALUES));
48 100 : if (i<n-1) PetscCall(MatSetValue(K,i,i+1,-kappa,INSERT_VALUES));
49 : }
50 :
51 1 : PetscCall(MatAssemblyBegin(K,MAT_FINAL_ASSEMBLY));
52 1 : PetscCall(MatAssemblyEnd(K,MAT_FINAL_ASSEMBLY));
53 :
54 : /* C is a tridiagonal */
55 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&C));
56 1 : PetscCall(MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,n,n));
57 1 : PetscCall(MatSetFromOptions(C));
58 :
59 1 : PetscCall(MatGetOwnershipRange(C,&Istart,&Iend));
60 101 : for (i=Istart;i<Iend;i++) {
61 100 : if (i>0) PetscCall(MatSetValue(C,i,i-1,-tau,INSERT_VALUES));
62 100 : PetscCall(MatSetValue(C,i,i,tau*3.0,INSERT_VALUES));
63 100 : if (i<n-1) PetscCall(MatSetValue(C,i,i+1,-tau,INSERT_VALUES));
64 : }
65 :
66 1 : PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY));
67 1 : PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY));
68 :
69 : /* M is a diagonal matrix */
70 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&M));
71 1 : PetscCall(MatSetSizes(M,PETSC_DECIDE,PETSC_DECIDE,n,n));
72 1 : PetscCall(MatSetFromOptions(M));
73 1 : PetscCall(MatGetOwnershipRange(M,&Istart,&Iend));
74 101 : for (i=Istart;i<Iend;i++) PetscCall(MatSetValue(M,i,i,mu,INSERT_VALUES));
75 1 : PetscCall(MatAssemblyBegin(M,MAT_FINAL_ASSEMBLY));
76 1 : PetscCall(MatAssemblyEnd(M,MAT_FINAL_ASSEMBLY));
77 :
78 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
79 : Create the eigensolver and solve the problem
80 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
81 :
82 1 : PetscCall(PEPCreate(PETSC_COMM_WORLD,&pep));
83 1 : A[0] = K; A[1] = C; A[2] = M;
84 1 : PetscCall(PEPSetOperators(pep,3,A));
85 1 : PetscCall(PEPSetProblemType(pep,PEP_HYPERBOLIC));
86 1 : PetscCall(PEPSetType(pep,PEPSTOAR));
87 :
88 : /*
89 : Set interval and other settings for spectrum slicing
90 : */
91 1 : int0 = -11.3;
92 1 : int1 = -9.5;
93 1 : PetscCall(PEPSetInterval(pep,int0,int1));
94 1 : PetscCall(PEPSetWhichEigenpairs(pep,PEP_ALL));
95 1 : PetscCall(PEPGetST(pep,&st));
96 1 : PetscCall(STSetType(st,STSINVERT));
97 1 : PetscCall(STGetKSP(st,&ksp));
98 1 : PetscCall(KSPSetType(ksp,KSPPREONLY));
99 1 : PetscCall(KSPGetPC(ksp,&pc));
100 1 : PetscCall(PCSetType(pc,PCCHOLESKY));
101 :
102 : /*
103 : Test interface functions of STOAR solver
104 : */
105 1 : PetscCall(PEPSTOARGetDetectZeros(pep,&detect));
106 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Detect zeros before changing = %d",(int)detect));
107 1 : PetscCall(PEPSTOARSetDetectZeros(pep,PETSC_TRUE));
108 1 : PetscCall(PEPSTOARGetDetectZeros(pep,&detect));
109 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)detect));
110 :
111 1 : PetscCall(PEPSTOARGetLocking(pep,&lock));
112 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Locking flag before changing = %d",(int)lock));
113 1 : PetscCall(PEPSTOARSetLocking(pep,PETSC_TRUE));
114 1 : PetscCall(PEPSTOARGetLocking(pep,&lock));
115 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)lock));
116 :
117 1 : PetscCall(PEPSTOARGetCheckEigenvalueType(pep,&checket));
118 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Check eigenvalue type flag before changing = %d",(int)checket));
119 1 : PetscCall(PEPSTOARSetCheckEigenvalueType(pep,PETSC_FALSE));
120 1 : PetscCall(PEPSTOARGetCheckEigenvalueType(pep,&checket));
121 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)checket));
122 :
123 1 : PetscCall(PEPSTOARGetDimensions(pep,&nev,&ncv,&mpd));
124 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Sub-solve dimensions before changing = [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]",nev,ncv,mpd));
125 1 : PetscCall(PEPSTOARSetDimensions(pep,30,60,60));
126 1 : PetscCall(PEPSTOARGetDimensions(pep,&nev,&ncv,&mpd));
127 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]\n",nev,ncv,mpd));
128 :
129 1 : PetscCall(PEPSetFromOptions(pep));
130 :
131 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 : Compute all eigenvalues in interval and display info
133 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134 :
135 1 : PetscCall(PEPSetUp(pep));
136 1 : PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
137 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Inertias (after setup):\n"));
138 3 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
139 1 : PetscCall(PetscFree(shifts));
140 1 : PetscCall(PetscFree(inertias));
141 :
142 1 : PetscCall(PEPSolve(pep));
143 1 : PetscCall(PEPGetDimensions(pep,&nev,NULL,NULL));
144 1 : PetscCall(PEPGetInterval(pep,&int0,&int1));
145 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1));
146 :
147 1 : if (showinertia) {
148 0 : PetscCall(PEPSTOARGetInertias(pep,&ns,&shifts,&inertias));
149 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Used %" PetscInt_FMT " shifts (inertia):\n",ns));
150 0 : for (i=0;i<ns;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
151 0 : PetscCall(PetscFree(shifts));
152 0 : PetscCall(PetscFree(inertias));
153 : }
154 :
155 1 : PetscCall(PEPErrorView(pep,PEP_ERROR_BACKWARD,NULL));
156 :
157 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
158 : Clean up
159 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
160 1 : PetscCall(PEPDestroy(&pep));
161 1 : PetscCall(MatDestroy(&M));
162 1 : PetscCall(MatDestroy(&C));
163 1 : PetscCall(MatDestroy(&K));
164 1 : PetscCall(SlepcFinalize());
165 : return 0;
166 : }
167 :
168 : /*TEST
169 :
170 : test:
171 : requires: !single
172 : args: -showinertia 0
173 :
174 : TEST*/
|
