如果我不使用打印子程序，NaN Cholesky 在 Fortran 中的分解

Posted 2023-03-11

技术标签:

【中文标题】如果我不使用打印子程序，NaN Cholesky 在 Fortran 中的分解【英文标题】：NaN Cholesky's decomposition in Fortran if I don't use a printing subroutine 【发布时间】：2017-11-04 22:29:06 【问题描述】：

我编写了一个函数来获得 Cholesky 对 REAL 矩阵的分解。问题是它在矩阵的某些位置返回 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

【问题讨论】：

哪些景点？请阅读How to Ask 和minimal reproducible example 我们需要一个确切的案例并查看完整的输入数据和它们的完整输出。我们需要能够复制您的计算。请注意，您在分配时依赖于 Fortran 2003 自动数组分配。某些编译器默认不启用此功能，尤其是旧版本的 Intel Fortran。你如何编译源代码？您没有显示整个代码，所以我不能 100% 确定，但据我所知，这段代码是非法的 - cholesky 有一个假定形状的虚拟参数，但没有接口调用点的作用域。请显示完整的代码，以便检查。也请使用隐式无。再次因为您没有显示整个程序我无法确定，但似乎程序中没有声明cholesky， @VladimirF，我已经添加了输出计算。我正在使用 gfortran 作为编译器。您应该使用-fcheck=all -Wall -g -fbacktrace 编译并再次运行。只要您怀疑有错误，就应该这样做。 【参考方案1】：

我相信你的问题出在这两行：

a=loadm(na)

a=matmul(transpose(a),a)

我没有loadm 的代码，但您的a 矩阵是4x3。然后您调用matmul 并将4x4 放入其中（因此使用未为其保留的内存，其他东西可以写入）。 cholesky 函数会将其视为 4x3，而不是 4x4。

更新：

它实际上在gfortran 中工作，它自动分配了变量。但它在pgfortran 中失败了。

问题出在：cholesky 的返回值没有正确定义。最简单的方法就是返回u：

FUNCTION cholesky(a) result (u)

取出cholesky数组的定义，不要在函数末尾赋值。

【讨论】：

我认为这也可能是问题所在，所以我在第二次编辑中进行了修复，将新变量分配给每个新的计算矩阵，但问题仍然存在。现在，at 和 am 未分配。我现在分配了，还是失败了发现并修复了另一个问题。稍后会补充回答。 @MSalmer 我认为问题主要有两种可能。一种是您的代码中仍然存在一些错误（从问题中的代码中看不到），以及另一个是您的 gfortran 版本有一些编译器错误，例如，对于具有可分配返回数组的函数。（您可以通过gfortran -v查看版本。）无论如何，我相信如果您准备一个最小可编译示例（不使用“dispmodule”）会非常有用，以便我们可以尝试代码。此外，如上所述，附加 -fcheck=all -Wall -g -fbacktrace 之类的选项可能会提供一些提示。【参考方案2】：

英特尔 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)   !<<<=====

【讨论】：

以上是关于如果我不使用打印子程序，NaN Cholesky 在 Fortran 中的分解的主要内容，如果未能解决你的问题，请参考以下文章

如何在 TensorFlow 中调试 NaN 值？

我不明白如何在 Java 中正确使用 NaN [重复]

线性代数笔记： Cholesky分解

使用 OpenMP 进行 Cholesky 分解

矩阵分解----Cholesky分解

SimplicialLLT 返回错误的cholesky 因子