Actual source code: ex7.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2011, Universitat Politecnica de Valencia, Spain

  6:    This file is part of SLEPc.
  7:       
  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:   "This example works for both real and complex numbers.\n\n"
 24:   "The command line options are:\n"
 25:   "  -f1 <filename>, where <filename> = matrix (A) file in PETSc binary form.\n"
 26:   "  -f2 <filename>, where <filename> = matrix (B) file in PETSc binary form.\n"
 27:   "  -evecs <filename>, output file to save computed eigenvectors.\n"
 28:   "  -ninitial <nini>, number of user-provided initial guesses.\n"
 29:   "  -finitial <filename>, where <filename> contains <nini> vectors (binary).\n"
 30:   "  -nconstr <ncon>, number of user-provided constraints.\n"
 31:   "  -fconstr <filename>, where <filename> contains <ncon> vectors (binary).\n\n";

 33: #include <slepceps.h>

 37: int main(int argc,char **argv)
 38: {
 39:   Mat            A,B;             /* matrices */
 40:   EPS            eps;             /* eigenproblem solver context */
 41:   const EPSType  type;
 42:   PetscReal      tol;
 43:   Vec            xr,xi,*is,*ds;
 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;

 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(PETSC_NULL,"-f1",filename,PETSC_MAX_PATH_LEN,&flg);
 58:   if (!flg) {
 59:     SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix A with the -f1 option.");
 60:   }

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

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

 85:   MatGetVecs(A,PETSC_NULL,&xr);
 86:   MatGetVecs(A,PETSC_NULL,&xi);

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

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

120:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
121:                 Create the eigensolver and set various options
122:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

134:   /* 
135:      If the user provided initial guesses or constraints, pass them here
136:   */
137:   EPSSetInitialSpace(eps,nini,is);
138:   EPSSetDeflationSpace(eps,ncon,ds);

140:   /*
141:      Set solver parameters at runtime
142:   */
143:   EPSSetFromOptions(eps);

145:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
146:                       Solve the eigensystem
147:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

149:   EPSSolve(eps);

151:   /*
152:      Optional: Get some information from the solver and display it
153:   */
154:   EPSGetIterationNumber(eps,&its);
155:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);
156:   EPSGetOperationCounters(eps,PETSC_NULL,PETSC_NULL,&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,PETSC_NULL,PETSC_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",tol,maxit);

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

169:   EPSPrintSolution(eps,PETSC_NULL);
170:   /*
171:      Save eigenvectors, if requested
172:   */
173:   PetscOptionsGetString(PETSC_NULL,"-evecs",filename,PETSC_MAX_PATH_LEN,&evecs);
174:   EPSGetConverged(eps,&nconv);
175:   if (nconv>0 && evecs) {
176:     PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_WRITE,&viewer);
177:     EPSIsHermitian(eps,&ishermitian);
178:     for (i=0;i<nconv;i++) {
179:       EPSGetEigenvector(eps,i,xr,xi);
180:       VecView(xr,viewer);
181:       if (!ishermitian) { VecView(xi,viewer); }
182:     }
183:     PetscViewerDestroy(&viewer);
184:   }
185: 
186:   /* 
187:      Free work space
188:   */
189:   EPSDestroy(&eps);
190:   MatDestroy(&A);
191:   MatDestroy(&B);
192:   VecDestroy(&xr);
193:   VecDestroy(&xi);
194:   if (nini>0) { VecDestroyVecs(nini,&is); }
195:   if (ncon>0) { VecDestroyVecs(ncon,&ds); }
196:   SlepcFinalize();
197:   return 0;
198: }