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 );
 34: PetscErrorCode MatLaplacian2D_GetDiagonal( Mat A, Vec diag );

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

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

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

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

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

 67:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 68:                 Create the eigensolver and set various options
 69:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

 82:   /*
 83:      Set solver parameters at runtime
 84:   */
 85:   EPSSetFromOptions(eps);

 87:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 88:                       Solve the eigensystem
 89:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

105:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
106:                     Display solution and clean up
107:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

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

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

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

170:   dd  = 4.0;
171:   dl  = -1.0;
172:   du  = -1.0;

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

182: /*
183:     Matrix-vector product subroutine for the 2D Laplacian.

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

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

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

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

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

232: PetscErrorCode MatLaplacian2D_GetDiagonal( Mat A, Vec diag )
233: {

237:   VecSet(diag,4.0);
238:   return(0);
239: }