Actual source code: ex36.c

slepc-3.11.2 2019-07-30
Report Typos and Errors
  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*/