Actual source code: ex17.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2010, Universidad 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            error, tol, re, im;
 38:   PetscScalar          kr, ki;
 40:   PetscInt             nev, maxit, i, its, nconv;
 41:   char                 filename[256];
 42:   PetscViewer          viewer;
 43:   PetscTruth           flg;

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

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

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

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

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

 74:   PetscOptionsGetString(PETSC_NULL,"-K",filename,256,&flg);
 75:   if (!flg) {
 76:     SETERRQ(1,"Must indicate a file name for matrix K with the -K option.");
 77:   }
 78:   PetscViewerBinaryOpen(PETSC_COMM_WORLD,filename,FILE_MODE_READ,&viewer);
 79:   MatLoad(viewer,MATAIJ,&K);
 80:   PetscViewerDestroy(viewer);

 82:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 83:                 Create the eigensolver and set various options
 84:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 86:   /* 
 87:      Create eigensolver context
 88:   */
 89:   QEPCreate(PETSC_COMM_WORLD,&qep);

 91:   /* 
 92:      Set matrices
 93:   */
 94:   QEPSetOperators(qep,M,C,K);

 96:   /*
 97:      Set solver parameters at runtime
 98:   */
 99:   QEPSetFromOptions(qep);

101:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
102:                       Solve the eigensystem
103:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

105:   QEPSolve(qep);
106:   QEPGetIterationNumber(qep, &its);
107:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

109:   /*
110:      Optional: Get some information from the solver and display it
111:   */
112:   QEPGetType(qep,&type);
113:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
114:   QEPGetDimensions(qep,&nev,PETSC_NULL,PETSC_NULL);
115:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
116:   QEPGetTolerances(qep,&tol,&maxit);
117:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
120:                     Display solution and clean up
121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

123:   /* 
124:      Get number of converged approximate eigenpairs
125:   */
126:   QEPGetConverged(qep,&nconv);
127:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);
128: 

130:   if (nconv>0) {
131:     /*
132:        Display eigenvalues and relative errors
133:     */
134:     PetscPrintf(PETSC_COMM_WORLD,
135:          "           k          ||(k^2M+Ck+K)x||/||kx||\n"
136:          "   ----------------- -------------------------\n" );

138:     for( i=0; i<nconv; i++ ) {
139:       /* 
140:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
141:         ki (imaginary part)
142:       */
143:       QEPGetEigenpair(qep,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
144:       /*
145:          Compute the relative error associated to each eigenpair
146:       */
147:       QEPComputeRelativeError(qep,i,&error);

149: #ifdef PETSC_USE_COMPLEX
150:       re = PetscRealPart(kr);
151:       im = PetscImaginaryPart(kr);
152: #else
153:       re = kr;
154:       im = ki;
155: #endif 
156:       if (im!=0.0) {
157:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j    %12g\n",re,im,error);
158:       } else {
159:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
160:       }
161:     }
162:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
163:   }
164: 
165:   /* 
166:      Free work space
167:   */
168:   QEPDestroy(qep);
169:   MatDestroy(M);
170:   MatDestroy(C);
171:   MatDestroy(K);
172:   SlepcFinalize();
173:   return 0;
174: }