      program mpigausselmination
      include 'mpif.h'
      integer ndims,xmax,ymax,nnodes
      logical reorder
      parameter(ndims=2,num_neigh=4,reorder=.true.)
      parameter(xmax=400,ymax=400)
      integer comm,comm_cart,dims(ndims),neigh_def(ndims),ierr
      integer own_position(ndims),rank_source,rank_dest
      logical periods(ndims)
      integer status(MPI_STATUS_SIZE)
      double precision,allocatable,dimension(:,:) :: A,B,OA
      double precision,allocatable,dimension(:)   :: X,sX,sXO
      double precision m
      data dims /ndims * 0/
      integer  myid,nitems,nlines
      double precision time_begin,time_end 

      comm=MPI_COMM_WORLD
      call MPI_init(ierr)
      call MPI_COMM_RANK(comm,myid,ierr)
      call MPI_COMM_SIZE(comm,nnodes,ierr)
     
      nlines = xmax/nnodes
      allocate(OA(xmax,ymax+1))
      allocate(A(nlines,ymax+1))
      allocate(B(nlines,ymax+1))
      allocate(X(nlines))
      allocate(sX(xmax))
      allocate(sXO(xmax)) 
      nitems = nlines*(ymax+1)
       
!initialize A matrix,which includes coefficient matrix and right-hand side
!at its ymax+1'th colum

     
      if (myid == nnodes -1) then
        call random_number(OA)
        OA = aint(30.0*OA+1)
        call random_number(sXO)
        sXO = aint(30.0*sXO)+1
        do i = 1,xmax
           OA(i,ymax+1) = sum(OA(i,1:ymax)*sXO)
        enddo  
        do np = 0,nnodes-2
           call MPI_SEND(OA(np*nlines+1:(np+1)*nlines,:),nitems,&
             MPI_DOUBLE_PRECISION,np,0,comm,ierr)
        enddo
        A = OA(myid*nlines+1:(myid+1)*nlines,:)

     else
        call MPI_RECV(A,nitems,MPI_DOUBLE_PRECISION,nnodes-1,0,comm,&
           status,ierr)
     endif
!start timer

    
     time_begin =  MPI_Wtime()

     do np = 0,nnodes-1
       lbase = np*nlines
       if (myid == np ) then
!solve local gauss elimination
         Do i =1,nlines
           Do j=i+1,nlines
             m = a(j,i+lbase)/a(i,i+lbase)
             a(j,i+lbase:ymax+1) = a(j,i+lbase:ymax+1) - m*a(i,i+lbase:ymax+1)
            enddo
         enddo
         B = A
       endif
       call MPI_BCAST(B,nitems,MPI_DOUBLE_PRECISION,np,comm,ierr)
       if (myid > np) then
          Do i=1,nlines
           do j=1,nlines
            m = a(j,i+lbase)/B(i,i+lbase)
            a(j,i+lbase:ymax+1) = a(j,i+lbase:ymax+1) - m*B(i,i+lbase:ymax+1)
          enddo
        enddo
       endif
      
     enddo
  
    
    ls = myid*nlines+1
    le = (myid+1)*nlines

    if (myid == nnodes-1) then
       sX(le) = A(nlines,ymax+1)/A(nlines,ymax)
       j = nlines-1
       do i=le-1,ls,-1
          sX(i) = (a(j,ymax+1)-sum(a(j,i+1:ymax)*sX(i+1:ymax)))/a(j,i) 
          j = j-1
       enddo 
       if (nnodes > 1) then
          call MPI_SEND(sX,ymax,MPI_DOUBLE_PRECISION,myid-1,0,comm,ierr)
       endif 
    
    endif
       
    do np = nnodes-2,0,-1
      if (myid == np) then
        call MPI_RECV(sX,ymax,MPI_DOUBLE_PRECISION,&
              myid+1,0,comm,status,ierr)
        j = nlines
        do i=le,ls,-1
          sX(i) = (a(j,ymax+1)-sum(a(j,i+1:ymax)*sX(i+1:ymax)))/a(j,i) 
          j = j-1
        enddo 
        if (myid > 0) then
           call MPI_SEND(sX,ymax,MPI_DOUBLE_PRECISION,myid-1,0,comm,ierr)
        endif
      endif
    enddo

    time_end =  MPI_Wtime()
    if (myid == 0) then
        print *,'Elapsed time for solving linear equation with dimension of'
        print *,'coefficient matrix',xmax,'by ',ymax,'is'
        print *,time_end-time_begin,' seconds'
    endif

    call MPI_FINALIZE(ierr)

      end
