[Fortran] Delphi Fortran DLL передача процедуры в качестве параметра

Уважаемые Гуру, очень нужна Ваша помощь.

Стоит задача переделать программу на Фортране в подпрограмму, получить DLL и далее обращаться к ней из Дельфи.

Получилась работающая версия:

      subroutine radau5solver(neqn,t,y,rtol,atol,h0)
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: radau5solver
      !DEC$ ATTRIBUTES ALIAS : "_radau5solver" :: radau5solver
      !DEC$ ATTRIBUTES VALUE :: neqn
      !DEC$ ATTRIBUTES REFERENCE :: t
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: rtol
      !DEC$ ATTRIBUTES REFERENCE :: atol
      !DEC$ ATTRIBUTES VALUE :: h0
 
      integer neqn
      double precision y(md),dy(md),t(0:md),
     +                 h0,rtol(md),atol(md)

type TDArrayR = array of Real;
...
procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayR;
                      y: TDArrayR;
                      rtol: TDArrayR;
                      atol: TDArrayR;
                      h0: Real
                      ); stdcall; external 'radau5solver.dll' name 'radau5solver';
...
radau5solver(neqn,t,y,rtol,atol,h0);

В бывшей программе на Фортране есть основная расчетная подпрограмма:

            call radau5(neqn,odef,t(i),y,t(i+1),h,
     +                  rtol,atol,itol,
     +                  odejac ,ijac,mljac,mujac,
     +                  odemas ,imas,mlmas,mumas,
     +                  solout,iout,
     +                  work,lwork,iwork,liwork,rpar,ipar,idid)

      subroutine odef(n,x,y,f,rpar,ipar)
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      integer ierr
      ierr = 0
      call feval(n,x,y,y,f,ierr,rpar,ipar)
      if (ierr.ne.0) then
         print *, 'RADAU5D: ERROR: ',
     +            'RADAU5 can not handle FEVAL IERR'
         stop
      endif
      return
      end
 
      subroutine feval(neqn,t,y,yprime,f,ierr,rpar,ipar)
      integer neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*)
 
      f(1) = y(2)
      f(2) = ((1-y(1)**2)*y(2)-y(1))/1.0d-6
      
      return
      end

Эти функции представляют собой систему дифференциальных уравнений. В моей основной задаче будут меняться и вид уравнений(от задачи к задаче), и коэффициенты(внутри одной задачи много раз). Оставлять реализацию в Фортране не хотелось бы.
Есть ли какие-нибудь идеи, как этого добиться?

Моя попытка выглядит следующим образом. Пробуем подменить подпрограмму odef из Дельфи.

      interface
      subroutine odef(n,x,y,f,rpar,ipar)
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      end subroutine odef
      end interface
 
c      subroutine odef(n,x,y,f,rpar,ipar)
c      integer n,ipar(*)
c      double precision x,y(n),f(n),rpar(*)
c      integer ierr
c      ierr = 0
c      call feval(n,x,y,y,f,ierr,rpar,ipar)
c      if (ierr.ne.0) then
c         print *, 'RADAU5D: ERROR: ',
c     +            'RADAU5 can not handle FEVAL IERR'
c         stop
c      endif
c      return
c      end

type TDArrayR = array of Real;
     Todef = procedure (
              n: Integer;
              x: Real;
              y: TDArrayR;
              f: TDArrayR;
              rpar: array of Real;
              ipar: array of Integer
              ); stdcall;
 
procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayR;
                      y: TDArrayR;
                      odef: Todef;
                      rtol: TDArrayR;
                      atol: TDArrayR;
                      h0: Real
                      ); stdcall; external 'radau5solver.dll' name 'radau5solver';
 
procedure odef(
              n: Integer;
              x: Real;
              y: TDArrayR;
              f: TDArrayR;
              rpar: array of Real;
              ipar: array of Integer
              ); stdcall;
begin                                  //сюда возвращает, ошибка при переходе к следующему шагу
  ShowMessage('OoOoo');
//  f[0] := y[1];
//  f[1] := ((1-sqr(y[0]))*y[1]-y[0])/0.000001;
end;
 
radau5solver(neqn,t,y,odef,rtol,atol,h0);

