      SUBROUTINE MATPDE( N, NX, NY, A, ROWPTR, COLIND, RHS )              
*     ..
*     .. Scalar Arguments ..
      INTEGER            N, NX, NY
*     .. 
*     .. Array Arguments .. 
      INTEGER            ROWPTR( * ), COLIND( * ) 
      DOUBLE PRECISION   A( * ), RHS( * )
*
*  Purpose
*  =======
*
*  Forms the five-point central differences operator for the elliptic 
*  PDE                    
*
*   -(P Ux)x -(Q Uy)y  +  R Ux + (R U)x  +  S Uy + (S U)y + T U = F
* 
*  where P, Q, R, S and T are the functions of x and y.  The domain 
*  is the unit square (0,1)x(0,1), and the boundary condition are 
*  Dirichlet.                                            
*
*  The matrix is a block tridiagonal matrix, there are NY blocks of
*  size NX by NX on the diagonal (each block is a tridiagonal matrix),
*  and then NY-1 blocks of size NX by NX on the sub- and super-block
*  diagonal positions.  
*
*  Important note: the matrix A is stored in compressed ROW format. 
*
*  Arguments
*  =========
*
*  N       (output) INTEGER 
*          The order of matrix A.  N = NX*NY 
*
*  NX      (input) INTEGER 
*          The number of mesh points in the X-axis
* 
*  NY      (input) INTEGER 
*          The number of mesh points in the Y-axis
*
*  A       (output) DOUBLE PRECISION array, dimension (NZ)  
*          The numerical values of the nonzero elements of matrix. 
*          where NZ is the number of nonzero elements in the generated
*          matrix. Specifically, NZ = NX*NY + (NX-1)*NY*2 + NX*(NY-1)*2. 
*
*  ROWPTR  (output) INTEGER array, dimension (N+1)  
*          The row start pointers.  
*
*  COLIND  (output) INTEGER array, dimension (NZ)  
*          the column indices.
*
*  RHS     (output) DOUBLE PRECISION array, dimension (N)  
*          The right hand side b of the linear system Ax = b. 
*
*
*  Internal parameters 
*  ===================
*
*  BETA, GAMMA are the coefficients of the elliptic operator,
*  BNDRY( 4 ) are the boundary Dirichlet conditions. These parameters
*  may be changed, which will result in different matrices and spectral
*  distribution. 
*                                                                       
*  ==================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, HALF
      PARAMETER	         ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 ) 
      DOUBLE PRECISION   BETA, GAMMA
      PARAMETER          ( BETA = 2.0D+1, GAMMA = 0.0D+0 )
*     .. 
*     .. Scalar variables .. 
      INTEGER		 I, J, IX, JY, INDEX, JD, IAI 
      IMPLICIT REAL*8    (A-H,O-Z) 
*
*     .. Array variables .. 
      DOUBLE PRECISION   BNDRY( 4 ) 
*     ..
*     .. External functions .. 
      DOUBLE PRECISION   F, PC, QC, RC, SC, U
      EXTERNAL           F, PC, QC, RC, SC, U
*     ..
*     .. Intrinsic functions .. 
      INTRINSIC	         MOD, DBLE
*                                                                       
*     set up initial boundary conditions 
*
      BNDRY( 1 ) = ZERO  
      BNDRY( 2 ) = ZERO   
      BNDRY( 3 ) = ZERO   
      BNDRY( 4 ) = ZERO   
*
*     set up other parameters 
*
      N = NX*NY 
      HX = ONE / (DBLE(NX) + ONE) 
      HY = ONE / (DBLE(NY) + ONE)
      HY2 = HY*HY                                                     
      RA = (HY/HX)**2                                                 
      RB  = HY2/HX                                                    
*                                                                       
*     FORM ROWPTR.                                                             
*
      ROWPTR(1) = 1                                                       
      DO 10 I = 1,N                                                     
         IAI = 1                                                      
         IF (I.GT.NX) 
     $      IAI = IAI+1                                      
         IF(MOD(I-1,NX).NE.0) 
     $     IAI = IAI+1                             
         IF(MOD(I,NX).NE.0) 
     $     IAI = IAI+1                               
         IF(I.LE.N-NX) 
     $     IAI = IAI+1                                    
         ROWPTR(I+1) = ROWPTR(I) + IAI                                        
 10   CONTINUE                                                        
