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[] = "Standard symmetric eigenproblem for the 3-D Laplacian built with the DM interface.\n\n"
12 : "Use -seed <k> to modify the random initial vector.\n"
13 : "Use -da_grid_x <nx> etc. to change the problem size.\n\n";
14 :
15 : #include <slepceps.h>
16 : #include <petscdmda.h>
17 : #include <petsctime.h>
18 :
19 0 : PetscErrorCode GetExactEigenvalues(PetscInt M,PetscInt N,PetscInt P,PetscInt nconv,PetscReal *exact)
20 : {
21 0 : PetscInt n,i,j,k,l;
22 0 : PetscReal *evals,ax,ay,az,sx,sy,sz;
23 :
24 0 : PetscFunctionBeginUser;
25 0 : ax = PETSC_PI/2/(M+1);
26 0 : ay = PETSC_PI/2/(N+1);
27 0 : az = PETSC_PI/2/(P+1);
28 0 : n = PetscCeilReal(PetscPowReal((PetscReal)nconv,0.33333)+1);
29 0 : PetscCall(PetscMalloc1(n*n*n,&evals));
30 : l = 0;
31 0 : for (i=1;i<=n;i++) {
32 0 : sx = PetscSinReal(ax*i);
33 0 : for (j=1;j<=n;j++) {
34 0 : sy = PetscSinReal(ay*j);
35 0 : for (k=1;k<=n;k++) {
36 0 : sz = PetscSinReal(az*k);
37 0 : evals[l++] = 4.0*(sx*sx+sy*sy+sz*sz);
38 : }
39 : }
40 : }
41 0 : PetscCall(PetscSortReal(n*n*n,evals));
42 0 : for (i=0;i<nconv;i++) exact[i] = evals[i];
43 0 : PetscCall(PetscFree(evals));
44 0 : PetscFunctionReturn(PETSC_SUCCESS);
45 : }
46 :
47 3 : PetscErrorCode FillMatrix(DM da,Mat A)
48 : {
49 3 : PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,idx;
50 3 : PetscScalar v[7];
51 3 : MatStencil row,col[7];
52 :
53 3 : PetscFunctionBeginUser;
54 3 : PetscCall(DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0));
55 3 : PetscCall(DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm));
56 :
57 33 : for (k=zs;k<zs+zm;k++) {
58 330 : for (j=ys;j<ys+ym;j++) {
59 3300 : for (i=xs;i<xs+xm;i++) {
60 3000 : row.i=i; row.j=j; row.k=k;
61 3000 : col[0].i=row.i; col[0].j=row.j; col[0].k=row.k;
62 3000 : v[0]=6.0;
63 3000 : idx=1;
64 3000 : if (k>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k-1; idx++; }
65 3000 : if (j>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j-1; col[idx].k=k; idx++; }
66 3000 : if (i>0) { v[idx]=-1.0; col[idx].i=i-1; col[idx].j=j; col[idx].k=k; idx++; }
67 3000 : if (i<mx-1) { v[idx]=-1.0; col[idx].i=i+1; col[idx].j=j; col[idx].k=k; idx++; }
68 3000 : if (j<my-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j+1; col[idx].k=k; idx++; }
69 3000 : if (k<mz-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k+1; idx++; }
70 3000 : PetscCall(MatSetValuesStencil(A,1,&row,idx,col,v,INSERT_VALUES));
71 : }
72 : }
73 : }
74 3 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
75 3 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
76 3 : PetscFunctionReturn(PETSC_SUCCESS);
77 : }
78 :
79 3 : int main(int argc,char **argv)
80 : {
81 3 : Mat A; /* operator matrix */
82 3 : EPS eps; /* eigenproblem solver context */
83 3 : EPSType type;
84 3 : DM da;
85 3 : Vec v0;
86 3 : PetscReal error,tol,re,im,*exact;
87 3 : PetscScalar kr,ki;
88 3 : PetscInt M,N,P,m,n,p,nev,maxit,i,its,nconv,seed;
89 3 : PetscLogDouble t1,t2,t3;
90 3 : PetscBool flg,terse;
91 3 : PetscRandom rctx;
92 :
93 3 : PetscFunctionBeginUser;
94 3 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
95 :
96 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n3-D Laplacian Eigenproblem\n\n"));
97 :
98 : /* show detailed info unless -terse option is given by user */
99 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
100 :
101 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102 : Compute the operator matrix that defines the eigensystem, Ax=kx
103 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 :
105 3 : PetscCall(DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
106 : DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,10,10,10,
107 : PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
108 : 1,1,NULL,NULL,NULL,&da));
109 3 : PetscCall(DMSetFromOptions(da));
110 3 : PetscCall(DMSetUp(da));
111 :
112 : /* print DM information */
113 3 : PetscCall(DMDAGetInfo(da,NULL,&M,&N,&P,&m,&n,&p,NULL,NULL,NULL,NULL,NULL,NULL));
114 3 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Grid partitioning: %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",m,n,p));
115 :
116 : /* create and fill the matrix */
117 3 : PetscCall(DMCreateMatrix(da,&A));
118 3 : PetscCall(FillMatrix(da,A));
119 :
120 : /* create random initial vector */
121 3 : seed = 1;
122 3 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL));
123 3 : PetscCheck(seed>=0,PETSC_COMM_WORLD,PETSC_ERR_USER_INPUT,"Seed must be >=0");
124 3 : PetscCall(MatCreateVecs(A,&v0,NULL));
125 3 : PetscCall(PetscRandomCreate(PETSC_COMM_WORLD,&rctx));
126 3 : PetscCall(PetscRandomSetFromOptions(rctx));
127 6 : for (i=0;i<seed;i++) { /* simulate different seeds in the random generator */
128 3 : PetscCall(VecSetRandom(v0,rctx));
129 : }
130 :
131 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132 : Create the eigensolver and set various options
133 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134 :
135 : /*
136 : Create eigensolver context
137 : */
138 3 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
139 :
140 : /*
141 : Set operators. In this case, it is a standard eigenvalue problem
142 : */
143 3 : PetscCall(EPSSetOperators(eps,A,NULL));
144 3 : PetscCall(EPSSetProblemType(eps,EPS_HEP));
145 :
146 : /*
147 : Set specific solver options
148 : */
149 3 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
150 3 : PetscCall(EPSSetTolerances(eps,1e-8,PETSC_DEFAULT));
151 3 : PetscCall(EPSSetInitialSpace(eps,1,&v0));
152 :
153 : /*
154 : Set solver parameters at runtime
155 : */
156 3 : PetscCall(EPSSetFromOptions(eps));
157 :
158 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159 : Solve the eigensystem
160 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161 :
162 3 : PetscCall(PetscTime(&t1));
163 3 : PetscCall(EPSSetUp(eps));
164 3 : PetscCall(PetscTime(&t2));
165 3 : PetscCall(EPSSolve(eps));
166 3 : PetscCall(PetscTime(&t3));
167 3 : if (!terse) {
168 0 : PetscCall(EPSGetIterationNumber(eps,&its));
169 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its));
170 :
171 : /*
172 : Optional: Get some information from the solver and display it
173 : */
174 0 : PetscCall(EPSGetType(eps,&type));
175 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
176 0 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
177 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
178 0 : PetscCall(EPSGetTolerances(eps,&tol,&maxit));
179 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit));
180 : }
181 :
182 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183 : Display solution and clean up
184 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185 :
186 3 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
187 : else {
188 : /*
189 : Get number of converged approximate eigenpairs
190 : */
191 0 : PetscCall(EPSGetConverged(eps,&nconv));
192 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %" PetscInt_FMT "\n\n",nconv));
193 :
194 0 : if (nconv>0) {
195 0 : PetscCall(PetscMalloc1(nconv,&exact));
196 0 : PetscCall(GetExactEigenvalues(M,N,P,nconv,exact));
197 : /*
198 : Display eigenvalues and relative errors
199 : */
200 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,
201 : " k ||Ax-kx||/||kx|| Eigenvalue Error \n"
202 : " ----------------- ------------------ ------------------\n"));
203 :
204 0 : for (i=0;i<nconv;i++) {
205 : /*
206 : Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
207 : ki (imaginary part)
208 : */
209 0 : PetscCall(EPSGetEigenpair(eps,i,&kr,&ki,NULL,NULL));
210 : /*
211 : Compute the relative error associated to each eigenpair
212 : */
213 0 : PetscCall(EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error));
214 :
215 : #if defined(PETSC_USE_COMPLEX)
216 : re = PetscRealPart(kr);
217 : im = PetscImaginaryPart(kr);
218 : #else
219 0 : re = kr;
220 0 : im = ki;
221 : #endif
222 0 : PetscCheck(im==0.0,PETSC_COMM_WORLD,PETSC_ERR_PLIB,"Eigenvalue should be real");
223 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," %12g %12g %12g\n",(double)re,(double)error,(double)PetscAbsReal(re-exact[i])));
224 : }
225 0 : PetscCall(PetscFree(exact));
226 0 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
227 : }
228 : }
229 :
230 : /*
231 : Show computing times
232 : */
233 3 : PetscCall(PetscOptionsHasName(NULL,NULL,"-showtimes",&flg));
234 3 : if (flg) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Elapsed time: %g (setup), %g (solve)\n",(double)(t2-t1),(double)(t3-t2)));
235 :
236 : /*
237 : Free work space
238 : */
239 3 : PetscCall(EPSDestroy(&eps));
240 3 : PetscCall(MatDestroy(&A));
241 3 : PetscCall(VecDestroy(&v0));
242 3 : PetscCall(PetscRandomDestroy(&rctx));
243 3 : PetscCall(DMDestroy(&da));
244 3 : PetscCall(SlepcFinalize());
245 : return 0;
246 : }
247 :
248 : /*TEST
249 :
250 : testset:
251 : args: -eps_nev 8 -terse
252 : requires: double
253 : output_file: output/ex19_1.out
254 : test:
255 : suffix: 1_krylovschur
256 : args: -eps_type krylovschur -eps_ncv 64
257 : test:
258 : suffix: 1_lobpcg
259 : args: -eps_type lobpcg -eps_tol 1e-7
260 : test:
261 : suffix: 1_blopex
262 : args: -eps_type blopex -eps_tol 1e-7 -eps_blopex_blocksize 4
263 : requires: blopex
264 :
265 : TEST*/
|
