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[] = "Test the RII solver with a user-provided KSP.\n\n"
12 : "This is a simplified version of ex20.\n"
13 : "The command line options are:\n"
14 : " -n <n>, where <n> = number of grid subdivisions.\n";
15 :
16 : /*
17 : Solve 1-D PDE
18 : -u'' = lambda*u
19 : on [0,1] subject to
20 : u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
21 : */
22 :
23 : #include <slepcnep.h>
24 :
25 : /*
26 : User-defined routines
27 : */
28 : PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
29 : PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
30 :
31 : /*
32 : User-defined application context
33 : */
34 : typedef struct {
35 : PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
36 : PetscReal h; /* mesh spacing */
37 : } ApplicationCtx;
38 :
39 1 : int main(int argc,char **argv)
40 : {
41 1 : NEP nep;
42 1 : KSP ksp;
43 1 : PC pc;
44 1 : Mat F,J;
45 1 : ApplicationCtx ctx;
46 1 : PetscInt n=128,lag,its;
47 1 : PetscBool terse,flg,cct,herm;
48 1 : PetscReal thres;
49 :
50 1 : PetscFunctionBeginUser;
51 1 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
52 1 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
53 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n));
54 1 : ctx.h = 1.0/(PetscReal)n;
55 1 : ctx.kappa = 1.0;
56 :
57 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
58 : Create a standalone KSP with appropriate settings
59 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
60 :
61 1 : PetscCall(KSPCreate(PETSC_COMM_WORLD,&ksp));
62 1 : PetscCall(KSPSetType(ksp,KSPBCGS));
63 1 : PetscCall(KSPGetPC(ksp,&pc));
64 1 : PetscCall(PCSetType(pc,PCBJACOBI));
65 1 : PetscCall(KSPSetFromOptions(ksp));
66 :
67 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
68 : Prepare nonlinear eigensolver context
69 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
70 :
71 1 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
72 :
73 : /* Create Function and Jacobian matrices; set evaluation routines */
74 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
75 1 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
76 1 : PetscCall(MatSetFromOptions(F));
77 1 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
78 1 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
79 1 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
80 :
81 1 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
82 1 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
83 1 : PetscCall(MatSetFromOptions(J));
84 1 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
85 1 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
86 1 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
87 :
88 1 : PetscCall(NEPSetType(nep,NEPRII));
89 1 : PetscCall(NEPRIISetKSP(nep,ksp));
90 1 : PetscCall(NEPRIISetMaximumIterations(nep,6));
91 1 : PetscCall(NEPRIISetConstCorrectionTol(nep,PETSC_TRUE));
92 1 : PetscCall(NEPRIISetHermitian(nep,PETSC_TRUE));
93 1 : PetscCall(NEPSetFromOptions(nep));
94 :
95 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
96 : Solve the eigensystem
97 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
98 :
99 1 : PetscCall(NEPSolve(nep));
100 1 : PetscCall(PetscObjectTypeCompare((PetscObject)nep,NEPRII,&flg));
101 1 : if (flg) {
102 1 : PetscCall(NEPRIIGetMaximumIterations(nep,&its));
103 1 : PetscCall(NEPRIIGetLagPreconditioner(nep,&lag));
104 1 : PetscCall(NEPRIIGetDeflationThreshold(nep,&thres));
105 1 : PetscCall(NEPRIIGetConstCorrectionTol(nep,&cct));
106 1 : PetscCall(NEPRIIGetHermitian(nep,&herm));
107 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Maximum inner iterations of RII is %" PetscInt_FMT "\n",its));
108 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Preconditioner rebuilt every %" PetscInt_FMT " iterations\n",lag));
109 1 : if (thres>0.0) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using deflation threshold=%g\n",(double)thres));
110 1 : if (cct) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Using a constant correction tolerance\n"));
111 1 : if (herm) PetscCall(PetscPrintf(PETSC_COMM_WORLD," Hermitian version of scalar equation\n"));
112 1 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n"));
113 : }
114 :
115 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116 : Display solution and clean up
117 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
118 :
119 : /* show detailed info unless -terse option is given by user */
120 1 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
121 1 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
122 : else {
123 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
124 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
125 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
126 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
127 : }
128 :
129 1 : PetscCall(NEPDestroy(&nep));
130 1 : PetscCall(KSPDestroy(&ksp));
131 1 : PetscCall(MatDestroy(&F));
132 1 : PetscCall(MatDestroy(&J));
133 1 : PetscCall(SlepcFinalize());
134 : return 0;
135 : }
136 :
137 : /* ------------------------------------------------------------------- */
138 : /*
139 : FormFunction - Computes Function matrix T(lambda)
140 :
141 : Input Parameters:
142 : . nep - the NEP context
143 : . lambda - the scalar argument
144 : . ctx - optional user-defined context, as set by NEPSetFunction()
145 :
146 : Output Parameters:
147 : . fun - Function matrix
148 : . B - optionally different preconditioning matrix
149 : */
150 22 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
151 : {
152 22 : ApplicationCtx *user = (ApplicationCtx*)ctx;
153 22 : PetscScalar A[3],c,d;
154 22 : PetscReal h;
155 22 : PetscInt i,n,j[3],Istart,Iend;
156 22 : PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
157 :
158 22 : PetscFunctionBeginUser;
159 : /*
160 : Compute Function entries and insert into matrix
161 : */
162 22 : PetscCall(MatGetSize(fun,&n,NULL));
163 22 : PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
164 22 : if (Istart==0) FirstBlock=PETSC_TRUE;
165 22 : if (Iend==n) LastBlock=PETSC_TRUE;
166 22 : h = user->h;
167 22 : c = user->kappa/(lambda-user->kappa);
168 22 : d = n;
169 :
170 : /*
171 : Interior grid points
172 : */
173 2794 : for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
174 2772 : j[0] = i-1; j[1] = i; j[2] = i+1;
175 2772 : A[0] = A[2] = -d-lambda*h/6.0; A[1] = 2.0*(d-lambda*h/3.0);
176 2772 : PetscCall(MatSetValues(fun,1,&i,3,j,A,INSERT_VALUES));
177 : }
178 :
179 : /*
180 : Boundary points
181 : */
182 22 : if (FirstBlock) {
183 22 : i = 0;
184 22 : j[0] = 0; j[1] = 1;
185 22 : A[0] = 2.0*(d-lambda*h/3.0); A[1] = -d-lambda*h/6.0;
186 22 : PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
187 : }
188 :
189 22 : if (LastBlock) {
190 22 : i = n-1;
191 22 : j[0] = n-2; j[1] = n-1;
192 22 : A[0] = -d-lambda*h/6.0; A[1] = d-lambda*h/3.0+lambda*c;
193 22 : PetscCall(MatSetValues(fun,1,&i,2,j,A,INSERT_VALUES));
194 : }
195 :
196 : /*
197 : Assemble matrix
198 : */
199 22 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
200 22 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
201 22 : if (fun != B) {
202 0 : PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
203 0 : PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
204 : }
205 22 : PetscFunctionReturn(PETSC_SUCCESS);
206 : }
207 :
208 : /* ------------------------------------------------------------------- */
209 : /*
210 : FormJacobian - Computes Jacobian matrix T'(lambda)
211 :
212 : Input Parameters:
213 : . nep - the NEP context
214 : . lambda - the scalar argument
215 : . ctx - optional user-defined context, as set by NEPSetJacobian()
216 :
217 : Output Parameters:
218 : . jac - Jacobian matrix
219 : . B - optionally different preconditioning matrix
220 : */
221 19 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
222 : {
223 19 : ApplicationCtx *user = (ApplicationCtx*)ctx;
224 19 : PetscScalar A[3],c;
225 19 : PetscReal h;
226 19 : PetscInt i,n,j[3],Istart,Iend;
227 19 : PetscBool FirstBlock=PETSC_FALSE,LastBlock=PETSC_FALSE;
228 :
229 19 : PetscFunctionBeginUser;
230 : /*
231 : Compute Jacobian entries and insert into matrix
232 : */
233 19 : PetscCall(MatGetSize(jac,&n,NULL));
234 19 : PetscCall(MatGetOwnershipRange(jac,&Istart,&Iend));
235 19 : if (Istart==0) FirstBlock=PETSC_TRUE;
236 19 : if (Iend==n) LastBlock=PETSC_TRUE;
237 19 : h = user->h;
238 19 : c = user->kappa/(lambda-user->kappa);
239 :
240 : /*
241 : Interior grid points
242 : */
243 2413 : for (i=(FirstBlock? Istart+1: Istart);i<(LastBlock? Iend-1: Iend);i++) {
244 2394 : j[0] = i-1; j[1] = i; j[2] = i+1;
245 2394 : A[0] = A[2] = -h/6.0; A[1] = -2.0*h/3.0;
246 2394 : PetscCall(MatSetValues(jac,1,&i,3,j,A,INSERT_VALUES));
247 : }
248 :
249 : /*
250 : Boundary points
251 : */
252 19 : if (FirstBlock) {
253 19 : i = 0;
254 19 : j[0] = 0; j[1] = 1;
255 19 : A[0] = -2.0*h/3.0; A[1] = -h/6.0;
256 19 : PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
257 : }
258 :
259 19 : if (LastBlock) {
260 19 : i = n-1;
261 19 : j[0] = n-2; j[1] = n-1;
262 19 : A[0] = -h/6.0; A[1] = -h/3.0-c*c;
263 19 : PetscCall(MatSetValues(jac,1,&i,2,j,A,INSERT_VALUES));
264 : }
265 :
266 : /*
267 : Assemble matrix
268 : */
269 19 : PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
270 19 : PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
271 19 : PetscFunctionReturn(PETSC_SUCCESS);
272 : }
273 :
274 : /*TEST
275 :
276 : test:
277 : suffix: 1
278 : args: -nep_target 21 -nep_rii_lag_preconditioner 2 -terse
279 : requires: !single
280 :
281 : TEST*/
|
