Actual source code: ex24.c
slepc-3.21.0 2024-03-30
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[] = "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,hasmat;
50: PetscFunctionBeginUser;
51: PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
53: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
54: PetscCall(PetscOptionsGetInt(NULL,NULL,"-m",&m,&flag));
55: if (!flag) m=n;
56: PetscCall(PetscOptionsGetReal(NULL,NULL,"-target",&target,NULL));
57: N = n*m;
58: 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));
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Compute the 5-point stencil Laplacian
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
65: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N));
66: PetscCall(MatSetFromOptions(A));
68: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
69: for (II=Istart;II<Iend;II++) {
70: i = II/n; j = II-i*n;
71: if (i>0) PetscCall(MatSetValue(A,II,II-n,-1.0,INSERT_VALUES));
72: if (i<m-1) PetscCall(MatSetValue(A,II,II+n,-1.0,INSERT_VALUES));
73: if (j>0) PetscCall(MatSetValue(A,II,II-1,-1.0,INSERT_VALUES));
74: if (j<n-1) PetscCall(MatSetValue(A,II,II+1,-1.0,INSERT_VALUES));
75: PetscCall(MatSetValue(A,II,II,4.0,INSERT_VALUES));
76: }
78: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
79: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
80: PetscCall(MatCreateVecs(A,&x,NULL));
81: PetscCall(MatGetLocalSize(A,&nloc,&mloc));
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create shell matrix to perform spectrum folding
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: PetscCall(PetscNew(&ctx));
87: ctx->A = A;
88: ctx->target = target;
89: PetscCall(VecDuplicate(x,&ctx->w));
91: PetscCall(MatCreateShell(PETSC_COMM_WORLD,nloc,mloc,N,N,ctx,&M));
92: PetscCall(MatShellSetOperation(M,MATOP_MULT,(void(*)(void))MatMult_Fold));
94: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
95: Create the eigensolver and set various options
96: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
99: PetscCall(EPSSetOperators(eps,M,NULL));
100: PetscCall(EPSSetProblemType(eps,EPS_HEP));
101: PetscCall(EPSSetWhichEigenpairs(eps,EPS_SMALLEST_REAL));
102: PetscCall(EPSSetFromOptions(eps));
104: PetscCall(PetscObjectTypeCompareAny((PetscObject)eps,&flag,EPSGD,EPSJD,EPSBLOPEX,EPSLOBPCG,EPSRQCG,""));
105: if (flag) {
106: /*
107: Build preconditioner specific for this application (diagonal of A^2)
108: */
109: PetscCall(MatGetRowSum(A,x));
110: PetscCall(VecScale(x,-1.0));
111: PetscCall(VecShift(x,20.0));
112: PetscCall(MatCreate(PETSC_COMM_WORLD,&P));
113: PetscCall(MatSetSizes(P,PETSC_DECIDE,PETSC_DECIDE,N,N));
114: PetscCall(MatSetFromOptions(P));
115: PetscCall(MatDiagonalSet(P,x,INSERT_VALUES));
116: /*
117: Set diagonal preconditioner
118: */
119: PetscCall(EPSGetST(eps,&st));
120: PetscCall(STSetType(st,STPRECOND));
121: PetscCall(STSetPreconditionerMat(st,P));
122: PetscCall(MatDestroy(&P));
123: PetscCall(STGetKSP(st,&ksp));
124: PetscCall(KSPGetPC(ksp,&pc));
125: PetscCall(PCSetType(pc,PCJACOBI));
126: PetscCall(STPrecondGetKSPHasMat(st,&hasmat));
127: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Preconditioned solver, hasmat=%s\n",hasmat?"true":"false"));
128: }
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Solve the eigensystem
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: PetscCall(EPSSolve(eps));
135: PetscCall(EPSGetType(eps,&type));
136: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
137: PetscCall(EPSGetTolerances(eps,&tol,NULL));
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Display solution and clean up
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: PetscCall(EPSGetConverged(eps,&nconv));
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %" PetscInt_FMT "\n\n",nconv));
145: if (nconv>0) {
146: PetscCall(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: PetscCall(EPSGetEigenpair(eps,i,NULL,NULL,x,NULL));
151: PetscCall(RayleighQuotient(A,x,&lambda));
152: PetscCall(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: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
160: if (!terse) {
161: PetscCall(PetscPrintf(PETSC_COMM_WORLD,
162: " k ||Ax-kx||\n"
163: " ----------------- ------------------\n"));
164: 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: for (i=0;i<nconv;i++) errok = (errok && error[i]<5.0*tol)? PETSC_TRUE: PETSC_FALSE;
168: 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: PetscCall(PetscPrintf(PETSC_COMM_WORLD," nconv=%" PetscInt_FMT " eigenvalues computed up to the required tolerance:",nconv));
171: for (i=0;i<nconv;i++) PetscCall(PetscPrintf(PETSC_COMM_WORLD," %.5f",(double)evals[i]));
172: }
173: }
174: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
175: PetscCall(PetscFree2(evals,error));
176: }
178: PetscCall(EPSDestroy(&eps));
179: PetscCall(MatDestroy(&A));
180: PetscCall(MatDestroy(&M));
181: PetscCall(VecDestroy(&ctx->w));
182: PetscCall(VecDestroy(&x));
183: PetscCall(PetscFree(ctx));
184: PetscCall(SlepcFinalize());
185: return 0;
186: }
188: /*
189: Matrix-vector product subroutine for the spectrum folding.
190: y <-- (A-t*I)^2*x
191: */
192: PetscErrorCode MatMult_Fold(Mat M,Vec x,Vec y)
193: {
194: CTX_FOLD *ctx;
195: PetscScalar sigma;
197: PetscFunctionBeginUser;
198: PetscCall(MatShellGetContext(M,&ctx));
199: sigma = -ctx->target;
200: PetscCall(MatMult(ctx->A,x,ctx->w));
201: PetscCall(VecAXPY(ctx->w,sigma,x));
202: PetscCall(MatMult(ctx->A,ctx->w,y));
203: PetscCall(VecAXPY(y,sigma,ctx->w));
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: /*
208: Computes the Rayleigh quotient of a vector x
209: r <-- x^T*A*x (assumes x has unit norm)
210: */
211: PetscErrorCode RayleighQuotient(Mat A,Vec x,PetscScalar *r)
212: {
213: Vec Ax;
215: PetscFunctionBeginUser;
216: PetscCall(VecDuplicate(x,&Ax));
217: PetscCall(MatMult(A,x,Ax));
218: PetscCall(VecDot(Ax,x,r));
219: PetscCall(VecDestroy(&Ax));
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: Computes the residual norm of an approximate eigenvector x, |A*x-lambda*x|
225: */
226: PetscErrorCode ComputeResidualNorm(Mat A,PetscScalar lambda,Vec x,PetscReal *r)
227: {
228: Vec Ax;
230: PetscFunctionBeginUser;
231: PetscCall(VecDuplicate(x,&Ax));
232: PetscCall(MatMult(A,x,Ax));
233: PetscCall(VecAXPY(Ax,-lambda,x));
234: PetscCall(VecNorm(Ax,NORM_2,r));
235: PetscCall(VecDestroy(&Ax));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*TEST
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
255: TEST*/