Actual source code: ex24.c
slepc-3.6.1 2015-09-03
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-2015, 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,its,i,j;
59: PetscReal error,re,target=2.1;
60: PetscScalar lambda;
61: PetscBool flag;
64: SlepcInitialize(&argc,&argv,(char*)0,help);
66: PetscOptionsGetInt(NULL,"-n",&n,NULL);
67: PetscOptionsGetInt(NULL,"-m",&m,&flag);
68: if (!flag) m=n;
69: PetscOptionsGetReal(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);
149: EPSGetIterationNumber(eps,&its);
150: PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
151: EPSGetType(eps,&type);
152: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
154: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
155: Display solution and clean up
156: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: EPSGetConverged(eps,&nconv);
159: PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %D\n\n",nconv);
160: if (nconv>0) {
161: /*
162: Display result
163: */
164: PetscPrintf(PETSC_COMM_WORLD,
165: " k ||Ax-kx||\n"
166: " ----------------- ------------------\n");
168: for (i=0;i<nconv;i++) {
169: /*
170: Get i-th eigenvector, compute eigenvalue approximation from
171: Rayleigh quotient and compute residual norm
172: */
173: EPSGetEigenpair(eps,i,NULL,NULL,x,NULL);
174: RayleighQuotient(A,x,&lambda);
175: ComputeResidualNorm(A,lambda,x,&error);
177: #if defined(PETSC_USE_COMPLEX)
178: re = PetscRealPart(lambda);
179: #else
180: re = lambda;
181: #endif
182: PetscPrintf(PETSC_COMM_WORLD," %12f %12.2g\n",(double)re,(double)error);
183: }
184: PetscPrintf(PETSC_COMM_WORLD,"\n");
185: }
187: EPSDestroy(&eps);
188: MatDestroy(&A);
189: MatDestroy(&M);
190: VecDestroy(&ctx->w);
191: VecDestroy(&x);
192: PetscFree(ctx);
193: SlepcFinalize();
194: return 0;
195: }
199: /*
200: Matrix-vector product subroutine for the spectrum folding.
201: y <-- (A-t*I)^2*x
202: */
203: PetscErrorCode MatMult_Fold(Mat M,Vec x,Vec y)
204: {
205: CTX_FOLD *ctx;
206: PetscScalar sigma;
210: MatShellGetContext(M,&ctx);
211: sigma = -ctx->target;
212: MatMult(ctx->A,x,ctx->w);
213: VecAXPY(ctx->w,sigma,x);
214: MatMult(ctx->A,ctx->w,y);
215: VecAXPY(y,sigma,ctx->w);
216: return(0);
217: }
221: /*
222: Computes the Rayleigh quotient of a vector x
223: r <-- x^T*A*x (assumes x has unit norm)
224: */
225: PetscErrorCode RayleighQuotient(Mat A,Vec x,PetscScalar *r)
226: {
227: Vec Ax;
231: VecDuplicate(x,&Ax);
232: MatMult(A,x,Ax);
233: VecDot(Ax,x,r);
234: VecDestroy(&Ax);
235: return(0);
236: }
240: /*
241: Computes the residual norm of an approximate eigenvector x, |A*x-lambda*x|
242: */
243: PetscErrorCode ComputeResidualNorm(Mat A,PetscScalar lambda,Vec x,PetscReal *r)
244: {
245: Vec Ax;
249: VecDuplicate(x,&Ax);
250: MatMult(A,x,Ax);
251: VecAXPY(Ax,-lambda,x);
252: VecNorm(Ax,NORM_2,r);
253: VecDestroy(&Ax);
254: return(0);
255: }