, Они требуют двойной точности.

аюсь вычислить обратную сложную матрицу с помощью ZGETRI, но даже если она выполняется без ошибок (info = 0), она не дает мне правильную обратную матрицу, и я абсолютно не знаю, откуда возникла ошибка.

PROGRAM solvelinear
implicit none
INTEGER                        :: i,j,info,lwork
INTEGER,dimension(3)        :: ipiv
COMPLEX(16), dimension(3,3) :: C,Cinv,M,LU
COMPLEX(16),allocatable :: work(:)

info=0
lwork=100
allocate(work(lwork))
ipiv=0
work=(0.d0,0.d0)

C(1,1)=(0.d0,-1.d0)
C(1,2)=(1.d0,5.d0)
C(1,3)=(2.d0,-2.d0)
C(2,1)=(4.d0,-1.d0)
C(2,2)=(2.d0,-3.d0)
C(2,3)=(-1.d0,2.d0)
C(3,1)=(1.d0,0.d0)
C(3,2)=(3.d0,-2.d0)
C(3,3)=(0.d0,1.d0)

write(*,*)"C = "
do i=1,3
   write(*,10)(C(i,j),j=1,3)
end do

!-- LU factorisation
LU=C
CALL ZGETRF(3,3,LU,3,ipiv,info)
write(*,*)'info = ',info
write(*,*)"LU = "
do i=1,3
   write(*,10)(LU(i,j),j=1,3)
end do

!-- Inversion of matrix C using the LU

Cinv=LU
CALL ZGETRI(3,Cinv,3,ipiv,work,lwork,info)
write(*,*)'info = ',info
write(*,*)"Cinv = "
do i=1,3
   write(*,10)(Cinv(i,j),j=1,3)
end do

!-- computation of C^-1 * C to check the inverse
M = matmul(Cinv,C)
write(*,*)"M = "
do i=1,3
   write(*,10)(M(i,j),j=1,3)
end do
          10 FORMAT(3('(',F20.10,',',F20.10,') '))

END PROGRAM solvelinear

Я собираю сifort (и мойLAPACK Библиотеки версии 3.7.1 также скомпилированы с ifort). Makefile:

#$Id: Makefile $
.SUFFIXES: .f90 .f .c .o
FC = ifort
FFLAGS = -g -check all -zmuldefs -i8
LIBS = -L/path/to/lapack-3.7.1 -llapack -ltmglib -lrefblas
MAIN = prog.o
EXEC = xx
all:  ${MAIN} Makefile
    ${FC} ${FFLAGS} -o ${EXEC} ${MAIN} ${LIBS}
.f.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c 
#$Id: Makefile $
.SUFFIXES: .f90 .f .c .o
FC = ifort
FFLAGS = -g -check all -zmuldefs -i8
LIBS = -L/path/to/lapack-3.7.1 -llapack -ltmglib -lrefblas
MAIN = prog.o
EXEC = xx
all:  ${MAIN} Makefile
    ${FC} ${FFLAGS} -o ${EXEC} ${MAIN} ${LIBS}
.f.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<
.f90.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<
lt; .f90.o: ${MODS} Makefile ${FC} ${FFLAGS} -c
#$Id: Makefile $
.SUFFIXES: .f90 .f .c .o
FC = ifort
FFLAGS = -g -check all -zmuldefs -i8
LIBS = -L/path/to/lapack-3.7.1 -llapack -ltmglib -lrefblas
MAIN = prog.o
EXEC = xx
all:  ${MAIN} Makefile
    ${FC} ${FFLAGS} -o ${EXEC} ${MAIN} ${LIBS}
.f.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<
.f90.o: ${MODS} Makefile
    ${FC} ${FFLAGS} -c $<
lt;

У меня нет ошибок при компиляции. Вот мой вывод:

 C = 
