c
c     Sample_Code_13
c     Solves a linear elliptic PDE using Chebyshev Collocation Method
c 
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION A(NX1*NY1,NX1*NY1),B(NX1*NY1)
      INTEGER IPVT(NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
      OPEN(11,FILE='Sample_Code_13_output.txt')
c
      PI    = 4.d0*ATAN(1.d0)
      x1    = 0.d0
      x2    = 1.d0
      y1    = 0.d0
      y2    = 1.d0
c
      call Chebyshev
c
      call MATRIX(A)
c
      call ludcmp(A,NX1*NY1,NX1*NY1,ipvt,dummy)
c
      do i = 1, NX1*NY1
         B(i) = 0.d0
      enddo
c
c     Implement boundary conditions
c
      do j = 1, NY1
         B((  1-1)*NY1+j) = 1.d0
         B((NX1-1)*NY1+j) = 1.d0
      enddo
c
      do i = 1, NX1
         B((i-1)*NY1+  1) = 1.d0
      enddo
c
      do i = 2, NX1-1
         B((i-1)*NY1+NY1) = 0.d0
      enddo
c
      call lubksb(A,NX1*NY1,NX1*NY1,ipvt,B)
c
      write(11,'(20x,17f10.4)') (X(i), i = 1, NX1)
      write(11,*)
      do j = 1, NY1
      write(11,'(f10.4,10x,17f10.4)') Y(j),(B((i-1)*NY1+j), i = 1, NX1)
      enddo
c
1000  stop 
      end
c
c
c
      SUBROUTINE MATRIX(A)
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION A(NX1*NY1,NX1*NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
c
      do i = 1, NX1*NY1
      do j = 1, NX1*NY1
        A(i,j) = 0.d0
      enddo
      enddo
c
c     Construct Residuals using Chebyshev Collocation
c
	do i = 2, NX1-1
        do j = 2, NY1-1
          do n = 1, NX1
            A((i-1)*NY1+j,(n-1)*NY1+j) = A((i-1)*NY1+j,(n-1)*NY1+j) + 
     !                                   DDX(i,n)
          enddo
          do m = 1, NY1
            A((i-1)*NY1+j,(i-1)*NY1+m) = A((i-1)*NY1+j,(i-1)*NY1+m) + 
     !                                   DDY(j,m)
          enddo
        enddo
      enddo
c
c     Implement boundary conditions
c
	do j = 1, NY1 
         A((  1-1)*NY1+j,(  1-1)*NY1+j) = 1.d0
         A((NX1-1)*NY1+j,(NX1-1)*NY1+j) = 1.d0
      enddo
c
	do i = 1, NX1 
         A((i-1)*NY1+  1,(i-1)*NY1+  1) = 1.d0
      enddo
c
	do i = 2, NX1-1 
         A((i-1)*NY1+NY1,(i-1)*NY1+NY1) = 1.d0
      enddo
c
      return   
      end 
c
c
c
      SUBROUTINE Chebyshev
      IMPLICIT DOUBLE PRECISION (a-h,o-z) 
      parameter (NX=16,NY=16,NX1=NX+1,NY1=NY+1)
      DOUBLE PRECISION XCHEB(NX1),CX(NX1),CBARX(NX1)
      DOUBLE PRECISION YCHEB(NY1),CY(NY1),CBARY(NY1)
      COMMON /CHEB/ X1,X2,Y1,Y2,PI,Alpha,Beta
      COMMON /TRANS/ X(NX1),Y(NY1)
      COMMON /COLLEC/ DX(NX1,NX1),DDX(NX1,NX1),DY(NY1,NY1),DDY(NY1,NY1)
c
c     Vector XCHEB (1,-1) contains the Gauss-Lobatto grid. 
      do i = 1, NX1 
        XCHEB(i) = dcos((i-1)*pi/NX)
        CBARX(i)  = 1.d0  
        CX(i)     = 1.d0  
      enddo
      do i = 1, NY1 
        YCHEB(i)  = dcos((i-1)*pi/NY)
        CBARY(i)  = 1.d0  
        CY(i)     = 1.d0  
      enddo
c
      CX(1)      = 2.d0  
      CBARX(1)   = 2.d0  
      CBARX(NX1) = 2.d0  
      CY(1)      = 2.d0  
      CBARY(1)   = 2.d0  
      CBARY(NY1) = 2.d0  
c
c     Vector X maps XCHEB to the interval (x1,x2)	==> X(i) = Alpha + Beta*XCHEB(i)
      Alpha1 = (x1+x2)/2.d0
      Beta1  = (x1-x2)/2.d0
      Alpha2 = (y1+y2)/2.d0
      Beta2  = (y1-y2)/2.d0
c
      do i = 1, NX1 
        X(i) = Alpha1 + Beta1*XCHEB(i)
      enddo
c
      do i = 1, NY1 
        Y(i) = Alpha2 + Beta2*YCHEB(i)
      enddo
c
      do i = 1, NX1
         ii = i-1
         do j = 1, NX1
            jj = j-1
            if (ii .ne. jj) then
              DX(i,j) = CBARX(i)/CBARX(j)*(-1)**(ii+jj)
     !                 /(XCHEB(i)-XCHEB(j))
            else
              if (ii .ne. 0 .and. ii .ne. NX) then
              DX(i,j) = -XCHEB(j)/2.d0/(1.d0-XCHEB(j)**2)
              endif
            endif
         enddo
      enddo
      DX( 1, 1) = (2.d0*dble(NX)**2+1.d0)/6.d0
      DX(NX1,NX1) =-(2.d0*dble(NX)**2+1.d0)/6.d0
c
      do i=1,NX1
        do j=1,NX1 
          DDX(i,j) = 0.d0
          do k=1,NX1
            DDX(i,j) = DDX(i,j) + DX(i,k)*DX(k,j)
          enddo
        enddo
      enddo
c
      do i = 1, NX1
         do j = 1, NX1
             DX(i,j) =  DX(i,j)/BETA1
            DDX(i,j) = DDX(i,j)/BETA1/BETA1
         enddo
      enddo
c
      do i = 1, NY1
         ii = i-1
         do j = 1, NY1
            jj = j-1
            if (ii .ne. jj) then
              DY(i,j) = CBARY(i)/CBARY(j)*(-1)**(ii+jj)
     !                 /(YCHEB(i)-YCHEB(j))
            else
              if (ii .ne. 0 .and. ii .ne. NY) then
              DY(i,j) = -YCHEB(j)/2.d0/(1.d0-YCHEB(j)**2)
              endif
            endif
         enddo
      enddo
      DY( 1, 1) = (2.d0*dble(NY)**2+1.d0)/6.d0
      DY(NY1,NY1) =-(2.d0*dble(NY)**2+1.d0)/6.d0
c
      do i=1,NY1
        do j=1,NY1 
          DDY(i,j) = 0.d0
          do k=1,NY1
            DDY(i,j) = DDY(i,j) + DY(i,k)*DY(k,j)
          enddo
        enddo
      enddo
c
      do i = 1, NY1
         do j = 1, NY1
             DY(i,j) =  DY(i,j)/BETA2
            DDY(i,j) = DDY(i,j)/BETA2/BETA2
         enddo
      enddo
c
      return   
      end 
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
      SUBROUTINE ludcmp(a,n,np,indx,d)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER n,np,indx(n),NMAX
      DOUBLE PRECISION d,a(np,np),TINY
      PARAMETER (NMAX=500,TINY=1.0e-20)
      INTEGER i,imax,j,k
      DOUBLE PRECISION aamax,dum,sum,vv(NMAX)
      d=1.
      do 12 i=1,n
        aamax=0.d0
        do 11 j=1,n
          if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j))
