- Code: Select all
program testInv
implicit none
double precision A(3,3), Ainv(3,3)
integer i
A(1,1) = 1d0
A(1,2) = 2d0
A(1,3) = 3d0
A(2,1) = 4d0
A(2,2) = 5d0
A(2,3) = 6d0
A(3,1) = 7d0
A(3,2) = 8d0
A(3,3) = 10d0
write(*,*) 'A = '
do i = 1,3
write(*,*) A(i,1), A(i,2), A(i,3)
end do
call invmat(3,A,Ainv)
write(*,*) 'Ainv = '
do i = 1,3
write(*,*) Ainv(i,1), Ainv(i,2), Ainv(i,3)
end do
end
This is lapackInv.for:
- Code: Select all
SUBROUTINE INVMAT(N,A,AINV)
C INVERSE OF A REAL*8 (NXN) MATRIX
C A = ORIGINAL MATRIX
C AINV = INVERSE OF A
IMPLICIT REAL*8 (A-H, O-Z)
DIMENSION A(N,N), AINV(N,N), ipiv(N), iwork(N)
C
AINV=A
c First do LU-factorization of A, the value of which is stored
c in AINV:
call dgetrf(N, N, AINV, N, ipiv, info_lu)
c Now AINV is the LU-factorization of A, and now the inverse can be
c computed and stored in AINV:
call dgetri(N, AINV, N, ipiv, iwork, N, info_inv)
RETURN
END
BTW, the INVMAT subroutine was modified from someone else's code, and I left the implicit variable declarations from it intact.
I compiled the code like this:
- Code: Select all
my_fort90_comp -O -c testInv.for
my_fort90_comp -O -c lapackInv.for
my_fort90_comp -o testInv testInv.o lapackInv.o /path/to/mylapack_compiled_with_my_fort90_comp.a /path/to/myblas_compiled_with_my_fort90_comp.a
where my_fort90_comp is g95, gfortran, or ifort, and /path/to/my{lapack,blas}_compiled_with_my_fort90_comp.a is the LAPACK installation compiled with the appropriate Fortran compiler, so the code compiled with g95 is linked to the LAPACK and BLAS libraries compiled with g95, the code compiled with gfortran is linked to the LAPACK and BLAS libraries compiled with gfortran, etc. I consistently used the reference BLAS.
The resulting binary worked with ifort 9.1, but segfaulted on g95 0.9 or gfortran 4.1.1. The platform is CentOS 4.5.
I also tried compiling and running the following code, which is all in one main program:
- Code: Select all
program testInv
implicit none
double precision A(3,3), Ainv(3,3)
integer ipiv(3), iwork(3), info_lu, info_inv, i
A(1,1) = 1d0
A(1,2) = 2d0
A(1,3) = 3d0
A(2,1) = 4d0
A(2,2) = 5d0
A(2,3) = 6d0
A(3,1) = 7d0
A(3,2) = 8d0
A(3,3) = 10d0
write(*,*) 'A = '
do i = 1,3
write(*,*) A(i,1), A(i,2), A(i,3)
end do
Ainv = A
call dgetrf(3, 3, Ainv, 3, ipiv, info_lu)
call dgetri(3, Ainv, 3, ipiv, iwork, 3, info_inv)
write(*,*) 'Ainv = '
do i = 1,3
write(*,*) Ainv(i,1), Ainv(i,2), Ainv(i,3)
end do
end
This segfaults with g95, but not with gfortran and ifort.
Why the difference? Am I just dealing with bugs in g95 and gfortran, or is there something that I'm not seeing in my code that ifort tolerates but the other two compilers don't?
