Получение значения NaN в приложении Fortran - PullRequest
Новая
       35

Получение значения NaN в приложении Fortran

1 голос
/ 16 декабря 2011

Я разрабатываю приложение на Фортране для численного решения краевой задачи для ODE второго порядка типа: -y '' + q (x) * y = r (x). В этом приложении я использую алгоритм исключения Гаусса, чтобы решить линейную систему уравнений и записать решение в файл. Но за вектор решения я получаю NaN. Почему это случилось? Вот некоторый код. 

      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

Где s (i, j) - матрица коэффициентов системы, а c (i) - вектор решения.

1 Ответ

10 голосов
/ 16 декабря 2011

Вы должны никогда написать свои собственные процедуры для исключения по Гауссу или аналогичные матричные операции.Вездесущие пакеты, такие как LAPACK, будут иметь более быстрые и точные версии, чем все, что вы, вероятно, будете кодировать самостоятельно;в LAPACK вы бы использовали комбинацию _getrf и _getrs для общих матриц, но если у вас есть полосчатые или симметричные матрицы, для них тоже есть специальные подпрограммы.Вам должно быть довольно легко найти и установить оптимизированный пакет линейной алгебры для вашей системы.(Ищите названия пакетов, такие как атлас, пламя или готоблас).

В приведенном выше коде должно быть call gaussion_ellimination(s, c, error) около начала вашей gaussian_solve подпрограммы, ваши temp_array (и temp и m и w) также должны иметь двойную точность, чтобы избежать потери точности из вашей матрицы двойной точности, проверка на точное равенство нулю с плавающей запятой - рискованная стратегия, и я бы проверил вашу входную матрицу - если естьпри любых линейных вырождениях вы получите все NaN (особенно, если это последний вектор строки, линейно вырожденный с любым из предыдущих).

Если это не решит проблему, вы можете использовать сигнальные NaN, чтобы найтитам, где впервые возникает проблема - Вынудите gfortran сначала остановить программу NaN - но вам действительно лучше просто использовать существующие пакеты для подобных вещей, которые были написаны людьми, которые изучали численное решениемножеств линейных уравнений по годам.

...