Actual source code: ex7.c

slepc-3.6.1 2015-09-03
Report Typos and Errors
  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2015, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.

  8:    SLEPc is free software: you can redistribute it and/or modify it under  the
  9:    terms of version 3 of the GNU Lesser General Public License as published by
 10:    the Free Software Foundation.

 12:    SLEPc  is  distributed in the hope that it will be useful, but WITHOUT  ANY
 13:    WARRANTY;  without even the implied warranty of MERCHANTABILITY or  FITNESS
 14:    FOR  A  PARTICULAR PURPOSE. See the GNU Lesser General Public  License  for
 15:    more details.

 17:    You  should have received a copy of the GNU Lesser General  Public  License
 18:    along with SLEPc. If not, see <http://www.gnu.org/licenses/>.
 19:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 20: */

 22: static char help[] = "Solves a generalized eigensystem Ax=kBx with matrices loaded from a file.\n"
 23:   "The command line options are:\n"
 24:   "  -f1 <filename> -f2 <filename>, PETSc binary files containing A and B.\n"
 25:   "  -evecs <filename>, output file to save computed eigenvectors.\n"
 26:   "  -ninitial <nini>, number of user-provided initial guesses.\n"
 27:   "  -finitial <filename>, binary file containing <nini> vectors.\n"
 28:   "  -nconstr <ncon>, number of user-provided constraints.\n"
 29:   "  -fconstr <filename>, binary file containing <ncon> vectors.\n\n";

 31: #include <slepceps.h>

 35: int main(int argc,char **argv)
 36: {
 37:   Mat            A,B;             /* matrices */
 38:   EPS            eps;             /* eigenproblem solver context */
 39:   ST             st;
 40:   KSP            ksp;
 41:   EPSType        type;
 42:   PetscReal      tol;
 43:   Vec            xr,xi,*Iv,*Cv;
 44:   PetscInt       nev,maxit,i,its,lits,nconv,nini=0,ncon=0;
 45:   char           filename[PETSC_MAX_PATH_LEN];
 46:   PetscViewer    viewer;
 47:   PetscBool      flg,evecs,ishermitian,terse;

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

 52:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 53:         Load the matrices that define the eigensystem, Ax=kBx
 54:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 56:   PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized eigenproblem stored in file.\n\n");
 57:   PetscOptionsGetString(NULL,"-f1",filename,PETSC_MAX_PATH_LEN,&flg);
 58:   if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix A with the -f1 option");

 60: #if defined(PETSC_USE_COMPLEX)
 61:   PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
 62: #else
 63:   PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
 64: #endif
 65:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 66:   MatCreate(PETSC_COMM_WORLD,&A);
 67:   MatSetFromOptions(A);
 68:   MatLoad(A,viewer);
 69:   PetscViewerDestroy(&viewer);

 71:   PetscOptionsGetString(NULL,"-f2",filename,PETSC_MAX_PATH_LEN,&flg);
 72:   if (flg) {
 73:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 74:     MatCreate(PETSC_COMM_WORLD,&B);
 75:     MatSetFromOptions(B);
 76:     MatLoad(B,viewer);
 77:     PetscViewerDestroy(&viewer);
 78:   } else {
 79:     PetscPrintf(PETSC_COMM_WORLD," Matrix B was not provided, setting B=I\n\n");
 80:     B = NULL;
 81:   }

 83:   MatCreateVecs(A,NULL,&xr);
 84:   MatCreateVecs(A,NULL,&xi);

 86:   /*
 87:      Read user constraints if available
 88:   */
 89:   PetscOptionsGetInt(NULL,"-nconstr",&ncon,&flg);
 90:   if (flg) {
 91:     if (ncon<=0) SETERRQ(PETSC_COMM_WORLD,1,"The number of constraints must be >0");
 92:     PetscOptionsGetString(NULL,"-fconstr",filename,PETSC_MAX_PATH_LEN,&flg);
 93:     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must specify the name of the file storing the constraints");
 94:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 95:     VecDuplicateVecs(xr,ncon,&Cv);
 96:     for (i=0;i<ncon;i++) {
 97:       VecLoad(Cv[i],viewer);
 98:     }
 99:     PetscViewerDestroy(&viewer);
100:   }

102:   /*
103:      Read initial guesses if available
104:   */
105:   PetscOptionsGetInt(NULL,"-ninitial",&nini,&flg);
106:   if (flg) {
107:     if (nini<=0) SETERRQ(PETSC_COMM_WORLD,1,"The number of initial vectors must be >0");
108:     PetscOptionsGetString(NULL,"-finitial",filename,PETSC_MAX_PATH_LEN,&flg);
109:     if (!flg) SETERRQ(PETSC_COMM_WORLD,1,"Must specify the name of the file containing the initial vectors");
110:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
111:     VecDuplicateVecs(xr,nini,&Iv);
112:     for (i=0;i<nini;i++) {
113:       VecLoad(Iv[i],viewer);
114:     }
115:     PetscViewerDestroy(&viewer);
116:   }

118:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119:                 Create the eigensolver and set various options
120:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

122:   /*
123:      Create eigensolver context
124:   */
125:   EPSCreate(PETSC_COMM_WORLD,&eps);

127:   /*
128:      Set operators. In this case, it is a generalized eigenvalue problem
129:   */
130:   EPSSetOperators(eps,A,B);

132:   /*
133:      If the user provided initial guesses or constraints, pass them here
134:   */
135:   EPSSetInitialSpace(eps,nini,Iv);
136:   EPSSetDeflationSpace(eps,ncon,Cv);

138:   /*
139:      Set solver parameters at runtime
140:   */
141:   EPSSetFromOptions(eps);

143:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
144:                       Solve the eigensystem
145:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

147:   EPSSolve(eps);

149:   /*
150:      Optional: Get some information from the solver and display it
151:   */
152:   EPSGetIterationNumber(eps,&its);
153:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
154:   EPSGetST(eps,&st);
155:   STGetKSP(st,&ksp);
156:   KSPGetTotalIterations(ksp,&lits);
157:   PetscPrintf(PETSC_COMM_WORLD," Number of linear iterations of the method: %D\n",lits);
158:   EPSGetType(eps,&type);
159:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
160:   EPSGetDimensions(eps,&nev,NULL,NULL);
161:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
162:   EPSGetTolerances(eps,&tol,&maxit);
163:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%D\n",(double)tol,maxit);

165:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
166:                     Display solution and clean up
167:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

169:   /*
170:      Show detailed info unless -terse option is given by user
171:    */
172:   PetscOptionsHasName(NULL,"-terse",&terse);
173:   if (terse) {
174:     EPSErrorView(eps,EPS_ERROR_RELATIVE,NULL);
175:   } else {
176:     PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_ASCII_INFO_DETAIL);
177:     EPSReasonView(eps,PETSC_VIEWER_STDOUT_WORLD);
178:     EPSErrorView(eps,EPS_ERROR_RELATIVE,PETSC_VIEWER_STDOUT_WORLD);
179:     PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD);
180:   }

182:   /*
183:      Save eigenvectors, if requested
184:   */
185:   PetscOptionsGetString(NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
186:   EPSGetConverged(eps,&nconv);
187:   if (nconv>0 && evecs) {
188:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
189:     EPSIsHermitian(eps,&ishermitian);
190:     for (i=0;i<nconv;i++) {
191:       EPSGetEigenvector(eps,i,xr,xi);
192:       VecView(xr,viewer);
193: #if !defined(PETSC_USE_COMPLEX)
194:       if (!ishermitian) { VecView(xi,viewer); }
195: #endif
196:     }
197:     PetscViewerDestroy(&viewer);
198:   }

200:   /*
201:      Free work space
202:   */
203:   EPSDestroy(&eps);
204:   MatDestroy(&A);
205:   MatDestroy(&B);
206:   VecDestroy(&xr);
207:   VecDestroy(&xi);
208:   if (nini > 0) {
209:     VecDestroyVecs(nini,&Iv);
210:   }
211:   if (ncon > 0) {
212:     VecDestroyVecs(ncon,&Cv);
213:   }
214:   SlepcFinalize();
215:   return 0;
216: }