program main
  real A(4,4),AA(4,4),A2(4,4),AB(4,4)
  data A/4.0, 4.0, 4.0, 4.0, 4.0, 3.0, 3.0, 3.0, 4.0, 3.0, 2.0, 2.0, 4.0, 3.0, 2.0, 1.0/
  data A2/1.0, 3.0, 2.0, 5.0, 4.0, 4.0, 3.0, 3.0, 5.0, 1.0, 1.0, 2.0, 1.0, 4.0, 3.0, 2.0/

  print *, 'testing inverse'
  print *, 'A=  ', A
  call inverse(4,A,AA)
  print *, 'AA= ', AA
  call inverse(4,AA,AB)
  print *, 'AB= ', AB
  print *, ' '
  print *, 'A2= ', A2
  call inverse(4,A2,AA)
  print *, 'AA= ', AA
  call inverse(4,AA,AB)
  print *, 'AB= ', AB
  print *, 'end'
end

subroutine inverse(n,A,AA)
  integer n
  real A(n,n), AA(n,n)
  real TEMP(n)
  integer ROW(n), COL(n)
  integer HOLD, I_PIVOT, J_PIVOT
  real PIVOT, ABS_PIVOT, NORM1

  NORM1=0.0

  do i=1,n
    do j=1,n
      AA(i,j)=A(i,j)
    end do
  end do

  do I=1,n
    do J=1,n
      if(abs(AA(I,J))> NORM1) NORM1=abs(AA(I,J))
    end do 
  end do 

  do K=1,n
    ROW(K)=K 
    COL(K)=K 
  end do 

  do K=1,n
    PIVOT=AA(ROW(K),COL(K))
    I_PIVOT=K 
    J_PIVOT=K 
    do I=K,N
      do J=K,N
        ABS_PIVOT=abs(PIVOT)
        if (abs(AA(ROW(I),COL(J)))> ABS_PIVOT) then
          I_PIVOT=I 
          J_PIVOT=J 
          PIVOT=AA(ROW(I),COL(J))
        end if 
      end do 
    end do 

    if (ABS_PIVOT < EPSILON(NORM1) * NORM1) then
      print *, 'matrix is singular'
      return
    end if 

    HOLD=ROW(K)
    ROW(K)=ROW(I_PIVOT)
    ROW(I_PIVOT)=HOLD 
    HOLD=COL(K)
    COL(K)=COL(J_PIVOT)
    COL(J_PIVOT)=HOLD 

    AA(ROW(K),COL(K))=1.0/PIVOT 
    do J=1,n
      if (J .ne. K) then
        AA(ROW(K),COL(J))=AA(ROW(K),COL(J))* AA(ROW(K),COL(K))
      end if 
    end do 

    do I=1,n
      if (K .ne. I) then
        do J=1,n
          if (K .ne. J) then
            AA(ROW(I),COL(J))=AA(ROW(I),COL(J))- AA(ROW(I),COL(K))* AA(ROW(K),COL(J))
          end if 
        end do 
        AA(ROW(I),COL(K))=-AA(ROW(I),COL(K))* AA(ROW(K),COL(K))
      end if 
    end do 
  end do 

  do J=1,n
    do I=1,n
      TEMP(COL(I))=AA(ROW(I),J)
    end do 
    do I=1,n
      AA(I,J)=TEMP(I)
    end do 
  end do 

  do I=1,n
    do J=1,n
      TEMP(ROW(J))=AA(I,COL(J))
    end do 
    do J=1,n
      AA(I,J)=TEMP(J)
    end do 
  end do 
  return
end