При запуске программы получается следующее: запускается подпрограмма, подгружается DLL, выполняются все строчки DLL до call odef(...), далее возврат в Дельфи на begin procedure odef, а на следующем шаге выскакивает ошибка EAccessViolation "Access violation at address XXX in module 'YYY'. Read of address 00000011".
Подскажите, пожалуйста, в чем причина и как заставить работать.
Почему ошибка при переходе к следующему шагу после входа в процедуру?

Так тоже не работает:

      interface
      subroutine odef(n,x,y,f,rpar,ipar)
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: odef
      !DEC$ ATTRIBUTES ALIAS : "_odef" :: odef
      !DEC$ ATTRIBUTES VALUE :: n
      !DEC$ ATTRIBUTES VALUE :: x
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: f
      !DEC$ ATTRIBUTES REFERENCE :: rpar
      !DEC$ ATTRIBUTES REFERENCE :: ipar      
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      end subroutine odef
      end interface

Я ничего не понимаю в Фортране, но кое-что слышал о Delphi.
Во-первых, нужно уточнить соответствия типов и их представлений.
Real в Delphi это устаревший 4-х байтовый тип и судя по описаниям в Фортран-коде как например:

Цитата gimlis @ 24.03.11, 08:26

double precision x,y(n),f(n),rpar(*)

в Delphi-программе ему должен соответствовать тип Double.
Кроме того, конструкция вида "array of тип" в Delphi может означать две вещи:
при описании нового типа означает динамический массив,
при указании в качестве типа формальных параметров подпрограмм - открытый массив.
Это специальные типы Delphi и в формировании интерфейса обмена между модулями из-под компиляторов других языков они не могут.

К сожалению, помочь по существу не могу из-за не знания Фортрана. Советую попросить модераторов отобразить тему в Delphi-раздел.

arj99

Спасибо за отклик.
Соответствия типов проверял: в Delphi 7 Real и Double по 8 байт. На всякий случай заменил на Double.
Динамический массив - фактически указатель на область памяти.

В первой части первого сообщения описана работающая программа с Real и динамическими массивами. К сожалению, надо копать глубже.

Цитата gimlis @ 24.03.11, 10:26

Соответствия типов проверял: в Delphi 7 Real и Double по 8 байт.

Да, прошу прощения, попутал с Single

Цитата gimlis @ 24.03.11, 10:26

Динамический массив - фактически указатель на область памяти.

Вы не до конца правы - это только верхушка айсберга, также динамические массивы Delphi содержат дополнительную информацию, как-то длина и счетчик ссылок.

arj99

Да я это знаю. Смещения -8 и -4 соответственно. Элементы располагаются далее. Фортрану передается указатель на массив и количество элементов (отдельно) + за типами данных требуется следить (за размером).

Одну ошибку нашел. В открытом массиве надо задавать тип передачи данных по ссылке через var, в отличии от динамических массивов, которые сами представляют из себя указатель.

Одна проблема решается так:

procedure odef(
              n: Integer;
              x: Double;
              y: TDArrayD;
              f: TDArrayD;
          var rpar: array of Double;
          var ipar: array of Integer
              ); stdcall;
begin
  ShowMessage('OoOoo');
  ShowMessage(IntToStr(n));
//  f[0] := y[1];
//  f[1] := ((1-sqr(y[0]))*y[1]-y[0])/0.000001;
end;

Теперь программа возвращается в процедуру, выкидывает окошко ОоОоо, далее выводит n, конец процедуры сопровождается invalid pointer operation.

Но еще момент: значение N передается в radau5solver из Delphi. При вызове подпрограммы из DLL N как и требуется равен 2 (передается верно), перед вызом call odef(...) в DLL N также равен 2, но в Delphi передается значение N равное 1637308.

Сначала хотелось бы разобраться с неверной передачей из Фортрана в Дельфи, а далее занимать проблемой указателей/памяти. Есть идеи?

Признаться честно я порядком путаюсь в вашей схеме взаимодействия частей программы. Если не лень - можно отдельным постом полный Delphi-код + аналог фортран-программы в псевдокоде с описанием и там, и там значений параметров в смысле которые должны быть просто параметрами расчета, а которые вернуться результатом выполнения.

Программа Дельфи:

unit TEST;
 
interface
 
uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, ExtCtrls, TeeProcs, TeEngine, Chart, Grids, Series;
 