*                                                                       
      J = 1                                                           
      DO 110 JY=1,NY                                                  
         YJ = HY*DBLE(JY)
         DO 100 IX=1,NX                                               
*
            XI = HX*DBLE(IX) 
            INDEX = IX + (JY-1)*NX                                    
            RHS(INDEX) = HY2*F(BETA,GAMMA,XI,YJ) 
            P12  = PC (XI + HALF*HX,YJ) 
            PM12 = PC (XI - HALF*HX,YJ)  
            Q12  = QC (XI,YJ + HALF*HY)
            QM12 = QC (XI,YJ - HALF*HY) 
            RIJ  = RC (BETA,XI,YJ)                                         
            R1   = RC (BETA,XI + HX,YJ)                                    
            RM1  = RC (BETA,XI - HX,YJ)                                    
            SIJ  = SC (GAMMA,XI,YJ)                                         
            S1   = SC (GAMMA,XI,YJ + HY)                                    
            SM1  = SC (GAMMA,XI,YJ - HY)                                    
*                                                                       
*           DIAGONAL.   
*
            A(J) = RA*(P12 + PM12)  +  Q12 + QM12  + HY2*TC(XI,YJ)    
            COLIND(J) = INDEX                                             
            JD = J                                                    
            J = J+1                                                   
*                                                                       
*           LOWEST BAND. 
*
            IF (JY.EQ.1) 
     $         GOTO 20                                      
            A(J) = - ( QM12 + HALF*HY*(SIJ+SM1) )
            COLIND(J) = INDEX - NX                                     
            J = J+1                                                
*                                                                       
*           SUB-DIAGONAL.  
*
 20         IF (IX.EQ.1) 
     $         GOTO 30                                      
            A(J) = - ( RA*PM12 + HALF*RB*(RIJ+RM1) ) 
            COLIND(J) = INDEX - 1                                      
            J = J+1                                                
*                                                                       
*           SUPER-DIAGONAL.  
*
 30         IF (IX.EQ.NX) 
     $         GOTO 40                                     
            A(J) = -RA*P12 + HALF*RB*(RIJ+R1) 
            COLIND(J) = INDEX + 1                                      
            J = J+1                                                
*                                                                       
*           HIGHEST BAND.  
*
 40         IF (JY.EQ.NY) 
     $         GOTO 50                                     
            A(J) = -Q12 + HALF*HY*(SIJ+S1) 
            COLIND(J) = INDEX + NX                                     
            J = J+1                                                
*                                                                       
 50         CONTINUE                                                  
*
*           BOUNDARY CONDITIONS.     
*
            IF (JY.NE.1) 
     $         GOTO 60                                      
            COEF = - ( QM12 + HALF*HY*(SIJ+SM1) )
            IF (BNDRY(1).EQ.0)                                     
     $         RHS(INDEX) = RHS(INDEX) - COEF*U(XI,YJ-HY)          
            IF (BNDRY(1).EQ.1)  
     $         A(JD) = A(JD) + COEF               
*                                                                       
 60         IF (IX.NE.1) GOTO 70                                     
               COEF = - ( RA*PM12 + HALF*RB*(RIJ+RM1) ) 
            IF (BNDRY(2).EQ.0)                                     
     $         RHS(INDEX) = RHS(INDEX) - COEF*U(XI-HX,YJ)          
            IF (BNDRY(2).EQ.1)  
     $         A(JD) = A(JD) + COEF               
*                                                                       
 70         IF (IX.NE.NX) GOTO 80                                     
               COEF = -RA*P12 + HALF*RB*(RIJ+R1) 
            IF (BNDRY(3).EQ.0)                                     
     $         RHS(INDEX) = RHS(INDEX) - COEF*U(XI+HX,YJ)          
            IF (BNDRY(3).EQ.1)  
     $         A(JD) = A(JD) + COEF               
*                                                                       
 80         IF (JY.NE.NY) 
     $         GOTO 100                                    
            COEF = -Q12 + HALF*HY*(SIJ+S1) 
            IF (BNDRY(4).EQ.0)                                     
     $         RHS(INDEX) = RHS(INDEX) - COEF*U(XI,YJ+HY)          
            IF (BNDRY(4).EQ.1)  
     $         A(JD) = A(JD) + COEF               
