Actual source code: ex36.c
slepc-3.17.1 2022-04-11
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*);
42: int main(int argc,char **argv)
43: {
44: Mat A; /* operator matrix */
45: EPS eps; /* eigenproblem solver context */
46: ST st; /* spectral transformation context */
47: MFN mfn; /* matrix function solver object to compute exp(A)*v */
48: FN f;
49: EPSType type;
50: PetscScalar alpha,beta,tau1,tau2,delta1,delta2,L,h;
51: PetscInt n=30,i,Istart,Iend,nev;
52: PetscBool isShell,terse;
54: SlepcInitialize(&argc,&argv,(char*)0,help);
55: #if defined(PETSC_HAVE_COMPLEX)
56: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
57: PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model with matrix exponential, n=%" PetscInt_FMT "\n\n",n);
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Generate the matrix
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: alpha = 2.0;
64: beta = 5.45;
65: delta1 = 0.008;
66: delta2 = 0.004;
67: L = 0.51302;
69: PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL);
70: PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
71: PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL);
72: PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL);
73: PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL);
75: h = 1.0 / (PetscReal)(n+1);
76: tau1 = delta1 / ((h*L)*(h*L));
77: tau2 = delta2 / ((h*L)*(h*L));
79: MatCreate(PETSC_COMM_WORLD,&A);
80: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n);
81: MatSetFromOptions(A);
82: MatSetUp(A);
84: MatGetOwnershipRange(A,&Istart,&Iend);
85: for (i=Istart;i<Iend;i++) {
86: if (i<n) { /* upper blocks */
87: if (i>0) MatSetValue(A,i,i-1,tau1,INSERT_VALUES);
88: if (i<n-1) MatSetValue(A,i,i+1,tau1,INSERT_VALUES);
89: MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES);
90: MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES);
91: } else { /* lower blocks */
92: if (i>n) MatSetValue(A,i,i-1,tau2,INSERT_VALUES);
93: if (i<2*n-1) MatSetValue(A,i,i+1,tau2,INSERT_VALUES);
94: MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES);
95: MatSetValue(A,i,i-n,-beta,INSERT_VALUES);
96: }
97: }
98: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
99: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
101: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102: Create the eigensolver and set various options
103: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
105: EPSCreate(PETSC_COMM_WORLD,&eps);
106: EPSSetOperators(eps,A,NULL);
107: EPSSetProblemType(eps,EPS_NHEP);
108: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
109: EPSGetST(eps,&st);
110: STSetType(st,STSHELL);
111: EPSSetFromOptions(eps);
113: /*
114: Initialize shell spectral transformation
115: */
116: PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
117: if (isShell) {
119: /* Create the MFN object to be used by the spectral transform */
120: MFNCreate(PETSC_COMM_WORLD,&mfn);
121: MFNSetOperator(mfn,A);
122: MFNGetFN(mfn,&f);
123: FNSetType(f,FNEXP);
124: FNSetScale(f,0.03,1.0); /* this can be set with -fn_scale */
125: MFNSetFromOptions(mfn);
127: /* Set callback functions */
128: STShellSetApply(st,STApply_Exp);
129: STShellSetBackTransform(st,STBackTransform_Exp);
130: STShellSetContext(st,mfn);
131: PetscObjectSetName((PetscObject)st,"STEXP");
132: }
134: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: Solve the eigensystem
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: EPSSolve(eps);
139: EPSGetType(eps,&type);
140: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
141: EPSGetDimensions(eps,&nev,NULL,NULL);
142: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev);
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Display solution and clean up
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
148: /* show detailed info unless -terse option is given by user */
149: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
150: if (terse) EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
151: else {
152: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
153: EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
154: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
155: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
156: }
157: EPSDestroy(&eps);
158: MatDestroy(&A);
159: if (isShell) MFNDestroy(&mfn);
160: #else
161: SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires C99 complex numbers");
162: #endif
163: SlepcFinalize();
164: return 0;
165: }
167: /* ------------------------------------------------------------------- */
168: /*
169: STBackTransform_Exp - Undoes the exp(A) transformation by taking logarithms.
171: Input Parameters:
172: + st - spectral transformation context
173: - n - number of eigenvalues to transform
175: Input/Output Parameters:
176: + eigr - pointer to real part of eigenvalues
177: - eigi - pointer to imaginary part of eigenvalues
178: */
179: PetscErrorCode STBackTransform_Exp(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
180: {
181: #if defined(PETSC_HAVE_COMPLEX)
182: PetscInt j;
183: MFN mfn;
184: FN fn;
185: PetscScalar tau,eta;
186: #if !defined(PETSC_USE_COMPLEX)
187: PetscComplex theta,lambda;
188: #endif
191: STShellGetContext(st,&mfn);
192: MFNGetFN(mfn,&fn);
193: FNGetScale(fn,&tau,&eta);
194: for (j=0;j<n;j++) {
195: #if defined(PETSC_USE_COMPLEX)
196: eigr[j] = PetscLogComplex(eigr[j]/eta)/tau;
197: #else
198: theta = PetscCMPLX(eigr[j],eigi[j])/eta;
199: lambda = PetscLogComplex(theta)/tau;
200: eigr[j] = PetscRealPartComplex(lambda);
201: eigi[j] = PetscImaginaryPartComplex(lambda);
202: #endif
203: }
204: PetscFunctionReturn(0);
205: #else
206: return 0;
207: #endif
208: }
210: /*
211: STApply_Exp - Applies the operator exp(tau*A) to a given vector using an MFN object.
213: Input Parameters:
214: + st - spectral transformation context
215: - x - input vector
217: Output Parameter:
218: . y - output vector
219: */
220: PetscErrorCode STApply_Exp(ST st,Vec x,Vec y)
221: {
222: MFN mfn;
225: STShellGetContext(st,&mfn);
226: MFNSolve(mfn,x,y);
227: PetscFunctionReturn(0);
228: }
230: /*TEST
232: testset:
233: args: -eps_nev 4 -mfn_ncv 16 -terse
234: requires: c99_complex !single
235: filter: sed -e "s/-2/+2/g"
236: output_file: output/ex36_1.out
237: test:
238: suffix: 1
239: requires: !__float128
240: test:
241: suffix: 1_quad
242: args: -eps_tol 1e-14
243: requires: __float128
245: test:
246: suffix: 2
247: args: -n 56 -eps_nev 2 -st_type sinvert -eps_target -390 -eps_target_magnitude -eps_type power
248: args: -eps_power_shift_type {{constant rayleigh}} -eps_two_sided {{0 1}} -eps_tol 1e-14 -terse
249: requires: c99_complex !single
250: filter: sed -e "s/[+-]0\.0*i//g"
252: test:
253: suffix: 3
254: 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
255: requires: c99_complex !single
256: 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/"
258: TEST*/