Actual source code: ex36.c

slepc-3.16.1 2021-11-17
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2021, 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;

 55:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 56: #if defined(PETSC_HAVE_COMPLEX)
 57:   PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 58:   PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model with matrix exponential, n=%D\n\n",n);

 60:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 61:         Generate the matrix
 62:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 64:   alpha  = 2.0;
 65:   beta   = 5.45;
 66:   delta1 = 0.008;
 67:   delta2 = 0.004;
 68:   L      = 0.51302;

 70:   PetscOptionsGetScalar(NULL,NULL,"-L",&L,NULL);
 71:   PetscOptionsGetScalar(NULL,NULL,"-alpha",&alpha,NULL);
 72:   PetscOptionsGetScalar(NULL,NULL,"-beta",&beta,NULL);
 73:   PetscOptionsGetScalar(NULL,NULL,"-delta1",&delta1,NULL);
 74:   PetscOptionsGetScalar(NULL,NULL,"-delta2",&delta2,NULL);

 76:   h = 1.0 / (PetscReal)(n+1);
 77:   tau1 = delta1 / ((h*L)*(h*L));
 78:   tau2 = delta2 / ((h*L)*(h*L));

 80:   MatCreate(PETSC_COMM_WORLD,&A);
 81:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n);
 82:   MatSetFromOptions(A);
 83:   MatSetUp(A);

 85:   MatGetOwnershipRange(A,&Istart,&Iend);
 86:   for (i=Istart;i<Iend;i++) {
 87:     if (i<n) {  /* upper blocks */
 88:       if (i>0) { MatSetValue(A,i,i-1,tau1,INSERT_VALUES); }
 89:       if (i<n-1) { MatSetValue(A,i,i+1,tau1,INSERT_VALUES); }
 90:       MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES);
 91:       MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES);
 92:     } else {  /* lower blocks */
 93:       if (i>n) { MatSetValue(A,i,i-1,tau2,INSERT_VALUES); }
 94:       if (i<2*n-1) { MatSetValue(A,i,i+1,tau2,INSERT_VALUES); }
 95:       MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES);
 96:       MatSetValue(A,i,i-n,-beta,INSERT_VALUES);
 97:     }
 98:   }
 99:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
100:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

102:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
103:                 Create the eigensolver and set various options
104:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

106:   EPSCreate(PETSC_COMM_WORLD,&eps);
107:   EPSSetOperators(eps,A,NULL);
108:   EPSSetProblemType(eps,EPS_NHEP);
109:   EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
110:   EPSGetST(eps,&st);
111:   STSetType(st,STSHELL);
112:   EPSSetFromOptions(eps);

114:   /*
115:      Initialize shell spectral transformation
116:   */
117:   PetscObjectTypeCompare((PetscObject)st,STSHELL,&isShell);
118:   if (isShell) {

120:     /* Create the MFN object to be used by the spectral transform */
121:     MFNCreate(PETSC_COMM_WORLD,&mfn);
122:     MFNSetOperator(mfn,A);
123:     MFNGetFN(mfn,&f);
124:     FNSetType(f,FNEXP);
125:     FNSetScale(f,0.03,1.0);  /* this can be set with -fn_scale */
126:     MFNSetFromOptions(mfn);

128:     /* Set callback functions */
129:     STShellSetApply(st,STApply_Exp);
130:     STShellSetBackTransform(st,STBackTransform_Exp);
131:     STShellSetContext(st,mfn);
132:     PetscObjectSetName((PetscObject)st,"STEXP");
133:   }

135:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
136:                       Solve the eigensystem
137:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

139:   EPSSolve(eps);
140:   EPSGetType(eps,&type);
141:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
142:   EPSGetDimensions(eps,&nev,NULL,NULL);
143:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146:                     Display solution and clean up
147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

149:   /* show detailed info unless -terse option is given by user */
150:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
151:   if (terse) {
152:     EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
153:   } else {
154:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
155:     EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
156:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
157:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
158:   }
159:   EPSDestroy(&eps);
160:   MatDestroy(&A);
161:   if (isShell) { MFNDestroy(&mfn); }
162: #else
163:   SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"This example requires C99 complex numbers");
164: #endif
165:   SlepcFinalize();
166:   return ierr;
167: }

169: /* ------------------------------------------------------------------- */
170: /*
171:    STBackTransform_Exp - Undoes the exp(A) transformation by taking logarithms.

173:    Input Parameters:
174: +  st - spectral transformation context
175: -  n  - number of eigenvalues to transform

177:    Input/Output Parameters:
178: +  eigr - pointer to real part of eigenvalues
179: -  eigi - pointer to imaginary part of eigenvalues
180: */
181: PetscErrorCode STBackTransform_Exp(ST st,PetscInt n,PetscScalar *eigr,PetscScalar *eigi)
182: {
183: #if defined(PETSC_HAVE_COMPLEX)
185:   PetscInt       j;
186:   MFN            mfn;
187:   FN             fn;
188:   PetscScalar    tau,eta;
189: #if !defined(PETSC_USE_COMPLEX)
190:   PetscComplex   theta,lambda;
191: #endif

194:   STShellGetContext(st,&mfn);
195:   MFNGetFN(mfn,&fn);
196:   FNGetScale(fn,&tau,&eta);
197:   for (j=0;j<n;j++) {
198: #if defined(PETSC_USE_COMPLEX)
199:     eigr[j] = PetscLogComplex(eigr[j]/eta)/tau;
200: #else
201:     theta   = PetscCMPLX(eigr[j],eigi[j])/eta;
202:     lambda  = PetscLogComplex(theta)/tau;
203:     eigr[j] = PetscRealPartComplex(lambda);
204:     eigi[j] = PetscImaginaryPartComplex(lambda);
205: #endif
206:   }
207:   return(0);
208: #else
209:   return 0;
210: #endif
211: }

213: /*
214:    STApply_Exp - Applies the operator exp(tau*A) to a given vector using an MFN object.

216:    Input Parameters:
217: +  st - spectral transformation context
218: -  x  - input vector

220:    Output Parameter:
221: .  y - output vector
222: */
223: PetscErrorCode STApply_Exp(ST st,Vec x,Vec y)
224: {
225:   MFN            mfn;

229:   STShellGetContext(st,&mfn);
230:   MFNSolve(mfn,x,y);
231:   return(0);
232: }

234: /*TEST

236:    testset:
237:       args: -eps_nev 4 -mfn_ncv 16 -terse
238:       requires: c99_complex !single
239:       filter: sed -e "s/-2/+2/g"
240:       output_file: output/ex36_1.out
241:       test:
242:          suffix: 1
243:          requires: !__float128
244:       test:
245:          suffix: 1_quad
246:          args: -eps_tol 1e-14
247:          requires: __float128

249:    test:
250:       suffix: 2
251:       args: -n 56 -eps_nev 2 -st_type sinvert -eps_target -390 -eps_target_magnitude -eps_type power
252:       args: -eps_power_shift_type {{constant rayleigh}} -eps_two_sided {{0 1}} -eps_tol 1e-14 -terse
253:       requires: c99_complex !single
254:       filter: sed -e "s/[+-]0\.0*i//g"

256:    test:
257:       suffix: 3
258:       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
259:       requires: c99_complex !single
260:       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/"

262: TEST*/