Actual source code: ex19.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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,PETSC_ERR_USER_INPUT,"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,PETSC_ERR_PLIB,"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
267:       test:
268:          suffix: 1_blopex
269:          args: -eps_type blopex -eps_tol 1e-7 -eps_blopex_blocksize 4
270:          requires: blopex

272: TEST*/