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