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[] = "Solves a problem associated to the Brusselator wave model in chemical reactions, illustrating the use of shell matrices.\n\n"
12 : "The command line options are:\n"
13 : " -n <n>, where <n> = block dimension of the 2x2 block matrix.\n"
14 : " -L <L>, where <L> = bifurcation parameter.\n"
15 : " -alpha <alpha>, -beta <beta>, -delta1 <delta1>, -delta2 <delta2>,\n"
16 : " where <alpha> <beta> <delta1> <delta2> = model parameters.\n\n";
17 :
18 : #include <slepceps.h>
19 :
20 : /*
21 : This example computes the eigenvalues with largest real part of the
22 : following matrix
23 :
24 : A = [ tau1*T+(beta-1)*I alpha^2*I
25 : -beta*I tau2*T-alpha^2*I ],
26 :
27 : where
28 :
29 : T = tridiag{1,-2,1}
30 : h = 1/(n+1)
31 : tau1 = delta1/(h*L)^2
32 : tau2 = delta2/(h*L)^2
33 : */
34 :
35 : /*
36 : Matrix operations
37 : */
38 : PetscErrorCode MatMult_Brussel(Mat,Vec,Vec);
39 : PetscErrorCode MatMultTranspose_Brussel(Mat,Vec,Vec);
40 : PetscErrorCode MatGetDiagonal_Brussel(Mat,Vec);
41 :
42 : typedef struct {
43 : Mat T;
44 : Vec x1,x2,y1,y2;
45 : PetscScalar alpha,beta,tau1,tau2,sigma;
46 : } CTX_BRUSSEL;
47 :
48 12 : int main(int argc,char **argv)
49 : {
50 12 : Mat A; /* eigenvalue problem matrix */
51 12 : EPS eps; /* eigenproblem solver context */
52 12 : EPSType type;
53 12 : PetscScalar delta1,delta2,L,h;
54 12 : PetscInt N=30,n,i,Istart,Iend,nev;
55 12 : CTX_BRUSSEL *ctx;
56 12 : PetscBool terse;
57 12 : PetscViewer viewer;
58 :
59 12 : PetscFunctionBeginUser;
60 12 : PetscCall(SlepcInitialize(&argc,&argv,(char*)0,help));
61 :
62 12 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&N,NULL));
63 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%" PetscInt_FMT "\n\n",N));
64 :
65 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
66 : Generate the matrix
67 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
68 :
69 : /*
70 : Create shell matrix context and set default parameters
71 : */
72 12 : PetscCall(PetscNew(&ctx));
73 12 : ctx->alpha = 2.0;
74 12 : ctx->beta = 5.45;
75 12 : delta1 = 0.008;
76 12 : delta2 = 0.004;
77 12 : L = 0.51302;
78 :
79 : /*
80 : Look the command line for user-provided parameters
81 : */
82 12 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
83 12 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&ctx->alpha,NULL));
84 12 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&ctx->beta,NULL));
85 12 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
86 12 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));
87 :
88 : /*
89 : Create matrix T
90 : */
91 12 : PetscCall(MatCreate(PETSC_COMM_WORLD,&ctx->T));
92 12 : PetscCall(MatSetSizes(ctx->T,PETSC_DECIDE,PETSC_DECIDE,N,N));
93 12 : PetscCall(MatSetFromOptions(ctx->T));
94 :
95 12 : PetscCall(MatGetOwnershipRange(ctx->T,&Istart,&Iend));
96 482 : for (i=Istart;i<Iend;i++) {
97 470 : if (i>0) PetscCall(MatSetValue(ctx->T,i,i-1,1.0,INSERT_VALUES));
98 470 : if (i<N-1) PetscCall(MatSetValue(ctx->T,i,i+1,1.0,INSERT_VALUES));
99 470 : PetscCall(MatSetValue(ctx->T,i,i,-2.0,INSERT_VALUES));
100 : }
101 12 : PetscCall(MatAssemblyBegin(ctx->T,MAT_FINAL_ASSEMBLY));
102 12 : PetscCall(MatAssemblyEnd(ctx->T,MAT_FINAL_ASSEMBLY));
103 12 : PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
104 :
105 : /*
106 : Fill the remaining information in the shell matrix context
107 : and create auxiliary vectors
108 : */
109 12 : h = 1.0 / (PetscReal)(N+1);
110 12 : ctx->tau1 = delta1 / ((h*L)*(h*L));
111 12 : ctx->tau2 = delta2 / ((h*L)*(h*L));
112 12 : ctx->sigma = 0.0;
113 12 : PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->x1));
114 12 : PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->x2));
115 12 : PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->y1));
116 12 : PetscCall(VecCreateMPIWithArray(PETSC_COMM_WORLD,1,n,PETSC_DECIDE,NULL,&ctx->y2));
117 :
118 : /*
119 : Create the shell matrix
120 : */
121 12 : PetscCall(MatCreateShell(PETSC_COMM_WORLD,2*n,2*n,2*N,2*N,(void*)ctx,&A));
122 12 : PetscCall(MatShellSetOperation(A,MATOP_MULT,(void(*)(void))MatMult_Brussel));
123 12 : PetscCall(MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_Brussel));
124 12 : PetscCall(MatShellSetOperation(A,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_Brussel));
125 :
126 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127 : Create the eigensolver and set various options
128 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129 :
130 : /*
131 : Create eigensolver context
132 : */
133 12 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
134 :
135 : /*
136 : Set operators. In this case, it is a standard eigenvalue problem
137 : */
138 12 : PetscCall(EPSSetOperators(eps,A,NULL));
139 12 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
140 :
141 : /*
142 : Ask for the rightmost eigenvalues
143 : */
144 12 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
145 :
146 : /*
147 : Set other solver options at runtime
148 : */
149 12 : PetscCall(EPSSetFromOptions(eps));
150 :
151 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152 : Solve the eigensystem
153 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
154 :
155 12 : PetscCall(EPSSolve(eps));
156 :
157 : /*
158 : Optional: Get some information from the solver and display it
159 : */
160 12 : PetscCall(EPSGetType(eps,&type));
161 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
162 12 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
163 12 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
164 :
165 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166 : Display solution and clean up
167 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
168 :
169 : /* show detailed info unless -terse option is given by user */
170 12 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
171 12 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
172 : else {
173 0 : PetscCall(PetscViewerASCIIGetStdout(PETSC_COMM_WORLD,&viewer));
174 0 : PetscCall(PetscViewerPushFormat(viewer,PETSC_VIEWER_ASCII_INFO_DETAIL));
175 0 : PetscCall(EPSConvergedReasonView(eps,viewer));
176 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,viewer));
177 0 : PetscCall(PetscViewerPopFormat(viewer));
178 : }
179 12 : PetscCall(EPSDestroy(&eps));
180 12 : PetscCall(MatDestroy(&A));
181 12 : PetscCall(MatDestroy(&ctx->T));
182 12 : PetscCall(VecDestroy(&ctx->x1));
183 12 : PetscCall(VecDestroy(&ctx->x2));
184 12 : PetscCall(VecDestroy(&ctx->y1));
185 12 : PetscCall(VecDestroy(&ctx->y2));
186 12 : PetscCall(PetscFree(ctx));
187 12 : PetscCall(SlepcFinalize());
188 : return 0;
189 : }
190 :
191 30698 : PetscErrorCode MatMult_Brussel(Mat A,Vec x,Vec y)
192 : {
193 30698 : PetscInt n;
194 30698 : const PetscScalar *px;
195 30698 : PetscScalar *py;
196 30698 : CTX_BRUSSEL *ctx;
197 :
198 30698 : PetscFunctionBeginUser;
199 30698 : PetscCall(MatShellGetContext(A,&ctx));
200 30698 : PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
201 30698 : PetscCall(VecGetArrayRead(x,&px));
202 30698 : PetscCall(VecGetArray(y,&py));
203 30698 : PetscCall(VecPlaceArray(ctx->x1,px));
204 30698 : PetscCall(VecPlaceArray(ctx->x2,px+n));
205 30698 : PetscCall(VecPlaceArray(ctx->y1,py));
206 30698 : PetscCall(VecPlaceArray(ctx->y2,py+n));
207 :
208 30698 : PetscCall(MatMult(ctx->T,ctx->x1,ctx->y1));
209 30698 : PetscCall(VecScale(ctx->y1,ctx->tau1));
210 30698 : PetscCall(VecAXPY(ctx->y1,ctx->beta-1.0+ctx->sigma,ctx->x1));
211 30698 : PetscCall(VecAXPY(ctx->y1,ctx->alpha*ctx->alpha,ctx->x2));
212 :
213 30698 : PetscCall(MatMult(ctx->T,ctx->x2,ctx->y2));
214 30698 : PetscCall(VecScale(ctx->y2,ctx->tau2));
215 30698 : PetscCall(VecAXPY(ctx->y2,-ctx->beta,ctx->x1));
216 30698 : PetscCall(VecAXPY(ctx->y2,-ctx->alpha*ctx->alpha+ctx->sigma,ctx->x2));
217 :
218 30698 : PetscCall(VecRestoreArrayRead(x,&px));
219 30698 : PetscCall(VecRestoreArray(y,&py));
220 30698 : PetscCall(VecResetArray(ctx->x1));
221 30698 : PetscCall(VecResetArray(ctx->x2));
222 30698 : PetscCall(VecResetArray(ctx->y1));
223 30698 : PetscCall(VecResetArray(ctx->y2));
224 30698 : PetscFunctionReturn(PETSC_SUCCESS);
225 : }
226 :
227 250 : PetscErrorCode MatMultTranspose_Brussel(Mat A,Vec x,Vec y)
228 : {
229 250 : PetscInt n;
230 250 : const PetscScalar *px;
231 250 : PetscScalar *py;
232 250 : CTX_BRUSSEL *ctx;
233 :
234 250 : PetscFunctionBeginUser;
235 250 : PetscCall(MatShellGetContext(A,&ctx));
236 250 : PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
237 250 : PetscCall(VecGetArrayRead(x,&px));
238 250 : PetscCall(VecGetArray(y,&py));
239 250 : PetscCall(VecPlaceArray(ctx->x1,px));
240 250 : PetscCall(VecPlaceArray(ctx->x2,px+n));
241 250 : PetscCall(VecPlaceArray(ctx->y1,py));
242 250 : PetscCall(VecPlaceArray(ctx->y2,py+n));
243 :
244 250 : PetscCall(MatMultTranspose(ctx->T,ctx->x1,ctx->y1));
245 250 : PetscCall(VecScale(ctx->y1,ctx->tau1));
246 250 : PetscCall(VecAXPY(ctx->y1,ctx->beta-1.0+ctx->sigma,ctx->x1));
247 250 : PetscCall(VecAXPY(ctx->y1,-ctx->beta,ctx->x2));
248 :
249 250 : PetscCall(MatMultTranspose(ctx->T,ctx->x2,ctx->y2));
250 250 : PetscCall(VecScale(ctx->y2,ctx->tau2));
251 250 : PetscCall(VecAXPY(ctx->y2,ctx->alpha*ctx->alpha,ctx->x1));
252 250 : PetscCall(VecAXPY(ctx->y2,-ctx->alpha*ctx->alpha+ctx->sigma,ctx->x2));
253 :
254 250 : PetscCall(VecRestoreArrayRead(x,&px));
255 250 : PetscCall(VecRestoreArray(y,&py));
256 250 : PetscCall(VecResetArray(ctx->x1));
257 250 : PetscCall(VecResetArray(ctx->x2));
258 250 : PetscCall(VecResetArray(ctx->y1));
259 250 : PetscCall(VecResetArray(ctx->y2));
260 250 : PetscFunctionReturn(PETSC_SUCCESS);
261 : }
262 :
263 1 : PetscErrorCode MatGetDiagonal_Brussel(Mat A,Vec diag)
264 : {
265 1 : Vec d1,d2;
266 1 : PetscInt n;
267 1 : PetscScalar *pd;
268 1 : MPI_Comm comm;
269 1 : CTX_BRUSSEL *ctx;
270 :
271 1 : PetscFunctionBeginUser;
272 1 : PetscCall(MatShellGetContext(A,&ctx));
273 1 : PetscCall(PetscObjectGetComm((PetscObject)A,&comm));
274 1 : PetscCall(MatGetLocalSize(ctx->T,&n,NULL));
275 1 : PetscCall(VecGetArray(diag,&pd));
276 1 : PetscCall(VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd,&d1));
277 1 : PetscCall(VecCreateMPIWithArray(comm,1,n,PETSC_DECIDE,pd+n,&d2));
278 :
279 1 : PetscCall(VecSet(d1,-2.0*ctx->tau1 + ctx->beta - 1.0 + ctx->sigma));
280 1 : PetscCall(VecSet(d2,-2.0*ctx->tau2 - ctx->alpha*ctx->alpha + ctx->sigma));
281 :
282 1 : PetscCall(VecDestroy(&d1));
283 1 : PetscCall(VecDestroy(&d2));
284 1 : PetscCall(VecRestoreArray(diag,&pd));
285 1 : PetscFunctionReturn(PETSC_SUCCESS);
286 : }
287 :
288 : /*TEST
289 :
290 : test:
291 : suffix: 1
292 : args: -n 50 -eps_nev 4 -eps_two_sided {{0 1}} -eps_type {{krylovschur lapack}} -terse
293 : requires: !single
294 : filter: grep -v method
295 :
296 : test:
297 : suffix: 2
298 : args: -eps_nev 8 -eps_max_it 300 -eps_target -28 -rg_type interval -rg_interval_endpoints -40,-20,-.1,.1 -terse
299 : requires: !single
300 :
301 : test:
302 : suffix: 3
303 : args: -n 50 -eps_nev 4 -eps_balance twoside -terse
304 : requires: double
305 : filter: grep -v method
306 : output_file: output/ex9_1.out
307 :
308 : test:
309 : suffix: 4
310 : args: -eps_smallest_imaginary -eps_ncv 24 -terse
311 : requires: !complex !single
312 :
313 : test:
314 : suffix: 4_complex
315 : args: -eps_smallest_imaginary -eps_ncv 24 -terse
316 : requires: complex !single
317 :
318 : test:
319 : suffix: 5
320 : args: -eps_nev 4 -eps_target_real -eps_target -3 -terse
321 : requires: !single
322 :
323 : test:
324 : suffix: 6
325 : args: -eps_nev 2 -eps_target_imaginary -eps_target 3i -terse
326 : requires: complex !single
327 :
328 : test:
329 : suffix: 7
330 : args: -n 40 -eps_nev 1 -eps_type arnoldi -eps_smallest_real -eps_refined -eps_ncv 42 -eps_max_it 300 -terse
331 : requires: double
332 :
333 : test:
334 : suffix: 8
335 : args: -eps_nev 2 -eps_target -30 -eps_type jd -st_matmode shell -eps_jd_fix 0.0001 -eps_jd_const_correction_tol 0 -terse
336 : requires: !single
337 : filter: sed -e "s/[+-]0\.0*i//g"
338 : timeoutfactor: 2
339 :
340 : test:
341 : suffix: 9
342 : args: -eps_largest_imaginary -eps_ncv 24 -terse
343 : requires: !single
344 :
345 : TEST*/
|
