1
votes

so reading the following question (Correct use of FORTRAN INTENT() for large arrays) I learned that defining a variable with intent(in) isn't enough, since when the variable is passed to another subroutine/function, it can be changed again. So how can I avoid this? In the original thread they talked about putting the subroutine into a module, but that doesn't help for me. For example I want to calculate the determinant of a matrix with a LU-factorization. Therefore I use the Lapack function zgetrf, but however this function alters my input matrix and the compiler don't displays any warnings. So what can I do?

module matHelper
    implicit none
    contains

    subroutine initMat(AA)
        real*8                      ::  u
        double complex, dimension(:,:), intent(inout)   ::  AA
        integer                     ::  row, col, counter

        counter = 1
        do row=1,size(AA,1)
            do col=1,size(AA,2)
                AA(row,col)=cmplx(counter ,0)
                counter=counter+1 
            end do
        end do

    end subroutine initMat

    !subroutine to write a Matrix to file
    !Input: AA      -   double complex matrix
    !       fid     -   integer file id
    !       fname   -   file name
    !       stat        -   integer status =replace[0] or old[1]
    subroutine writeMat(AA,fid, fname, stat)
        integer                     ::  fid, stat
        character(len=*)                ::  fname
        double complex, dimension(:,:), intent(in)  ::  AA
        integer                     ::  row, col
        character (len=64)                ::  fmtString

        !opening file with given options
        if(fid  /= 0) then
            if(stat == 0) then
                open(unit=fid, file=fname, status='replace', &
                    action='write')
            else if(stat ==1) then
                open(unit=fid, file=fname, status='old', &
                    action='write')
            else
                print*, 'Error while trying to open file with Id', fid
                return
            end if
        end if

        !initializing matrix print format
        write(fmtString,'(I0)') size(aa,2)
        fmtString = '('// trim(fmtString) //'("{",ES10.3, ",", 1X, ES10.3,"}",:,1X))'
        !write(*,*) fmtString

        !writing matrix to file by iterating through each row
        do row=1,size(aa,1)
            write(fid,fmt = fmtString) AA(row,:)
        enddo
        write(fid,*) ''
    end subroutine writeMat



    !function to calculate the determinant of the input
    !Input: AA              -   double complex matrix
    !Output determinantMat  -   double complex, 
    !                           0 if AA not a square matrix
    function determinantMat(AA)
        double complex, dimension(:,:), intent(in)  ::  AA
        double complex              ::  determinantMat
        integer, dimension(min(size(AA,1),size(AA,2)))&
                                    ::  ipiv
        integer                     ::  ii, info

        !check if not square matrix, then set determinant to 0
        if(size(AA,1)/= size(AA,2)) then
            determinantMat = 0
            return
        end if

        !compute LU facotirzation with LAPACK function
        call zgetrf(size(AA,1),size(AA,2), AA,size(AA,1), ipiv,info)

        if(info /= 0) then
            determinantMat = cmplx(0.D0, 0.D0)
            return
        end if
        determinantMat = cmplx(1.D0, 0.D0)
        !determinant of triangular matrix is product of diagonal elements
        do ii=1,size(AA,1)
            if(ipiv(ii) /= ii) then
                !a permutation was done, so a factor of -1 
                determinantMat = -determinantMat *AA(ii,ii)
            else
                !no permutation, so no -1 
                determinantMat = determinantMat*AA(ii,ii)
            end if      
        end do

    end function determinantMat

end module matHelper
!***********************************************************************


!module which stores matrix elements, dimension, trace, determinant

program test
    use matHelper
    implicit none
    double complex, dimension(:,:), allocatable ::  AA, BB
    integer                                 ::  n, fid

    fid  = 0;

    allocate(AA(3,3))
    call initMat(AA)
    call writeMat(AA,0,' ', 0)
    print*, 'Determinante: ',determinantMat(AA) !changes AA
    call writeMat(AA,0, ' ', 0)
end program test

PS: I am using the ifort compiler v15.0.3 20150407

1
is it so hard to read the lapack docs and be aware when the routines you use modify their arguments?agentp
sure it is not that hard, but I was curious and wanted to know whether there was a way to avoid this behavior so I wouldn't have to use a temporary arrayv.tralala
its actually not clear what you are asking. I took it to be "how can I get the compiler to alert me if an external routine changes an intent(in) variable". There seems no way around making a copy, but simply using parentheses (AA) in the sub call should force a copy to be made.agentp
Wow, thank you for the tip. What is it actually doing when calling the lapack function with (AA) in the argument? Why is that making a copy?v.tralala
By my testing the parenthesis force a copy to be made. Honestly I was waiting for some standards expert to chime in on whether that can be counted on by standard.agentp

1 Answers

2
votes

I do not have ifort at home, but you may want to try compiling with '-check interfaces' and maybe with '-ipo'. You may need the path to 'zgetrf' for the '-check interfaces' to work, and if that is not source then it may not help. If you declare 'function determinantMat' as 'PURE FUNCTION determinantMat' then I am pretty sure it would complain because 'zgetrf' is not known to be PURE nor ELEMENTAL. Try ^this stuff^ first.

If LAPACK has a module, then zgetrf could be known to be, or not be, PURE/ELEMENTAL. https://software.intel.com/en-us/articles/blas-and-lapack-fortran95-mod-files

I would suggest you add to your compile line:

-check interfaces -ipo

During initial build I like (Take it out for speed once it works):

-check all -warn all

Making a temporary array is one way around it. (I have not compiled this, so it is only a conceptual exemplar.)

PURE FUNCTION determinantMat(AA)
USE LAPACK95                       !--New Line--!
IMPLICIT NONE                      !--New Line--!
double complex, dimension(:,:)  , intent(IN   )  ::  AA
double complex                                   ::  determinantMat !<- output
!--internals--
integer, dimension(min(size(AA,1),size(AA,2)))   ::  ipiv
!!--Next line is new--
double complex, dimension(size(AA,1),size(AA,2)) ::  AA_Temp  !!<- I have no idea if this will work, you may need an allocatable??
integer                                          ::  ii, info

!check if not square matrix, then set determinant to 0
if(size(AA,1)/= size(AA,2)) then
    determinantMat = 0
    return
end if

!compute LU factorization with LAPACK function
!!--Next line is new--
AA_Temp = AA  !--Initialise AA_Temp to be the same as AA--!
call zgetrf(size(AA_temp,1),size(AA_Temp,2), AA_Temp,size(AA_Temp,1), ipiv,info)

if(info /= 0) then
    determinantMat = cmplx(0.D0, 0.D0)
    return
end if

determinantMat = cmplx(1.D0, 0.D0)
!determinant of triangular matrix is product of diagonal elements
do ii=1,size(AA_Temp,1)
    if(ipiv(ii) /= ii) then
        !a permutation was done, so a factor of -1 
        determinantMat = -determinantMat *AA_Temp(ii,ii)
    else
        !no permutation, so no -1 
        determinantMat = determinantMat*AA_Temp(ii,ii)
    end if      
end do

end function determinantMat

With the 'USE LAPACK95' you probably do not need PURE, but if you wanted it to be PURE then you want to explicitly say so.