Actual source code: ex36.c
1: /*
2: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
3: SLEPc - Scalable Library for Eigenvalue Problem Computations
4: Copyright (c) 2002-, 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[] = "Use the matrix exponential to compute rightmost eigenvalues.\n\n"
12: "Same problem as ex9.c but with explicitly created matrix. 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>
19: #include <slepcmfn.h>
21: /*
22: This example computes the eigenvalues with largest real part of the
23: following matrix
25: A = [ tau1*T+(beta-1)*I alpha^2*I
26: -beta*I tau2*T-alpha^2*I ],
28: where
30: T = tridiag{1,-2,1}
31: h = 1/(n+1)
32: tau1 = delta1/(h*L)^2
33: tau2 = delta2/(h*L)^2
35: but it builds A explicitly, as opposed to ex9.c
36: */
38: /* Routines for shell spectral transformation */
39: PetscErrorCode STApply_Exp(ST,Vec,Vec);
40: PetscErrorCode STBackTransform_Exp(ST,PetscInt,PetscScalar*,PetscScalar*);
41: PetscErrorCode STDestroy_Exp(ST);
43: int main(int argc,char **argv)
44: {
45: Mat A; /* operator matrix */
46: EPS eps; /* eigenproblem solver context */
47: ST st; /* spectral transformation context */
48: MFN mfn; /* matrix function solver object to compute exp(A)*v */
49: FN f;
50: EPSType type;
51: PetscScalar alpha,beta,tau1,tau2,delta1,delta2,L,h;
52: PetscInt n=30,i,Istart,Iend,nev;
53: PetscBool isShell,terse;
55: PetscFunctionBeginUser;
56: PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
57: #if defined(PETSC_HAVE_COMPLEX)
58: PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
59: PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model with matrix exponential, n=%" PetscInt_FMT "\n\n",n));
61: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
62: Generate the matrix
63: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
65: alpha = 2.0;
66: beta = 5.45;
67: delta1 = 0.008;
68: delta2 = 0.004;
69: L = 0.51302;
71: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
72: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
73: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL));
74: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
75: PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));
77: h = 1.0 / (PetscReal)(n+1);
78: tau1 = delta1 / ((h*L)*(h*L));
79: tau2 = delta2 / ((h*L)*(h*L));
81: PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
82: PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n));
83: PetscCall(MatSetFromOptions(A));
85: PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
86: for (i=Istart;i<Iend;i++) {
87: if (i<n) { /* upper blocks */
88: if (i>0) PetscCall(MatSetValue(A,i,i-1,tau1,INSERT_VALUES));
89: if (i<n-1) PetscCall(MatSetValue(A,i,i+1,tau1,INSERT_VALUES));
90: PetscCall(MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES));
91: PetscCall(MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES));
92: } else { /* lower blocks */
93: if (i>n) PetscCall(MatSetValue(A,i,i-1,tau2,INSERT_VALUES));
94: if (i<2*n-1) PetscCall(MatSetValue(A,i,i+1,tau2,INSERT_VALUES));
95: PetscCall(MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES));
96: PetscCall(MatSetValue(A,i,i-n,-beta,INSERT_VALUES));
97: }
98: }
99: PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
100: PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
102: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103: Create the eigensolver and set various options
104: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
106: PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
107: PetscCall(EPSSetOperators(eps,A,NULL));
108: PetscCall(EPSSetProblemType(eps,EPS_NHEP));
109: PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
110: PetscCall(EPSGetST(eps,&st));
111: PetscCall(STSetType(st,STSHELL));
112: PetscCall(EPSSetFromOptions(eps));
114: /*
115: Initialize shell spectral transformation
116: */
117: PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell));
118: if (isShell) {
120: /* Create the MFN object to be used by the spectral transform */
121: PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
122: PetscCall(MFNSetOperator(mfn,A));
123: PetscCall(MFNGetFN(mfn,&f));
124: PetscCall(FNSetType(f,FNEXP));
125: PetscCall(FNSetScale(f,0.03,1.0)); /* this can be set with -fn_scale */
126: PetscCall(MFNSetFromOptions(mfn));
128: /* Set callback functions */
129: PetscCall(STShellSetApply(st,STApply_Exp));
130: PetscCall(STShellSetBackTransform(st,STBackTransform_Exp));
131: PetscCall(STShellSetContext(st,mfn));
132: PetscCall(STShellSetDestroy(st,STDestroy_Exp));
133: PetscCall(PetscObjectSetName((PetscObject)st,"STEXP"));
134: }
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Solve the eigensystem
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: PetscCall(EPSSolve(eps));
141: PetscCall(EPSGetType(eps,&type));
142: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
143: PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
144: PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Display solution and clean up
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /* show detailed info unless -terse option is given by user */
151: PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
152: if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
153: else {
154: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
155: PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
156: PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
157: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
158: }
159: PetscCall(EPSDestroy(&eps));
160: PetscCall(MatDestroy(&A));
161: #else
162: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires C99 complex numbers");
163: #endif
164: PetscCall(SlepcFinalize());
165: return 0;
166: }
168: /* ------------------------------------------------------------------- */
169: /*
170: STBackTransform_Exp - Undoes the exp(A) transformation by taking logarithms.
172: Input Parameters:
173: + st - spectral transformation context
174: - n - number of eigenvalues to transform
176: Input/Output Parameters:
177: + eigr - pointer to real part of eigenvalues
178: - eigi - pointer to imaginary part of eigenvalues
179: */
180: PetscErrorCode STBackTransform_Exp(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
181: {
182: #if defined(PETSC_HAVE_COMPLEX)
183: PetscInt j;
184: MFN mfn;
185: FN fn;
186: PetscScalar tau,eta;
187: #if !defined(PETSC_USE_COMPLEX)
188: PetscComplex theta,lambda;
189: #endif
191: PetscFunctionBeginUser;
192: PetscCall(STShellGetContext(st,&mfn));
193: PetscCall(MFNGetFN(mfn,&fn));
194: PetscCall(FNGetScale(fn,&tau,&eta));
195: for (j=0;j<n;j++) {
196: #if defined(PETSC_USE_COMPLEX)
197: eigr[j] = PetscLogComplex(eigr[j]/eta)/tau;
198: #else
199: theta = PetscCMPLX(eigr[j],eigi[j])/eta;
200: lambda = PetscLogComplex(theta)/tau;
201: eigr[j] = PetscRealPartComplex(lambda);
202: eigi[j] = PetscImaginaryPartComplex(lambda);
203: #endif
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: #else
207: return 0;
208: #endif
209: }
211: /*
212: STApply_Exp - Applies the operator exp(tau*A) to a given vector using an MFN object.
214: Input Parameters:
215: + st - spectral transformation context
216: - x - input vector
218: Output Parameter:
219: . y - output vector
220: */
221: PetscErrorCode STApply_Exp(ST st,Vec x,Vec y)
222: {
223: MFN mfn;
225: PetscFunctionBeginUser;
226: PetscCall(STShellGetContext(st,&mfn));
227: PetscCall(MFNSolve(mfn,x,y));
228: PetscFunctionReturn(PETSC_SUCCESS);
229: }
231: /*
232: STDestroy_Exp - This routine destroys the shell ST context.
234: Input Parameter:
235: . st - spectral transformation context
236: */
237: PetscErrorCode STDestroy_Exp(ST st)
238: {
239: MFN mfn;
241: PetscFunctionBeginUser;
242: PetscCall(STShellGetContext(st,&mfn));
243: PetscCall(MFNDestroy(&mfn));
244: PetscFunctionReturn(PETSC_SUCCESS);
245: }
247: /*TEST
249: testset:
250: args: -eps_nev 4 -mfn_ncv 16 -terse
251: requires: c99_complex !single
252: filter: sed -e "s/-2/+2/g"
253: output_file: output/ex36_1.out
254: test:
255: suffix: 1
256: requires: !__float128
257: test:
258: suffix: 1_quad
259: args: -eps_tol 1e-14
260: requires: __float128
262: test:
263: suffix: 2
264: args: -n 56 -eps_nev 2 -st_type sinvert -eps_target -390 -eps_target_magnitude -eps_type power
265: args: -eps_power_shift_type {{constant rayleigh}} -eps_two_sided {{0 1}} -eps_tol 1e-14 -terse
266: requires: c99_complex !single
267: filter: sed -e "s/[+-]0\.0*i//g"
269: test:
270: suffix: 3
271: args: -n 100 -st_type sinvert -eps_type ciss -rg_type ellipse -rg_ellipse_center 0 -rg_ellipse_radius 6 -eps_all -eps_tol 1e-6 -terse
272: requires: c99_complex !single
273: filter: sed -e "s/-3.37036-3.55528i, -3.37036+3.55528i/-3.37036+3.55528i, -3.37036-3.55528i/" -e "s/-1.79853-3.03216i, -1.79853+3.03216i/-1.79853+3.03216i, -1.79853-3.03216i/" -e "s/-0.67471-2.52856i, -0.67471+2.52856i/-0.67471+2.52856i, -0.67471-2.52856i/" -e "s/0.00002-2.13950i, 0.00002+2.13950i/0.00002+2.13950i, 0.00002-2.13950i/"
275: TEST*/