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[] = "Simple 1-D nonlinear eigenproblem (matrix-free version).\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions\n\n";
14 :
15 : /*
16 : Solve 1-D PDE
17 : -u'' = lambda*u
18 : on [0,1] subject to
19 : u(0)=0, u'(1)=u(1)*lambda*kappa/(kappa-lambda)
20 : */
21 :
22 : #include <slepcnep.h>
23 :
24 : /*
25 : User-defined routines
26 : */
27 : PetscErrorCode FormFunction(NEP,PetscScalar,Mat,Mat,void*);
28 : PetscErrorCode FormJacobian(NEP,PetscScalar,Mat,void*);
29 :
30 : /*
31 : Matrix operations and context
32 : */
33 : PetscErrorCode MatMult_Fun(Mat,Vec,Vec);
34 : PetscErrorCode MatGetDiagonal_Fun(Mat,Vec);
35 : PetscErrorCode MatDestroy_Fun(Mat);
36 : PetscErrorCode MatDuplicate_Fun(Mat,MatDuplicateOption,Mat*);
37 : PetscErrorCode MatMult_Jac(Mat,Vec,Vec);
38 : PetscErrorCode MatGetDiagonal_Jac(Mat,Vec);
39 : PetscErrorCode MatDestroy_Jac(Mat);
40 :
41 : typedef struct {
42 : PetscScalar lambda,kappa;
43 : PetscReal h;
44 : PetscMPIInt next,prev;
45 : } MatCtx;
46 :
47 : /*
48 : User-defined application context
49 : */
50 : typedef struct {
51 : PetscScalar kappa; /* ratio between stiffness of spring and attached mass */
52 : PetscReal h; /* mesh spacing */
53 : } ApplicationCtx;
54 :
55 6 : int main(int argc,char **argv)
56 : {
57 6 : NEP nep; /* nonlinear eigensolver context */
58 6 : Mat F,J; /* Function and Jacobian matrices */
59 6 : ApplicationCtx ctx; /* user-defined context */
60 6 : MatCtx *ctxF,*ctxJ; /* contexts for shell matrices */
61 6 : PetscInt n=128,nev;
62 6 : KSP ksp;
63 6 : PC pc;
64 6 : PetscMPIInt rank,size;
65 6 : PetscBool terse;
66 :
67 6 : PetscFunctionBeginUser;
68 6 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
69 6 : PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size));
70 6 : PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank));
71 6 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
72 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Nonlinear Eigenproblem, n=%" PetscInt_FMT "\n\n",n));
73 6 : ctx.h = 1.0/(PetscReal)n;
74 6 : ctx.kappa = 1.0;
75 :
76 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
77 : Create nonlinear eigensolver context
78 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79 :
80 6 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
81 :
82 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83 : Create matrix data structure; set Function evaluation routine
84 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
85 :
86 6 : PetscCall(PetscNew(&ctxF));
87 6 : ctxF->h = ctx.h;
88 6 : ctxF->kappa = ctx.kappa;
89 6 : ctxF->next = rank==size-1? MPI_PROC_NULL: rank+1;
90 6 : ctxF->prev = rank==0? MPI_PROC_NULL: rank-1;
91 :
92 6 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctxF,&F));
93 6 : PetscCall(MatShellSetOperation(F,MATOP_MULT,(void(*)(void))MatMult_Fun));
94 6 : PetscCall(MatShellSetOperation(F,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun));
95 6 : PetscCall(MatShellSetOperation(F,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun));
96 6 : PetscCall(MatShellSetOperation(F,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun));
97 :
98 : /*
99 : Set Function matrix data structure and default Function evaluation
100 : routine
101 : */
102 6 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,NULL));
103 :
104 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105 : Create matrix data structure; set Jacobian evaluation routine
106 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107 :
108 6 : PetscCall(PetscNew(&ctxJ));
109 6 : ctxJ->h = ctx.h;
110 6 : ctxJ->kappa = ctx.kappa;
111 6 : ctxJ->next = rank==size-1? MPI_PROC_NULL: rank+1;
112 6 : ctxJ->prev = rank==0? MPI_PROC_NULL: rank-1;
113 :
114 6 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,(void*)ctxJ,&J));
115 6 : PetscCall(MatShellSetOperation(J,MATOP_MULT,(void(*)(void))MatMult_Jac));
116 6 : PetscCall(MatShellSetOperation(J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Jac));
117 6 : PetscCall(MatShellSetOperation(J,MATOP_DESTROY,(void(*)(void))MatDestroy_Jac));
118 :
119 : /*
120 : Set Jacobian matrix data structure and default Jacobian evaluation
121 : routine
122 : */
123 6 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,NULL));
124 :
125 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126 : Customize nonlinear solver; set runtime options
127 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128 :
129 6 : PetscCall(NEPSetType(nep,NEPRII));
130 6 : PetscCall(NEPRIISetLagPreconditioner(nep,0));
131 6 : PetscCall(NEPRIIGetKSP(nep,&ksp));
132 6 : PetscCall(KSPSetType(ksp,KSPBCGS));
133 6 : PetscCall(KSPGetPC(ksp,&pc));
134 6 : PetscCall(PCSetType(pc,PCJACOBI));
135 :
136 : /*
137 : Set solver parameters at runtime
138 : */
139 6 : PetscCall(NEPSetFromOptions(nep));
140 :
141 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
142 : Solve the eigensystem
143 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144 :
145 6 : PetscCall(NEPSolve(nep));
146 6 : PetscCall(NEPGetDimensions(nep,&nev,NULL,NULL));
147 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
148 :
149 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150 : Display solution and clean up
151 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152 :
153 : /* show detailed info unless -terse option is given by user */
154 6 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
155 6 : if (terse) PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
156 : else {
157 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
158 0 : PetscCall(NEPConvergedReasonView(nep,PETSC_VIEWER_STDOUT_WORLD));
159 0 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
160 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
161 : }
162 6 : PetscCall(NEPDestroy(&nep));
163 6 : PetscCall(MatDestroy(&F));
164 6 : PetscCall(MatDestroy(&J));
165 6 : PetscCall(SlepcFinalize());
166 : return 0;
167 : }
168 :
169 : /* ------------------------------------------------------------------- */
170 : /*
171 : FormFunction - Computes Function matrix T(lambda)
172 :
173 : Input Parameters:
174 : . nep - the NEP context
175 : . lambda - real part of the scalar argument
176 : . ctx - optional user-defined context, as set by NEPSetFunction()
177 :
178 : Output Parameters:
179 : . fun - Function matrix
180 : . B - optionally different preconditioning matrix
181 : */
182 90 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
183 : {
184 90 : MatCtx *ctxF;
185 :
186 90 : PetscFunctionBeginUser;
187 90 : PetscCall(MatShellGetContext(fun,&ctxF));
188 90 : ctxF->lambda = lambda;
189 90 : PetscFunctionReturn(PETSC_SUCCESS);
190 : }
191 :
192 : /* ------------------------------------------------------------------- */
193 : /*
194 : FormJacobian - Computes Jacobian matrix T'(lambda)
195 :
196 : Input Parameters:
197 : . nep - the NEP context
198 : . lambda - real part of the scalar argument
199 : . ctx - optional user-defined context, as set by NEPSetJacobian()
200 :
201 : Output Parameters:
202 : . jac - Jacobian matrix
203 : */
204 75 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
205 : {
206 75 : MatCtx *ctxJ;
207 :
208 75 : PetscFunctionBeginUser;
209 75 : PetscCall(MatShellGetContext(jac,&ctxJ));
210 75 : ctxJ->lambda = lambda;
211 75 : PetscFunctionReturn(PETSC_SUCCESS);
212 : }
213 :
214 : /* ------------------------------------------------------------------- */
215 74547 : PetscErrorCode MatMult_Fun(Mat A,Vec x,Vec y)
216 : {
217 74547 : MatCtx *ctx;
218 74547 : PetscInt i,n,N;
219 74547 : const PetscScalar *px;
220 74547 : PetscScalar *py,c,d,de,oe,upper=0.0,lower=0.0;
221 74547 : PetscReal h;
222 74547 : MPI_Comm comm;
223 :
224 74547 : PetscFunctionBeginUser;
225 74547 : PetscCall(MatShellGetContext(A,&ctx));
226 74547 : PetscCall(VecGetArrayRead(x,&px));
227 74547 : PetscCall(VecGetArray(y,&py));
228 74547 : PetscCall(VecGetSize(x,&N));
229 74547 : PetscCall(VecGetLocalSize(x,&n));
230 :
231 74547 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
232 74547 : PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,ctx->prev,0,&lower,1,MPIU_SCALAR,ctx->next,0,comm,MPI_STATUS_IGNORE));
233 74547 : PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,ctx->next,0,&upper,1,MPIU_SCALAR,ctx->prev,0,comm,MPI_STATUS_IGNORE));
234 :
235 74547 : h = ctx->h;
236 74547 : c = ctx->kappa/(ctx->lambda-ctx->kappa);
237 74547 : d = N;
238 74547 : de = 2.0*(d-ctx->lambda*h/3.0); /* diagonal entry */
239 74547 : oe = -d-ctx->lambda*h/6.0; /* offdiagonal entry */
240 74547 : py[0] = oe*upper + de*px[0] + oe*px[1];
241 6426189 : for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
242 74547 : if (ctx->next==MPI_PROC_NULL) de = d-ctx->lambda*h/3.0+ctx->lambda*c; /* diagonal entry of last row */
243 74547 : py[n-1] = oe*px[n-2] + de*px[n-1] + oe*lower;
244 :
245 74547 : PetscCall(VecRestoreArrayRead(x,&px));
246 74547 : PetscCall(VecRestoreArray(y,&py));
247 74547 : PetscFunctionReturn(PETSC_SUCCESS);
248 : }
249 :
250 : /* ------------------------------------------------------------------- */
251 6 : PetscErrorCode MatGetDiagonal_Fun(Mat A,Vec diag)
252 : {
253 6 : MatCtx *ctx;
254 6 : PetscInt n,N;
255 6 : PetscScalar *pd,c,d;
256 6 : PetscReal h;
257 :
258 6 : PetscFunctionBeginUser;
259 6 : PetscCall(MatShellGetContext(A,&ctx));
260 6 : PetscCall(VecGetSize(diag,&N));
261 6 : PetscCall(VecGetLocalSize(diag,&n));
262 6 : h = ctx->h;
263 6 : c = ctx->kappa/(ctx->lambda-ctx->kappa);
264 6 : d = N;
265 6 : PetscCall(VecSet(diag,2.0*(d-ctx->lambda*h/3.0)));
266 6 : PetscCall(VecGetArray(diag,&pd));
267 6 : pd[n-1] = d-ctx->lambda*h/3.0+ctx->lambda*c;
268 6 : PetscCall(VecRestoreArray(diag,&pd));
269 6 : PetscFunctionReturn(PETSC_SUCCESS);
270 : }
271 :
272 : /* ------------------------------------------------------------------- */
273 9 : PetscErrorCode MatDestroy_Fun(Mat A)
274 : {
275 9 : MatCtx *ctx;
276 :
277 9 : PetscFunctionBegin;
278 9 : PetscCall(MatShellGetContext(A,&ctx));
279 9 : PetscCall(PetscFree(ctx));
280 9 : PetscFunctionReturn(PETSC_SUCCESS);
281 : }
282 :
283 : /* ------------------------------------------------------------------- */
284 3 : PetscErrorCode MatDuplicate_Fun(Mat A,MatDuplicateOption op,Mat *B)
285 : {
286 3 : MatCtx *actx,*bctx;
287 3 : PetscInt m,n,M,N;
288 3 : MPI_Comm comm;
289 :
290 3 : PetscFunctionBegin;
291 3 : PetscCall(MatShellGetContext(A,&actx));
292 3 : PetscCall(MatGetSize(A,&M,&N));
293 3 : PetscCall(MatGetLocalSize(A,&m,&n));
294 :
295 3 : PetscCall(PetscNew(&bctx));
296 3 : bctx->h = actx->h;
297 3 : bctx->kappa = actx->kappa;
298 3 : bctx->lambda = actx->lambda;
299 3 : bctx->next = actx->next;
300 3 : bctx->prev = actx->prev;
301 :
302 3 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
303 3 : PetscCall(MatCreateShell(comm,m,n,M,N,(void*)bctx,B));
304 3 : PetscCall(MatShellSetOperation(*B,MATOP_MULT,(void(*)(void))MatMult_Fun));
305 3 : PetscCall(MatShellSetOperation(*B,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Fun));
306 3 : PetscCall(MatShellSetOperation(*B,MATOP_DESTROY,(void(*)(void))MatDestroy_Fun));
307 3 : PetscCall(MatShellSetOperation(*B,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_Fun));
308 3 : PetscFunctionReturn(PETSC_SUCCESS);
309 : }
310 :
311 : /* ------------------------------------------------------------------- */
312 210 : PetscErrorCode MatMult_Jac(Mat A,Vec x,Vec y)
313 : {
314 210 : MatCtx *ctx;
315 210 : PetscInt i,n;
316 210 : const PetscScalar *px;
317 210 : PetscScalar *py,c,de,oe,upper=0.0,lower=0.0;
318 210 : PetscReal h;
319 210 : MPI_Comm comm;
320 :
321 210 : PetscFunctionBeginUser;
322 210 : PetscCall(MatShellGetContext(A,&ctx));
323 210 : PetscCall(VecGetArrayRead(x,&px));
324 210 : PetscCall(VecGetArray(y,&py));
325 210 : PetscCall(VecGetLocalSize(x,&n));
326 :
327 210 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
328 210 : PetscCallMPI(MPI_Sendrecv(px,1,MPIU_SCALAR,ctx->prev,0,&lower,1,MPIU_SCALAR,ctx->next,0,comm,MPI_STATUS_IGNORE));
329 210 : PetscCallMPI(MPI_Sendrecv(px+n-1,1,MPIU_SCALAR,ctx->next,0,&upper,1,MPIU_SCALAR,ctx->prev,0,comm,MPI_STATUS_IGNORE));
330 :
331 210 : h = ctx->h;
332 210 : c = ctx->kappa/(ctx->lambda-ctx->kappa);
333 210 : de = -2.0*h/3.0; /* diagonal entry */
334 210 : oe = -h/6.0; /* offdiagonal entry */
335 210 : py[0] = oe*upper + de*px[0] + oe*px[1];
336 18094 : for (i=1;i<n-1;i++) py[i] = oe*px[i-1] +de*px[i] + oe*px[i+1];
337 210 : if (ctx->next==MPI_PROC_NULL) de = -h/3.0-c*c; /* diagonal entry of last row */
338 210 : py[n-1] = oe*px[n-2] + de*px[n-1] + oe*lower;
339 :
340 210 : PetscCall(VecRestoreArrayRead(x,&px));
341 210 : PetscCall(VecRestoreArray(y,&py));
342 210 : PetscFunctionReturn(PETSC_SUCCESS);
343 : }
344 :
345 : /* ------------------------------------------------------------------- */
346 0 : PetscErrorCode MatGetDiagonal_Jac(Mat A,Vec diag)
347 : {
348 0 : MatCtx *ctx;
349 0 : PetscInt n;
350 0 : PetscScalar *pd,c;
351 0 : PetscReal h;
352 :
353 0 : PetscFunctionBeginUser;
354 0 : PetscCall(MatShellGetContext(A,&ctx));
355 0 : PetscCall(VecGetLocalSize(diag,&n));
356 0 : h = ctx->h;
357 0 : c = ctx->kappa/(ctx->lambda-ctx->kappa);
358 0 : PetscCall(VecSet(diag,-2.0*h/3.0));
359 0 : PetscCall(VecGetArray(diag,&pd));
360 0 : pd[n-1] = -h/3.0-c*c;
361 0 : PetscCall(VecRestoreArray(diag,&pd));
362 0 : PetscFunctionReturn(PETSC_SUCCESS);
363 : }
364 :
365 : /* ------------------------------------------------------------------- */
366 6 : PetscErrorCode MatDestroy_Jac(Mat A)
367 : {
368 6 : MatCtx *ctx;
369 :
370 6 : PetscFunctionBegin;
371 6 : PetscCall(MatShellGetContext(A,&ctx));
372 6 : PetscCall(PetscFree(ctx));
373 6 : PetscFunctionReturn(PETSC_SUCCESS);
374 : }
375 :
376 : /*TEST
377 :
378 : testset:
379 : nsize: {{1 2}}
380 : args: -terse
381 : requires: !single
382 : output_file: output/ex21_1.out
383 : filter: sed -e "s/[+-]0\.0*i//g" -e "s/+0i//g"
384 : test:
385 : suffix: 1_rii
386 : args: -nep_type rii -nep_target 4
387 : test:
388 : suffix: 1_slp
389 : args: -nep_type slp -nep_slp_pc_type jacobi -nep_slp_ksp_type bcgs -nep_target 10
390 :
391 : TEST*/
|
