I programmed a function to obtain Cholesky's decomposition of a REAL matrix. The problem is that it returns NaN values in some spots of the matrix. Here is the code:
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
And the program code is:
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" is a module (Jonasson, 2009) for pretty-printing matrices (you can take a look about it here: https://notendur.hi.is/jonasson/greinar/dispmodule-report.pdf ). I discovered that if I use the subroutine DISP in this part:
a=matmul(transpose(a),a)
CALL DISP('ATA =',a,ADVANCE='DOUBLE')
the problem of NaN disappears!(notice that in this last part, there is not "!" before "CALL DISP").
Why is this happening? I am going to use the Cholesky decomposition function for another function, and I can't call DISP all the time to avoid problems. How can I fix this?
Edit: Here's the numerical result without compiling the program with the DISP subroutine:
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
Here is the numerical result when I do compile the DISP subroutine:
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
EDIT 2: I assigned new variables to every new computed matrix, but the problem remains.
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
UPDATE: Minimum compilable example. I took out the subroutine DISP() with its respective module. I also took out the module and function that made the comprobation of the symmetry of the matrix. Moreover, I took out the function that I wrote for reading the matrix from a file (I wrote the same matrix directly on the program) and the function for printing the matrix. The problem gets worse, now there are no numbers, not even the first row as before, all the matrix becomes NaN after the use of cholesky(). The gfortran version is 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
Intel Fortran identified that you are using s incorrectly in Cholesky. The value of s is not defined when you are first using it.
Thus, the result can be anything.
Probably you wanted to set it to 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
The reported code is:
s = 0 !you probably wanted
DO i=2,n,1
DO k=1,i-1,1
s=s+(u(k,i)**2) !<<<=====