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