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 Krylov-Schur.\n\n"
12 : "This is based on ex12.c. 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 3 : int main(int argc,char **argv)
19 : {
20 3 : Mat A,B; /* matrices */
21 3 : Mat As,Bs; /* matrices distributed in subcommunicators */
22 3 : Mat Au; /* matrix used to modify A on subcommunicators */
23 3 : EPS eps; /* eigenproblem solver context */
24 3 : ST st; /* spectral transformation context */
25 3 : KSP ksp;
26 3 : PC pc;
27 3 : Vec v;
28 3 : PetscMPIInt size,rank;
29 3 : PetscInt N,n=35,m,Istart,Iend,II,nev,ncv,mpd,i,j,k,*inertias,npart,nval,nloc,nlocs,mlocs;
30 3 : PetscBool flag,showinertia=PETSC_TRUE,lock,detect;
31 3 : PetscReal int0,int1,*shifts,keep,*subint,*evals;
32 3 : PetscScalar lambda;
33 3 : char vlist[4000];
34 :
35 3 : PetscFunctionBeginUser;
36 3 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
37 3 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
38 3 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
39 :
40 3 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-showinertia",&showinertia,NULL));
41 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
42 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
43 3 : if (!flag) m=n;
44 3 : N = n*m;
45 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum-slicing test, N=%" PetscInt_FMT " (%" PetscInt_FMT "x%" PetscInt_FMT " grid)\n\n",N,n,m));
46 :
47 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
48 : Compute the matrices that define the eigensystem, Ax=kBx
49 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
50 :
51 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
52 3 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
53 3 : PetscCall(MatSetFromOptions(A));
54 :
55 3 : PetscCall(MatCreate(PETSC_COMM_WORLD,&B));
56 3 : PetscCall(MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N));
57 3 : PetscCall(MatSetFromOptions(B));
58 :
59 3 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
60 2453 : for (II=Istart;II<Iend;II++) {
61 2450 : i = II/n; j = II-i*n;
62 2450 : if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
63 2450 : if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
64 2450 : if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
65 2450 : if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
66 2450 : PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
67 2450 : PetscCall(MatSetValue(B,II,II,2.0,INSERT_VALUES));
68 : }
69 3 : if (Istart==0) {
70 2 : PetscCall(MatSetValue(B,0,0,6.0,INSERT_VALUES));
71 2 : PetscCall(MatSetValue(B,0,1,-1.0,INSERT_VALUES));
72 2 : PetscCall(MatSetValue(B,1,0,-1.0,INSERT_VALUES));
73 2 : PetscCall(MatSetValue(B,1,1,1.0,INSERT_VALUES));
74 : }
75 :
76 3 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
77 3 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
78 3 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
79 3 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
80 :
81 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
82 : Create the eigensolver and set various options
83 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
84 :
85 3 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
86 3 : PetscCall(EPSSetOperators(eps,A,B));
87 3 : PetscCall(EPSSetProblemType(eps,EPS_GHEP));
88 3 : PetscCall(EPSSetType(eps,EPSKRYLOVSCHUR));
89 :
90 : /*
91 : Set interval and other settings for spectrum slicing
92 : */
93 3 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_ALL));
94 3 : int0 = 1.1; int1 = 1.3;
95 3 : PetscCall(EPSSetInterval(eps,int0,int1));
96 3 : PetscCall(EPSGetST(eps,&st));
97 3 : PetscCall(STSetType(st,STSINVERT));
98 3 : if (size>1) PetscCall(EPSKrylovSchurSetPartitions(eps,size));
99 3 : PetscCall(EPSKrylovSchurGetKSP(eps,&ksp));
100 3 : PetscCall(KSPGetPC(ksp,&pc));
101 3 : PetscCall(KSPSetType(ksp,KSPPREONLY));
102 3 : PetscCall(PCSetType(pc,PCCHOLESKY));
103 :
104 : /*
105 : Test interface functions of Krylov-Schur solver
106 : */
107 3 : PetscCall(EPSKrylovSchurGetRestart(eps,&keep));
108 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Restart parameter before changing = %g",(double)keep));
109 3 : PetscCall(EPSKrylovSchurSetRestart(eps,0.4));
110 3 : PetscCall(EPSKrylovSchurGetRestart(eps,&keep));
111 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %g\n",(double)keep));
112 :
113 3 : PetscCall(EPSKrylovSchurGetDetectZeros(eps,&detect));
114 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Detect zeros before changing = %d",(int)detect));
115 3 : PetscCall(EPSKrylovSchurSetDetectZeros(eps,PETSC_TRUE));
116 3 : PetscCall(EPSKrylovSchurGetDetectZeros(eps,&detect));
117 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)detect));
118 :
119 3 : PetscCall(EPSKrylovSchurGetLocking(eps,&lock));
120 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Locking flag before changing = %d",(int)lock));
121 3 : PetscCall(EPSKrylovSchurSetLocking(eps,PETSC_FALSE));
122 3 : PetscCall(EPSKrylovSchurGetLocking(eps,&lock));
123 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to %d\n",(int)lock));
124 :
125 3 : PetscCall(EPSKrylovSchurGetDimensions(eps,&nev,&ncv,&mpd));
126 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Sub-solve dimensions before changing = [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]",nev,ncv,mpd));
127 3 : PetscCall(EPSKrylovSchurSetDimensions(eps,30,60,60));
128 3 : PetscCall(EPSKrylovSchurGetDimensions(eps,&nev,&ncv,&mpd));
129 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," ... changed to [%" PetscInt_FMT ",%" PetscInt_FMT ",%" PetscInt_FMT "]\n",nev,ncv,mpd));
130 :
131 3 : if (size>1) {
132 2 : PetscCall(EPSKrylovSchurGetPartitions(eps,&npart));
133 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using %" PetscInt_FMT " partitions\n",npart));
134 :
135 2 : PetscCall(PetscMalloc1(npart+1,&subint));
136 2 : subint[0] = int0;
137 2 : subint[npart] = int1;
138 4 : for (i=1;i<npart;i++) subint[i] = int0+i*(int1-int0)/npart;
139 2 : PetscCall(EPSKrylovSchurSetSubintervals(eps,subint));
140 2 : PetscCall(PetscFree(subint));
141 2 : PetscCall(EPSKrylovSchurGetSubintervals(eps,&subint));
142 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using sub-interval separations = "));
143 4 : for (i=1;i<npart;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %g",(double)subint[i]));
144 2 : PetscCall(PetscFree(subint));
145 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
146 : }
147 :
148 3 : PetscCall(EPSSetFromOptions(eps));
149 :
150 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151 : Compute all eigenvalues in interval and display info
152 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153 :
154 3 : PetscCall(EPSSetUp(eps));
155 3 : PetscCall(EPSKrylovSchurGetInertias(eps,&k,&shifts,&inertias));
156 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Inertias after EPSSetUp:\n"));
157 11 : for (i=0;i<k;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
158 3 : PetscCall(PetscFree(shifts));
159 3 : PetscCall(PetscFree(inertias));
160 :
161 3 : PetscCall(EPSSolve(eps));
162 3 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
163 3 : PetscCall(EPSGetInterval(eps,&int0,&int1));
164 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1));
165 :
166 3 : if (showinertia) {
167 0 : PetscCall(EPSKrylovSchurGetInertias(eps,&k,&shifts,&inertias));
168 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Used %" PetscInt_FMT " shifts (inertia):\n",k));
169 0 : for (i=0;i<k;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," .. %g (%" PetscInt_FMT ")\n",(double)shifts[i],inertias[i]));
170 0 : PetscCall(PetscFree(shifts));
171 0 : PetscCall(PetscFree(inertias));
172 : }
173 :
174 3 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
175 :
176 3 : if (size>1) {
177 2 : PetscCall(EPSKrylovSchurGetSubcommInfo(eps,&k,&nval,&v));
178 2 : PetscCall(PetscMalloc1(nval,&evals));
179 60 : for (i=0;i<nval;i++) {
180 58 : PetscCall(EPSKrylovSchurGetSubcommPairs(eps,i,&lambda,v));
181 58 : evals[i] = PetscRealPart(lambda);
182 : }
183 2 : PetscCall(PetscFormatRealArray(vlist,sizeof(vlist),"%f",nval,evals));
184 2 : PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD," Process %d has worked in sub-interval %" PetscInt_FMT ", containing %" PetscInt_FMT " eigenvalues: %s\n",(int)rank,k,nval,vlist));
185 2 : PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT));
186 2 : PetscCall(VecDestroy(&v));
187 2 : PetscCall(PetscFree(evals));
188 :
189 2 : PetscCall(EPSKrylovSchurGetSubcommMats(eps,&As,&Bs));
190 2 : PetscCall(MatGetLocalSize(A,&nloc,NULL));
191 2 : PetscCall(MatGetLocalSize(As,&nlocs,&mlocs));
192 2 : PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD," Process %d owns %" PetscInt_FMT " rows of the global matrices, and %" PetscInt_FMT " rows in the subcommunicator\n",(int)rank,nloc,nlocs));
193 2 : PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT));
194 :
195 : /* modify A on subcommunicators */
196 2 : PetscCall(MatCreate(PetscObjectComm((PetscObject)As),&Au));
197 2 : PetscCall(MatSetSizes(Au,nlocs,mlocs,N,N));
198 2 : PetscCall(MatSetFromOptions(Au));
199 2 : PetscCall(MatGetOwnershipRange(Au,&Istart,&Iend));
200 2452 : for (II=Istart;II<Iend;II++) PetscCall(MatSetValue(Au,II,II,0.5,INSERT_VALUES));
201 2 : PetscCall(MatAssemblyBegin(Au,MAT_FINAL_ASSEMBLY));
202 2 : PetscCall(MatAssemblyEnd(Au,MAT_FINAL_ASSEMBLY));
203 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Updating internal matrices\n"));
204 2 : PetscCall(EPSKrylovSchurUpdateSubcommMats(eps,1.1,-5.0,Au,1.0,0.0,NULL,DIFFERENT_NONZERO_PATTERN,PETSC_TRUE));
205 2 : PetscCall(MatDestroy(&Au));
206 2 : PetscCall(EPSSolve(eps));
207 2 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
208 2 : PetscCall(EPSGetInterval(eps,&int0,&int1));
209 2 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," After update, found %" PetscInt_FMT " eigenvalues in interval [%g,%g]\n",nev,(double)int0,(double)int1));
210 : }
211 3 : PetscCall(EPSDestroy(&eps));
212 3 : PetscCall(MatDestroy(&A));
213 3 : PetscCall(MatDestroy(&B));
214 3 : PetscCall(SlepcFinalize());
215 : return 0;
216 : }
217 :
218 : /*TEST
219 :
220 : test:
221 : suffix: 1
222 : nsize: 2
223 : args: -showinertia 0 -log_exclude eps,st,rg,bv,ds
224 : requires: !single
225 :
226 : test:
227 : suffix: 2
228 : nsize: 1
229 : args: -showinertia 0 -log_exclude eps,st,rg,bv,ds
230 : requires: !single
231 :
232 : TEST*/
|
