Actual source code: ex13.c

  1: /*
  2:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  3:    SLEPc - Scalable Library for Eigenvalue Problem Computations
  4:    Copyright (c) 2002-2009, 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:   const EPSType  type;
 40:   PetscReal            error, tol, re, im;
 41:   PetscScalar          kr, ki;
 43:   PetscInt             N, n=10, m, Istart, Iend, II, J, nev, maxit, i, j, its, nconv, nulldim=0;
 44:   PetscScalar          v;
 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==PETSC_FALSE ) 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:     v = -1.0; i = II/n; j = II-i*n;
 71:     if(i>0) { J=II-n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 72:     if(i<m-1) { J=II+n; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 73:     if(j>0) { J=II-1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 74:     if(j<n-1) { J=II+1; MatSetValues(A,1,&II,1,&J,&v,INSERT_VALUES); }
 75:     v=4.0; MatSetValues(A,1,&II,1,&II,&v,INSERT_VALUES);
 76:     if (II>=nulldim) { v=4.0; MatSetValues(B,1,&II,1,&II,&v,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:      Set solver parameters at runtime
101:   */
102:   EPSSetFromOptions(eps);

104:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
105:                       Solve the eigensystem
106:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

108:   EPSSolve(eps);
109:   EPSGetIterationNumber(eps, &its);
110:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

112:   /*
113:      Optional: Get some information from the solver and display it
114:   */
115:   EPSGetType(eps,&type);
116:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
117:   EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
118:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
119:   EPSGetTolerances(eps,&tol,&maxit);
120:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

122:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
123:                     Display solution and clean up
124:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

126:   /* 
127:      Get number of converged approximate eigenpairs
128:   */
129:   EPSGetConverged(eps,&nconv);
130:   PetscPrintf(PETSC_COMM_WORLD," Number of converged approximate eigenpairs: %d\n\n",nconv);
131: 

133:   if (nconv>0) {
134:     /*
135:        Display eigenvalues and relative errors
136:     */
137:     PetscPrintf(PETSC_COMM_WORLD,
138:          "           k          ||Ax-kBx||/||kx||\n"
139:          "   ----------------- ------------------\n" );

141:     for( i=0; i<nconv; i++ ) {
142:       /* 
143:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
144:         ki (imaginary part)
145:       */
146:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
147:       /*
148:          Compute the relative error associated to each eigenpair
149:       */
150:       EPSComputeRelativeError(eps,i,&error);

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