Fortran - Moving from Lahey/IMSL to ifc/pgf90

Author

Moving from Lahey/IMSL to ifc/pgf90

Jose L G?mez Dans

2004-07-19, 8:56 am

Hi!
A number of codes have fallen into my hands. These codes were prepared
using lahey and the IMSL libraries on Linux. The version of Lahey that
we have is old, and a pain to install on newer Linux installations, so
we want to use the code with other Linux f90 compilers that we have
around. While the code is mostly well-behaved and standards-compliant,
the calls to libraries aren't, as they assume the IMSL overloads
(i.e., A .x. B to multiply two matrices, FFT(a) to carry out an
FFT...). There are only a few different calls, but they are all over
the place.

So, I could investigate how to re-write these calls, and link against
LAPACK (and FFTPACK/FFTW, as I believe these are the only libraries
used), but I was wondering whether anyone has had any experiences on
this (specifically, I gather that LAPACK/FFTW with the Intel Compiler
is a bit of a black art to set-up!).

Many thanks
Jose

Jason Nielsen

2004-07-19, 8:56 am

On Mon, 19 Jul 2004, Jose L G?mez Dans wrote:

> Hi!
> A number of codes have fallen into my hands. These codes were prepared
> using lahey and the IMSL libraries on Linux. The version of Lahey that
> we have is old, and a pain to install on newer Linux installations, so
> we want to use the code with other Linux f90 compilers that we have
> around. While the code is mostly well-behaved and standards-compliant,
> the calls to libraries aren't, as they assume the IMSL overloads
> (i.e., A .x. B to multiply two matrices, FFT(a) to carry out an
> FFT...). There are only a few different calls, but they are all over
> the place.
>
> So, I could investigate how to re-write these calls, and link against
> LAPACK (and FFTPACK/FFTW, as I believe these are the only libraries
> used), but I was wondering whether anyone has had any experiences on
> this (specifically, I gather that LAPACK/FFTW with the Intel Compiler
> is a bit of a black art to set-up!).
>
> Many thanks
> Jose
>

Hi Jose,

Well here is my suggestion... as always someone else on this list will
likely give better advice ;-)! I would make a module interface so that
you have to change the code you have inherited to the minimum... for
example something like this:

module example
use types ! defines the precision of dp
implicit none

interface operator (.x.)
module procedure mm, vm, mv, vv
end interface

contains

function mm(A,B) result(out)
implicit none
real(dp), dimension(:,:), intent(in):: A, B
real(dp), dimension(size(A,1),size(B,2)):: out
call dgemm_f95(A,B,out)
end function mm

function mv(A,B) result(out)
implicit none
real(dp), dimension(:,:), intent(in):: A
real(dp), dimension(:), intent(in):: B
real(dp), dimension(size(A,1)):: out
call dgemv_f95(A,B,out)
end function mv

function vm(A,B) result(out)
implicit none
real(dp), dimension(:,:), intent(in):: B
real(dp), dimension(:), intent(in):: A
real(dp), dimension(size(B,2)):: out
call dgevm_f95(B,A,out)
end function vm

function vv(A,B) result(out)
implicit none
real(dp), dimension(:), intent(in):: B
real(dp), dimension(:), intent(in):: A
real(dp):: out
out=dot_product(A,B)
end function vv

SUBROUTINE dgemm_f95(a,b,c)
! .. Use Statements ..
! .. Parameters ..
REAL (dp), PARAMETER :: zero = 0.0_dp
REAL (dp), PARAMETER :: one = 1.0_dp
! .. Array Arguments ..
REAL (dp), INTENT (IN) :: a(:,:), b(:,:)
REAL (dp), INTENT (OUT) :: c(:,:)
! .. Local Scalar ..
INTEGER :: k, m, n
! .. External Procedures ..
EXTERNAL dgemm
m = SIZE(c,1)
n = SIZE(c,2)
k = SIZE(a,2)
CALL dgemm('N','N',m,n,k,one,a,m,b,k,zero,c,m
)
END SUBROUTINE dgemm_f95

SUBROUTINE dgemv_f95(a,b,c)
! .. Use Statements ..
! .. Parameters ..
REAL (dp), PARAMETER :: one = 1.0_dp
REAL (dp), PARAMETER :: zero = 0.0_dp
! .. Array Arguments ..
REAL (dp), INTENT (IN) :: a(:,:), b(:)
REAL (dp), INTENT (OUT) :: c(:)
! .. Local Scalar ..
INTEGER :: m, n
! .. External Procedures ..
EXTERNAL dgemv
m = SIZE(a,1)
n = SIZE(a,2)
CALL dgemv('N',m,n,one,a,m,b,1,zero,c,1)
END SUBROUTINE dgemv_f95

SUBROUTINE dgevm_f95(a,b,c)
! .. Use Statements ..
! .. Parameters ..
REAL (dp), PARAMETER :: one = 1.0_dp
REAL (dp), PARAMETER :: zero = 0.0_dp
! .. Array Arguments ..
REAL (dp), INTENT (IN) :: a(:,:), b(:)
REAL (dp), INTENT (OUT) :: c(:)
! .. Local Scalar ..
INTEGER :: m, n
! .. External Procedures ..
EXTERNAL dgemv
m = SIZE(a,1)
n = SIZE(a,2)
CALL dgemv('T',m,n,one,a,m,b,1,zero,c,1)
END SUBROUTINE dgevm_f95

end module example

So in the code you inherited you won't have to change any of the A.x.B
calls as long use you 'use' the module example.... at the same time
leveraging an optimized BLAS such as ATLAS. You can do similar things
with any of the fft libraries so that they match up to calls to IMSL fft
routines. Good luck.

Cheers,
Jason

Sponsored Links