Actual source code: ex19.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:       
  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>

 31: PetscErrorCode GetExactEigenvalues(PetscInt M,PetscInt N,PetscInt P,PetscInt nconv,PetscReal *exact)
 32: {
 33:   PetscInt       n,i,j,k,l;
 34:   PetscReal      *evals,ax,ay,az,sx,sy,sz;

 38:   ax = PETSC_PI/2/(M+1);
 39:   ay = PETSC_PI/2/(N+1);
 40:   az = PETSC_PI/2/(P+1);
 41:   n = ceil(pow(nconv,0.33333)+1);
 42:   PetscMalloc(n*n*n*sizeof(PetscReal),&evals);
 43:   l = 0;
 44:   for (i=1;i<=n;i++) {
 45:     sx = sin(ax*i);
 46:     for (j=1;j<=n;j++) {
 47:       sy = sin(ay*j);
 48:       for (k=1;k<=n;k++) {
 49:         sz = sin(az*k);
 50:         evals[l++] = 4.0*(sx*sx+sy*sy+sz*sz);
 51:       }
 52:     }
 53:   }
 54:   PetscSortReal(n*n*n,evals);
 55:   for (i=0;i<nconv;i++) exact[i] = evals[i];
 56:   PetscFree(evals);
 57:   return(0);
 58: }

 62: PetscErrorCode FillMatrix(DM da,Mat A)
 63: {
 65:   PetscInt       i,j,k,mx,my,mz,xm,ym,zm,xs,ys,zs,idx;
 66:   PetscScalar    v[7];
 67:   MatStencil     row,col[7];

 70:   DMDAGetInfo(da,0,&mx,&my,&mz,0,0,0,0,0,0,0,0,0);
 71:   DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);

 73:   for (k=zs;k<zs+zm;k++) {
 74:     for (j=ys;j<ys+ym;j++) {
 75:       for (i=xs;i<xs+xm;i++) {
 76:         row.i=i; row.j=j; row.k=k;
 77:         col[0].i=row.i; col[0].j=row.j; col[0].k=row.k;
 78:         v[0]=6.0;
 79:         idx=1;
 80:         if (k>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k-1; idx++; }
 81:         if (j>0) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j-1; col[idx].k=k; idx++; }
 82:         if (i>0) { v[idx]=-1.0; col[idx].i=i-1; col[idx].j=j; col[idx].k=k; idx++; }
 83:         if (i<mx-1) { v[idx]=-1.0; col[idx].i=i+1; col[idx].j=j; col[idx].k=k; idx++; }
 84:         if (j<my-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j+1; col[idx].k=k; idx++; }
 85:         if (k<mz-1) { v[idx]=-1.0; col[idx].i=i; col[idx].j=j; col[idx].k=k+1; idx++; }
 86:         MatSetValuesStencil(A,1,&row,idx,col,v,INSERT_VALUES);
 87:       }
 88:     }
 89:   }
 90:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 91:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 92:   return(0);
 93: }

 97: int main(int argc,char **argv)
 98: {
 99:   Mat            A;               /* operator matrix */
100:   EPS            eps;             /* eigenproblem solver context */
101:   const EPSType  type;
102:   DM             da;
103:   Vec            v0;
104:   PetscReal      error,tol,re,im,*exact;
105:   PetscScalar    kr,ki;
106:   PetscInt       M,N,P,m,n,p,nev,maxit,i,its,nconv,seed;
107:   PetscLogDouble t1,t2,t3;
108:   PetscBool      flg;
109:   PetscRandom    rctx;

112:   SlepcInitialize(&argc,&argv,(char*)0,help);

114:   PetscPrintf(PETSC_COMM_WORLD,"\n3-D Laplacian Eigenproblem\n\n");

116:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
117:      Compute the operator matrix that defines the eigensystem, Ax=kx
118:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

120:   DMDACreate3d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,DMDA_BOUNDARY_NONE,
121:                       DMDA_BOUNDARY_NONE,DMDA_STENCIL_STAR,-10,-10,-10,
122:                       PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,
123:                       1,1,PETSC_NULL,PETSC_NULL,PETSC_NULL,&da);

125:   /* print DM information */
126:   DMDAGetInfo(da,PETSC_NULL,&M,&N,&P,&m,&n,&p,
127:                      PETSC_NULL,PETSC_NULL,PETSC_NULL,PETSC_NULL,
128:                      PETSC_NULL,PETSC_NULL);
129:   PetscPrintf(PETSC_COMM_WORLD," Grid partitioning: %D %D %D\n",m,n,p);

