Actual source code: ex3.c

  2: static char help[] = "Solves the same eigenproblem as in example ex2, but using a shell matrix. "
  3:   "The problem is a standard symmetric eigenproblem corresponding to the 2-D Laplacian operator.\n\n"
  4:   "The command line options are:\n"
  5:   "  -n <n>, where <n> = number of grid subdivisions in both x and y dimensions.\n\n";

 7:  #include slepceps.h
  8: #include "petscblaslapack.h"

 10: /* 
 11:    User-defined routines
 12: */
 13: PetscErrorCode MatLaplacian2D_Mult( Mat A, Vec x, Vec y );

 17: int main( int argc, char **argv )
 18: {
 19:   Mat         A;               /* operator matrix */
 20:   EPS         eps;             /* eigenproblem solver context */
 21:   EPSType     type;
 22:   PetscReal   error, tol, re, im;
 23:   PetscScalar kr, ki;
 24:   PetscMPIInt size;
 26:   PetscInt    N, n=10;
 27:   int         nev, maxit, i, its, nconv;

 29:   SlepcInitialize(&argc,&argv,(char*)0,help);
 30:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 31:   if (size != 1) SETERRQ(1,"This is a uniprocessor example only!");

 33:   PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);
 34:   N = n*n;
 35:   PetscPrintf(PETSC_COMM_WORLD,"\n2-D Laplacian Eigenproblem (matrix-free version), N=%d (%dx%d grid)\n\n",N,n,n);

 37:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 38:      Compute the operator matrix that defines the eigensystem, Ax=kx
 39:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 41:   MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,&n,&A);
 42:   MatSetFromOptions(A);
 43:   MatShellSetOperation(A,MATOP_MULT,(void(*)())MatLaplacian2D_Mult);
 44:   MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)())MatLaplacian2D_Mult);

 46:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 47:                 Create the eigensolver and set various options
 48:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 50:   /* 
 51:      Create eigensolver context
 52:   */
 53:   EPSCreate(PETSC_COMM_WORLD,&eps);

 55:   /* 
 56:      Set operators. In this case, it is a standard eigenvalue problem
 57:   */
 58:   EPSSetOperators(eps,A,PETSC_NULL);
 59:   EPSSetProblemType(eps,EPS_HEP);

 61:   /*
 62:      Set solver parameters at runtime
 63:   */
 64:   EPSSetFromOptions(eps);

 66:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 67:                       Solve the eigensystem
 68:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 70:   EPSSolve(eps);
 71:   EPSGetIterationNumber(eps, &its);
 72:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

 74:   /*
 75:      Optional: Get some information from the solver and display it
 76:   */
 77:   EPSGetType(eps,&type);
 78:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 79:   EPSGetDimensions(eps,&nev,PETSC_NULL);
 80:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
 81:   EPSGetTolerances(eps,&tol,&maxit);
 82:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 85:                     Display solution and clean up
 86:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 88:   /* 
 89:      Get number of converged approximate eigenpairs
 90:   */
 91:   EPSGetConverged(eps,&nconv);
 92:   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);
 93: 

 95:   if (nconv>0) {
 96:     /*
 97:        Display eigenvalues and relative errors
 98:     */
 99:     PetscPrintf(PETSC_COMM_WORLD,
100:          "           k          ||Ax-kx||/||kx||\n"
101:          "   ----------------- ------------------\n" );

103:     for( i=0; i<nconv; i++ ) {
104:       /* 
105:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
106:         ki (imaginary part)
107:       */
108:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
109:       /*
110:          Compute the relative error associated to each eigenpair
111:       */
112:       EPSComputeRelativeError(eps,i,&error);

114: #ifdef PETSC_USE_COMPLEX
115:       re = PetscRealPart(kr);
116:       im = PetscImaginaryPart(kr);
117: #else
118:       re = kr;
119:       im = ki;
120: #endif 
121:       if (im!=0.0) {
122:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
123:       } else {
124:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
125:       }
126:     }
127:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
128:   }
129: 
130:   /* 
131:      Free work space
132:   */
133:   EPSDestroy(eps);
134:   MatDestroy(A);
135:   SlepcFinalize();
136:   return 0;
137: }

139: /*
140:     Compute the matrix vector multiplication y<---T*x where T is a nx by nx
141:     tridiagonal matrix with DD on the diagonal, DL on the subdiagonal, and 
142:     DU on the superdiagonal.
143:  */
144: static void tv( int nx, PetscScalar *x, PetscScalar *y )
145: {
146:   PetscScalar dd, dl, du;
147:   int         j;

149:   dd  = 4.0;
150:   dl  = -1.0;
151:   du  = -1.0;

153:   y[0] =  dd*x[0] + du*x[1];
154:   for( j=1; j<nx-1; j++ )
155:     y[j] = dl*x[j-1] + dd*x[j] + du*x[j+1];
156:   y[nx-1] = dl*x[nx-2] + dd*x[nx-1];
157: }

161: /*
162:     Matrix-vector product subroutine for the 2D Laplacian.

164:     The matrix used is the 2 dimensional discrete Laplacian on unit square with
165:     zero Dirichlet boundary condition.
166:  
167:     Computes y <-- A*x, where A is the block tridiagonal matrix
168:  
169:                  | T -I          | 
170:                  |-I  T -I       |
171:              A = |   -I  T       |
172:                  |        ...  -I|
173:                  |           -I T|
174:  
175:     The subroutine TV is called to compute y<--T*x.
176:  */
177: PetscErrorCode MatLaplacian2D_Mult( Mat A, Vec x, Vec y )
178: {
179:   void           *ctx;
181:   int            nx, lo, j, one=1;
182:   PetscScalar    *px, *py, dmone=-1.0;
183: 
185:   MatShellGetContext( A, &ctx );
186:   nx = *(int *)ctx;
187:   VecGetArray( x, &px );
188:   VecGetArray( y, &py );

190:   tv( nx, &px[0], &py[0] );
191:   BLASaxpy_( &nx, &dmone, &px[nx], &one, &py[0], &one );

193:   for( j=2; j<nx; j++ ) {
194:     lo = (j-1)*nx;
195:     tv( nx, &px[lo], &py[lo]);
196:     BLASaxpy_( &nx, &dmone, &px[lo-nx], &one, &py[lo], &one );
197:     BLASaxpy_( &nx, &dmone, &px[lo+nx], &one, &py[lo], &one );
198:   }

200:   lo = (nx-1)*nx;
201:   tv( nx, &px[lo], &py[lo]);
202:   BLASaxpy_( &nx, &dmone, &px[lo-nx], &one, &py[lo], &one );

204:   VecRestoreArray( x, &px );
205:   VecRestoreArray( y, &py );
206:   return(0);
207: }