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