Actual source code: ex3.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[] = "Solves the same eigenproblem as in example ex2, but using a shell matrix. "
 23:   "The problem is a standard symmetric eigenproblem corresponding to the 2-D Laplacian operator.\n\n"
 24:   "The command line options are:\n"
 25:   "  -n <n>, where <n> = number of grid subdivisions in both x and y dimensions.\n\n";

 27:  #include slepceps.h
 28: #include "petscblaslapack.h"

 30: /* 
 31:    User-defined routines
 32: */
 33: PetscErrorCode MatLaplacian2D_Mult( Mat A, Vec x, Vec y );

 37: int main( int argc, char **argv )
 38: {
 39:   Mat            A;               /* operator matrix */
 40:   EPS            eps;             /* eigenproblem solver context */
 41:   const EPSType  type;
 42:   PetscReal      error, tol, re, im;
 43:   PetscScalar    kr, ki;
 44:   PetscMPIInt    size;
 46:   PetscInt       N, n=10, nev, maxit, i, its, nconv;

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

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

 56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 57:      Compute the operator matrix that defines the eigensystem, Ax=kx
 58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 60:   MatCreateShell(PETSC_COMM_WORLD,N,N,N,N,&n,&A);
 61:   MatSetFromOptions(A);
 62:   MatShellSetOperation(A,MATOP_MULT,(void(*)())MatLaplacian2D_Mult);
 63:   MatShellSetOperation(A,MATOP_MULT_TRANSPOSE,(void(*)())MatLaplacian2D_Mult);

 65:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 66:                 Create the eigensolver and set various options
 67:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 69:   /* 
 70:      Create eigensolver context
 71:   */
 72:   EPSCreate(PETSC_COMM_WORLD,&eps);

 74:   /* 
 75:      Set operators. In this case, it is a standard eigenvalue problem
 76:   */
 77:   EPSSetOperators(eps,A,PETSC_NULL);
 78:   EPSSetProblemType(eps,EPS_HEP);

 80:   /*
 81:      Set solver parameters at runtime
 82:   */
 83:   EPSSetFromOptions(eps);

 85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 86:                       Solve the eigensystem
 87:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

 89:   EPSSolve(eps);
 90:   EPSGetIterationNumber(eps, &its);
 91:   PetscPrintf(PETSC_COMM_WORLD," Number of iterations of the method: %d\n",its);

 93:   /*
 94:      Optional: Get some information from the solver and display it
 95:   */
 96:   EPSGetType(eps,&type);
 97:   PetscPrintf(PETSC_COMM_WORLD," Solution method: %s\n\n",type);
 98:   EPSGetDimensions(eps,&nev,PETSC_NULL,PETSC_NULL);
 99:   PetscPrintf(PETSC_COMM_WORLD," Number of requested eigenvalues: %d\n",nev);
100:   EPSGetTolerances(eps,&tol,&maxit);
101:   PetscPrintf(PETSC_COMM_WORLD," Stopping condition: tol=%.4g, maxit=%d\n",tol,maxit);

103:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
104:                     Display solution and clean up
105:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

107:   /* 
108:      Get number of converged approximate eigenpairs
109:   */
110:   EPSGetConverged(eps,&nconv);
111:   PetscPrintf(PETSC_COMM_WORLD," Number of converged eigenpairs: %d\n\n",nconv);
112: 

114:   if (nconv>0) {
115:     /*
116:        Display eigenvalues and relative errors
117:     */
118:     PetscPrintf(PETSC_COMM_WORLD,
119:          "           k          ||Ax-kx||/||kx||\n"
120:          "   ----------------- ------------------\n" );

122:     for( i=0; i<nconv; i++ ) {
123:       /* 
124:         Get converged eigenpairs: i-th eigenvalue is stored in kr (real part) and
125:         ki (imaginary part)
126:       */
127:       EPSGetEigenpair(eps,i,&kr,&ki,PETSC_NULL,PETSC_NULL);
128:       /*
129:          Compute the relative error associated to each eigenpair
130:       */
131:       EPSComputeRelativeError(eps,i,&error);

133: #ifdef PETSC_USE_COMPLEX
134:       re = PetscRealPart(kr);
135:       im = PetscImaginaryPart(kr);
136: #else
137:       re = kr;
138:       im = ki;
139: #endif 
140:       if (im!=0.0) {
141:         PetscPrintf(PETSC_COMM_WORLD," %9f%+9f j %12g\n",re,im,error);
142:       } else {
143:         PetscPrintf(PETSC_COMM_WORLD,"   %12f       %12g\n",re,error);
144:       }
145:     }
146:     PetscPrintf(PETSC_COMM_WORLD,"\n" );
147:   }
148: 
149:   /* 
150:      Free work space
151:   */
152:   EPSDestroy(eps);
153:   MatDestroy(A);
154:   SlepcFinalize();
155:   return 0;
156: }

158: /*
159:     Compute the matrix vector multiplication y<---T*x where T is a nx by nx
160:     tridiagonal matrix with DD on the diagonal, DL on the subdiagonal, and 
161:     DU on the superdiagonal.
162:  */
163: static void tv( int nx, PetscScalar *x, PetscScalar *y )
164: {
165:   PetscScalar dd, dl, du;
166:   int         j;

168:   dd  = 4.0;
169:   dl  = -1.0;
170:   du  = -1.0;

172:   y[0] =  dd*x[0] + du*x[1];
173:   for( j=1; j<nx-1; j++ )
174:     y[j] = dl*x[j-1] + dd*x[j] + du*x[j+1];
175:   y[nx-1] = dl*x[nx-2] + dd*x[nx-1];
176: }

180: /*
181:     Matrix-vector product subroutine for the 2D Laplacian.

183:     The matrix used is the 2 dimensional discrete Laplacian on unit square with
184:     zero Dirichlet boundary condition.
185:  
186:     Computes y <-- A*x, where A is the block tridiagonal matrix
187:  
188:                  | T -I          | 
189:                  |-I  T -I       |
190:              A = |   -I  T       |
191:                  |        ...  -I|
192:                  |           -I T|
193:  
194:     The subroutine TV is called to compute y<--T*x.
195:  */
196: PetscErrorCode MatLaplacian2D_Mult( Mat A, Vec x, Vec y )
197: {
198:   void           *ctx;
200:   int            nx, lo, j, one=1;
201:   PetscScalar    *px, *py, dmone=-1.0;
202: 
204:   MatShellGetContext( A, &ctx );
205:   nx = *(int *)ctx;
206:   VecGetArray( x, &px );
207:   VecGetArray( y, &py );

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

212:   for( j=2; j<nx; j++ ) {
213:     lo = (j-1)*nx;
214:     tv( nx, &px[lo], &py[lo]);
215:     BLASaxpy_( &nx, &dmone, &px[lo-nx], &one, &py[lo], &one );
216:     BLASaxpy_( &nx, &dmone, &px[lo+nx], &one, &py[lo], &one );
217:   }

219:   lo = (nx-1)*nx;
220:   tv( nx, &px[lo], &py[lo]);
221:   BLASaxpy_( &nx, &dmone, &px[lo-nx], &one, &py[lo], &one );

223:   VecRestoreArray( x, &px );
224:   VecRestoreArray( y, &py );
225:   return(0);
226: }