2016-11-13 4 views
1

次の質問を読む(Correct use of FORTRAN INTENT() for large arrays)インテント(in)で変数を定義するだけでは不十分だということを知った。その変数が別のサブルーチン/再び変わった。だから私はこれを避けることができますか?元のスレッドでは、サブルーチンをモジュールに入れることについて話しましたが、それは私には役に立ちません。たとえば、LU分解を使って行列の行列式を計算したいとします。したがって、Lapack関数zgetrfを使用しますが、この関数は入力行列を変更し、コンパイラは警告を表示しません。じゃあどうすればいい?インテント(in)で変数を変更できないようにする

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:私は使用していますのifortコンパイラv15.0.3 20150407

+3

lapackのドキュメントを読むのが難しく、使用するルーチンが引数を変更するときは気になるのですか? – agentp

+0

それは難しいとは限りませんが、私は興味があり、この動作を回避する方法があるかどうかを知りたがっているので、一時的な配列を使用する必要はありません –

+1

実際にあなたが求めているものは明確ではありません。私は、 "外部ルーチンが' intent(in) '変数を変更した場合、コンパイラに警告を出す方法を教えてくれました。コピーを作成する方法はないようですが、サブ呼び出しで '(AA)'だけを使用するだけで、コピーが強制的に作成されます。 – agentp

答えて

2

私は自宅でのifortを持っていませんが、「-checkインタフェース」で、多分「-ipoでコンパイルしようとする場合があります' '-check interfaces'が動作するためには、 'zgetrf'へのパスが必要な場合があります。ソースがない場合は役立たない可能性があります。 'function determinantMat'を 'PURE FUNCTION determinantMat'と宣言すると、 'zgetrf'が純粋でもエレメンタリーでもないことが分かっているので、私はそれが本当に不平を言うと確信しています。試してみてください。

LAPACKにモジュールがある場合、zgetrfはPURE/ELEMENTALであるかどうかを知ることができます。私はあなたがコンパイル行に追加することをお勧めhttps://software.intel.com/en-us/articles/blas-and-lapack-fortran95-mod-files

-check interfaces -ipo 

初期構築時には、私は(それが動作後の速度のためにそれを取り出し)好き:

-check all -warn all 

は、一時的な配列である作りますその周りの一方通行。 (私はこれをコンパイルしていないので、それは単に概念的な模範である。)「を使用LAPACK95」で

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 

をあなたはおそらくPURE必要はありませんが、あなたはそれがPUREになりたいならば、あなたは明示的に言いたいですそう。

+0

コピーをうまく動かしても、チェックインターフェイスはうまくいきません。 'ifort:コマンドラインエラー:オプション '-check''のための認識できないキーワード' interfaces 'これで私はコンパイルオプションを変更しましたが、コードを変更しませんでした –

+0

残念ですが' -warn interfaces '私の記憶は厄介なので、あなたはまた、マニュアルを精査する必要があります。 – Holmz

+0

はい。 'AA 'を' intent(in) 'と宣言するだけでは不十分です。関数も純粋なものとして宣言しなければなりません。しかし、それが純粋であると宣言された場合、私は警告オプションを必要とせず、それでも私にメッセージを与えるでしょう。 –

関連する問題