type
  TForm1 = class(TForm)
    Chart1: TChart;
    Series1: TLineSeries;
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;
 
var
  Form1: TForm1;
 
type TDArrayD = array of Double;
     Todef = procedure (
              n: Integer;
              x: Double;
              y: TDArrayD;
              f: TDArrayD;
          var rpar: array of Double;
          var ipar: array of Integer
              ); stdcall;
 
procedure odef(
              n: Integer;
              x: Double;
              y: TDArrayD;
              f: TDArrayD;
          var rpar: array of Double;
          var ipar: array of Integer
              ); stdcall;
 
procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayD;
                      y: TDArrayD;
                      odef: Todef;
                      rtol: TDArrayD;
                      atol: TDArrayD;
                      h0: Double
                      ); stdcall;
 
implementation
 
{$R *.dfm}
 
procedure odef(
              n: Integer;
              x: Double;
              y: TDArrayD;
              f: TDArrayD;
          var rpar: array of Double;
          var ipar: array of Integer
              ); stdcall;
begin
  ShowMessage('OoOoo');
  ShowMessage(IntToStr(n));
//  f[0] := y[1];
//  f[1] := ((1-sqr(y[0]))*y[1]-y[0])/0.000001;
end;
 
procedure radau5solver(
                      neqn: Integer;
                      t: TDArrayD;
                      y: TDArrayD;
                      odef: Todef;
                      rtol: TDArrayD;
                      atol: TDArrayD;
                      h0: Double
                      ); stdcall; external 'radau5solver.dll' name 'radau5solver';
 
procedure TForm1.FormCreate(Sender: TObject);
var md, neqn: Integer;
    h0: Double;
    t, y, rtol, atol: TDArrayD;
    Fil: TextFile;
    K, i: Integer;
    Outs: array of array of Double;
begin
    md := 32;
    neqn := 2;
    SetLength(t,md+1);
    SetLength(y,md);
    SetLength(rtol,md);
    SetLength(atol,md);
    h0 := 0.000001;
    rtol[0] := 0.000001;
    atol[0] := 0.000000001;
    t[0] := 0;
    t[1] := 2;
    y[0] := 2;
    y[1] := 0;
    radau5solver(neqn,t,y,odef,rtol,atol,h0);
 
    AssignFile(Fil, 'vdpolRADAU5.txt');
    {$I-}                           // îòêëþ÷åíèå êîíòðîëÿ îøèáîê ââîäà-âûâîäà
    Reset(Fil);                     // îòêðûòèå ôàéëà äëÿ ÷òåíèÿ
    {$I+}                           // âêëþ÷åíèå êîíòðîëÿ îøèáîê ââîäà-âûâîäà
    if IOResult <> 0 then           // åñëè åñòü îøèáêà îòêðûòèÿ, òî
      begin
      ShowMessage('Îøèáêà îòêðûòèÿ ôàéëà');
      Exit;                         // âûõîä èç ïðîöåäóðû ïðè îøèáêå îòêðûòèÿ ôàéëà
      end;
 
    K := 0;
    while not Eof(Fil) do
    begin
    SetLength(Outs,3,K+1);
    Readln(Fil, Outs[0,K], Outs[1,K], Outs[2,K]);
    Inc(K);
    end;
    CloseFile(Fil);
 
    for i := 0 to K-1 do
        begin
        Series1.AddXY(Outs[0,i],Outs[1,i]);
        end;
 
end;
 
end.

