Ошибка QR-разложения Фортрана - PullRequest
0 голосов
/ 22 мая 2018

У меня проблема с методом разложения QR.Я использую подпрограмму dgeqrf для декомпозиции, но в компиляторе ошибок нет, но после этого возникает проблема.Я не нашел, где ошибка.Другой вопрос: A = Q * R => если матрица имеет ноль, может ли разложение быть нулевым или потерять ранг.

program decomposition

!CONTAINS
!subroutine Qrdecomposition(A_mat, R)
real,dimension(2,2)   :: A_mat    !real,dimension(2,2),intent(inout)   
:: A_mat
real,dimension(2,2)   :: R        !real,dimension(2,2),intent(out)     
:: R
real,dimension(2,2)                  :: A
integer                              :: M,N,LDA,LWORK,INFO
real,allocatable, dimension(:,:)     :: TAU
real,allocatable, dimension(:,:)     :: WORK
external   dgeqrf
M=2
N=2
LDA=2
LWORK=2
INFO=0
A_mat(1,1)=4
A_mat(1,2)=1
A_mat(2,1)=3
A_mat(2,2)=1
A=A_mat

call dgeqrf(M,N,A,TAU,WORK,LWORK,INFO)
R=A
print *,R,WORK,LWORK

!end subroutine Qrdecomposition
end program decomposition

1 Ответ

0 голосов
/ 23 мая 2018

Я вижу три ошибки в вашем коде:

1) Вы забыли аргумент LDA для dgeqrf,

2) TAU и WORK должны быть явно выделены,

3) Все массивы должны быть объявлены с двойной точностью, чтобы соответствовать интерфейсу dgeqrf:

program decomposition

!CONTAINS
!subroutine Qrdecomposition(A_mat, R)
! Note: using '8' for the kind parameter is not the best style but I'm doing it here for brevity.
real(8),dimension(2,2)   :: A_mat    !real,dimension(2,2),intent(inout)
real(8),dimension(2,2)   :: R        !real,dimension(2,2),intent(out)
real(8),dimension(2,2)                  :: A
integer                              :: M,N,LDA,LWORK,INFO
real(8),allocatable, dimension(:,:)     :: TAU
real(8),allocatable, dimension(:,:)     :: WORK
external   dgeqrf
M=2
N=2
LDA=2
LWORK=2
INFO=0
A_mat(1,1)=4
A_mat(1,2)=1
A_mat(2,1)=3
A_mat(2,2)=1
A=A_mat

allocate(tau(M,N), work(M,N))
call dgeqrf(M,N,A,LDA,TAU,WORK,LWORK,INFO)
R=A
print *,R,WORK,LWORK

!end subroutine Qrdecomposition
end program decomposition

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

РЕДАКТИРОВАТЬ Точка 3 была отмечена Ройгвибом, см. ниже.

...