2

私は Fortran でプログラミングしており、Lapack パッケージの DGETRI マトリックス インバーターを使用しようとしています。

http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html

しかし、非常に奇妙なことに、それは私のすべての変数をいじっているようです。この非常に単純な例では、DGETRI が A を含まないにもかかわらず、DGETRI が適用されると、プログラムの開始時に初期化された行列 A が変化します。

誰が何が起こっているのか教えてもらえますか? ありがとう！

PROGRAM solvelinear

REAL(8), dimension(2,2)     :: A,Ainv
REAL(8),allocatable         :: work(:)
INTEGER                     :: info,lwork,i
INTEGER,dimension(2)        :: ipiv

info=0
lwork=10000
allocate(work(lwork))
work=0
ipiv=0

A = reshape((/1,-1,3,3/),(/2,2/))
Ainv = reshape((/1,-1,3,3/),(/2,2/))

CALL DGETRI(2,Ainv,2,Ipiv,work,lwork,info)

print*,"A = "
do i=1,2
  print*,A(i,:)
end do

END PROGRAM solve linear

これは出力です：

 A =
   1.0000000000000000        0.0000000000000000
  -1.0000000000000000       0.33333333333333331
4

1 に答える 1

4

DGETRI を呼び出す前に、LU 分解を計算する必要があります。ルーチンの double バージョンを使用しているため、LU を計算する必要がありますDGETRF(ZGETRFは複雑なバージョンです)。

次のコードは、正しい値をコンパイルして返します

PROGRAM solvelinear
implicit none
REAL(8), dimension(3,3)     :: A,Ainv,M,LU
REAL(8),allocatable              :: work(:)
INTEGER                        :: info,lwork
INTEGER,dimension(3)        :: ipiv
INTEGER                        :: i

info=0
lwork=100
allocate(work(lwork))
work=0
ipiv=0

A = reshape((/1,-1,3,3,1,1,1,1,1,1/),(/3,3/))

!-- LU factorisation
LU = A
CALL DGETRF(3,3,LU,3,ipiv,info)

!-- Inversion of matrix A using the LU
Ainv=LU
CALL DGETRI(3,Ainv,3,Ipiv,work,lwork,info)

!-- computation of A^-1 * A to check the inverse
M = matmul(Ainv,A)

print*,"M = "
do i=1,3
  print*,M(i,:)
end do

END PROGRAM solvelinear

出力

M = 
1.0000000000000000        0.0000000000000000        0.0000000000000000     
0.0000000000000000        1.0000000000000000       -5.5511151231257827E-017
-4.4408920985006262E-016  -2.2204460492503131E-016   1.0000000000000000

ちなみに、元のコードは、gfortran 4.8.2 でコンパイルすると、予想される A の変更されていない値を返します。

于 2014-10-21T08:46:19.640 に答える