Actual source code: ex19.c

slepc-3.7.3 2016-09-29
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2016, 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,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,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 ierr;
259: }