Actual source code: ex36.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[] = "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: */
39: /* Routines for shell spectral transformation */
40: PetscErrorCode STApply_Exp(ST,Vec,Vec);
41: PetscErrorCode STBackTransform_Exp(ST,PetscInt,PetscScalar*,PetscScalar*);
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;
56: SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
57: #if defined(PETSC_HAVE_COMPLEX)
58: PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
59: PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model with matrix exponential, n=%D\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: PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL);
72: PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
73: PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL);
74: PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL);
75: 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: MatCreate(PETSC_COMM_WORLD,&A);
82: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n);
83: MatSetFromOptions(A);
84: MatSetUp(A);
86: MatGetOwnershipRange(A,&Istart,&Iend);
87: for (i=Istart;i<Iend;i++) {
88: if (i<n) { /* upper blocks */
89: if (i>0) { MatSetValue(A,i,i-1,tau1,INSERT_VALUES); }
90: if (i<n-1) { MatSetValue(A,i,i+1,tau1,INSERT_VALUES); }
91: MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES);
92: MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES);
93: } else { /* lower blocks */
94: if (i>n) { MatSetValue(A,i,i-1,tau2,INSERT_VALUES); }
95: if (i<2*n-1) { MatSetValue(A,i,i+1,tau2,INSERT_VALUES); }
96: MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES);
97: MatSetValue(A,i,i-n,-beta,INSERT_VALUES);
98: }
99: }
100: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
101: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create the eigensolver and set various options
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: EPSCreate(PETSC_COMM_WORLD,&eps);
108: EPSSetOperators(eps,A,NULL);
109: EPSSetProblemType(eps,EPS_NHEP);
110: EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
111: EPSGetST(eps,&st);
112: STSetType(st,STSHELL);
113: EPSSetFromOptions(eps);
115: /*
116: Initialize shell spectral transformation
117: */
118: PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
119: if (isShell) {
121: /* Create the MFN object to be used by the spectral transform */
122: MFNCreate(PETSC_COMM_WORLD,&mfn);
123: MFNSetOperator(mfn,A);
124: MFNGetFN(mfn,&f);
125: FNSetType(f,FNEXP);
126: FNSetScale(f,0.03,1.0); /* this can be set with -fn_scale */
127: MFNSetFromOptions(mfn);
129: /* Set callback functions */
130: STShellSetApply(st,STApply_Exp);
131: STShellSetBackTransform(st,STBackTransform_Exp);
132: STShellSetContext(st,mfn);
133: PetscObjectSetName((PetscObject)st,"STEXP");
134: }
136: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
137: Solve the eigensystem
138: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140: EPSSolve(eps);
141: EPSGetType(eps,&type);
142: PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
143: EPSGetDimensions(eps,&nev,NULL,NULL);
144: PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Display solution and clean up
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /* show detailed info unless -terse option is given by user */
151: PetscOptionsHasName(NULL,NULL,"-terse",&terse);
152: if (terse) {
153: EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
154: } else {
155: PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
156: EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
157: EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
158: PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
159: }
160: EPSDestroy(&eps);
161: MatDestroy(&A);
162: if (isShell) { MFNDestroy(&mfn); }
163: #else
164: SETERRQ(PETSC_COMM_WORLD,1,"This examples requires C99 complex numbers");
165: #endif
166: SlepcFinalize();
167: return ierr;
168: }
170: /* ------------------------------------------------------------------- */
171: /*
172: STBackTransform_Exp - Undoes the exp(A) transformation by taking logarithms.
174: Input Parameters:
175: + st - spectral transformation context
176: - n - number of eigenvalues to transform
178: Input/Output Parameters:
179: + eigr - pointer to real part of eigenvalues
180: - eigi - pointer to imaginary part of eigenvalues
181: */
182: PetscErrorCode STBackTransform_Exp(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
183: {
184: #if defined(PETSC_HAVE_COMPLEX)
186: PetscInt j;
187: MFN mfn;
188: FN fn;
189: PetscScalar tau,eta;
190: #if !defined(PETSC_USE_COMPLEX)
191: PetscComplex theta,lambda;
192: #endif
195: STShellGetContext(st,(void**)&mfn);
196: MFNGetFN(mfn,&fn);
197: FNGetScale(fn,&tau,&eta);
198: for (j=0;j<n;j++) {
199: #if defined(PETSC_USE_COMPLEX)
200: eigr[j] = PetscLogComplex(eigr[j]/eta)/tau;
201: #else
202: theta = (eigr[j]+eigi[j]*PETSC_i)/eta;
203: lambda = PetscLogComplex(theta)/tau;
204: eigr[j] = PetscRealPartComplex(lambda);
205: eigi[j] = PetscImaginaryPartComplex(lambda);
206: #endif
207: }
208: return(0);
209: #else
210: return 0;
211: #endif
212: }
214: /*
215: STApply_Exp - Applies the operator exp(tau*A) to a given vector using an MFN object.
217: Input Parameters:
218: + st - spectral transformation context
219: - x - input vector
221: Output Parameter:
222: . y - output vector
223: */
224: PetscErrorCode STApply_Exp(ST st,Vec x,Vec y)
225: {
226: MFN mfn;
230: STShellGetContext(st,(void**)&mfn);
231: MFNSolve(mfn,x,y);
232: return(0);
233: }
235: /*TEST
237: testset:
238: args: -eps_nev 4 -mfn_ncv 16 -terse
239: requires: c99_complex !single
240: filter: sed -e "s/-2/+2/g"
241: output_file: output/ex36_1.out
242: test:
243: suffix: 1
244: requires: !__float128
245: test:
246: suffix: 1_quad
247: args: -eps_tol 1e-14
248: requires: __float128
250: test:
251: suffix: 2
252: args: -n 56 -eps_nev 2 -st_type sinvert -eps_target -390 -eps_target_magnitude -eps_type power
253: args: -eps_power_shift_type {{constant rayleigh}} -eps_two_sided {{0 1}} -eps_tol 1e-14 -terse
254: requires: c99_complex !single
255: filter: sed -e "s/[+-]0.0*i//"
256: output_file: output/ex36_2.out
258: TEST*/