program mpisteep
** mpi control declarations  
      include 'mpif.h'
      integer myrank, master /0/, ierr, serr, rerr, tag 
      integer status(MPI_STATUS_SIZE)  
      real done, donetest, donetop     
** other declarations      
      integer count, maxcount 
      real xa(0:4,4),x(5),x0(4),ya(0:4),y0,y,dy(4),ndy,udy(4)
      real delta,eps,h,m     
** mpi data initalizations    
      done = 0.0
      donetest = 1.0
      donetop = 2.0
      tag = 1
*  master node parameter and data initalizations       
      delta = 0.01
      h = .75
      eps = 10**(-6)          
      maxcount = 100
      do j = 1,4
         xa(0,j) = 1
         x0(j) = xa(0,j)
      enddo         
      y0 = 1000
      y  = 0
      count = 0   
** mpi driver
      call mpi_init(ierr)
      call mpi_comm_rank(MPI_COMM_WORLD, myrank, status)
      if (myrank .EQ. master) then
         goto 200 
      else
         goto 300
      endif
** mpi cleanup
100   call mpi_finalize(ierr)
      goto 400
** master node loop
200   do while (done .LT. donetest)   
* set seek directions        
        do i = 1,4
           do j = 1,4
              xa(i,j)= xa(0,j)
           enddo
           xa(i,i)= xa(i,i)+delta 
        enddo  
* compute values in seek directions and then gradient
        do i = 0,4
           do j = 1,4
              x(j)=xa(i,j)
              enddo
           x(5) = done
           call mpi_send(x,5,MPI_REAL,i+1,tag,MPI_COMM_WORLD,serr)    
           enddo  
        do i = 0,4
          call mpi_recv(ya(i),1,MPI_REAL,i+1,tag,MPI_COMM_WORLD,status,rerr)
          enddo
        do i = 1,4
          dy(i) = (ya(i)-ya(0))/delta
          enddo 
        ndy = vnorm(dy)
        do i = 1,4
           udy(i) = dy(i)/ndy
           enddo 
* update seek point
        if (ndy .GT. 1) then     
           m=1
         else 
          m = ndy
         endif 
        do j = 1,4
          xa(0,j) = xa(0,j)-m*h*udy(j)
          enddo 
        y0 = y
        y = ya(0)
        count = count+1
* check if done        
        if ((eps .GT. abs(y-y0)) .OR. (count .EQ. maxcount)) then 
          done = donetop
          do i = 0,4
            x(5) = done
            call mpi_send(x,5,MPI_REAL,i+1,tag,MPI_COMM_WORLD,serr)    
          enddo
          endif
      enddo  
      print *, "x values: ", x(1),x(2),x(3),x(4)
      print *, "y value: ", y 
      goto 100             
** slave node program loop
300   do while (done .LT. donetest)
        call mpi_recv(x,5,MPI_REAL,0,tag,MPI_COMM_WORLD,status,rerr)         
        if (x(5) .LT. donetest) then
          y = func(x)
          call mpi_send(y,1,MPI_REAL,0,tag,MPI_COMM_WORLD,serr)     
        else 
           done = donetop
        endif    
      enddo    
      goto 100          
400   end
*
** functions and subroutines
*
      real function func(x) 
      real x(4)
      func = x(1)**4+x(2)**2+x(3)**2+x(4)**4
      return    
      end
*  
      real function vnorm(x) 
      real x(4)
      vnorm = (x(1)**2+x(2)**2+x(3)**2+x(4)**2)**(0.5)
      return    
      end