      include 'ELIM.f'
      PROGRAM LINEQ
c     LINEEAR SYSTEM WITH 3-DIAGONAL MATRIX
c     WITH HELP OF EXTERNAL SUBROUTINE ELIM.f
c     USING DISC

c     DECLARATIONS

      integer N,p

c     ON AND FROM THE CONSOL
    1 write(*,*)'PROGRAM LINEQ'
      write(*,*)'MATRIX DIMENSION N PLEASE'
      read (*,*) N
      
c     CALLING SUBROUTINE SUBLINEQ()
      call SUBLINEQ(N)
c     WHAT TO DO AT THE END?
      write(*,*)' 0 - END'
      write(*,*)' 1 - REPEAT'
      read (*,*) p
      if(p.eq.0)then
      goto 100
      endif
      goto 1
  100 stop 100
      end
      
      subroutine SUBLINEQ(N)
c     DECLARATIONS
c     Ax=f
c     A IS GIVEN AS: 
c     UNDER DIAGONAL,U
c     DIAGONAL D, 
c     AND UPPER DIAGONAL P.
c     WARNING: DIMENSIONS FROM ZERO TO N!
c
c     |d0,p0,0 ,0 ,0 ||x0|    |f0|
c     |u1,d1,p1,0 ,0 ||x1|    |f1|
c     |0 ,u2,d2,p2,0 ||x2| =  |f2|
c     |0 ,0 ,u3,d3,p3||x3|    |f3|
c     |0 ,0 ,0 ,u4,d4||x4|    |f4|
c

      integer N,j,k
      real*8 A,f,x,MAG,O,E
      dimension U(0:N),D(0:N),P(0:N),f(0:N),x(0:N)
      dimension MAG(0:N,0:4)
c     CONSTANTS
      O=real(0)
      E=real(1)

c     ELIM(N,MAG,X):
c     MAG(k,0) U(k)
c     MAG(k,1) D(k)
c     MAG(k,2) P(k)
c     MAG(k,3) O
c     MAG(k,4) f(k)

c     READING DATA (CONSOL) 
      write(*,*)' U(0)-U(N) PLEASE'
      do k=0,N,1
      read (*,*) U(k)
      enddo
      write(*,*)' D(0)-D(N) PLEASE'
      do k=0,N,1
      read (*,*) D(k)
      enddo
      write(*,*)' P(0)-P(N) PLEASE'
      do k=0,N,1
      read (*,*) P(k)
      enddo
      write(*,*)' f(0)-f(N) PLEASE'
      do k=0,N,1
      read (*,*) f(k)
      enddo
      write(*,*)'END, <ENTER> TO CONTINUE'
      read (*,*)
      
c     WRITTING INPUT DATA TO DISC
c     U
      open(unit=10,file='U')
      rewind(10)
      do k=0,N,1
      write(10,*)'k=',k,'U=', U(k)
      enddo
      close(10)
      
c     D
      open(unit=11,file='D')
      rewind(11)
      do k=0,N,1
      write(11,*)'k=',k,'D=', D(k)
      enddo
      close(11)
      
c     P
      open(unit=12,file='P')
      rewind(12)
      do k=0,N,1
      write(12,*)'k=',k,'P=', P(k)
      enddo
      close(12)
      
c     f
      open(unit=13,file='f')
      rewind(13)
      do k=0,N,1
      write(13,*)'k=',k,'f=', f(k)
      enddo
      close(13)
      

c     LOADING MAG
      do k=0,N,1
      if(k.eq.0)then
      MAG(k,0)=O
      MAG(k,1)=D(k)
      MAG(k,2)=P(k)
      MAG(k,3)=O
      MAG(k,4)=f(k)
      elseif(k.eq.N)then
      MAG(k,0)=U(k)
      MAG(k,1)=D(k)
      MAG(k,2)=O
      MAG(k,3)=O
      MAG(k,4)=f(k)
      else
      MAG(k,0)=U(k)
      MAG(k,1)=D(k)
      MAG(k,2)=P(k)
      MAG(k,3)=O
      MAG(k,4)=f(k)
      endif
      enddo

c     CALLING ELIM.f
      call ELIM(N,MAG,x)

c     OUTPUT x ON DISC
      open(unit=20,file='x')
      rewind(20)
      do k=0,N,1
      write(20,*)'k=',k,'x=',x(k)
      enddo
      close(20)

c     GOING BACK
      return
      end