Actual source code: ex19.c
slepc-3.6.3 2016-03-29
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain
6: This file is part of SLEPc.
8: SLEPc is free software: you can redistribute it and/or modify it under the
9: terms of version 3 of the GNU Lesser General Public License as published by
10: the Free Software Foundation.
12: SLEPc is distributed in the hope that it will be useful, but WITHOUT ANY
13: WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14: FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15: more details.
17: You should have received a copy of the GNU Lesser General Public License
18: along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
19: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
20: */
22: static char help[] = "Standard symmetric eigenproblem for the 3-D Laplacian built with the DM interface.\n\n"
23: "Use -seed <k> to modify the random initial vector.\n"
24: "Use -da_grid_x <nx> etc. to change the problem size.\n\n";
26: #include <slepceps.h>
27: #include <petscdmda.h>
28: #include <petsctime.h>
32: PetscErrorCode GetExactEigenvalues(PetscInt M,PetscInt N,PetscInt P,PetscInt nconv,PetscReal *exact)
33: {
34: PetscInt n,i,j,k,l;
35: PetscReal *evals,ax,ay,az,sx,sy,sz;
39: ax = PETSC_PI/2/(M+1);
40: ay = PETSC_PI/2/(N+1);
41: az = PETSC_PI/2/(P+1);
42: n = PetscCeilReal(PetscPowReal((PetscReal)nconv,0.33333)+1);
43: PetscMalloc1(n*n*n,&evals);
44: l = 0;
45: for (i=1;i<=n;i++) {
46: sx = PetscSinReal(ax*i);
47: for (j=1;j<=n;j++) {
48: sy = PetscSinReal(ay*j);
49: for (k=1;k<=n;k++) {
50: sz = PetscSinReal(az*k);
51: evals[l++] = 4.0*(sx*sx+sy*sy+sz*sz);
52: }
53: }
54: }
55: PetscSortReal(n*n*n,evals);
56: for (i=0;i<nconv;i++) exact[i] = evals[i];
57: PetscFree(evals);
58: return(0);
59: }
63: PetscErrorCode FillMatrix(DM da,Mat A)
64: {
66: PetscInt i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,idx;
67: PetscScalar v[7];
68: MatStencil row,col[7];
71: DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
72: DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);
74: for (k=zs;k<zs+zm;k++) {
75: for (j=ys;j<ys+ym;j++) {
76: for (i=xs;i<xs+xm;i++) {
77: row.i=i; row.j=j; row.k=k;
78: col[0].i=row.i; col[0].j=row.j; col[0].k=row.k;
79: v[0]=6.0;
80: idx=1;
81: if (k>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k-1; idx++; }
82: if (j>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j-1; col[idx].k=k; idx++; }
83: if (i>0) { v[idx]=-1.0; col[idx].i=i-1; col[idx].j=j; col[idx].k=k; idx++; }
84: if (i<mx-1) { v[idx]=-1.0; col[idx].i=i+1; col[idx].j=j; col[idx].k=k; idx++; }
85: if (j<my-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j+1; col[idx].k=k; idx++; }
86: if (k<mz-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k+1; idx++; }
87: MatSetValuesStencil(A,1,&row,idx,col,v,INSERT_VALUES);
88: }
89: }
90: }
91: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
92: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
93: return(0);
94: }
98: int main(int argc,char **argv)
99: {
100: Mat A; /* operator matrix */
101: EPS eps; /* eigenproblem solver context */
102: EPSType type;
103: DM da;
104: Vec v0;
105: PetscReal error,tol,re,im,*exact;
106: PetscScalar kr,ki;
107: PetscInt M,N,P,m,n,p,nev,maxit,i,its,nconv,seed;
108: PetscLogDouble t1,t2,t3;
109: PetscBool flg;
110: PetscRandom rctx;
113: SlepcInitialize(&argc,&argv,(char*)0,help);
115: PetscPrintf(PETSC_COMM_WORLD,"\n3-D Laplacian Eigenproblem\n\n");
117: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118: Compute the operator matrix that defines the eigensystem, Ax=kx
119: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121: DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,
122: DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,-10,-10,-10,
123: PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
124: 1,1,NULL,NULL,NULL,&da);
126: /* print DM information */
127: DMDAGetInfo(da,NULL,&M,&N,&P,&m,&n,&p,NULL,NULL,NULL,NULL,NULL,NULL);
128: PetscPrintf(PETSC_COMM_WORLD," Grid partitioning: %D %D %D\n",m,n,p);
130: /* create and fill the matrix */
131: DMCreateMatrix(da,&A);
132: FillMatrix(da,A);
134: /* create random initial vector */
135: seed = 1;
136: PetscOptionsGetInt(NULL,"-seed",&seed,NULL);
137: if (seed<0) SETERRQ(PETSC_COMM_WORLD,1,"Seed must be >=0");
138: MatCreateVecs(A,&v0,NULL);
139: PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
140: PetscRandomSetFromOptions(rctx);
141: for (i=0;i<seed;i++) { /* simulate different seeds in the random generator */
142: VecSetRandom(v0,rctx);
143: }
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Create the eigensolver and set various options
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: /*
150: Create eigensolver context
151: */
152: EPSCreate(PETSC_COMM_WORLD,&eps);
154: /*
155: Set operators. In this case, it is a standard eigenvalue problem
156: */
157: EPSSetOperators(eps,A,NULL);
158: EPSSetProblemType(eps,EPS_HEP);
160: /*
161: Set specific solver options
162: */
163: EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
164: EPSSetTolerances(eps,1e-8,PETSC_DEFAULT);
165: EPSSetInitialSpace(eps,1,&v0);
167: /*
168: Set solver parameters at runtime
169: */
170: EPSSetFromOptions(eps);
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Solve the eigensystem
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
176: PetscTime(&t1);
177: EPSSetUp(eps);
178: PetscTime(&t2);
179: EPSSolve(eps);
180: PetscTime(&t3);
181: EPSGetIterationNumber(eps,&its);
182: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
184: /*
185: Optional: Get some information from the solver and display it
186: */
187: EPSGetType(eps,&type);
188: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
189: EPSGetDimensions(eps,&nev,NULL,NULL);
190: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
191: EPSGetTolerances(eps,&tol,&maxit);
192: PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);
194: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195: Display solution and clean up
196: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198: /*
199: Get number of converged approximate eigenpairs
200: */
201: EPSGetConverged(eps,&nconv);
202: PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);
204: if (nconv>0) {
205: PetscMalloc1(nconv,&exact);
206: GetExactEigenvalues(M,N,P,nconv,exact);
207: /*
208: Display eigenvalues and relative errors
209: */
210: PetscPrintf(PETSC_COMM_WORLD,
211: " k ||Ax-kx||/||kx|| Eigenvalue Error \n"
212: " ----------------- ------------------ ------------------\n");
214: for (i=0;i<nconv;i++) {
215: /*
216: Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
217: ki (imaginary part)
218: */
219: EPSGetEigenpair(eps,i,&kr,&ki,NULL,NULL);
220: /*
221: Compute the relative error associated to each eigenpair
222: */
223: EPSComputeError(eps,i,EPS_ERROR_RELATIVE,&error);
225: #if defined(PETSC_USE_COMPLEX)
226: re = PetscRealPart(kr);
227: im = PetscImaginaryPart(kr);
228: #else
229: re = kr;
230: im = ki;
231: #endif
232: if (im!=0.0) SETERRQ(PETSC_COMM_WORLD,1,"Eigenvalue should be real");
233: else {
234: PetscPrintf(PETSC_COMM_WORLD," %12g %12g %12g\n",(double)re,(double)error,(double)PetscAbsReal(re-exact[i]));
235: }
236: }
237: PetscFree(exact);
238: PetscPrintf(PETSC_COMM_WORLD,"\n");
239: }
241: /*
242: Show computing times
243: */
244: PetscOptionsHasName(NULL,"-showtimes",&flg);
245: if (flg) {
246: PetscPrintf(PETSC_COMM_WORLD," Elapsed time: %g (setup), %g (solve)\n",(double)(t2-t1),(double)(t3-t2));
247: }
249: /*
250: Free work space
251: */
252: EPSDestroy(&eps);
253: MatDestroy(&A);
254: VecDestroy(&v0);
255: PetscRandomDestroy(&rctx);
256: DMDestroy(&da);
257: SlepcFinalize();
258: return 0;
259: }