!----------------------------------------------------------------- ! This program compute the inverse of a general square matrix ! using BLAS/LAPACK. ! The matrix has to be given from the input file which should ! have the first line the # of rows followd by the matrix entries ! example of the inpute file ! 3 ! 1.0 2.0 4.0 ! 4.9 8.9 1.0 ! 3.0 1.9 9.0 ! Jayanti Prasad, Sept 20, 2011 ! Comments & Feedback : prasad.jayanti@gmail.com !--------------------------------------------------------------- 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,str numarg = iargc ( ) if (numarg .lt. 1) then write(*,*)"ERROR ! Plase give the input file prefix" write(*,*)"./compute_inverse <input file>" stop end if CALL GETARG(1, STR) read (STR,'(A)')inputfile 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 ) write(*,*)"-----------------A------------------------------------------------" do i = 1, n write(*,'(6(E12.4))')(A(i,j),j=1,n) end do write(*,*)"-----------------A^-1 ---------------------------------------------" do i = 1, n write(*,'(6(E12.4))')(B(i,j),j=1,n) end do C = MATMUL(A,B) write(*,*)"checking---" write(*,*)"----------------A^-1 X A -------------------------------------------" do i = 1, n write(*,'(6(F12.4))')(C(i,j),j=1,n) end do write(*,'(6(F12.6))')(sqrt(B(I,i)),i=1,n) end program inverse