вектор функций фортран

Question

Недавно я наткнулся на термин «вектор функций», о котором раньше никогда не слышал.Мне нужно определить «вектор функций», чтобы решить систему нелинейных алгебраических уравнений.У Числовых Рецептов есть рутина, названная "тритоном", которая делает это.В C ++, Numeric Recipes определяет новый класс векторов функций для меня, поэтому все, что мне нужно сделать, это вызвать nr3.h, используя предоставленную библиотеку.Кто-нибудь знает, как это сделать в Fortran90 / 95/03/08 (у меня работает любой из этих стандартов)?Я хотел бы знать, потому что я не очень хорошо разбираюсь в C ++ и поэтому предпочел бы работать в Fortran.

Код для C ++ можно найти здесь: http://numerical.recipes/forum/showthread.php?t=1703&highlight=403

Обратите внимание на код функции «VecDoub y (3)».Это не является родным для C ++ (я не думаю).Однако из-за способности в C ++ определять новые классы VecDoub определен в C ++.Есть ли способ сделать это в Фортране?

Спасибо.

Answer 1

Как сказал @Peter в комментарии, в этом "векторе функций" нет ничего особенного.VecDoub немного больше, чем регулярный вектор двойников, а vecfunc на самом деле является функцией векторов.Тогда newt - это просто процедура, которая принимает вектор и функцию векторов.

Ближайший эквивалент в фортране будет выглядеть примерно так:

program vector_functions
  implicit none
  logical :: check
  integer :: i
  double precision :: x(3) = [0, 5, 7], y(3)

  print *, 'Initial guess: ', x
  print *, 'Calling newt...'
  call newt(x, check, vecfunc)
  if(check) then
    print *, 'Check is true (convergence to a local minimum).'
    print *, 'Try another initial guess.'
  else
    print *, 'Check is false (a \"normal\" return).'
  end if
  print *, 'Solution from newt: ', x
  y = vecfunc(x)
  print *, 'Value of the solution vector: ', y

contains
  function vecfunc(x) result(y)
    double precision, intent(in) :: x(:)
    double precision :: y(size(x))
    y(1) = x(1) * x(2) * x(3) - 6
    y(2) = x(1) * x(1) * x(2) + x(2) * x(2) * x(3) + x(3) * x(3) * x(1) - 23
    y(3) = exp(x(1) + x(2) + x(3)) - 403    
  end
  subroutine newt(x, check, f)
    double precision, intent(inout) :: x(:)
    logical, intent(out) :: check
    interface
      function f(a)
        double precision, intent(in) :: a(:)
        double precision :: f(size(a))
      end
    end interface
    ! custom implementation of newt here...
  end
end