大家好,
我想向大家请教一个关于求逆矩阵的问题。我的矩阵是52行,52列的复数矩阵。我想求它的逆矩阵,我想用MKL函数库里的ZGETRF和ZGEFTI函数来求。为此我先做了一个测试,分别各自求了一个2行2列和一个3行3列的复数矩阵的逆矩阵。我把我的代码贴在下面了。 [Fortran] 纯文本查看 复制代码 PROGRAM TE
IMPLICIT NONE
INTEGER :: i, j, dime !Ingeter for loop and order for complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: a(:,:) !Original complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: b(:,:) !Inversed complex matrix
DOUBLE COMPLEX, ALLOCATABLE :: c(:,:) !Multiplication of the original andinversed complex matrices
!Variables used in LU decomposition subroutine (ZGETRF)
INTEGER :: nu_r !Number ofrow in target matrix
INTEGER :: nu_c !Number ofcolumn in target matrix
INTEGER :: lu_lda !Leadingdimension of target matrix
INTEGER, ALLOCATABLE :: lu_ipiv(:) !Arraywith dimension of .GE. MIN(nu_r,nu_c)
INTEGER :: lu_info !Judgement(0-successful;<0-error and infor = number, then
!'1-number'th arguement hasillegal value;>0 and info =
!number - LU matrix(number,number) is zero)
!Variables used in inverse subroutine (ZGETRI)
INTEGER :: nu_o !Order ofLU matrix
DOUBLE COMPLEX, ALLOCATABLE :: in_work(:) !Workspacearray with dimension of .GE. MAX(1,in_lwork)
INTEGER :: in_lda !Leadingdimension of LU decomposed matrix .GE. MAX(1,nu_o)
INTEGER :: in_lwork !Size ofin_work array with value .GE. nu_o
INTEGER :: in_info !Judgement(0-successful;<0-error and infor = number, then
!'1-number'th arguement hasillegal value;>0 and info =
!number; then, LU matrix(number,number) is zero and inversion
!cannot be accomplished)
dime = 3 !Thisparameter is used to control the size of the complex matrix
ALLOCATE (a(dime,dime))
ALLOCATE (b(dime,dime))
ALLOCATE (c(dime,dime))
DO i = 1, dime, 1
DO j = 1, dime, 1
a(i,j) =CMPLX(DBLE(i), DBLE(j))
END DO
END DO
b = a
OPEN (UNIT=3, FILE='te_in.dat', STATUS='UNKNOWN')
WRITE (UNIT=3, FMT='(A1)') 'a'
DO i = 1, dime, 1
DO j = 1, dime, 1
WRITE (UNIT=3,FMT=*) a(i,j)
END DO
END DO
WRITE (UNIT=3, FMT=*)
!Initialising parameters for LU decomposition subroutine(ZGETRF)
nu_r = dime
nu_c = dime
lu_lda = dime
ALLOCATE (lu_ipiv(dime))
lu_ipiv = 0
!
!Initialising parameters for inverse subroutine (ZGETRI)
nu_o = dime
in_lwork = dime
ALLOCATE (in_work(in_lwork))
in_work = (0.0d0, 0.0d0)
in_lda = dime
!
CALL ZGETRF(nu_r,nu_c,b,lu_lda,lu_ipiv,lu_info)
CALL ZGETRI(nu_o,b,in_lda,lu_ipiv,in_work,in_lwork,in_info)
WRITE (UNIT=3, FMT='(A1)') 'b'
DO i = 1, dime, 1
DO j = 1, dime, 1
WRITE (UNIT=3,FMT=*) b(i,j)
END DO
END DO
WRITE (UNIT=3, FMT=*)
c = MATMUL(a,b)
WRITE (UNIT=3, FMT='(A1)') 'c'
DO i = 1, dime, 1
DO j = 1, dime, 1
WRITE (UNIT=3,FMT=*) c(i,j)
END DO
END DO
WRITE (UNIT=3, FMT=*)
DEALLOCATE (lu_ipiv)
DEALLOCATE (a)
DEALLOCATE (b)
DEALLOCATE (c)
CLOSE (UNIT=3)
STOP
END PROGRAM TE 当测试2行2列的复数矩阵时,结果是正确的,因为原矩阵和它的逆矩阵乘积是一个单位矩阵。结果如下。 a (1.00000000000000,1.00000000000000) (1.00000000000000,2.00000000000000) (2.00000000000000,1.00000000000000) (2.00000000000000,2.00000000000000)
b (-2.00000000000000,2.00000000000000) (2.00000000000000,-1.00000000000000) (1.00000000000000,-2.00000000000000) (-1.00000000000000,1.00000000000000)
c (1.00000000000000,-2.220446049250313E-016) (0.000000000000000E+000,2.220446049250313E-016) (0.000000000000000E+000,4.440892098500626E-016) (1.00000000000000,0.000000000000000E+000)
但是当测试3行3列的矩阵时,结果却是错误的,原矩阵和它的逆矩阵乘积不是单位矩阵。结果如下。 a (1.00000000000000,1.00000000000000) (1.00000000000000,2.00000000000000) (1.00000000000000,3.00000000000000) (2.00000000000000,1.00000000000000) (2.00000000000000,2.00000000000000) (2.00000000000000,3.00000000000000) (3.00000000000000,1.00000000000000) (3.00000000000000,2.00000000000000) (3.00000000000000,3.00000000000000)
b (723205779577744.,262983919846453.) (-1.446411559155488E+015,-525967839692903.) (723205779577744.,262983919846451.) (-1.446411559155489E+015,-525967839692904.) (2.892823118310976E+015,1.051935679385808E+015) (-1.446411559155487E+015,-525967839692903.) (723205779577745.,262983919846452.) (-1.446411559155488E+015,-525967839692905.) (723205779577743.,262983919846452.)
c (0.437500000000000,0.000000000000000E+000) (1.25000000000000,0.000000000000000E+000) (-0.250000000000000,0.000000000000000E+000) (0.125000000000000,0.000000000000000E+000) (0.500000000000000,0.000000000000000E+000) (0.750000000000000,0.000000000000000E+000) (-0.500000000000000,0.000000000000000E+000) (0.000000000000000E+000,-1.00000000000000) (1.50000000000000,0.500000000000000)
我有测试了更高阶的复数矩阵,发现结果也是错误的(原矩阵和它的逆矩阵乘积不是单位矩阵)。我用来编译FORTRAN程序的脚本文件如下。 #!/bin/bash module switch PrgEnv-cray PrgEnv-intel #module load intel #----------------------------------------------------------# ifort -O3 te.f90 \ -Wl,--start-group${MKLROOT}/lib/intel64/libmkl_intel_lp64.a \ ${MKLROOT}/lib/intel64/libmkl_core.a \ ${MKLROOT}/lib/intel64/libmkl_sequential.a \ -Wl,--end-group-lpthread -lm –ldl
我没有找出其中的错误。能否麻烦大家帮我看下,我的代码中哪里出了问题呢?要如何改正呢?
另外,有没有专门计算复数矩阵的逆矩阵的程序,或算法推荐给我呢?多谢大家的帮助啦。 |