Actual source code: ex7.c

slepc-3.17.1 2022-04-11
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-, 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[] = "Solves a generalized eigensystem Ax=kBx with matrices loaded from a file.\n"
 12:   "The command line options are:\n"
 13:   "  -f1 <filename> -f2 <filename>, PETSc binary files containing A and B.\n"
 14:   "  -evecs <filename>, output file to save computed eigenvectors.\n"
 15:   "  -ninitial <nini>, number of user-provided initial guesses.\n"
 16:   "  -finitial <filename>, binary file containing <nini> vectors.\n"
 17:   "  -nconstr <ncon>, number of user-provided constraints.\n"
 18:   "  -fconstr <filename>, binary file containing <ncon> vectors.\n\n";

 20: #include <slepceps.h>

 22: int main(int argc,char **argv)
 23: {
 24:   Mat            A,B;             /* matrices */
 25:   EPS            eps;             /* eigenproblem solver context */
 26:   ST             st;
 27:   KSP            ksp;
 28:   EPSType        type;
 29:   PetscReal      tol;
 30:   Vec            xr,xi,*Iv,*Cv;
 31:   PetscInt       nev,maxit,i,its,lits,nconv,nini=0,ncon=0;
 32:   char           filename[PETSC_MAX_PATH_LEN];
 33:   PetscViewer    viewer;
 34:   PetscBool      flg,evecs,ishermitian,terse;

 36:   SlepcInitialize(&argc,&argv,(char*)0,help);

 38:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 39:         Load the matrices that define the eigensystem, Ax=kBx
 40:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 42:   PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n");
 43:   PetscOptionsGetString(NULL,NULL,"-f1",filename,sizeof(filename),&flg);

 46: #if defined(PETSC_USE_COMPLEX)
 47:   PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
 48: #else
 49:   PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
 50: #endif
 51:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 52:   MatCreate(PETSC_COMM_WORLD,&A);
 53:   MatSetFromOptions(A);
 54:   MatLoad(A,viewer);
 55:   PetscViewerDestroy(&viewer);

 57:   PetscOptionsGetString(NULL,NULL,"-f2",filename,sizeof(filename),&flg);
 58:   if (flg) {
 59:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 60:     MatCreate(PETSC_COMM_WORLD,&B);
 61:     MatSetFromOptions(B);
 62:     MatLoad(B,viewer);
 63:     PetscViewerDestroy(&viewer);
 64:   } else {
 65:     PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
 66:     B = NULL;
 67:   }

 69:   MatCreateVecs(A,NULL,&xr);
 70:   MatCreateVecs(A,NULL,&xi);

 72:   /*
 73:      Read user constraints if available
 74:   */
 75:   PetscOptionsGetInt(NULL,NULL,"-nconstr",&ncon,&flg);
 76:   if (flg) {
 78:     PetscOptionsGetString(NULL,NULL,"-fconstr",filename,sizeof(filename),&flg);
 80:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 81:     VecDuplicateVecs(xr,ncon,&Cv);
 82:     for (i=0;i<ncon;i++) VecLoad(Cv[i],viewer);
 83:     PetscViewerDestroy(&viewer);
 84:   }

 86:   /*
 87:      Read initial guesses if available
 88:   */
 89:   PetscOptionsGetInt(NULL,NULL,"-ninitial",&nini,&flg);
 90:   if (flg) {
 92:     PetscOptionsGetString(NULL,NULL,"-finitial",filename,sizeof(filename),&flg);
 94:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 95:     VecDuplicateVecs(xr,nini,&Iv);
 96:     for (i=0;i<nini;i++) VecLoad(Iv[i],viewer);
 97:     PetscViewerDestroy(&viewer);
 98:   }

100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101:                 Create the eigensolver and set various options
102:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

104:   /*
105:      Create eigensolver context
106:   */
107:   EPSCreate(PETSC_COMM_WORLD,&eps);

109:   /*
110:      Set operators. In this case, it is a generalized eigenvalue problem
111:   */
112:   EPSSetOperators(eps,A,B);

114:   /*
115:      If the user provided initial guesses or constraints, pass them here
116:   */
117:   EPSSetInitialSpace(eps,nini,Iv);
118:   EPSSetDeflationSpace(eps,ncon,Cv);

120:   /*
121:      Set solver parameters at runtime
122:   */
123:   EPSSetFromOptions(eps);

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:                       Solve the eigensystem
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

129:   EPSSolve(eps);

131:   /*
132:      Optional: Get some information from the solver and display it
133:   */
134:   EPSGetIterationNumber(eps,&its);
135:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %" PetscInt_FMT "\n",its);
136:   EPSGetST(eps,&st);
137:   STGetKSP(st,&ksp);
138:   KSPGetTotalIterations(ksp,&lits);
139:   PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %" PetscInt_FMT "\n",lits);
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: %" PetscInt_FMT "\n",nev);
144:   EPSGetTolerances(eps,&tol,&maxit);
145:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%" PetscInt_FMT "\n",(double)tol,maxit);

147:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
148:                     Display solution and clean up
149:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

151:   /*
152:      Show detailed info unless -terse option is given by user
153:    */
154:   PetscOptionsHasName(NULL,NULL,"-terse",&terse);
155:   if (terse) EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
156:   else {
157:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
158:     EPSConvergedReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
159:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
160:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
161:   }

163:   /*
164:      Save eigenvectors, if requested
165:   */
166:   PetscOptionsGetString(NULL,NULL,"-evecs",filename,sizeof(filename),&evecs);
167:   EPSGetConverged(eps,&nconv);
168:   if (nconv>0 && evecs) {
169:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
170:     EPSIsHermitian(eps,&ishermitian);
171:     for (i=0;i<nconv;i++) {
172:       EPSGetEigenvector(eps,i,xr,xi);
173:       VecView(xr,viewer);
174: #if !defined(PETSC_USE_COMPLEX)
175:       if (!ishermitian) VecView(xi,viewer);
176: #endif
177:     }
178:     PetscViewerDestroy(&viewer);
179:   }

181:   /*
182:      Free work space
183:   */
184:   EPSDestroy(&eps);
185:   MatDestroy(&A);
186:   MatDestroy(&B);
187:   VecDestroy(&xr);
188:   VecDestroy(&xi);
189:   if (nini > 0) VecDestroyVecs(nini,&Iv);
190:   if (ncon > 0) VecDestroyVecs(ncon,&Cv);
191:   SlepcFinalize();
192:   return 0;
193: }

195: /*TEST

197:    test:
198:       suffix: 1
199:       args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62b.petsc -eps_nev 4 -terse
200:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

202:    test:
203:       suffix: ciss_1
204:       args: -f1 ${DATAFILESPATH}/matrices/complex/mhd1280a.petsc -f2 ${DATAFILESPATH}/matrices/complex/mhd1280b.petsc -eps_type ciss -eps_ciss_usest 0 -eps_ciss_quadrule chebyshev -rg_type ring -rg_ring_center 0 -rg_ring_radius .49 -rg_ring_width 0.2 -rg_ring_startangle .25 -rg_ring_endangle .49 -terse -eps_max_it 1
205:       requires: double complex datafilespath !defined(PETSC_USE_64BIT_INDICES)
206:       timeoutfactor: 2

208:    test:
209:       suffix: 3  # test problem (A,A)
210:       args: -f1 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -f2 ${SLEPC_DIR}/share/slepc/datafiles/matrices/bfw62a.petsc -eps_nev 4 -terse
211:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

213: TEST*/