В Фортране 6000 строчек кода. Попытаюсь показать (уж как смогу), что к чему.

      subroutine radau5solver(neqn,t,y,odef,rtol,atol,h0)
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: radau5solver
      !DEC$ ATTRIBUTES ALIAS : "_radau5solver" :: radau5solver
      !DEC$ ATTRIBUTES VALUE :: neqn
      !DEC$ ATTRIBUTES REFERENCE :: t
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: odef
      !DEC$ ATTRIBUTES REFERENCE :: rtol
      !DEC$ ATTRIBUTES REFERENCE :: atol
      !DEC$ ATTRIBUTES VALUE :: h0
 
      interface
      subroutine odef(n,x,y,f,rpar,ipar)
      !DEC$ ATTRIBUTES STDCALL, DLLEXPORT :: odef
      !DEC$ ATTRIBUTES ALIAS : "_odef" :: odef
      !DEC$ ATTRIBUTES VALUE :: n
      !DEC$ ATTRIBUTES VALUE :: x
      !DEC$ ATTRIBUTES REFERENCE :: y
      !DEC$ ATTRIBUTES REFERENCE :: f
      !DEC$ ATTRIBUTES REFERENCE :: rpar
      !DEC$ ATTRIBUTES REFERENCE :: ipar      
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      end subroutine odef
      end interface
 
      integer md
      parameter (md=32)
      integer lwork, liwork
      parameter (lwork=4*md*md+12*md+20,liwork=3*md+20)
 
      integer neqn,ndisc,mljac,mujac,mlmas,mumas,ind(md),
     +        iwork(20+4*md),ipar(md+1),idid
      double precision y(md),dy(md),t(0:md),
     +                 h0,rtol(md),atol(md),
     +                 work(lwork),rpar(1)
 
..... тело subroutine radau5solver
 
            call radau5(neqn,[B]odef[/B],t(i),y,t(i+1),h,             --> вызов процедуры интегрирования
     +                  rtol,atol,itol,
     +                  odejac ,ijac,mljac,mujac,
     +                  odemas ,imas,mlmas,mumas,
     +                  solout,iout,
     +                  work,lwork,iwork,liwork,rpar,ipar,idid)
 
..... тело subroutine radau5solver
 
      return
      end
 
 
 
      SUBROUTINE RADAU5(N,FCN,X,Y,XEND,H,                            --> процедура интегрирование
     &                  RTOL,ATOL,ITOL,
     &                  JAC ,IJAC,MLJAC,MUJAC,
     &                  MAS ,IMAS,MLMAS,MUMAS,
     &                  SOLOUT,IOUT,
     &                  WORK,LWORK,IWORK,LIWORK,RPAR,IPAR,IDID)
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION Y(N),ATOL(*),RTOL(*),WORK(LWORK),IWORK(LIWORK)
      DIMENSION RPAR(*),IPAR(*)
      LOGICAL IMPLCT,JBAND,ARRET,STARTN,PRED
      EXTERNAL FCN,JAC,MAS,SOLOUT
 
..... тело subroutine radau5solver
 
      CALL FCN(N,X,Y,Y0,RPAR,IPAR)                                  --> вызов процедуры подсчета значений функций
 
..... тело subroutine radau5solver
 
      RETURN
      END

odef в виде подпрограммы на Фортране (в исходном варианте вызывается) выглядит так:

      subroutine odef(n,x,y,f,rpar,ipar)
      integer n,ipar(*)
      double precision x,y(n),f(n),rpar(*)
      integer ierr
      ierr = 0
      call feval(n,x,y,y,f,ierr,rpar,ipar)
      if (ierr.ne.0) then
         print *, 'RADAU5D: ERROR: ',
     +            'RADAU5 can not handle FEVAL IERR'
         stop
      endif
      return
      end
 
      subroutine feval(neqn,t,y,yprime,f,ierr,rpar,ipar)
      integer neqn,ierr,ipar(*)
      double precision t,y(neqn),yprime(neqn),f(neqn),rpar(*)
 
      f(1) = y(2)
      f(2) = ((1-y(1)**2)*y(2)-y(1))/1.0d-6
      
      return
      end

Если что-то добавить надо, пишите.

Не очень понял, что именно ты делаешь, но ты учитываешь, что в фортране аргументы, если не указано иное, передаются по ссылке? Даже константы (для них в этом случае создается временная переменная)
Не знаю точно, что означают строки, начинающиеся с "!DEF$", но похоже описывают способ передачи параметров.
Описание массива y(md) означает массив с нумерацией элементов от 1 до md
Описание t(0:md) нумерация от 0 до md

Всем спасибо за помощь. Проблема решена. В Фортране надо указывать interface (или возможно только атрибут STDCALL) при передаче ссылки во внутреннюю подпрограмму в качестве фактического параметра (в моем случае FCN).

Вообще-то еще следовало бы учитывать, что массивы в Паскале нумеруются с 0, а в Фортране - с 1.
Естественно, в обоих случаях должно быть одно и то же соглашение о вызове (например, stdcall).
Ну и передача по ссылке по умолчанию в Фортране и по значению по умолчанию в Паскале.