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 28 : PetscErrorCode BuildSplitMatrices(PetscInt n,PetscReal a,Mat *Id,Mat *A,Mat *B)
33 : {
34 28 : PetscInt i,Istart,Iend;
35 28 : PetscReal h,xi;
36 28 : PetscScalar b;
37 :
38 28 : PetscFunctionBeginUser;
39 28 : h = PETSC_PI/(PetscReal)(n+1);
40 :
41 : /* Identity matrix */
42 28 : PetscCall(MatCreateConstantDiagonal(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,n,n,1.0,Id));
43 28 : PetscCall(MatSetOption(*Id,MAT_HERMITIAN,PETSC_TRUE));
44 :
45 : /* A = 1/h^2*tridiag(1,-2,1) + a*I */
46 28 : PetscCall(MatCreate(PETSC_COMM_WORLD,A));
47 28 : PetscCall(MatSetSizes(*A,PETSC_DECIDE,PETSC_DECIDE,n,n));
48 28 : PetscCall(MatSetFromOptions(*A));
49 28 : PetscCall(MatGetOwnershipRange(*A,&Istart,&Iend));
50 3100 : for (i=Istart;i<Iend;i++) {
51 3072 : if (i>0) PetscCall(MatSetValue(*A,i,i-1,1.0/(h*h),INSERT_VALUES));
52 3072 : if (i<n-1) PetscCall(MatSetValue(*A,i,i+1,1.0/(h*h),INSERT_VALUES));
53 3072 : PetscCall(MatSetValue(*A,i,i,-2.0/(h*h)+a,INSERT_VALUES));
54 : }
55 28 : PetscCall(MatAssemblyBegin(*A,MAT_FINAL_ASSEMBLY));
56 28 : PetscCall(MatAssemblyEnd(*A,MAT_FINAL_ASSEMBLY));
57 28 : PetscCall(MatSetOption(*A,MAT_HERMITIAN,PETSC_TRUE));
58 :
59 : /* B = diag(b(xi)) */
60 28 : PetscCall(MatCreate(PETSC_COMM_WORLD,B));
61 28 : PetscCall(MatSetSizes(*B,PETSC_DECIDE,PETSC_DECIDE,n,n));
62 28 : PetscCall(MatSetFromOptions(*B));
63 28 : PetscCall(MatGetOwnershipRange(*B,&Istart,&Iend));
64 3100 : for (i=Istart;i<Iend;i++) {
65 3072 : xi = (i+1)*h;
66 3072 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
67 3072 : PetscCall(MatSetValue(*B,i,i,b,INSERT_VALUES));
68 : }
69 28 : PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY));
70 28 : PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY));
71 28 : PetscCall(MatSetOption(*B,MAT_HERMITIAN,PETSC_TRUE));
72 28 : PetscFunctionReturn(PETSC_SUCCESS);
73 : }
74 :
75 : /*
76 : Compute Function matrix T(lambda)
77 : */
78 362 : PetscErrorCode FormFunction(NEP nep,PetscScalar lambda,Mat fun,Mat B,void *ctx)
79 : {
80 362 : ApplicationCtx *user = (ApplicationCtx*)ctx;
81 362 : PetscInt i,n,Istart,Iend;
82 362 : PetscReal h,xi;
83 362 : PetscScalar b;
84 :
85 362 : PetscFunctionBeginUser;
86 362 : PetscCall(MatGetSize(fun,&n,NULL));
87 362 : h = PETSC_PI/(PetscReal)(n+1);
88 362 : PetscCall(MatGetOwnershipRange(fun,&Istart,&Iend));
89 55914 : for (i=Istart;i<Iend;i++) {
90 55552 : if (i>0) PetscCall(MatSetValue(fun,i,i-1,1.0/(h*h),INSERT_VALUES));
91 55552 : if (i<n-1) PetscCall(MatSetValue(fun,i,i+1,1.0/(h*h),INSERT_VALUES));
92 55552 : xi = (i+1)*h;
93 55552 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
94 55552 : PetscCall(MatSetValue(fun,i,i,-lambda-2.0/(h*h)+user->a+PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES));
95 : }
96 362 : PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY));
97 362 : PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY));
98 362 : if (fun != B) {
99 0 : PetscCall(MatAssemblyBegin(fun,MAT_FINAL_ASSEMBLY));
100 0 : PetscCall(MatAssemblyEnd(fun,MAT_FINAL_ASSEMBLY));
101 : }
102 362 : PetscFunctionReturn(PETSC_SUCCESS);
103 : }
104 :
105 : /*
106 : Compute Jacobian matrix T'(lambda)
107 : */
108 116 : PetscErrorCode FormJacobian(NEP nep,PetscScalar lambda,Mat jac,void *ctx)
109 : {
110 116 : ApplicationCtx *user = (ApplicationCtx*)ctx;
111 116 : PetscInt i,n,Istart,Iend;
112 116 : PetscReal h,xi;
113 116 : PetscScalar b;
114 :
115 116 : PetscFunctionBeginUser;
116 116 : PetscCall(MatGetSize(jac,&n,NULL));
117 116 : h = PETSC_PI/(PetscReal)(n+1);
118 116 : PetscCall(MatGetOwnershipRange(jac,&Istart,&Iend));
119 11764 : for (i=Istart;i<Iend;i++) {
120 11648 : xi = (i+1)*h;
121 11648 : b = -4.1+xi*(1.0-PetscExpReal(xi-PETSC_PI));
122 11648 : PetscCall(MatSetValue(jac,i,i,-1.0-user->tau*PetscExpScalar(-user->tau*lambda)*b,INSERT_VALUES));
123 : }
124 116 : PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY));
125 116 : PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY));
126 116 : PetscFunctionReturn(PETSC_SUCCESS);
127 : }
128 :
129 19 : int main(int argc,char **argv)
130 : {
131 19 : NEP nep; /* nonlinear eigensolver context */
132 19 : Mat Id,A,B,J,F; /* problem matrices */
133 19 : FN f1,f2,f3; /* functions to define the nonlinear operator */
134 19 : ApplicationCtx ctx; /* user-defined context */
135 19 : Mat mats[3];
136 19 : FN funs[3];
137 19 : PetscScalar coeffs[2];
138 19 : PetscInt n=128;
139 19 : PetscReal tau=0.001,a=20;
140 19 : PetscBool split=PETSC_TRUE;
141 :
142 19 : PetscFunctionBeginUser;
143 19 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
144 19 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
145 19 : PetscCall(PetscOptionsGetReal(NULL,NULL,"-tau",&tau,NULL));
146 19 : PetscCall(PetscOptionsGetBool(NULL,NULL,"-split",&split,NULL));
147 19 : 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 19 : PetscCall(NEPCreate(PETSC_COMM_WORLD,&nep));
154 19 : PetscCall(NEPSetTolerances(nep,1e-9,PETSC_CURRENT));
155 :
156 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
157 : First solve
158 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
159 :
160 19 : if (split) {
161 14 : PetscCall(BuildSplitMatrices(n,a,&Id,&A,&B));
162 : /* f1=-lambda */
163 14 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f1));
164 14 : PetscCall(FNSetType(f1,FNRATIONAL));
165 14 : coeffs[0] = -1.0; coeffs[1] = 0.0;
166 14 : PetscCall(FNRationalSetNumerator(f1,2,coeffs));
167 : /* f2=1.0 */
168 14 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f2));
169 14 : PetscCall(FNSetType(f2,FNRATIONAL));
170 14 : coeffs[0] = 1.0;
171 14 : PetscCall(FNRationalSetNumerator(f2,1,coeffs));
172 : /* f3=exp(-tau*lambda) */
173 14 : PetscCall(FNCreate(PETSC_COMM_WORLD,&f3));
174 14 : PetscCall(FNSetType(f3,FNEXP));
175 14 : PetscCall(FNSetScale(f3,-tau,1.0));
176 14 : mats[0] = A; funs[0] = f2;
177 14 : mats[1] = Id; funs[1] = f1;
178 14 : mats[2] = B; funs[2] = f3;
179 14 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
180 : } else {
181 : /* callback form */
182 5 : ctx.tau = tau;
183 5 : ctx.a = a;
184 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
185 5 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
186 5 : PetscCall(MatSetFromOptions(F));
187 5 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
188 5 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
189 5 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
190 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
191 5 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
192 5 : PetscCall(MatSetFromOptions(J));
193 5 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
194 5 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
195 5 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
196 : }
197 :
198 : /* Set solver parameters at runtime */
199 19 : PetscCall(NEPSetFromOptions(nep));
200 :
201 : /* Solve the eigensystem */
202 19 : PetscCall(NEPSolve(nep));
203 19 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
204 :
205 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206 : Second solve, with problem matrices of size 2*n
207 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208 :
209 19 : n *= 2;
210 19 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\n1-D Delay Eigenproblem, n=%" PetscInt_FMT ", tau=%g\n\n",n,(double)tau));
211 19 : if (split) {
212 14 : PetscCall(MatDestroy(&Id));
213 14 : PetscCall(MatDestroy(&A));
214 14 : PetscCall(MatDestroy(&B));
215 14 : PetscCall(BuildSplitMatrices(n,a,&Id,&A,&B));
216 14 : mats[0] = A;
217 14 : mats[1] = Id;
218 14 : mats[2] = B;
219 14 : PetscCall(NEPSetSplitOperator(nep,3,mats,funs,SUBSET_NONZERO_PATTERN));
220 : } else {
221 : /* callback form */
222 5 : PetscCall(MatDestroy(&F));
223 5 : PetscCall(MatDestroy(&J));
224 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&F));
225 5 : PetscCall(MatSetSizes(F,PETSC_DECIDE,PETSC_DECIDE,n,n));
226 5 : PetscCall(MatSetFromOptions(F));
227 5 : PetscCall(MatSeqAIJSetPreallocation(F,3,NULL));
228 5 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
229 5 : PetscCall(NEPSetFunction(nep,F,F,FormFunction,&ctx));
230 5 : PetscCall(MatCreate(PETSC_COMM_WORLD,&J));
231 5 : PetscCall(MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n));
232 5 : PetscCall(MatSetFromOptions(J));
233 5 : PetscCall(MatSeqAIJSetPreallocation(J,3,NULL));
234 5 : PetscCall(MatMPIAIJSetPreallocation(F,3,NULL,1,NULL));
235 5 : PetscCall(NEPSetJacobian(nep,J,FormJacobian,&ctx));
236 : }
237 :
238 : /* Solve the eigensystem */
239 19 : PetscCall(NEPSolve(nep));
240 19 : PetscCall(NEPErrorView(nep,NEP_ERROR_RELATIVE,NULL));
241 :
242 19 : PetscCall(NEPDestroy(&nep));
243 19 : if (split) {
244 14 : PetscCall(MatDestroy(&Id));
245 14 : PetscCall(MatDestroy(&A));
246 14 : PetscCall(MatDestroy(&B));
247 14 : PetscCall(FNDestroy(&f1));
248 14 : PetscCall(FNDestroy(&f2));
249 14 : PetscCall(FNDestroy(&f3));
250 : } else {
251 5 : PetscCall(MatDestroy(&F));
252 5 : PetscCall(MatDestroy(&J));
253 : }
254 19 : 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*/
|
