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[] = "Tests multiple calls to NEPSolve() with different matrix size.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = number of grid subdivisions.\n"
14 : " -tau <tau>, where <tau> is the delay parameter.\n"
15 : " -split <0/1>, to select the split form in the problem definition (enabled by default).\n";
16 :
17 : /* Based on ex22.c (delay) */
18 :
19 : #include <slepcnep.h>
20 :
21 : /*
22 : User-defined application context
23 : */
24 : typedef struct {
25 : PetscScalar tau;
26 : PetscReal a;
27 : } ApplicationCtx;
28 :
29 : /*
30 : Create problem matrices in split form
31 : */
32 30 : PetscErrorCode BuildSplitMatrices(PetscInt n,PetscReal a,Mat *Id,Mat *A,Mat *B)
33 : {
34 30 : PetscInt i,Istart,Iend;
35 30 : PetscReal h,xi;
36 30 : PetscScalar b;
37 :
38 30 : PetscFunctionBeginUser;
39 30 : h = PETSC_PI/(PetscReal)(n+1);
40 :
41 : /* Identity matrix */
42 30 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,Id));
43 30 : PetscCall(MatSetOption(*Id,MAT_HERMITIAN,PETSC_TRUE));
44 :
45 : /* A = 1/h^2*tridiag(1,-2,1) + a*I */
46 30 : PetscCall(MatCreate(PETSC_COMM_WORLD,A));
47 30 : PetscCall(MatSetSizes(*A,PETSC_DECIDE,PETSC_DECIDE,n,n));
48 30 : PetscCall(MatSetFromOptions(*A));
49 30 : PetscCall(MatGetOwnershipRange(*A,&Istart,&Iend));
50 3486 : for (i=Istart;i<Iend;i++) {
51 3456 : if (i>0) PetscCall(MatSetValue(*A,i,i-1,1.0/(h*h),INSERT_VALUES));
52 3456 : if (i<n-1) PetscCall(MatSetValue(*A,i,i+1,1.0/(h*h),INSERT_VALUES));
53 3456 : PetscCall(MatSetValue(*A,i,i,-2.0/(h*h)+a,INSERT_VALUES));
54 : }
55 30 : PetscCall(MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY));
56 30 : PetscCall(MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY));
57 30 : PetscCall(MatSetOption(*A,MAT_HERMITIAN,PETSC_TRUE));
58 :
59 : /* B = diag(b(xi)) */
60 30 : PetscCall(MatCreate(PETSC_COMM_WORLD,B));
61 30 : PetscCall(MatSetSizes(*B,PETSC_DECIDE,PETSC_DECIDE,n,n));
62 30 : PetscCall(MatSetFromOptions(*B));
63 30 : PetscCall(MatGetOwnershipRange(*B,&Istart,&Iend));
64 3486 : for (i=Istart;i<Iend;i++) {
65 3456 : xi = (i+1)*h;
66 3456 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
67 3456 : PetscCall(MatSetValue(*B,i,i,b,INSERT_VALUES));
68 : }
69 30 : PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY));
70 30 : PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY));
71 30 : PetscCall(MatSetOption(*B,MAT_HERMITIAN,PETSC_TRUE));
72 30 : PetscFunctionReturn(PETSC_SUCCESS);
73 : }
74 :
75 : /*
76 : Compute Function matrix T(lambda)
77 : */
78 576 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
79 : {
80 576 : ApplicationCtx *user = (ApplicationCtx*)ctx;
81 576 : PetscInt i,n,Istart,Iend;
82 576 : PetscReal h,xi;
83 576 : PetscScalar b;
84 :
85 576 : PetscFunctionBeginUser;
86 576 : PetscCall(MatGetSize(fun,&n,NULL));
87 576 : h = PETSC_PI/(PetscReal)(n+1);
88 576 : PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
89 99008 : for (i=Istart;i<Iend;i++) {
90 98432 : if (i>0) PetscCall(MatSetValue(fun,i,i-1,1.0/(h*h),INSERT_VALUES));
91 98432 : if (i<n-1) PetscCall(MatSetValue(fun,i,i+1,1.0/(h*h),INSERT_VALUES));
92 98432 : xi = (i+1)*h;
93 98432 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
94 98432 : PetscCall(MatSetValue(fun,i,i,-lambda-2.0/(h*h)+user->a+PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES));
95 : }
96 576 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
97 576 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
98 576 : if (fun != B) {
99 0 : PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
100 0 : PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
101 : }
102 576 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 :
105 : /*
106 : Compute Jacobian matrix T'(lambda)
107 : */
108 190 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
109 : {
110 190 : ApplicationCtx *user = (ApplicationCtx*)ctx;
111 190 : PetscInt i,n,Istart,Iend;
112 190 : PetscReal h,xi;
113 190 : PetscScalar b;
114 :
115 190 : PetscFunctionBeginUser;
116 190 : PetscCall(MatGetSize(jac,&n,NULL));
117 190 : h = PETSC_PI/(PetscReal)(n+1);
118 190 : PetscCall(MatGetOwnershipRange(jac,&Istart,&Iend));
119 26558 : for (i=Istart;i<Iend;i++) {
120 26368 : xi = (i+1)*h;
121 26368 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
122 26368 : PetscCall(MatSetValue(jac,i,i,-1.0-user->tau*PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES));
123 : }
124 190 : PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
125 190 : PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
126 190 : PetscFunctionReturn(PETSC_SUCCESS);
127 : }
128 :
129 21 : int main(int argc,char **argv)
130 : {
131 21 : NEP nep; /* nonlinear eigensolver context */
132 21 : Mat Id,A,B,J,F; /* problem matrices */
133 21 : FN f1,f2,f3; /* functions to define the nonlinear operator */
134 21 : ApplicationCtx ctx; /* user-defined context */
135 21 : Mat mats[3];
136 21 : FN funs[3];
137 21 : PetscScalar coeffs[2];
138 21 : PetscInt n=128;
139 21 : PetscReal tau=0.001,a=20;
140 21 : PetscBool split=PETSC_TRUE;
141 :
142 21 : PetscFunctionBeginUser;
143 21 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
144 21 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
145 21 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL));
146 21 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
147 21 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau));
148 :
149 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150 : Create nonlinear eigensolver and set options
151 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152 :
153 21 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
154 21 : PetscCall(NEPSetTolerances(nep,1e-9,PETSC_DEFAULT));
155 :
156 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157 : First solve
158 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159 :
160 21 : if (split) {
161 15 : PetscCall(BuildSplitMatrices(n,a,&Id,&A,&B));
162 : /* f1=-lambda */
163 15 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f1));
164 15 : PetscCall(FNSetType(f1,FNRATIONAL));
165 15 : coeffs[0] = -1.0; coeffs[1] = 0.0;
166 15 : PetscCall(FNRationalSetNumerator(f1,2,coeffs));
167 : /* f2=1.0 */
168 15 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f2));
169 15 : PetscCall(FNSetType(f2,FNRATIONAL));
170 15 : coeffs[0] = 1.0;
171 15 : PetscCall(FNRationalSetNumerator(f2,1,coeffs));
172 : /* f3=exp(-tau*lambda) */
173 15 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f3));
174 15 : PetscCall(FNSetType(f3,FNEXP));
175 15 : PetscCall(FNSetScale(f3,-tau,1.0));
176 15 : mats[0] = A; funs[0] = f2;
177 15 : mats[1] = Id; funs[1] = f1;
178 15 : mats[2] = B; funs[2] = f3;
179 15 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
180 : } else {
181 : /* callback form */
182 6 : ctx.tau = tau;
183 6 : ctx.a = a;
184 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
185 6 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
186 6 : PetscCall(MatSetFromOptions(F));
187 6 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
188 6 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
189 6 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
190 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
191 6 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
192 6 : PetscCall(MatSetFromOptions(J));
193 6 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
194 6 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
195 6 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
196 : }
197 :
198 : /* Set solver parameters at runtime */
199 21 : PetscCall(NEPSetFromOptions(nep));
200 :
201 : /* Solve the eigensystem */
202 21 : PetscCall(NEPSolve(nep));
203 21 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
204 :
205 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206 : Second solve, with problem matrices of size 2*n
207 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208 :
209 21 : n *= 2;
210 21 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau));
211 21 : if (split) {
212 15 : PetscCall(MatDestroy(&Id));
213 15 : PetscCall(MatDestroy(&A));
214 15 : PetscCall(MatDestroy(&B));
215 15 : PetscCall(BuildSplitMatrices(n,a,&Id,&A,&B));
216 15 : mats[0] = A;
217 15 : mats[1] = Id;
218 15 : mats[2] = B;
219 15 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
220 : } else {
221 : /* callback form */
222 6 : PetscCall(MatDestroy(&F));
223 6 : PetscCall(MatDestroy(&J));
224 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
225 6 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
226 6 : PetscCall(MatSetFromOptions(F));
227 6 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
228 6 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
229 6 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
230 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
231 6 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
232 6 : PetscCall(MatSetFromOptions(J));
233 6 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
234 6 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
235 6 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
236 : }
237 :
238 : /* Solve the eigensystem */
239 21 : PetscCall(NEPSolve(nep));
240 21 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
241 :
242 21 : PetscCall(NEPDestroy(&nep));
243 21 : if (split) {
244 15 : PetscCall(MatDestroy(&Id));
245 15 : PetscCall(MatDestroy(&A));
246 15 : PetscCall(MatDestroy(&B));
247 15 : PetscCall(FNDestroy(&f1));
248 15 : PetscCall(FNDestroy(&f2));
249 15 : PetscCall(FNDestroy(&f3));
250 : } else {
251 6 : PetscCall(MatDestroy(&F));
252 6 : PetscCall(MatDestroy(&J));
253 : }
254 21 : PetscCall(SlepcFinalize());
255 : return 0;
256 : }
257 :
258 : /*TEST
259 :
260 : testset:
261 : nsize: 2
262 : requires: !single
263 : output_file: output/test10_1.out
264 : test:
265 : suffix: 1
266 : args: -nep_type narnoldi -nep_target 0.55
267 : test:
268 : suffix: 1_rii
269 : args: -nep_type rii -nep_target 0.55 -nep_rii_hermitian -split {{0 1}}
270 : test:
271 : suffix: 1_narnoldi
272 : args: -nep_type narnoldi -nep_target 0.55 -nep_narnoldi_lag_preconditioner 2
273 : test:
274 : suffix: 1_slp
275 : args: -nep_type slp -nep_slp_st_pc_type redundant -split {{0 1}}
276 : test:
277 : suffix: 1_interpol
278 : args: -nep_type interpol -rg_type interval -rg_interval_endpoints .5,1,-.1,.1 -nep_target .7 -nep_interpol_st_pc_type redundant
279 : test:
280 : suffix: 1_narnoldi_sync
281 : args: -nep_type narnoldi -ds_parallel synchronized
282 :
283 : testset:
284 : args: -nep_nev 2 -rg_type interval -rg_interval_endpoints .5,15,-.1,.1 -nep_target .7
285 : requires: !single
286 : output_file: output/test10_2.out
287 : filter: sed -e "s/[+-]0\.0*i//g"
288 : test:
289 : suffix: 2_interpol
290 : args: -nep_type interpol -nep_interpol_pep_type jd -nep_interpol_st_pc_type sor
291 : test:
292 : suffix: 2_nleigs
293 : args: -nep_type nleigs -split {{0 1}}
294 : requires: complex
295 : test:
296 : suffix: 2_nleigs_real
297 : args: -nep_type nleigs -rg_interval_endpoints .5,15 -split {{0 1}}
298 : requires: !complex
299 :
300 : test:
301 : suffix: 3
302 : requires: complex !single
303 : args: -nep_type ciss -rg_type ellipse -rg_ellipse_center 10 -rg_ellipse_radius 9.5 -rg_ellipse_vscale 0.1 -split {{0 1}}
304 : timeoutfactor: 2
305 :
306 : TEST*/
|