(       0.00000,      -1.00000) (       1.00000,       5.00000) (       2.00000,      -2.00000) 
(       4.00000,      -1.00000) (       2.00000,      -3.00000) (      -1.00000,       2.00000) 
(       1.00000,       0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 info =                      0
 LU = 
(       4.00000,       0.00000) (       2.00000,  120470.58824) (       2.00000,      -2.00000) 
(       0.00000,       0.00000) (28003147.29412,      -3.00000) (      -1.00000,       2.00000) 
(       1.00000,       0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 info =                      0
 Cinv = 
(       0.00000,       0.00000) (      -0.00000,      -0.00000) (       2.00000,      -2.00000) 
(      -0.00000,       0.00000) (      -0.00000,      -3.00000) (      -1.00000,       2.00000) 
(      -0.00000,      -0.00000) (       3.00000,      -2.00000) (       0.00000,       1.00000) 
 M = 
(       2.00000,      -2.00000) (       2.00000,     -10.00000) (       2.00000,       2.00000) 
(      -4.00000,     -10.00000) (      -8.00000,       2.00000) (       4.00000,       2.00000) 
(      10.00000,     -10.00000) (       2.00000,     -10.00000) (      -0.00000,       8.00000) 

И М должен быть личностью, если я не ошибаюсь.

 Vladimir F23 окт. 2017 г., 20:00
Пожалуйста, не ставьте теги в заголовок в скобках, как вы это сделали. StackOverflow имеет богатую систему тегов, поместите теги ниже вопроса, где они принадлежат. Если они являются неотъемлемой частью заголовка, они могут оставаться там (например, «в Фортране» или «в функции Фортрана» или аналогично), но не просто размещать там список, список тегов находится ниже вопроса.

Ответы на вопрос(1)

Решение Вопроса

такими какREAL(4) или жеCOMPLEX(16).

Во-первых, это некрасиво и не переносимо.

Во-вторых, это может сбивать с толку для сложных переменных.

Здесь вы определяете свои переменные какCOMPLEX(16), ноZGETRIи все остальноеZ рутины, ожидаетCOMPLEX*16, ЭтоНЕ тот же самый.

COMPLEX*16 это нестандартная запись для комплексных чисел сREAL*8 компоненты.REAL*8 нестандартная запись для 8-байтовых действительных чисел, которые обычно эквивалентныDOUBLE PRECISION.

COMPLEX(16) комплексное число с двумяREAL(16) компоненты, при условии, что такой вид существует. В компиляторах, которые предоставляютREAL(16) это реальная четверная точность, а не двойная точность.

Таким образом, вы эффективно передаете 32-байтовые комплексные переменные вместо 16-байтовых комплексных переменных.

Существует достаточно ресурсов, чтобы научиться правильно использовать виды Фортрана. Вы можете начать с

integer, parameter :: dp = kind(1.d0)

а также

real(dp) :: x
complex(dp) :: c
 Vladimir F24 окт. 2017 г., 14:04
@B_runo Не совсем. Компилятор имел четверную точность, но четверная точность не совместима сZ рутины какZGETRF, Они требуют двойной точности.
 roygvib23 окт. 2017 г., 20:29
Мне интересно, может ли «лучший» подход использовать сложный * 16, а не сложный (kind (1.0d0)), чтобы соответствовать типу, используемому в оригинальном источникеnetlib.org/lapack/explore-3.1.1-html/zgetrf.f.html , но я не уверен, как это определяется в справочных (или последних или mkl или ...) библиотеках на практике.
 Vladimir F23 окт. 2017 г., 20:30
Я упоминаю об этом в ответе, но я не буду рекомендовать этот способ, потому что это не (стандартный) Фортран.
 Vladimir F23 окт. 2017 г., 21:28
@roygvib Основная причина для COMPLEX * 16 заключается в том, что не было стандартного DOUBLE COMPLEX, но DGETRF использует DOUBLE PRECISION, поэтому я бы рекомендовал использоватьkind(1.d0).
 roygvib23 окт. 2017 г., 21:40
Спасибо за вашу информацию, да, я согласен, что вид (1.d0) практически лучший (и я им пользуюсь). Меня беспокоит то, что если для действительного по умолчанию задано значение 64-битное (с помощью параметров компилятора или чего-то еще), возможны проблемы в будущем ... (В этом смысле сложный (real64) из iso_fortran_env может быть безопаснее). Но до этого Лапак сам использует 1.0D0 и сложный * 16 (или действительный * 8) смешанным образом, настолько хаотично ... lol (Но это просто крайне незначительная проблема, поэтому я также считаю сложным (вид (0 .d0)) будет наиболее практичным.)

Ваш ответ на вопрос