Receiving NaN value in Fortran app

Question

I'm developing an Fortran application for numerically solving Boundary Value Problem for second order ODE of the type: -y''+q(x)*y=r(x). In this application I use Gauss-ellimination algorithm to solve the linear system of equations and write the solution in file. But for the solution vector I receive NaN. Why is that happen? Here is some code.

      subroutine gaussian_solve(s, c, error)

    double precision, dimension(:,:), intent(in out) :: s
    double precision, dimension(:),   intent(in out) :: c
         integer        :: error

    if(error == 0) then
        call back_substitution(s, c)
    end if

     end subroutine gaussian_solve
!=========================================================================================

!================= Subroutine gaussian_ellimination ===============================
      subroutine gaussion_ellimination(s, c, error)

    double precision, dimension(:,:), intent(in out) :: s
    double precision, dimension(:),   intent(in out) :: c
    integer,            intent(out)  :: error

    real, dimension(size(s, 1)) :: temp_array
    integer, dimension(1)       :: ksave
    integer                     :: i, j, k, n
    real                        :: temp, m

    n = size(s, 1)

    if(n == 0) then
        error = -1 
        return
    end if

    if(n /= size(s, 2)) then
        error = -2 
        return
    end if

    if(n /= size(s, 2)) then 
        error = -3 
        return
    end if

    error = 0
    do i = 1, n-1
        ksave = maxloc(abs(s(i:n, i)))
        k = ksave(1) + i - 1
        if(s(k, i) == 0) then
            error = -4
            return
        end if

        if(k /= i) then
            temp_array = s(i, :)
            s(i, :) = s(k, :)
            s(k, :) = temp_array
            temp = c(i)
            c(i) = c(k)
            c(k) = temp
        end if

        do j = i + 1, n
            m = s(j, i)/s(i, i)
            s(j, :) = s(j, :) - m*s(i, :)
            c(j) = c(j) - m*c(i)
        end do
    end do

     end subroutine gaussion_ellimination
     !==========================================================================================

   !================= Subroutine back_substitution ========================================
     subroutine back_substitution(s, c)

    double precision, dimension(:,:), intent(in) :: s
    double precision, dimension(:),   intent(in out) :: c

    real    :: w
    integer :: i, j, n

    n = size(c)

    do i = n, 1, -1
        w = c(i)
        do j = i + 1, n
            w = w - s(i, j)*c(j)
        end do
        c(i) = w/s(i, i)
    end do

      end subroutine back_substitution

Where s(i, j) is the matrix of coefficients of the system and c(i) is the solution vector.

Answer 1

You should never write your own routines for gaussian elimination or similar matrix operations. Ubiquitous packages like LAPACK will have faster, more accurate versions than anything you're likely to code yourself; in LAPACK you'd use a combination of _getrf and _getrs for general matricies, but if you have banded or symettric matricies there are special routines for those, too. You should find it quite easy to find and install an optimized linear algebra package for your system. (Look for package names like atlas or flame or gotoblas).

In the code above, there should presumably be a call gaussion_ellimination(s, c, error) near the start of your gaussian_solve routine, your temp_array (and temp and m and w) should also be double precision to avoid losing precision from your double precision matrix, checking for exact equality to floating point zero is a risky strategy, and I'd check your input matrix - if there are any linear degeneracies you will get all NaNs (particularly if it's the last row vector which is linearly degenerate with any of the earlier ones).

If that doesn't fix the problem, you can use signaling NaNs to find out where the problem first arises - Force gfortran to stop program at first NaN - but you're really better off just using existing packages for stuff like this, which have been written by people who have studied the numerical solution of sets of linear equations for years.