m***@gmail.com
2016-10-19 22:43:13 UTC
Hello Everyone !
First of all, I am glad to see this group for help in FORTRAN. And I hope to get valuable suggestions and help here.
I am working in FORTRAN 90 where I need to calculate Inverse of a 7x7 Matrix accurately. I have used a subroutine in main program which augmenting the matrix with Identity matrix and then perform calculations for Inverse determination by pivoting. But the problem is that most of the time I do not get the correct answer. I verify my inverse using fact that : matrix A multiplied by its Inverse = Identity Matrix . But it do not gives Identity matrix when I use the Inverse calculated by the subroutine. The code is attached at end (named: MatInv.f90) .
I am working in Ubuntu 16.04 LTS. I have installed FORTRAN 90 and I compile my programs by command " gfortran porgram.f90 -o program" in terminal.
I searched on INTERNET that LAPACK is a library which can efficiently calculates inverse of large matrices. So following instruction in a video on YOUTUBE, I installed LAPACK by commands in terminal as follows:
1) tar zxvf lapack-3.6.0.gz
2) cd lapack-3.6.0
◇ gfortran
cp make.inc.example make.inc
4) make blaslib
5) make lapacklib
6) sudo ln -s $HOME/lapack-3.6.0/librefblas.a /usr/local/lib/libblas.a
7) sudo ln -s $HOME/lapack-3.6.0/liblapack.a /usr/local/lib/liblapack.a
(on checking in /usr/local/bin these libraries are present there)
Now I want to use the program which uses LAPACK to find inverse of a large matrix but I do not know how to compile the code using these libraries. After searching on INTERNET and using file matrix_inverse.f90 (attached at end named: inverse_mat.f90) when I use command :
gfortran my_program.f90 -llapack -lblas
It says:
(.text+0x20): undefined reference to `main'collect2: error: ld returned 1 exit status
I am completely stuck here. I do not know certain things as follows:
(I). Whether I need to use LAPACK or not to find Inverse of a 7x7 matrix
(II). Why my program (without LAPACK) do not give correct result for inverse. If the inverse is accurate it should satisfy condition:
matrix A multiplied by its Inverse = Identity Matrix
but it do not satisfy this condition.
(III). If I use LAPACK, then how to accuratly use it to write and compile a program for inverse
(IV). What is the reason of the error mentioned above
(V). How to check whether I have installed LAPACK correctly or not
(VI). How I can install and use LAPACK, step by step instruction.
(VII). Finally, without knowing answers of all above questions, how can I find ACCURATE Inverse of a 7x7 martrix in FORTRAN 90?
I shall be highly grateful for any suggestions and help.
Best Regards,
Mirza Younis Baig
####################################################################
File:
{code}
SUBROUTINE MATINV(MATSUM,INVERSE)
IMPLICIT NONE
Double Precision, INTENT(IN),DIMENSION(100,100) :: MATSUM
Double Precision, INTENT(OUT),DIMENSION(100,100) :: INVERSE
INTEGER :: I,J,NROWS,NCOLS,M,N,K,p,q
REAL :: summ,xmult,xmult2,xmult3
Double Precision, DIMENSION(100,100) :: A, Ainv, Ident
NROWS=7
NCOLS=7
Do I=1,NROWS
DO J=1,NCOLS
A(I,J)=MATSUM(I,J)
END Do
END DO
PRINT*,''
PRINT*,''
PRINT*,' ---------------------------------------------------'
!Echoing Input Matrix
PRINT*,'>>>>>>>>Input Matrix<<<<<<<<'
DO I=1,NROWS
Print '(7(1x,F9.4))', (A(I,J),J=1,NCOLS)
END DO
PRINT*,' ---------------------------------------------------'
!Defining Identity Matrix
Do I=1,NROWS
Do J=1,NCOLS
IF(I.EQ.J) THEN
Ident(I,J)=1.0
ELSE
Ident(I,J)=0.0
END IF
END DO
END DO
!Row Operations
! PRINT*,'>>>>>>>>Row Operations<<<<<<<<'
DO I=1,NROWS
DO J=1,NCOLS
IF(I.EQ.J) THEN
xmult=A(I,J)
DO K=1,NCOLS
A(I,K)=A(I,K)/xmult
Ident(I,K)=Ident(I,K)/xmult
END DO
DO M=I+1,NROWS
xmult2=A(M,J)
DO N=1,NCOLS
A(M,N)=A(M,N)-xmult2*A(I,N)
Ident(M,N)=Ident(M,N)-xmult2*Ident(I,N)
END DO
END DO
IF(I.GT.1) then
DO p=I-1,1,-1
xmult3=A(p,J)
Do q=1,NCOLS
A(p,q)=A(p,q)-xmult3*A(I,q)
Ident(p,q)=Ident(p,q)-xmult3*Ident(I,q)
END DO
END DO
END IF
END IF
END DO
! PRINT*,I,'>>>>'
! PRINT*,''
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(A(K,J),J=1,NCOLS)
! END DO
! PRINT*,'======================iii=========================='
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(Ident(K,J),J=1,NCOLS)
! END DO
! PRINT*,'==================================================='
END DO
! PRINT*,'>>>>>>>>Matrix Inverse<<<<<<<<'
! DO I=1,NROWS
! Print '(10(1x,F10.4))',(Ident(I,J),J=1,NCOLS)
! END DO
! Back Substitution
DO J=1,NCOLS
DO I=NROWS,1,-1
IF(I.EQ.NROWS) THEN
Ainv(I,J)=Ident(I,J)/A(I,I)
ELSE
summ=0.0
DO K=I+1,NROWS
summ=summ+(A(I,K)*Ainv(K,J))
END DO
Ainv(I,J)=(Ident(I,J)-summ)/A(I,I)
END IF
END DO
END DO
Do I=1,NROWS
DO J=1,NCOLS
INVERSE(I,J)=Ainv(I,J)
END Do
END DO
RETURN
END SUBROUTINE MATINV
{/code}
####################################################################
File:inverse_mat.f90
{code}
! Returns the inverse of a matrix calculated by finding the LU
! decomposition. Depends on LAPACK.
function inv(A) result(Ainv)
real(dp), dimension(:,:), intent(in) :: A
real(dp), dimension(size(A,1),size(A,2)) :: Ainv
real(dp), dimension(size(A,1)) :: work ! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info
! External procedures defined in LAPACK
external DGETRF
external DGETRI
! Store A in Ainv to prevent it from being overwritten by LAPACK
Ainv = A
n = size(A,1)
! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
call DGETRF(n, n, Ainv, n, ipiv, info)
if (info /= 0) then
stop 'Matrix is numerically singular!'
end if
! DGETRI computes the inverse of a matrix using the LU factorization
! computed by DGETRF.
call DGETRI(n, Ainv, n, ipiv, work, n, info)
if (info /= 0) then
stop 'Matrix inversion failed!'
end if
end function inv
{/code}
First of all, I am glad to see this group for help in FORTRAN. And I hope to get valuable suggestions and help here.
I am working in FORTRAN 90 where I need to calculate Inverse of a 7x7 Matrix accurately. I have used a subroutine in main program which augmenting the matrix with Identity matrix and then perform calculations for Inverse determination by pivoting. But the problem is that most of the time I do not get the correct answer. I verify my inverse using fact that : matrix A multiplied by its Inverse = Identity Matrix . But it do not gives Identity matrix when I use the Inverse calculated by the subroutine. The code is attached at end (named: MatInv.f90) .
I am working in Ubuntu 16.04 LTS. I have installed FORTRAN 90 and I compile my programs by command " gfortran porgram.f90 -o program" in terminal.
I searched on INTERNET that LAPACK is a library which can efficiently calculates inverse of large matrices. So following instruction in a video on YOUTUBE, I installed LAPACK by commands in terminal as follows:
1) tar zxvf lapack-3.6.0.gz
2) cd lapack-3.6.0
◇ gfortran
cp make.inc.example make.inc
4) make blaslib
5) make lapacklib
6) sudo ln -s $HOME/lapack-3.6.0/librefblas.a /usr/local/lib/libblas.a
7) sudo ln -s $HOME/lapack-3.6.0/liblapack.a /usr/local/lib/liblapack.a
(on checking in /usr/local/bin these libraries are present there)
Now I want to use the program which uses LAPACK to find inverse of a large matrix but I do not know how to compile the code using these libraries. After searching on INTERNET and using file matrix_inverse.f90 (attached at end named: inverse_mat.f90) when I use command :
gfortran my_program.f90 -llapack -lblas
It says:
(.text+0x20): undefined reference to `main'collect2: error: ld returned 1 exit status
I am completely stuck here. I do not know certain things as follows:
(I). Whether I need to use LAPACK or not to find Inverse of a 7x7 matrix
(II). Why my program (without LAPACK) do not give correct result for inverse. If the inverse is accurate it should satisfy condition:
matrix A multiplied by its Inverse = Identity Matrix
but it do not satisfy this condition.
(III). If I use LAPACK, then how to accuratly use it to write and compile a program for inverse
(IV). What is the reason of the error mentioned above
(V). How to check whether I have installed LAPACK correctly or not
(VI). How I can install and use LAPACK, step by step instruction.
(VII). Finally, without knowing answers of all above questions, how can I find ACCURATE Inverse of a 7x7 martrix in FORTRAN 90?
I shall be highly grateful for any suggestions and help.
Best Regards,
Mirza Younis Baig
####################################################################
File:
{code}
SUBROUTINE MATINV(MATSUM,INVERSE)
IMPLICIT NONE
Double Precision, INTENT(IN),DIMENSION(100,100) :: MATSUM
Double Precision, INTENT(OUT),DIMENSION(100,100) :: INVERSE
INTEGER :: I,J,NROWS,NCOLS,M,N,K,p,q
REAL :: summ,xmult,xmult2,xmult3
Double Precision, DIMENSION(100,100) :: A, Ainv, Ident
NROWS=7
NCOLS=7
Do I=1,NROWS
DO J=1,NCOLS
A(I,J)=MATSUM(I,J)
END Do
END DO
PRINT*,''
PRINT*,''
PRINT*,' ---------------------------------------------------'
!Echoing Input Matrix
PRINT*,'>>>>>>>>Input Matrix<<<<<<<<'
DO I=1,NROWS
Print '(7(1x,F9.4))', (A(I,J),J=1,NCOLS)
END DO
PRINT*,' ---------------------------------------------------'
!Defining Identity Matrix
Do I=1,NROWS
Do J=1,NCOLS
IF(I.EQ.J) THEN
Ident(I,J)=1.0
ELSE
Ident(I,J)=0.0
END IF
END DO
END DO
!Row Operations
! PRINT*,'>>>>>>>>Row Operations<<<<<<<<'
DO I=1,NROWS
DO J=1,NCOLS
IF(I.EQ.J) THEN
xmult=A(I,J)
DO K=1,NCOLS
A(I,K)=A(I,K)/xmult
Ident(I,K)=Ident(I,K)/xmult
END DO
DO M=I+1,NROWS
xmult2=A(M,J)
DO N=1,NCOLS
A(M,N)=A(M,N)-xmult2*A(I,N)
Ident(M,N)=Ident(M,N)-xmult2*Ident(I,N)
END DO
END DO
IF(I.GT.1) then
DO p=I-1,1,-1
xmult3=A(p,J)
Do q=1,NCOLS
A(p,q)=A(p,q)-xmult3*A(I,q)
Ident(p,q)=Ident(p,q)-xmult3*Ident(I,q)
END DO
END DO
END IF
END IF
END DO
! PRINT*,I,'>>>>'
! PRINT*,''
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(A(K,J),J=1,NCOLS)
! END DO
! PRINT*,'======================iii=========================='
! DO K=1,NROWS
! Print '(12(1x,F10.4))',(Ident(K,J),J=1,NCOLS)
! END DO
! PRINT*,'==================================================='
END DO
! PRINT*,'>>>>>>>>Matrix Inverse<<<<<<<<'
! DO I=1,NROWS
! Print '(10(1x,F10.4))',(Ident(I,J),J=1,NCOLS)
! END DO
! Back Substitution
DO J=1,NCOLS
DO I=NROWS,1,-1
IF(I.EQ.NROWS) THEN
Ainv(I,J)=Ident(I,J)/A(I,I)
ELSE
summ=0.0
DO K=I+1,NROWS
summ=summ+(A(I,K)*Ainv(K,J))
END DO
Ainv(I,J)=(Ident(I,J)-summ)/A(I,I)
END IF
END DO
END DO
Do I=1,NROWS
DO J=1,NCOLS
INVERSE(I,J)=Ainv(I,J)
END Do
END DO
RETURN
END SUBROUTINE MATINV
{/code}
####################################################################
File:inverse_mat.f90
{code}
! Returns the inverse of a matrix calculated by finding the LU
! decomposition. Depends on LAPACK.
function inv(A) result(Ainv)
real(dp), dimension(:,:), intent(in) :: A
real(dp), dimension(size(A,1),size(A,2)) :: Ainv
real(dp), dimension(size(A,1)) :: work ! work array for LAPACK
integer, dimension(size(A,1)) :: ipiv ! pivot indices
integer :: n, info
! External procedures defined in LAPACK
external DGETRF
external DGETRI
! Store A in Ainv to prevent it from being overwritten by LAPACK
Ainv = A
n = size(A,1)
! DGETRF computes an LU factorization of a general M-by-N matrix A
! using partial pivoting with row interchanges.
call DGETRF(n, n, Ainv, n, ipiv, info)
if (info /= 0) then
stop 'Matrix is numerically singular!'
end if
! DGETRI computes the inverse of a matrix using the LU factorization
! computed by DGETRF.
call DGETRI(n, Ainv, n, ipiv, work, n, info)
if (info /= 0) then
stop 'Matrix inversion failed!'
end if
end function inv
{/code}