11      continue
        if (aamax.eq.0.d0) pause 'singular matrix in ludcmp'
        vv(i)=1./aamax
12    continue
      do 19 j=1,n
        do 14 i=1,j-1
          sum=a(i,j)
          do 13 k=1,i-1
            sum=sum-a(i,k)*a(k,j)
13        continue
          a(i,j)=sum
14      continue
        aamax=0.d0
        do 16 i=j,n

          sum=a(i,j)
          do 15 k=1,j-1
            sum=sum-a(i,k)*a(k,j)
15        continue
          a(i,j)=sum
          dum=vv(i)*abs(sum)
          if (dum.ge.aamax) then
            imax=i
            aamax=dum
          endif
16      continue
        if (j.ne.imax)then
          do 17 k=1,n
            dum=a(imax,k)
            a(imax,k)=a(j,k)
            a(j,k)=dum
17        continue
          d=-d
          vv(imax)=vv(j)
        endif
        indx(j)=imax
        if(a(j,j).eq.0.d0)a(j,j)=TINY
        if(j.ne.n)then
          dum=1./a(j,j)

          do 18 i=j+1,n
            a(i,j)=a(i,j)*dum
18        continue
        endif
19    continue
      return
      END
c
      SUBROUTINE lubksb(a,n,np,indx,b)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      INTEGER n,np,indx(n)
      DOUBLE PRECISION a(np,np),b(n)
      INTEGER i,ii,j,ll
      DOUBLE PRECISION sum
      ii=0
      do 12 i=1,n
        ll=indx(i)
        sum=b(ll)
        b(ll)=b(i)
        if (ii.ne.0)then
          do 11 j=ii,i-1
            sum=sum-a(i,j)*b(j)
11        continue
        else if (sum.ne.0.d0) then
          ii=i
        endif
        b(i)=sum
12    continue
      do 14 i=n,1,-1
        sum=b(i)
        do 13 j=i+1,n
          sum=sum-a(i,j)*b(j)
13      continue
        b(i)=sum/a(i,i)
14    continue
      return
      END