      PROGRAM OSCHDD
c     PROGRAM TESTOWY DLA SCHUR_DWUDIAGONAL L=>R (ALPHA>0)
      implicit none
      include 'mpif.h'
      integer lp,s,M,R,size,rank,ierr

      call MPI_Init(ierr)
      call MPI_Comm_Size(MPI_Comm_World,size,ierr)
      call MPI_Comm_Rank(MPI_Comm_World,rank,ierr)
      lp=size
       s=rank

      open(unit=10,file='SCHDATA')
      rewind(10)
      read(10,*) R
      read(10,*) M
      close(10)
      

     

      call SUBTEST(lp,s,R,M)
      write(*,*)'GOTOWE s=',s
      call MPI_Finalize(ierr)
      stop 100
      end

      subroutine SUBTEST(lp,s,R,M)
      implicit none
      include 'mpif.h'
      integer lp,s,R,M,ierr
      integer j,k
      real*8 O,E,F,D,A,SF,SD,SA,X,XX

      dimension F(0:M),D(0:M),A(0:M),XX(0:M)
      dimension SF(0:R),SD(0:R),SA(0:R),X(0:R)
      

c     X WYNIK W PROCESORZE S, XX SKOMASOWANY WYNIK X WPROCERORZE s=0

c     STALE
      O=real(0)
      E=real(1)
c     CZYTANIE DANYCH CD
      open(unit=11,file='F')
      rewind(11)
      do k=0,M,1
      read(11,*)F(k)
      enddo
      close(11)
      
      open(unit=12,file='D')
      rewind(12)
      do k=0,M,1
      read(12,*) D(k)
      enddo
      close(12)
      
      open(unit=13, file='A')
      rewind(13)
      do k=0,M,1
      read(13,*) A(k)
      enddo
      close(13)           
      


c     TWORZENIE DANYCH DLA PODPROGRAMU SCHDD
      do k=0,R,1
      SF(k)=F(s*(R+1)+k)
      SD(k)=D(s*(R+1)+k)
      SA(k)=A(s*(R+1)+k)
      enddo

      
      call SCHDD(lp,s,R,M,SF,SD,SA,X)
      call MPI_Barrier(MPI_Comm_World,ierr)
      
c     KOMASACJA X DO XX W PROC. s=0
      call MPI_Gather
     $(X,R+1,MPI_Real8,XX,R+1,MPI_Real8,0,MPI_Comm_World,ierr)

c     OUTPUT XX
      if(s.eq.0)then
      open(unit=15,file='XX')
      rewind(15)
      do k=0,M,1
      write(15,*)'k=',k,' XX=',XX(k)
      enddo
      close(15)
      endif

c     SUBTEST OUT
      return
      end

      subroutine SCHDD(lp,s,R,M,SF,SD,SA,X)

      implicit none
      include 'mpif.h'
      integer lp,s,R,M,j,k,l,ierr
      real*8 SF,SD,SA,X,TX,W,WR,O,E,C,RHS,Z,PA,PAKA

      dimension SF(0:R),SD(0:R),SA(0:R),X(0:R),TX(0:R),W(0:R)
      dimension Z(0:(lp-1)),PA(0:2),PAKA(0:(3*lp-1))
      dimension RHS(0:(lp-1)),C(0:(lp-1)),WR(0:(lp-1))

c     STALE
      O=real(0)
      E=real(1)

c     PIERWSZE PRZYBLIZENIE ~X I W
      TX(0)=SF(0)/SD(0)
       W(0)=    E/SD(0)
      do k=1,R,1
      TX(k)=(SF(k)-SA(k)*TX(k-1))/SD(k)
       W(k)=        -SA(k)*W(k-1)/SD(k) 
      enddo
       

c     UKŁAD SCHURA: PRAWA STRONA  WSPOLCZYNNIKI C I W(R)
c     PAKOWANIE DO WYSYLKI
     
      
c     PRAWA STRONA
      PA(0)=TX(R)
c     WSPOLCZYNNIKI
      if(s.eq.0)then
      PA(1)=O
      else
      PA(1)=SA(0) 
      endif
      PA(2)=W(R)
c     ALLGATHER DLA PRAWYCH STRON I WSPOLCZYNNIKOW PA==>PAKA
      call MPI_Allgather
     $(PA,3,MPI_Real8,PAKA,3,MPI_Real8,MPI_Comm_World,ierr)

c     ROZPAKOWYWANIE PAKI
    
      do j=0,3*lp-1,1
      k=j/3
      l=mod(j,3)
      if(l.eq.0)then
      RHS(k)=PAKA(j)
      elseif(l.eq.1)then
      C(k)=PAKA(j)
      else
      WR(k)=PAKA(j)
      endif
      enddo
      

c     ROZWIAZYWANIE UKŁADU SCHURA W KAZDYM PROCESORZE
      Z(0)=RHS(0)
      do k=1,lp-1,1
      Z(k)=RHS(k)-C(k)*WR(k)*Z(k-1)      
      enddo

c     POPRAWKI: X W PROCESORZE NR s
      do k=0,R,1
      if(s.eq.0)then
      X(k)=TX(k)
      else
      X(k)=TX(k)-C(s)*W(k)*Z(s-1)
      endif
      enddo


c     SCHDD OUT
      return
      end

