Я хочу вычислить инверсию сложной матрицы с помощью кода FORTRAN и размер матрицы 52 на 52. Я использую функции ZGETRF и ZGETRI в библиотеке MKL для запуска теста. Если размер комплексной матрицы равен 2 на 2, эти две функции хорошо работают, потому что умножение исходной и обратной матриц является единичной матрицей;однако я считаю, что если размер комплексных матриц равен 3 на 3, результаты будут неправильными, поскольку умножение исходной и обратной матриц больше не является единичной матрицей.
Я вставляю свой тестовый код здесь, чтобывы можете сослаться на него.
Я запускаю тестовую программу со сценарием оболочки на кластере, и этот скрипт также вставляется сюда.
PROGRAM TE
IMPLICIT NONE
INTEGER :: i, j, dime
DOUBLE COMPLEX, ALLOCATABLE :: a(:,:)
DOUBLE COMPLEX, ALLOCATABLE :: b(:,:)
DOUBLE COMPLEX, ALLOCATABLE :: c(:,:)
!Variables used in LU decomposition subroutine (ZGETRF)
INTEGER :: nu_r !Number of row in target matrix
INTEGER :: nu_c !Number of column in target matrix
INTEGER :: lu_lda !Leading dimension of target matrix
INTEGER, ALLOCATABLE :: lu_ipiv(:) !Array with dimension of .GE. MIN(nu_r,nu_c)
INTEGER :: lu_info !Judgement (0-successful;<0-error and infor = number, then
!'1-number'th arguement has illegal value;>0 and info =
!number - LU matrix(number,number) is zero)
!
!Variables used in inverse subroutine (ZGETRI)
INTEGER :: nu_o !Order of LU matrix
DOUBLE COMPLEX, ALLOCATABLE :: in_work(:) !Workspace array with dimension of .GE. MAX(1,in_lwork)
INTEGER :: in_lda !Leading dimension of LU decomposed matrix .GE. MAX(1,nu_o)
INTEGER :: in_lwork !Size of in_work array with value .GE. nu_o
INTEGER :: in_info !Judgement (0-successful;<0-error and infor = number, then
!'1-number'th arguement has illegal value;>0 and info =
!number; then, LU matrix(number,number) is zero and inversion
!cannot be accomplished)
dime = 3
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)
Я также тестирую свой код с другими сложными матрицами более высокого порядка, и результаты все неверны. Я вставляю нижеприведенный скрипт оболочки, который я использую для отправки своей работы для запуска этой тестовой программы.
#!/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
Кто-нибудь может сказать, что не так с моим кодом? Что-то не так с функциями ZGETRF и ZGETRI? Кто-нибудь, пожалуйста, дайте мне несколько советов о том, как разобраться? Заранее большое спасибо.