Actual source code: ex24.c
slepc-3.11.2 2019-07-30
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2019, 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[] = "Spectrum folding for a standard symmetric eigenproblem.\n\n"
12: "The problem matrix is the 2-D Laplacian.\n\n"
13: "The command line options are:\n"
14: " -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
15: " -m <m>, where <m> = number of grid subdivisions in y dimension.\n";
17: #include <slepceps.h>
19: /*
20: User context for spectrum folding
21: */
22: typedef struct {
23: Mat A;
24: Vec w;
25: PetscReal target;
26: } CTX_FOLD;
28: /*
29: Auxiliary routines
30: */
31: PetscErrorCode MatMult_Fold(Mat,Vec,Vec);
32: PetscErrorCode RayleighQuotient(Mat,Vec,PetscScalar*);
33: PetscErrorCode ComputeResidualNorm(Mat,PetscScalar,Vec,PetscReal*);
35: int main(int argc,char **argv)
36: {
37: Mat A,M,P; /* problem matrix, shell matrix and preconditioner */
38: Vec x; /* eigenvector */
39: EPS eps; /* eigenproblem solver context */
40: ST st; /* spectral transformation */
41: KSP ksp;
42: PC pc;
43: EPSType type;
44: CTX_FOLD *ctx;
45: PetscInt nconv,N,n=10,m,nloc,mloc,Istart,Iend,II,i,j;
46: PetscReal *error,*evals,target=0.0,tol;
47: PetscScalar lambda;
48: PetscBool flag,terse,errok;
51: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
53: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
54: PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag);
55: if (!flag) m=n;
56: PetscOptionsGetReal(NULL,NULL,"-target",&target,NULL);
57: N = n*m;
58: PetscPrintf(PETSC_COMM_WORLD,"\nSpectrum Folding, N=%D (%Dx%D grid) target=%f\n\n",N,n,m,(double)target);
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Compute the 5-point stencil Laplacian
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: MatCreate(PETSC_COMM_WORLD,&A);
65: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
66: MatSetFromOptions(A);
67: MatSetUp(A);
69: MatGetOwnershipRange(A,&Istart,&Iend);
70: for (II=Istart;II<Iend;II++) {
71: i = II/n; j = II-i*n;
72: if (i>0) { MatSetValue(A,II,II-n,-1.0,INSERT_VALUES); }
73: if (i<m-1) { MatSetValue(A,II,II+n,-1.0,INSERT_VALUES); }
74: if (j>0) { MatSetValue(A,II,II-1,-1.0,INSERT_VALUES); }
75: if (j<n-1) { MatSetValue(A,II,II+1,-1.0,INSERT_VALUES); }
76: MatSetValue(A,II,II,4.0,INSERT_VALUES);
77: }
79: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
80: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
81: MatCreateVecs(A,&x,NULL);
82: MatGetLocalSize(A,&nloc,&mloc);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create shell matrix to perform spectrum folding
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: PetscNew(&ctx);
88: ctx->A = A;
89: ctx->target = target;
90: VecDuplicate(x,&ctx->w);
92: MatCreateShell(PETSC_COMM_WORLD,nloc,mloc,N,N,ctx,&M);
93: MatShellSetOperation(M,MATOP_MULT,(void(*)(void))MatMult_Fold);
95: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96: Create the eigensolver and set various options
97: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
99: EPSCreate(PETSC_COMM_WORLD,&eps);
100: EPSSetOperators(eps,M,NULL);
101: EPSSetProblemType(eps,EPS_HEP);
102: EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL);
103: EPSSetFromOptions(eps);
105: PetscObjectTypeCompareAny((PetscObject)eps,&flag,EPSBLOPEX,EPSRQCG,"");
106: if (flag) {
107: /*
108: Build preconditioner specific for this application (diagonal of A^2)
109: */
110: MatGetRowSum(A,x);
111: VecScale(x,-1.0);
112: VecShift(x,20.0);
113: MatCreate(PETSC_COMM_WORLD,&P);
114: MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,N);
115: MatSetFromOptions(P);
116: MatSetUp(P);
117: MatDiagonalSet(P,x,INSERT_VALUES);
118: /*
119: Set diagonal preconditioner
120: */
121: EPSGetST(eps,&st);
122: STSetType(st,STPRECOND);
123: STPrecondSetMatForPC(st,P);
124: MatDestroy(&P);
125: STGetKSP(st,&ksp);
126: KSPGetPC(ksp,&pc);
127: PCSetType(pc,PCJACOBI);
128: }
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Solve the eigensystem
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: EPSSolve(eps);
135: EPSGetType(eps,&type);
136: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
137: EPSGetTolerances(eps,&tol,NULL);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Display solution and clean up
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: EPSGetConverged(eps,&nconv);
144: PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
145: if (nconv>0) {
146: PetscMalloc2(nconv,&evals,nconv,&error);
147: for (i=0;i<nconv;i++) {
148: /* Get i-th eigenvector, compute eigenvalue approximation from
149: Rayleigh quotient and compute residual norm */
150: EPSGetEigenpair(eps,i,NULL,NULL,x,NULL);
151: RayleighQuotient(A,x,&lambda);
152: ComputeResidualNorm(A,lambda,x,&error[i]);
153: #if defined(PETSC_USE_COMPLEX)
154: evals[i] = PetscRealPart(lambda);
155: #else
156: evals[i] = lambda;
157: #endif
158: }
159: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
160: if (!terse) {
161: PetscPrintf(PETSC_COMM_WORLD,
162: " k ||Ax-kx||\n"
163: " ----------------- ------------------\n");
164: for (i=0;i<nconv;i++) {
165: PetscPrintf(PETSC_COMM_WORLD," %12f %12.2g\n",(double)evals[i],(double)error[i]);
166: }
167: } else {
168: errok = PETSC_TRUE;
169: for (i=0;i<nconv;i++) errok = (errok && error[i]<5.0*tol)? PETSC_TRUE: PETSC_FALSE;
170: if (!errok) {
171: PetscPrintf(PETSC_COMM_WORLD," Problem: some of the first %D relative errors are higher than the tolerance\n\n",nconv);
172: } else {
173: PetscPrintf(PETSC_COMM_WORLD," nconv=%D eigenvalues computed up to the required tolerance:",nconv);
174: for (i=0;i<nconv;i++) {
175: PetscPrintf(PETSC_COMM_WORLD," %.5f",(double)evals[i]);
176: }
177: }
178: }
179: PetscPrintf(PETSC_COMM_WORLD,"\n");
180: PetscFree2(evals,error);
181: }
183: EPSDestroy(&eps);
184: MatDestroy(&A);
185: MatDestroy(&M);
186: VecDestroy(&ctx->w);
187: VecDestroy(&x);
188: PetscFree(ctx);
189: SlepcFinalize();
190: return ierr;
191: }
193: /*
194: Matrix-vector product subroutine for the spectrum folding.
195: y <-- (A-t*I)^2*x
196: */
197: PetscErrorCode MatMult_Fold(Mat M,Vec x,Vec y)
198: {
199: CTX_FOLD *ctx;
200: PetscScalar sigma;
204: MatShellGetContext(M,&ctx);
205: sigma = -ctx->target;
206: MatMult(ctx->A,x,ctx->w);
207: VecAXPY(ctx->w,sigma,x);
208: MatMult(ctx->A,ctx->w,y);
209: VecAXPY(y,sigma,ctx->w);
210: return(0);
211: }
213: /*
214: Computes the Rayleigh quotient of a vector x
215: r <-- x^T*A*x (assumes x has unit norm)
216: */
217: PetscErrorCode RayleighQuotient(Mat A,Vec x,PetscScalar *r)
218: {
219: Vec Ax;
223: VecDuplicate(x,&Ax);
224: MatMult(A,x,Ax);
225: VecDot(Ax,x,r);
226: VecDestroy(&Ax);
227: return(0);
228: }
230: /*
231: Computes the residual norm of an approximate eigenvector x, |A*x-lambda*x|
232: */
233: PetscErrorCode ComputeResidualNorm(Mat A,PetscScalar lambda,Vec x,PetscReal *r)
234: {
235: Vec Ax;
239: VecDuplicate(x,&Ax);
240: MatMult(A,x,Ax);
241: VecAXPY(Ax,-lambda,x);
242: VecNorm(Ax,NORM_2,r);
243: VecDestroy(&Ax);
244: return(0);
245: }
247: /*TEST
249: test:
250: suffix: 1
251: args: -n 15 -eps_nev 1 -eps_ncv 12 -eps_max_it 1000 -eps_tol 1e-5 -terse
253: TEST*/