program inverse 
 implicit none 
 double precision, dimension(:,:),allocatable  :: A, B,C 
 integer, dimension(:), allocatable  :: IPIV
 double precision, dimension(:),allocatable ::  WORK
 integer :: i,j,ndim,INFO,LDA,LWORK,N,numarg   
 character (len=240)inputfile,outputfile,str 

 numarg = iargc()

  if (numarg .lt. 2) then
     write(*,*)"ERROR ! Plase give use as follows"
     write(*,*)"./compute_inverse  <input file> <output file>"
     stop
  end if

  CALL GETARG(1, STR)
  read (STR,'(A)')inputfile 

  CALL GETARG(2, STR)
  read (STR,'(A)')outputfile 
      
 open(11,file=inputfile)
 read(11,*)ndim 
 LDA=ndim
 LWORK=ndim
 N= ndim
 allocate(A(ndim,ndim))
 allocate(B(ndim,ndim))
 allocate(C(ndim,ndim))

 do i =  1, ndim 
  read(11,*)(A(i,j),j=1,ndim)
  end do   
 close(11)

! Copy A into B for inverse 
 B=A

 allocate(IPIV(ndim))
 allocate(WORK(ndim))

 call  DGETRF(N, N, B, LDA, IPIV, INFO)
 call  DGETRI(N, B, LDA, IPIV, WORK, LWORK, INFO )

  open(11,file=outputfile)
  write(11,*) n  
  do i = 1, n
     write(11,'(100(E12.4))')(B(i,j),j=1,n)
  end do
  close(11) 
  
  C = MATMUL(A,B) 
  
  write(*,*)"checking---" 
  write(*,*)"----------------A^-1 X A -------------------------------------------"
  do i = 1, n
     write(*,'(100(F12.4))')(C(i,j),j=1,n)
  end do

end program inverse