Actual source code: ex44.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[] = "Compute rightmost eigenvalues with Lyapunov inverse iteration.\n\n"
 12:   "Loads matrix from a file or builds the same problem as ex36.c (with fixed parameters).\n\n"
 13:   "The command line options are:\n"
 14:   "  -file <filename>, where <filename> = matrix file in PETSc binary form.\n"
 15:   "  -shift <sigma>, shift to make the matrix stable.\n"
 16:   "  -n <n>, block dimension of the 2x2 block matrix (if matrix is generated).\n\n";

 18: #include <slepceps.h>

 20: int main(int argc,char **argv)
 21: {
 22:   Mat            A;               /* operator matrix */
 23:   EPS            eps;             /* eigenproblem solver context */
 24:   EPSType        type;
 25:   PetscScalar    alpha,beta,tau1,tau2,delta1,delta2,L,h,sigma=0.0;
 26:   PetscInt       n=30,i,Istart,Iend,nev;
 27:   char           filename[PETSC_MAX_PATH_LEN];
 28:   PetscViewer    viewer;
 29:   PetscBool      flg,terse;

 32:   SlepcInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;

 34:   PetscOptionsGetString(NULL,NULL,"-file",filename,sizeof(filename),&flg);
 35:   if (flg) {

 37:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 38:                         Load the matrix from file
 39:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 40:     PetscPrintf(PETSC_COMM_WORLD,"\nEigenproblem stored in file.\n\n");
 41: #if defined(PETSC_USE_COMPLEX)
 42:     PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrix from a binary file...\n");
 43: #else
 44:     PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrix from a binary file...\n");
 45: #endif
 46:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 47:     MatCreate(PETSC_COMM_WORLD,&A);
 48:     MatSetFromOptions(A);
 49:     MatLoad(A,viewer);
 50:     PetscViewerDestroy(&viewer);

 52:   } else {

 54:     /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 55:           Generate Brusselator matrix
 56:        - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 57:     PetscOptionsGetInt(NULL,NULL,"-n",&n,NULL);
 58:     PetscPrintf(PETSC_COMM_WORLD,"\nBrusselator wave model, n=%D\n\n",n);

 60:     alpha  = 2.0;
 61:     beta   = 5.45;
 62:     delta1 = 0.008;
 63:     delta2 = 0.004;
 64:     L      = 0.51302;

 66:     h = 1.0 / (PetscReal)(n+1);
 67:     tau1 = delta1 / ((h*L)*(h*L));
 68:     tau2 = delta2 / ((h*L)*(h*L));

 70:     MatCreate(PETSC_COMM_WORLD,&A);
 71:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,2*n,2*n);
 72:     MatSetFromOptions(A);
 73:     MatSetUp(A);

 75:     MatGetOwnershipRange(A,&Istart,&Iend);
 76:     for (i=Istart;i<Iend;i++) {
 77:       if (i<n) {  /* upper blocks */
 78:         if (i>0) { MatSetValue(A,i,i-1,tau1,INSERT_VALUES); }
 79:         if (i<n-1) { MatSetValue(A,i,i+1,tau1,INSERT_VALUES); }
 80:         MatSetValue(A,i,i,-2.0*tau1+beta-1.0,INSERT_VALUES);
 81:         MatSetValue(A,i,i+n,alpha*alpha,INSERT_VALUES);
 82:       } else {  /* lower blocks */
 83:         if (i>n) { MatSetValue(A,i,i-1,tau2,INSERT_VALUES); }
 84:         if (i<2*n-1) { MatSetValue(A,i,i+1,tau2,INSERT_VALUES); }
 85:         MatSetValue(A,i,i,-2.0*tau2-alpha*alpha,INSERT_VALUES);
 86:         MatSetValue(A,i,i-n,-beta,INSERT_VALUES);
 87:       }
 88:     }
 89:     MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 90:     MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 91:   }

 93:   /* Shift the matrix to make it stable, A-sigma*I */
 94:   PetscOptionsGetScalar(NULL,NULL,"-shift",&sigma,NULL);
 95:   MatShift(A,-sigma);

 97:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 98:                 Create the eigensolver and set various options
 99:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

101:   EPSCreate(PETSC_COMM_WORLD,&eps);
102:   EPSSetOperators(eps,A,NULL);
103:   EPSSetProblemType(eps,EPS_NHEP);
104:   EPSSetWhichEigenpairs(eps,EPS_LARGEST_REAL);
105:   EPSSetFromOptions(eps);

107:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
108:                       Solve the eigensystem
109:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

111:   EPSSolve(eps);
112:   EPSGetType(eps,&type);
113:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
114:   EPSGetDimensions(eps,&nev,NULL,NULL);
115:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);

117:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
118:                     Display solution and clean up
119:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

121:   /* show detailed info unless -terse option is given by user */
122:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
123:   if (terse) {
124:     EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
125:   } else {
126:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
127:     EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
128:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
129:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
130:   }
131:   EPSDestroy(&eps);
132:   MatDestroy(&A);
133:   SlepcFinalize();
134:   return ierr;
135: }

137: /*TEST

139:    testset:
140:       args: -eps_nev 6 -shift 0.1 -eps_type {{krylovschur lyapii}} -eps_tol 1e-7 -terse
141:       requires: double
142:       filter: grep -v method | sed -e "s/-0.09981-2.13938i, -0.09981+2.13938i/-0.09981+2.13938i, -0.09981-2.13938i/" | sed -e "s/-0.77192-2.52712i, -0.77192+2.52712i/-0.77192+2.52712i, -0.77192-2.52712i/" | sed -e "s/-1.88445-3.02666i, -1.88445+3.02666i/-1.88445+3.02666i, -1.88445-3.02666i/"
143:       output_file: output/ex44_1.out
144:       test:
145:          suffix: 1
146:       test:
147:          suffix: 2
148:          args: -eps_lyapii_ranks 8,20 -options_left no

150:    test:
151:       suffix: 2_feast
152:       args: -eps_type feast -eps_interval -103,-90 -terse
153:       requires: feast !single

155: TEST*/