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