Fortran subroutine NaN issue when multiple outputs is assigned to the same variable in the main program

I tend to use temporary variables to “ignore” some subroutine outputs. Since Fortran doesn’t have a command like tilde (~) in Matlab, I have to get the outputs fron the subroutine but I assign them to the same temporary variable with the right size. It has been my preference to make things look cleaner, no practical purposes.

So far for years I had no problems with it; today, I am getting NaN. I isolated the problem to a matrix inverse operation in the “r_calc” subroutine but I am not sure what is happening. If I replace that with a transpose (which actually gives the same result for this matrix), the problem disappears.

My question is, is it bad practive to assign multiple subroutine outputs to the same variable in the main code? Why am I having an issue with that when I use matrix inverse? I appreciate any help.

Below is a minimum working example:

program example
implicit none

real(kind=8) :: temp1(3,3)
real(kind=8) :: temp2(3,3)
real(kind=8) :: temp3(3,3)
real(kind=8) :: temp4(3,3)
real(kind=8) :: temp5(3,3)
real(kind=8) :: r_PN_U(3)

call r_calc(4.2d0, &
    0d0, 0d0, &
    0d0, 0d0, &
    0d0, 0d0, &
    0d0,      &
    0d0, 0d0, &
    0d0,      &
    0d0, 0d0, &
    0d0, 0d0, &
    r_PN_U, &
    temp1, temp1, temp1, temp1, temp1)

print *, r_PN_U ! gives NaN in mat_inv_3 is used in r_calc subroutine

call r_calc(4.2d0, &
    0d0, 0d0, &
    0d0, 0d0, &
    0d0, 0d0, &
    0d0,      &
    0d0, 0d0, &
    0d0,      &
    0d0, 0d0, &
    0d0, 0d0, &
    r_PN_U, &
    temp1, temp2, temp3, temp4, temp5)

print *, r_PN_U ! gives correct values even with mat_inv_3

end program example

subroutine r_calc(rin, &
e, ed, &
f, fd, &
v, vd, &
vp, &
w, wd, &
wp, &
theta, thetad, &
y_pos, z_pos, &
r_PN_U, &
C_DU_theta, C_DU_beta, C_DU_zeta, C_DU, C_UD)
implicit none

real(kind=8), intent(in) :: rin

real(kind=8), intent(in) :: e,     ed
real(kind=8), intent(in) :: f,     fd
real(kind=8), intent(in) :: v,     vd
real(kind=8), intent(in) :: vp
real(kind=8), intent(in) :: w,     wd
real(kind=8), intent(in) :: wp
real(kind=8), intent(in) :: theta, thetad

real(kind=8), intent(in) :: y_pos, z_pos

real(kind=8), intent(out) :: r_PN_U(3)

real(kind=8), intent(out) :: C_DU_theta(3,3)
real(kind=8), intent(out) :: C_DU_beta (3,3)
real(kind=8), intent(out) :: C_DU_zeta (3,3)
real(kind=8), intent(out) :: C_DU      (3,3)
real(kind=8), intent(out) :: C_UD      (3,3)

real(kind=8) :: beta, zeta

beta  = -atan(wp) ! [rad], flap down angle
zeta  =  atan(vp) ! [rad], lead angle

call angle2dcm(theta, beta, zeta, C_DU_theta, C_DU_beta, C_DU_zeta, C_DU)

call mat_inv_3(C_DU, C_UD) ! results in NaN in r_PN_U output
! C_UD = transpose(C_DU) ! gives the same result as inverse, eliminates the NaN issue

r_PN_U = [rin+e+f, v, w] + matmul(C_UD, [0d0, y_pos, z_pos])

end subroutine r_calc

subroutine angle2dcm(phi, theta, psi, C_phi, C_theta, C_psi, C_out)
implicit none

! Calculates the direction cosine matrix in psi - theta - phi (3 - 2 - 1) order
! Difference from "angle2dcm" subroutine is the extra outputs

real(kind=8), intent(in)  :: phi, theta, psi

real(kind=8), intent(out) :: C_psi(3,3), C_theta(3,3), C_phi(3,3), C_out(3,3)

C_phi(1,1:3) = [1d0,       0d0,      0d0]
C_phi(2,1:3) = [0d0,  cos(phi), sin(phi)]
C_phi(3,1:3) = [0d0, -sin(phi), cos(phi)]

C_theta(1,1:3) = [cos(theta), 0d0, -sin(theta)]
C_theta(2,1:3) = [       0d0, 1d0,         0d0]
C_theta(3,1:3) = [sin(theta), 0d0,  cos(theta)]

C_psi(1,1:3) = [ cos(psi),  sin(psi), 0d0]
C_psi(2,1:3) = [-sin(psi),  cos(psi), 0d0]
C_psi(3,1:3) = [      0d0,       0d0, 1d0]

C_out = matmul(C_phi, matmul(C_theta,C_psi)) ! psi - theta - phi (3 - 2 - 1) order

end subroutine angle2dcm

subroutine mat_inv_3(A, B)
implicit none

real(kind=8), intent(in)  :: A(3,3)
real(kind=8), intent(out) :: B(3,3)

real(kind=8) :: det

det = 1d0/(A(1,1)*A(2,2)*A(3,3) - A(1,1)*A(2,3)*A(3,2)&
  - A(1,2)*A(2,1)*A(3,3) + A(1,2)*A(2,3)*A(3,1)&
  + A(1,3)*A(2,1)*A(3,2) - A(1,3)*A(2,2)*A(3,1))

B(1,1) = +det * (A(2,2)*A(3,3) - A(2,3)*A(3,2))
B(2,1) = -det * (A(2,1)*A(3,3) - A(2,3)*A(3,1))
B(3,1) = +det * (A(2,1)*A(3,2) - A(2,2)*A(3,1))
B(1,2) = -det * (A(1,2)*A(3,3) - A(1,3)*A(3,2))
B(2,2) = +det * (A(1,1)*A(3,3) - A(1,3)*A(3,1))
B(3,2) = -det * (A(1,1)*A(3,2) - A(1,2)*A(3,1))
B(1,3) = +det * (A(1,2)*A(2,3) - A(1,3)*A(2,2))
B(2,3) = -det * (A(1,1)*A(2,3) - A(1,3)*A(2,1))
B(3,3) = +det * (A(1,1)*A(2,2) - A(1,2)*A(2,1))

end subroutine

Go to Source
Author: Seyhan Gul