如果我不使用打印子例程，则 NaN Cholesky 在 Fortran 中的分解
NaN Cholesky's decomposition in Fortran if I don't use a printing subroutine
我编写了一个函数来获取 REAL 矩阵的 Cholesky 分解。问题是它在矩阵的某些位置 returns NaN 值。这是代码：
MODULE funciones2
IMPLICIT NONE
CONTAINS
FUNCTION cholesky(a)
USE ::comprob !Module to test whether or not a matrix is square and symmetric
IMPLICIT NONE
REAL, ALLOCATABLE :: cholesky(:,:),u(:,:)
REAL :: a(:,:),s
INTEGER :: i,j,k,n
LOGICAL :: comp
comp=symmat(a)
IF (comp) THEN
n=size(a,1)
ALLOCATE (u(n,n))
u=0.0
u(1,1)=sqrt(a(1,1))
DO j=2,n,1
u(1,j)=a(1,j)/u(1,1)
END DO
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2)
END DO
u(i,i)=sqrt(a(i,i)-s)
s=0
DO j=i+1,n,1
DO k=1,i-1,1
s=s+u(k,i)*u(k,j)
END DO
u(i,j)=(1/u(i,i))*(a(i,j)-s)
s=0
END DO
END DO
cholesky=u
ELSE
print*,"CHOLESKY: The matrix is not symmetric."
ENDIF
ENDFUNCTION
ENDMODULE
程序代码为：
PROGRAM pchol
USE :: funciones2
USE :: dispmodule
IMPLICIT NONE
REAL, ALLOCATABLE :: a(:,:),af(:,:)
CHARACTER(LEN=20) :: na
na='B.txt'
a=loadm(na)
CALL DISP('A =',a,ADVANCE='DOUBLE')
a=matmul(transpose(a),a)
!CALL DISP('ATA =',a,ADVANCE='DOUBLE')
af=cholesky(a)
CALL DISP('Ach ',af,ADVANCE='DOUBLE')
af=transpose(af)
CALL DISP('AchT ',af,ADVANCE='DOUBLE')
ENDPROGRAM
"dispmodule" 是一个用于漂亮打印矩阵的模块（Jonasson，2009）（您可以在这里查看：https://notendur.hi.is/jonasson/greinar/dispmodule-report.pdf）。我发现如果我在这部分使用子程序 DISP:
a=matmul(transpose(a),a)
CALL DISP('ATA =',a,ADVANCE='DOUBLE')
NaN的问题消失了！（注意这最后一部分，"CALL DISP"前没有"!"）
为什么会这样？我准备把Cholesky分解函数用于另一个函数，不能一直调用DISP，以免出问题。我该如何解决这个问题？
编辑：
这是没有使用 DISP 子程序编译程序的数值结果：
A =3.00000 2.00000 4.00000 5.00000
2.00000 -3.00000 1.00000 -2.00000
1.00000 1.00000 2.00000 0.00000
Ach 3.74166 0.26726 4.27618 2.93987
0.00000 NaN NaN NaN
0.00000 0.00000 NaN NaN
0.00000 0.00000 0.00000 NaN
AchT 3.74166 0.00000 0.00000 0.00000
0.26726 NaN 0.00000 0.00000
4.27618 NaN NaN 0.00000
2.93987 NaN NaN NaN
下面是我编译 DISP 子程序时的数值结果：
A =3.00000 2.00000 4.00000 5.00000
2.00000 -3.00000 1.00000 -2.00000
1.00000 1.00000 2.00000 0.00000
ATA =14.0000 1.0000 16.0000 11.0000
1.0000 14.0000 7.0000 16.0000
16.0000 7.0000 21.0000 18.0000
11.0000 16.0000 18.0000 29.0000
Ach 3.74166 0.26726 4.27618 2.93987
0.00000 3.73210 1.56940 4.07660
0.00000 0.00000 0.50128 -1.93350
0.00000 0.00000 0.00000 0.00648
AchT 3.74166 0.00000 0.00000 0.00000
0.26726 3.73210 0.00000 0.00000
4.27618 1.56940 0.50128 0.00000
2.93987 4.07660 -1.93350 0.00648
编辑 2：我为每个新计算的矩阵分配了新变量，但问题仍然存在。
PROGRAM pchol
USE :: funciones2
USE :: dispmodule
IMPLICIT NONE
REAL, ALLOCATABLE :: a(:,:),af(:,:),at(:,:),am(:,:),am1(:,:)
CHARACTER(LEN=20) :: na
INTEGER :: m,n
na='B.txt'
a=loadm(na)
!CALL DISP('A =',a,ADVANCE='DOUBLE')
m=size(a,1)
n=size(a,2)
ALLOCATE (at(n,m))
ALLOCATE (am(n,n))
ALLOCATE (af(n,n))
at=transpose(a)
am=matmul(at,a)
am1=am
!CALL DISP('ATA =',am1,ADVANCE='DOUBLE')
af=cholesky(am)
CALL DISP('Ach ',af,ADVANCE='DOUBLE')
af=transpose(af)
CALL DISP('AchT ',af,ADVANCE='DOUBLE')
ENDPROGRAM
更新：最小可编译示例。我取出了子程序 DISP() 及其各自的模块。我还取出了对矩阵对称性进行比较的模块和函数。此外，我还取出了我写的从文件中读取矩阵的函数（我直接在程序上写了相同的矩阵）和打印矩阵的函数。问题变得更糟，现在没有数字，甚至没有像以前那样的第一行，使用 cholesky() 后所有矩阵都变成 NaN。 gfortran 版本是 5.3.0
PROGRAM pcholstack
IMPLICIT NONE
REAL, ALLOCATABLE :: af(:,:),at(:,:),am(:,:),am1(:,:)
REAL :: a(3,4)
INTEGER :: m,n,i,j
a(1,1)=3
a(2,1)=2
a(3,1)=1
a(1,2)=2
a(2,2)=-3
a(3,2)=1
a(1,3)=4
a(2,3)=1
a(3,3)=2
a(1,4)=5
a(2,4)=-2
a(3,4)=1
m=size(a,1)
n=size(a,2)
do i=1,m,1
write(*,1017),(a(i,j),j=1,n)
end do
1017 format (10f15.2)
write(*,*)
ALLOCATE (at(n,m))
ALLOCATE (am(n,n))
ALLOCATE (af(n,n))
at=transpose(a)
am=matmul(at,a)
am1=am
write(*,*) am1(1,1)
write(*,*)
do i=1,n,1
write(*,1018),(at(i,j),j=1,m)
end do
1018 format (10f15.2)
write(*,*)
do i=1,n,1
write(*,1019),(am(i,j),j=1,n)
end do
1019 format (10f15.2)
write(*,*)
af=cholesky(am)
!af=transpose(af)
do i=1,n,1
write(*,1016),(af(i,j),j=1,n)
end do
1016 format (10f15.2)
CONTAINS
FUNCTION cholesky(a)
IMPLICIT NONE
REAL, ALLOCATABLE :: u(:,:),cholesky(:,:)
REAL :: a(:,:),s
INTEGER :: i,j,k,n
n=size(a,1)
ALLOCATE (u(n,n))
u=0.0
u(1,1)=sqrt(a(1,1))
DO j=2,n,1
u(1,j)=a(1,j)/u(1,1)
END DO
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2)
END DO
u(i,i)=sqrt(a(i,i)-s)
s=0
DO j=i+1,n,1
DO k=1,i-1,1
s=s+u(k,i)*u(k,j)
END DO
u(i,j)=(1/u(i,i))*(a(i,j)-s)
s=0
END DO
END DO
cholesky=u
ENDFUNCTION
ENDPROGRAM
我认为您的问题出在这两行中：
a=loadm(na)
a=matmul(transpose(a),a)
我没有 loadm 的代码，但您的 a 矩阵是 4x3。然后你调用 matmul 并将一个 4x4 放入其中（因此使用未为其保留的内存，其他东西可以写入）。 cholesky 函数会将其视为 4x3，而不是 4x4。
更新：
它实际上在gfortran中起作用，它自动分配了变量。但它在 pgfortran.
中失败了
问题所在：cholesky 的 return 值未正确定义。最简单的做法就是 return u:
FUNCTION cholesky(a) result (u)
取出cholesky数组的定义，函数末尾不要赋值
英特尔 Fortran 识别出您在 Cholesky 中使用 s 不正确。 s 的值在您第一次使用时未定义。
因此，结果可以是任何结果。
可能您想将其设置为 0。
> ifort cholesky.f90 -g -check -traceback -standard-semantics
> ./a.out
3.00 2.00 4.00 5.00
2.00 -3.00 1.00 -2.00
1.00 1.00 2.00 1.00
14.00000
3.00 2.00 1.00
2.00 -3.00 1.00
4.00 1.00 2.00
5.00 -2.00 1.00
14.00 1.00 16.00 12.00
1.00 14.00 7.00 17.00
16.00 7.00 21.00 20.00
12.00 17.00 20.00 30.00
forrtl: severe (194): Run-Time Check Failure. The variable 'cholesky$S$_1' is being used in 'cholesky.f90(82,9)' without being defined
Image PC Routine Line Source
a.out 0000000000407DB7 pcholstack_IP_cho 82 cholesky.f90
a.out 0000000000405B47 MAIN__ 54 cholesky.f90
a.out 000000000040261E Unknown Unknown Unknown
libc.so.6 00007F41B73496E5 Unknown Unknown Unknown
a.out 0000000000402529 Unknown Unknown Unknown
报告的代码是：
s = 0 !you probably wanted
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2) !<<<=====
我编写了一个函数来获取 REAL 矩阵的 Cholesky 分解。问题是它在矩阵的某些位置 returns NaN 值。这是代码：
MODULE funciones2
IMPLICIT NONE
CONTAINS
FUNCTION cholesky(a)
USE ::comprob !Module to test whether or not a matrix is square and symmetric
IMPLICIT NONE
REAL, ALLOCATABLE :: cholesky(:,:),u(:,:)
REAL :: a(:,:),s
INTEGER :: i,j,k,n
LOGICAL :: comp
comp=symmat(a)
IF (comp) THEN
n=size(a,1)
ALLOCATE (u(n,n))
u=0.0
u(1,1)=sqrt(a(1,1))
DO j=2,n,1
u(1,j)=a(1,j)/u(1,1)
END DO
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2)
END DO
u(i,i)=sqrt(a(i,i)-s)
s=0
DO j=i+1,n,1
DO k=1,i-1,1
s=s+u(k,i)*u(k,j)
END DO
u(i,j)=(1/u(i,i))*(a(i,j)-s)
s=0
END DO
END DO
cholesky=u
ELSE
print*,"CHOLESKY: The matrix is not symmetric."
ENDIF
ENDFUNCTION
ENDMODULE
程序代码为：
PROGRAM pchol
USE :: funciones2
USE :: dispmodule
IMPLICIT NONE
REAL, ALLOCATABLE :: a(:,:),af(:,:)
CHARACTER(LEN=20) :: na
na='B.txt'
a=loadm(na)
CALL DISP('A =',a,ADVANCE='DOUBLE')
a=matmul(transpose(a),a)
!CALL DISP('ATA =',a,ADVANCE='DOUBLE')
af=cholesky(a)
CALL DISP('Ach ',af,ADVANCE='DOUBLE')
af=transpose(af)
CALL DISP('AchT ',af,ADVANCE='DOUBLE')
ENDPROGRAM
"dispmodule" 是一个用于漂亮打印矩阵的模块（Jonasson，2009）（您可以在这里查看：https://notendur.hi.is/jonasson/greinar/dispmodule-report.pdf）。我发现如果我在这部分使用子程序 DISP:
a=matmul(transpose(a),a)
CALL DISP('ATA =',a,ADVANCE='DOUBLE')
NaN的问题消失了！（注意这最后一部分，"CALL DISP"前没有"!"）
为什么会这样？我准备把Cholesky分解函数用于另一个函数，不能一直调用DISP，以免出问题。我该如何解决这个问题？
编辑： 这是没有使用 DISP 子程序编译程序的数值结果：
A =3.00000 2.00000 4.00000 5.00000
2.00000 -3.00000 1.00000 -2.00000
1.00000 1.00000 2.00000 0.00000
Ach 3.74166 0.26726 4.27618 2.93987
0.00000 NaN NaN NaN
0.00000 0.00000 NaN NaN
0.00000 0.00000 0.00000 NaN
AchT 3.74166 0.00000 0.00000 0.00000
0.26726 NaN 0.00000 0.00000
4.27618 NaN NaN 0.00000
2.93987 NaN NaN NaN
下面是我编译 DISP 子程序时的数值结果：
A =3.00000 2.00000 4.00000 5.00000
2.00000 -3.00000 1.00000 -2.00000
1.00000 1.00000 2.00000 0.00000
ATA =14.0000 1.0000 16.0000 11.0000
1.0000 14.0000 7.0000 16.0000
16.0000 7.0000 21.0000 18.0000
11.0000 16.0000 18.0000 29.0000
Ach 3.74166 0.26726 4.27618 2.93987
0.00000 3.73210 1.56940 4.07660
0.00000 0.00000 0.50128 -1.93350
0.00000 0.00000 0.00000 0.00648
AchT 3.74166 0.00000 0.00000 0.00000
0.26726 3.73210 0.00000 0.00000
4.27618 1.56940 0.50128 0.00000
2.93987 4.07660 -1.93350 0.00648
编辑 2：我为每个新计算的矩阵分配了新变量，但问题仍然存在。
PROGRAM pchol
USE :: funciones2
USE :: dispmodule
IMPLICIT NONE
REAL, ALLOCATABLE :: a(:,:),af(:,:),at(:,:),am(:,:),am1(:,:)
CHARACTER(LEN=20) :: na
INTEGER :: m,n
na='B.txt'
a=loadm(na)
!CALL DISP('A =',a,ADVANCE='DOUBLE')
m=size(a,1)
n=size(a,2)
ALLOCATE (at(n,m))
ALLOCATE (am(n,n))
ALLOCATE (af(n,n))
at=transpose(a)
am=matmul(at,a)
am1=am
!CALL DISP('ATA =',am1,ADVANCE='DOUBLE')
af=cholesky(am)
CALL DISP('Ach ',af,ADVANCE='DOUBLE')
af=transpose(af)
CALL DISP('AchT ',af,ADVANCE='DOUBLE')
ENDPROGRAM
更新：最小可编译示例。我取出了子程序 DISP() 及其各自的模块。我还取出了对矩阵对称性进行比较的模块和函数。此外，我还取出了我写的从文件中读取矩阵的函数（我直接在程序上写了相同的矩阵）和打印矩阵的函数。问题变得更糟，现在没有数字，甚至没有像以前那样的第一行，使用 cholesky() 后所有矩阵都变成 NaN。 gfortran 版本是 5.3.0
PROGRAM pcholstack
IMPLICIT NONE
REAL, ALLOCATABLE :: af(:,:),at(:,:),am(:,:),am1(:,:)
REAL :: a(3,4)
INTEGER :: m,n,i,j
a(1,1)=3
a(2,1)=2
a(3,1)=1
a(1,2)=2
a(2,2)=-3
a(3,2)=1
a(1,3)=4
a(2,3)=1
a(3,3)=2
a(1,4)=5
a(2,4)=-2
a(3,4)=1
m=size(a,1)
n=size(a,2)
do i=1,m,1
write(*,1017),(a(i,j),j=1,n)
end do
1017 format (10f15.2)
write(*,*)
ALLOCATE (at(n,m))
ALLOCATE (am(n,n))
ALLOCATE (af(n,n))
at=transpose(a)
am=matmul(at,a)
am1=am
write(*,*) am1(1,1)
write(*,*)
do i=1,n,1
write(*,1018),(at(i,j),j=1,m)
end do
1018 format (10f15.2)
write(*,*)
do i=1,n,1
write(*,1019),(am(i,j),j=1,n)
end do
1019 format (10f15.2)
write(*,*)
af=cholesky(am)
!af=transpose(af)
do i=1,n,1
write(*,1016),(af(i,j),j=1,n)
end do
1016 format (10f15.2)
CONTAINS
FUNCTION cholesky(a)
IMPLICIT NONE
REAL, ALLOCATABLE :: u(:,:),cholesky(:,:)
REAL :: a(:,:),s
INTEGER :: i,j,k,n
n=size(a,1)
ALLOCATE (u(n,n))
u=0.0
u(1,1)=sqrt(a(1,1))
DO j=2,n,1
u(1,j)=a(1,j)/u(1,1)
END DO
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2)
END DO
u(i,i)=sqrt(a(i,i)-s)
s=0
DO j=i+1,n,1
DO k=1,i-1,1
s=s+u(k,i)*u(k,j)
END DO
u(i,j)=(1/u(i,i))*(a(i,j)-s)
s=0
END DO
END DO
cholesky=u
ENDFUNCTION
ENDPROGRAM
我认为您的问题出在这两行中：
a=loadm(na)
a=matmul(transpose(a),a)
我没有 loadm 的代码，但您的 a 矩阵是 4x3。然后你调用 matmul 并将一个 4x4 放入其中（因此使用未为其保留的内存，其他东西可以写入）。 cholesky 函数会将其视为 4x3，而不是 4x4。
更新：
它实际上在gfortran中起作用，它自动分配了变量。但它在 pgfortran.
问题所在：cholesky 的 return 值未正确定义。最简单的做法就是 return u:
FUNCTION cholesky(a) result (u)
取出cholesky数组的定义，函数末尾不要赋值
英特尔 Fortran 识别出您在 Cholesky 中使用 s 不正确。 s 的值在您第一次使用时未定义。
因此，结果可以是任何结果。
可能您想将其设置为 0。
> ifort cholesky.f90 -g -check -traceback -standard-semantics
> ./a.out
3.00 2.00 4.00 5.00
2.00 -3.00 1.00 -2.00
1.00 1.00 2.00 1.00
14.00000
3.00 2.00 1.00
2.00 -3.00 1.00
4.00 1.00 2.00
5.00 -2.00 1.00
14.00 1.00 16.00 12.00
1.00 14.00 7.00 17.00
16.00 7.00 21.00 20.00
12.00 17.00 20.00 30.00
forrtl: severe (194): Run-Time Check Failure. The variable 'cholesky$S$_1' is being used in 'cholesky.f90(82,9)' without being defined
Image PC Routine Line Source
a.out 0000000000407DB7 pcholstack_IP_cho 82 cholesky.f90
a.out 0000000000405B47 MAIN__ 54 cholesky.f90
a.out 000000000040261E Unknown Unknown Unknown
libc.so.6 00007F41B73496E5 Unknown Unknown Unknown
a.out 0000000000402529 Unknown Unknown Unknown
报告的代码是：
s = 0 !you probably wanted
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2) !<<<=====