131:   /* create and fill the matrix */
132:   DMGetMatrix(da,MATAIJ,&A);
133:   FillMatrix(da,A);

135:   /* create random initial vector */
136:   seed = 1;
137:   PetscOptionsGetInt(PETSC_NULL,"-seed",&seed,PETSC_NULL);
138:   if (seed<0) SETERRQ(PETSC_COMM_WORLD,1,"Seed must be >=0");
139:   MatGetVecs(A,&v0,PETSC_NULL);
140:   PetscRandomCreate(PETSC_COMM_WORLD,&rctx);
141:   PetscRandomSetFromOptions(rctx);
142:   for (i=0;i<seed;i++) {   /* simulate different seeds in the random generator */
143:     VecSetRandom(v0,rctx);
144:   }

146:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
147:                 Create the eigensolver and set various options
148:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

150:   /* 
151:      Create eigensolver context
152:   */
153:   EPSCreate(PETSC_COMM_WORLD,&eps);

155:   /* 
156:      Set operators. In this case, it is a standard eigenvalue problem
157:   */
158:   EPSSetOperators(eps,A,PETSC_NULL);
159:   EPSSetProblemType(eps,EPS_HEP);

161:   /*
162:      Set specific solver options
163:   */
164:   EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
165:   EPSSetTolerances(eps,1e-8,100);
166:   EPSSetInitialSpace(eps,1,&v0);

168:   /*
169:      Set solver parameters at runtime
170:   */
171:   EPSSetFromOptions(eps);

173:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
174:                       Solve the eigensystem
175:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

177:   PetscGetTime(&t1);
178:   EPSSetUp(eps);
179:   PetscGetTime(&t2);
180:   EPSSolve(eps);
181:   PetscGetTime(&t3);
182:   EPSGetIterationNumber(eps,&its);
183:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);

185:   /*
186:      Optional: Get some information from the solver and display it
187:   */
188:   EPSGetType(eps,&type);
189:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
190:   EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
191:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
192:   EPSGetTolerances(eps,&tol,&maxit);
193:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4G, maxit=%D\n",tol,maxit);

195:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
196:                     Display solution and clean up
197:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

199:   /* 
200:      Get number of converged approximate eigenpairs
201:   */
202:   EPSGetConverged(eps,&nconv);
203:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %D\n\n",nconv);

205:   if (nconv>0) {
206:     PetscMalloc(nconv*sizeof(PetscReal),&exact);
207:     GetExactEigenvalues(M,N,P,nconv,exact);
208:     /*
209:        Display eigenvalues and relative errors
210:     */
211:     PetscPrintf(PETSC_COMM_WORLD,
212:          "           k          ||Ax-kx||/||kx||   Eigenvalue Error \n"
213:          "   ----------------- ------------------ ------------------\n");

215:     for (i=0;i<nconv;i++) {
216:       /* 
217:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
218:         ki (imaginary part)
219:       */
220:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
221:       /*
222:          Compute the relative error associated to each eigenpair
223:       */
224:       EPSComputeRelativeError(eps,i,&error);

226: #if defined(PETSC_USE_COMPLEX)
227:       re = PetscRealPart(kr);
228:       im = PetscImaginaryPart(kr);
229: #else
230:       re = kr;
231:       im = ki;
232: #endif 
233:       if (im!=0.0) SETERRQ(PETSC_COMM_WORLD,1,"Eigenvalue should be real!");
234:       else {
235:         PetscPrintf(PETSC_COMM_WORLD,"   %12G       %12G        %12G\n",re,error,PetscAbsReal(re-exact[i]));
236:       }
237:     }
238:     PetscFree(exact);
239:     PetscPrintf(PETSC_COMM_WORLD,"\n");
240:   }
241: 
242:   /* 
243:      Show computing times
244:   */
245:   PetscOptionsHasName(PETSC_NULL,"-showtimes",&flg);
246:   if (flg) {
247:     PetscPrintf(PETSC_COMM_WORLD," Elapsed time: %G (setup), %G (solve)\n",t2-t1,t3-t2);
248:   }

250:   /* 
251:      Free work space
252:   */
253:   EPSDestroy(&eps);
254:   MatDestroy(&A);
255:   VecDestroy(&v0);
256:   PetscRandomDestroy(&rctx);
257:   DMDestroy(&da);
258:   SlepcFinalize();
259:   return 0;
260: }