Actual source code: ex17.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 quadratic eigenproblem (l^2*M + l*C + K)*x = 0 with matrices loaded from a file.\n\n"
 23:   "The command line options are:\n"
 24:   "  -M <filename>, where <filename> = matrix (M) file in PETSc binary form.\n"
 25:   "  -C <filename>, where <filename> = matrix (C) file in PETSc binary form.\n"
 26:   "  -K <filename>, where <filename> = matrix (K) file in PETSc binary form.\n\n";

 28: #include <slepcqep.h>

 32: int main(int argc,char **argv)
 33: {
 34:   Mat            M,C,K;           /* problem matrices */
 35:   QEP            qep;             /* quadratic eigenproblem solver context */
 36:   const QEPType  type;
 37:   PetscReal      tol;
 38:   PetscInt       nev,maxit,its;
 39:   char           filename[PETSC_MAX_PATH_LEN];
 40:   PetscViewer    viewer;
 41:   PetscBool      flg;

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

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 47:         Load the matrices that define the quadratic eigenproblem
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   PetscPrintf(PETSC_COMM_WORLD,"\nQuadratic eigenproblem stored in file.\n\n");
 51: #if defined(PETSC_USE_COMPLEX)
 52:   PetscPrintf(PETSC_COMM_WORLD," Reading COMPLEX matrices from binary files...\n");
 53: #else
 54:   PetscPrintf(PETSC_COMM_WORLD," Reading REAL matrices from binary files...\n");
 55: #endif

 57:   PetscOptionsGetString(PETSC_NULL,"-M",filename,PETSC_MAX_PATH_LEN,&flg);
 58:   if (!flg) {
 59:     SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix M with the -M option.");
 60:   }
 61:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 62:   MatCreate(PETSC_COMM_WORLD,&M);
 63:   MatSetFromOptions(M);
 64:   MatLoad(M,viewer);
 65:   PetscViewerDestroy(&viewer);

 67:   PetscOptionsGetString(PETSC_NULL,"-C",filename,PETSC_MAX_PATH_LEN,&flg);
 68:   if (!flg) {
 69:     SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix C with the -C option.");
 70:   }
 71:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 72:   MatCreate(PETSC_COMM_WORLD,&C);
 73:   MatSetFromOptions(C);
 74:   MatLoad(C,viewer);
 75:   PetscViewerDestroy(&viewer);

 77:   PetscOptionsGetString(PETSC_NULL,"-K",filename,PETSC_MAX_PATH_LEN,&flg);
 78:   if (!flg) {
 79:     SETERRQ(PETSC_COMM_WORLD,1,"Must indicate a file name for matrix K with the -K option.");
 80:   }
 81:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 82:   MatCreate(PETSC_COMM_WORLD,&K);
 83:   MatSetFromOptions(K);
 84:   MatLoad(K,viewer);
 85:   PetscViewerDestroy(&viewer);

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 88:                 Create the eigensolver and set various options
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 91:   /* 
 92:      Create eigensolver context
 93:   */
 94:   QEPCreate(PETSC_COMM_WORLD,&qep);

 96:   /* 
 97:      Set matrices
 98:   */
 99:   QEPSetOperators(qep,M,C,K);

101:   /*
102:      Set solver parameters at runtime
103:   */
104:   QEPSetFromOptions(qep);

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

110:   QEPSolve(qep);
111:   QEPGetIterationNumber(qep,&its);
112:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %D\n",its);

114:   /*
115:      Optional: Get some information from the solver and display it
116:   */
117:   QEPGetType(qep,&type);
118:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
119:   QEPGetDimensions(qep,&nev,PETSC_NULL,PETSC_NULL);
120:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %D\n",nev);
121:   QEPGetTolerances(qep,&tol,&maxit);
122:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4G, maxit=%D\n",tol,maxit);

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
125:                     Display solution and clean up
126:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

128:   QEPPrintSolution(qep,PETSC_NULL);
129:   QEPDestroy(&qep);
130:   MatDestroy(&M);
131:   MatDestroy(&C);
132:   MatDestroy(&K);
133:   SlepcFinalize();
134:   return 0;
135: }