Actual source code: ex13.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[] = "Generalized Symmetric eigenproblem.\n\n"
 23:   "The problem is Ax = lambda Bx, with:\n"
 24:   "   A = Laplacian operator in 2-D\n"
 25:   "   B = diagonal matrix with all values equal to 4 except nulldim zeros\n\n"
 26:   "The command line options are:\n"
 27:   "  -n <n>, where <n> = number of grid subdivisions in x dimension.\n"
 28:   "  -m <m>, where <m> = number of grid subdivisions in y dimension.\n"
 29:   "  -nulldim <k>, where <k> = dimension of the nullspace of B.\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;                  /* spectral transformation context */
 40:   const EPSType  type;
 41:   PetscReal            error, tol, re, im;
 42:   PetscScalar          kr, ki;
 44:   PetscInt             N, n=10, m, Istart, Iend, II, nev, maxit, i, j, its, nconv, nulldim=0;
 45:   PetscTruth           flag;

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

 49:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 50:   PetscOptionsGetInt(PETSC_NULL,"-m",&m,&flag);
 51:   if(!flag) m=n;
 52:   N = n*m;
 53:   PetscOptionsGetInt(PETSC_NULL,"-nulldim",&nulldim,PETSC_NULL);
 54:   PetscPrintf(PETSC_COMM_WORLD,"\nGeneralized Symmetric Eigenproblem, N=%d (%dx%d grid), null(B)=%d\n\n",N,n,m,nulldim);

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 57:      Compute the matrices that define the eigensystem, Ax=kBx
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   MatCreate(PETSC_COMM_WORLD,&A);
 61:   MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,N,N);
 62:   MatSetFromOptions(A);
 63: 
 64:   MatCreate(PETSC_COMM_WORLD,&B);
 65:   MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,N,N);
 66:   MatSetFromOptions(B);

 68:   MatGetOwnershipRange(A,&Istart,&Iend);
 69:   for( II=Istart; II<Iend; II++ ) {
 70:     i = II/n; j = II-i*n;
 71:     if(i>0) { MatSetValue(A,II,II-n,-1.0,INSERT_VALUES); }
 72:     if(i<m-1) { MatSetValue(A,II,II+n,-1.0,INSERT_VALUES); }
 73:     if(j>0) { MatSetValue(A,II,II-1,-1.0,INSERT_VALUES); }
 74:     if(j<n-1) { MatSetValue(A,II,II+1,-1.0,INSERT_VALUES); }
 75:     MatSetValue(A,II,II,4.0,INSERT_VALUES);
 76:     if (II>=nulldim) { MatSetValue(B,II,II,4.0,INSERT_VALUES); }
 77:   }

 79:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
 80:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
 81:   MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);
 82:   MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 85:                 Create the eigensolver and set various options
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   /* 
 89:      Create eigensolver context
 90:   */
 91:   EPSCreate(PETSC_COMM_WORLD,&eps);

 93:   /* 
 94:      Set operators. In this case, it is a generalized eigenvalue problem
 95:   */
 96:   EPSSetOperators(eps,A,B);
 97:   EPSSetProblemType(eps,EPS_GHEP);

 99:   /* 
100:      Use shift-and-invert to avoid solving linear systems with a singular B
101:      in case nulldim>0
102:   */
103:   EPSGetST(eps,&st);
104:   EPSSetTarget(eps,0.0);
105:   EPSSetWhichEigenpairs(eps,EPS_TARGET_MAGNITUDE);
106:   STSetType(st,STSINVERT);

108:   /*
109:      Set solver parameters at runtime
110:   */
111:   EPSSetFromOptions(eps);

113:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
114:                       Solve the eigensystem
115:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

117:   EPSSolve(eps);
118:   EPSGetIterationNumber(eps, &its);
119:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

121:   /*
122:      Optional: Get some information from the solver and display it
123:   */
124:   EPSGetType(eps,&type);
125:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
126:   EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
127:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
128:   EPSGetTolerances(eps,&tol,&maxit);
129:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

131:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
132:                     Display solution and clean up
133:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

135:   /* 
136:      Get number of converged approximate eigenpairs
137:   */
138:   EPSGetConverged(eps,&nconv);
139:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);
140: 

142:   if (nconv>0) {
143:     /*
144:        Display eigenvalues and relative errors
145:     */
146:     PetscPrintf(PETSC_COMM_WORLD,
147:          "           k          ||Ax-kBx||/||kx||\n"
148:          "   ----------------- ------------------\n" );

150:     for( i=0; i<nconv; i++ ) {
151:       /* 
152:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
153:         ki (imaginary part)
154:       */
155:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
156:       /*
157:          Compute the relative error associated to each eigenpair
158:       */
159:       EPSComputeRelativeError(eps,i,&error);

161: #ifdef PETSC_USE_COMPLEX
162:       re = PetscRealPart(kr);
163:       im = PetscImaginaryPart(kr);
164: #else
165:       re = kr;
166:       im = ki;
167: #endif 
168:       if (im!=0.0) {
169:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
170:       } else {
171:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
172:       }
173:     }
174:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
175:   }
176: 
177:   /* 
178:      Free work space
179:   */
180:   EPSDestroy(eps);
181:   MatDestroy(A);
182:   MatDestroy(B);
183:   SlepcFinalize();
184:   return 0;
185: }