*                                                                       
 100     CONTINUE                                                     
*
 110  CONTINUE                                                        
*                                                                       
      RETURN                                                          
*
*     End of MATPDE
*
      END                                                             
*                                                                       
      DOUBLE PRECISION FUNCTION  PC (X,Y)
      DOUBLE PRECISION   X, Y 
      INTRINSIC          EXP
      PC = EXP(-X*Y)                                                  
      RETURN                                                          
      END                                                             
*                                                                       
      DOUBLE PRECISION FUNCTION  QC (X,Y)
      DOUBLE PRECISION   X, Y 
      INTRINSIC          EXP
      QC = EXP(X*Y)                                                   
      RETURN                                                          
      END                                                             
*
      DOUBLE PRECISION FUNCTION TC (X,Y) 
      DOUBLE PRECISION   X, Y 
      DOUBLE PRECISION   ONE
      PARAMETER          ( ONE = 1.0D+0 ) 
      TC = ONE / (ONE + X + Y)
      RETURN                                                          
      END                                                             
*
      DOUBLE PRECISION FUNCTION  U (X,Y) 
      DOUBLE PRECISION   X, Y
      DOUBLE PRECISION   PI
      PARAMETER          (PI = 3.14159265358979D+0)
      INTRINSIC          SIN, EXP
      U = X * EXP(X*Y) * SIN(PI*X) * SIN(PI*Y)                        
      RETURN                                                          
      END                                                             
*                                                                       
      DOUBLE PRECISION FUNCTION  RC (BETA,X,Y) 
      DOUBLE PRECISION   BETA, X, Y 
      RC = BETA * (X+Y)                                           
      RETURN                                                          
      END                                                             
*                                                                       
      DOUBLE PRECISION FUNCTION  SC (GAMMA,X,Y)
      DOUBLE PRECISION   GAMMA, X, Y 
      SC = GAMMA * (X+Y)                                           
      RETURN                                                          
      END                                                             
*                                                                       
      DOUBLE PRECISION FUNCTION F (BETA, GAMMA, X,Y) 
      DOUBLE PRECISION   BETA, GAMMA, X, Y
*
*
      DOUBLE PRECISION   PI
      PARAMETER          (PI = 3.14159265358979D+0) 
      DOUBLE PRECISION   ONE, TWO
      PARAMETER          (ONE = 1.0D+0, TWO = 2.0D+0) 
*
      DOUBLE PRECISION   PXY, PXXY, QXY, QYXY, CXY, RXY, RXXY,
     $                   SXY, SYXY, EXY, SX, CX, SY, CY, A, AX,
     $                   AXX, AY, AYY
*
      DOUBLE PRECISION   PC, QC, RC, SC
      EXTERNAL	         PC, QC, RC, SC
      INTRINSIC          COS, SIN, EXP
*                                                                       
      PXY  = PC(X,Y)                                                  
      PXXY = -Y * EXP(-X*Y)                                           
      QXY  = QC(X,Y)                                                  
      QYXY = X * EXP(X*Y)                                             
      CXY  = ONE / ( ONE + X + Y)                                        
      RXY  = RC(BETA,X,Y)                                                  
      RXXY = BETA 
      SXY  = SC(GAMMA,X,Y)                                                  
      SYXY = GAMMA
*                                                                       
      EXY = EXP(X*Y)                                                  
      SX  = SIN(PI*X)                                                 
      CX  = COS(PI*X)                                                 
      SY  = SIN(PI*Y)                                                 
      CY  = COS(PI*Y)                                                
*                                                                       
      A   = X * EXY * SX * SY                                         
      AX  = EXY * (PI * X * CX + (ONE + X * Y) * SX) * SY 
      AXX = EXY * (PI * (-PI * X * SX + (X * Y + TWO) * CX) + Y * SX)
     $        * SY + Y * AX                                             
      AY  = EXY * (PI * CY + X * SY) * X * SX                         
      AYY = X * (EXY * PI * (-PI * SY + X * CY) * SX + AY)            
*                                                                       
      F = - (PXY*AXX + QXY*AYY)  +                                    
     $      (TWO*RXY - PXXY) * AX  +  (TWO*SXY - QYXY) * AY  + 
     *      (RXXY + SYXY + CXY) * A                                     
*
      RETURN                                                          
      END                                                             
