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