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[] = "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";
17 :
18 : #include <slepceps.h>
19 : #include <slepcmfn.h>
20 :
21 : /*
22 : This example computes the eigenvalues with largest real part of the
23 : following matrix
24 :
25 : A = [ tau1*T+(beta-1)*I alpha^2*I
26 : -beta*I tau2*T-alpha^2*I ],
27 :
28 : where
29 :
30 : T = tridiag{1,-2,1}
31 : h = 1/(n+1)
32 : tau1 = delta1/(h*L)^2
33 : tau2 = delta2/(h*L)^2
34 :
35 : but it builds A explicitly, as opposed to ex9.c
36 : */
37 :
38 : /* Routines for shell spectral transformation */
39 : PetscErrorCode STApply_Exp(ST,Vec,Vec);
40 : PetscErrorCode STBackTransform_Exp(ST,PetscInt,PetscScalar*,PetscScalar*);
41 :
42 6 : int main(int argc,char **argv)
43 : {
44 6 : Mat A; /* operator matrix */
45 6 : EPS eps; /* eigenproblem solver context */
46 6 : ST st; /* spectral transformation context */
47 6 : MFN mfn; /* matrix function solver object to compute exp(A)*v */
48 6 : FN f;
49 6 : EPSType type;
50 6 : PetscScalar alpha,beta,tau1,tau2,delta1,delta2,L,h;
51 6 : PetscInt n=30,i,Istart,Iend,nev;
52 6 : PetscBool isShell,terse;
53 :
54 6 : PetscFunctionBeginUser;
55 6 : PetscCall(SlepcInitialize(&argc,&argv,NULL,help));
56 : #if defined(PETSC_HAVE_COMPLEX)
57 6 : PetscCall(PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL));
58 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model with matrix exponential, n=%" PetscInt_FMT "\n\n",n));
59 :
60 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61 : Generate the matrix
62 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63 :
64 6 : alpha = 2.0;
65 6 : beta = 5.45;
66 6 : delta1 = 0.008;
67 6 : delta2 = 0.004;
68 6 : L = 0.51302;
69 :
70 6 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL));
71 6 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL));
72 6 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL));
73 6 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL));
74 6 : PetscCall(PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL));
75 :
76 6 : h = 1.0 / (PetscReal)(n+1);
77 6 : tau1 = delta1 / ((h*L)*(h*L));
78 6 : tau2 = delta2 / ((h*L)*(h*L));
79 :
80 6 : PetscCall(MatCreate(PETSC_COMM_WORLD,&A));
81 6 : PetscCall(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n));
82 6 : PetscCall(MatSetFromOptions(A));
83 :
84 6 : PetscCall(MatGetOwnershipRange(A,&Istart,&Iend));
85 714 : for (i=Istart;i<Iend;i++) {
86 708 : if (i<n) { /* upper blocks */
87 354 : if (i>0) PetscCall(MatSetValue(A,i,i-1,tau1,INSERT_VALUES));
88 354 : if (i<n-1) PetscCall(MatSetValue(A,i,i+1,tau1,INSERT_VALUES));
89 354 : PetscCall(MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES));
90 354 : PetscCall(MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES));
91 : } else { /* lower blocks */
92 354 : if (i>n) PetscCall(MatSetValue(A,i,i-1,tau2,INSERT_VALUES));
93 354 : if (i<2*n-1) PetscCall(MatSetValue(A,i,i+1,tau2,INSERT_VALUES));
94 354 : PetscCall(MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES));
95 708 : PetscCall(MatSetValue(A,i,i-n,-beta,INSERT_VALUES));
96 : }
97 : }
98 6 : PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY));
99 6 : PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY));
100 :
101 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
102 : Create the eigensolver and set various options
103 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
104 :
105 6 : PetscCall(EPSCreate(PETSC_COMM_WORLD,&eps));
106 6 : PetscCall(EPSSetOperators(eps,A,NULL));
107 6 : PetscCall(EPSSetProblemType(eps,EPS_NHEP));
108 6 : PetscCall(EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL));
109 6 : PetscCall(EPSGetST(eps,&st));
110 6 : PetscCall(STSetType(st,STSHELL));
111 6 : PetscCall(EPSSetFromOptions(eps));
112 :
113 : /*
114 : Initialize shell spectral transformation
115 : */
116 6 : PetscCall(PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell));
117 6 : if (isShell) {
118 :
119 : /* Create the MFN object to be used by the spectral transform */
120 1 : PetscCall(MFNCreate(PETSC_COMM_WORLD,&mfn));
121 1 : PetscCall(MFNSetOperator(mfn,A));
122 1 : PetscCall(MFNGetFN(mfn,&f));
123 1 : PetscCall(FNSetType(f,FNEXP));
124 1 : PetscCall(FNSetScale(f,0.03,1.0)); /* this can be set with -fn_scale */
125 1 : PetscCall(MFNSetFromOptions(mfn));
126 :
127 : /* Set callback functions */
128 1 : PetscCall(STShellSetApply(st,STApply_Exp));
129 1 : PetscCall(STShellSetBackTransform(st,STBackTransform_Exp));
130 1 : PetscCall(STShellSetContext(st,mfn));
131 1 : PetscCall(PetscObjectSetName((PetscObject)st,"STEXP"));
132 : }
133 :
134 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135 : Solve the eigensystem
136 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
137 :
138 6 : PetscCall(EPSSolve(eps));
139 6 : PetscCall(EPSGetType(eps,&type));
140 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type));
141 6 : PetscCall(EPSGetDimensions(eps,&nev,NULL,NULL));
142 6 : PetscCall(PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %" PetscInt_FMT "\n",nev));
143 :
144 : /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145 : Display solution and clean up
146 : - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147 :
148 : /* show detailed info unless -terse option is given by user */
149 6 : PetscCall(PetscOptionsHasName(NULL,NULL,"-terse",&terse));
150 6 : if (terse) PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL));
151 : else {
152 0 : PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL));
153 0 : PetscCall(EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD));
154 0 : PetscCall(EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD));
155 0 : PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
156 : }
157 6 : PetscCall(EPSDestroy(&eps));
158 6 : PetscCall(MatDestroy(&A));
159 6 : if (isShell) PetscCall(MFNDestroy(&mfn));
160 : #else
161 : SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires C99 complex numbers");
162 : #endif
163 6 : PetscCall(SlepcFinalize());
164 : return 0;
165 : }
166 :
167 : /* ------------------------------------------------------------------- */
168 : /*
169 : STBackTransform_Exp - Undoes the exp(A) transformation by taking logarithms.
170 :
171 : Input Parameters:
172 : + st - spectral transformation context
173 : - n - number of eigenvalues to transform
174 :
175 : Input/Output Parameters:
176 : + eigr - pointer to real part of eigenvalues
177 : - eigi - pointer to imaginary part of eigenvalues
178 : */
179 357 : PetscErrorCode STBackTransform_Exp(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
180 : {
181 : #if defined(PETSC_HAVE_COMPLEX)
182 357 : PetscInt j;
183 357 : MFN mfn;
184 357 : FN fn;
185 357 : PetscScalar tau,eta;
186 : #if !defined(PETSC_USE_COMPLEX)
187 357 : PetscComplex theta,lambda;
188 : #endif
189 :
190 357 : PetscFunctionBeginUser;
191 357 : PetscCall(STShellGetContext(st,&mfn));
192 357 : PetscCall(MFNGetFN(mfn,&fn));
193 357 : PetscCall(FNGetScale(fn,&tau,&eta));
194 1075 : for (j=0;j<n;j++) {
195 : #if defined(PETSC_USE_COMPLEX)
196 : eigr[j] = PetscLogComplex(eigr[j]/eta)/tau;
197 : #else
198 718 : theta = PetscCMPLX(eigr[j],eigi[j])/eta;
199 718 : lambda = PetscLogComplex(theta)/tau;
200 718 : eigr[j] = PetscRealPartComplex(lambda);
201 718 : eigi[j] = PetscImaginaryPartComplex(lambda);
202 : #endif
203 : }
204 357 : PetscFunctionReturn(PETSC_SUCCESS);
205 : #else
206 : return 0;
207 : #endif
208 : }
209 :
210 : /*
211 : STApply_Exp - Applies the operator exp(tau*A) to a given vector using an MFN object.
212 :
213 : Input Parameters:
214 : + st - spectral transformation context
215 : - x - input vector
216 :
217 : Output Parameter:
218 : . y - output vector
219 : */
220 56 : PetscErrorCode STApply_Exp(ST st,Vec x,Vec y)
221 : {
222 56 : MFN mfn;
223 :
224 56 : PetscFunctionBeginUser;
225 56 : PetscCall(STShellGetContext(st,&mfn));
226 56 : PetscCall(MFNSolve(mfn,x,y));
227 56 : PetscFunctionReturn(PETSC_SUCCESS);
228 : }
229 :
230 : /*TEST
231 :
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
244 :
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"
251 :
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/"
257 :
258 : TEST*/
|
