program jordan
 implicit none
 integer,parameter:: n = 3
 double precision A(n,n), B(n)
 integer i,j,k

 data A /0.0d0, 2.0d0, 1.0d0, 1.0d0, 3.0d0, 1.0d0, 1.0d0, 1.0d0, 3.0d0/
 data B /2.0d0, 6.0d0, 5.0d0/

 call simeq( A, B, n )

 write(*,'(i1,f4.1)') (k,B(k),k=1,n)
 write(*,'(3f4.1)') ((A(i,j),j=1,n),i=1,n)

end program jordan

subroutine simeq( A, B, n )
 implicit none
 integer n,i,j,k,ip,ipivot
 double precision A(n,n), B(n),amax,atmp,btmp

 do k = 1, n
    ipivot = k
    amax = dabs ( a(k,k) )
    do ip = k + 1, n
       if( amax .lt. dabs ( a(ip,k) ) ) then
          ipivot = ip
          amax = dabs ( a(ip,k) )
       end if
    end do

    if( ipivot .ne. k ) then
       do j = 1, n
          atmp = a(k,j)
          a(k,j) = a(ipivot,j)
          a(ipivot,j) = atmp
       end do
       btmp = b(k)
       b(k) = b(ipivot)
       b(ipivot) = btmp
    end if

    do j = k + 1, n
       a(k,j) = a(k,j) / a(k,k) 
    end do
    b(k) = b(k) / a(k,k) 

    do i = 1, n
      if ( i .ne. k ) then
         do j = k + 1, n
           a(i,j) = a(i,j) - a(i,k) * a(k,j) 
         end do
         b(i) = b(i) - a(i,k) * b(k)
      end if
    end do

 end do

end subroutine