library EnDll;

{ Important note about DLL memory management: ShareMem must be the
  first unit in your library's USES clause AND your project's (select
  Project-View Source) USES clause if your DLL exports any procedures or
  functions that pass strings as parameters or function results. This
  applies to all strings passed to and from your DLL--even those that
  are nested in records and classes. ShareMem is the interface unit to
  the BORLNDMM.DLL shared memory manager, which must be deployed along
  with your DLL. To avoid using BORLNDMM.DLL, pass string information
  using PChar or ShortString parameters. }

uses
  SysUtils,
  Classes,
  alap;

{$R *.RES}
Const Tiny=1.0E-20;

Type
        frissito=Record
                LU: array of array of RealType;         //LU felbontas matrixa
                BG: array of array of RealType;         //a Bartels-Golub frissites
                OBL:array of integer;                   //a bazis, amire LU felbontas keszul
                MBL:array of integer;                   //OBL bazis az aktualis U-nak megfelelo sorrendben
                LUP:array of integer;                   //az LU felbontas kozben torteno perturbaciok
                End;

var     valt:frissito;
        size:integer;

Procedure InterChangeRows(R1,R2:integer);
        var j:integer; w:RealType;
        Begin
        for j:=1 To size Do
                Begin
                        w:=valt.LU[R1,j];
                        valt.LU[R1,j]:=valt.LU[R2,j];
                        valt.LU[R2,j]:=w;
                End;
End;

Procedure ChangeRows(R1,R2:integer);     //LU-m sorcsereje, az LU-nak csak az U-hoz tartozo elemeit csereli
        var j:integer; w:RealType;
        Begin
        for j:=R1 To size Do
                Begin
                        w:=valt.LU[R1,j];
                        valt.LU[R1,j]:=valt.LU[R2,j];
                        valt.LU[R2,j]:=w;
                End;
End;

procedure LU_Init(m:integer;bl:array of integer); stdcall;
        var i,j,k:integer;
        IMax:Integer;
        {Dum,}Big{,S}:Extended;


        Begin   //---------------------------------

        SetLength(valt.LU,m);
        SetLength(valt.LUP,m-1);
        SetLength(valt.BG,m);

        for i:=1 to m do
                Begin
                SetLength(valt.LU[i],m);
                SetLength(valt.BG[i],m);
                if i<>m then
                        Begin
                        valt.LUP[i]:=-1;
                        End;
                End;

        SetLength(valt.OBL,m);
        SetLength(valt.MBL,m);
        size:=m;

        for i:=1 to m do
                Begin
                valt.OBL[i]:=bl[i];
                valt.MBL[i]:=bl[i];
                End;

//LU matrix feltoltese az aktualis bazisoszlopokkal
        for i:=1 to m do
                for j:=1 to m do
                        Begin
                        valt.LU[i,j]:=OriginalProblem.Matrix[j,valt.OBL[i]];
                        if i=j then
                                valt.BG[i,j]:=1
                        else
                                valt.BG[i,j]:=0;
                        End;

// LU felbontas elkeszitese-----------------------A. verzio
     {   For j:=1 To m Do
                Begin
                        For i:=1 To j-1 Do
                                Begin
                                        S:=valt.LU[i,j];
                                        For k:=1 To i-1 Do
                                                S:=S-(valt.LU[i,k]*valt.LU[k,j]);
                                        valt.LU[i,j]:=S;
                                End;    //for i

                        Big:=0.0;
                        IMax:=j;

                        For i:=j To m Do
                                Begin
                                        S:=valt.LU[i,j];
                                        For k:=1 To j-1 Do
                                                S:=S-(valt.LU[i,k]*valt.LU[k,j]);
                                        valt.LU[i,j]:=S;
                                        Dum:=Abs(S);
                                        If Dum>=Big Then
                                                Begin
                                                        Big:=Dum;
                                                        IMax:=i;
                                                End;
                                End;    //for i

                        If (j<>IMax) Then
                                Begin
                                        InterChangeRows(j,IMax);
                                        valt.LUP[j]:=IMax;
                                End;
                        If valt.LU[j,j]=0.0 Then
                                valt.LU[j,j]:=Tiny;
                        If j<>m Then
                                Begin
                                        Dum:=1/valt.LU[j,j];
                                        For i:=j+1 To m Do
                                                valt.LU[i,j]:=valt.LU[i,j]*Dum;
                                End;
                End;    //for j
                    }
// LU felbontas kesz---------------------------A. verzio

// LU felbontas elkeszitese-----------------------B. verzio

        IMax:=0;
        for i:=1 to m-1 do
                begin
                        if valt.LU[i,i]=0 then
                                begin
                                        Big:=0.0;
                                        for j:=i+1 to m do
        //legnagyobb abszolút értékű elem megkeresése az adott oszlopban
                                                if Abs(valt.LU[j,i])>Big then
                                                        begin
                                                                Big:= Abs(valt.LU[j,i]);
                                                                IMax:=j;
                                                        end;
                                        if Big<>0.0 then
                                                begin
                                                        ChangeRows(i,IMax);
                                                        valt.LUP[i]:=IMax;
                                                end;

                                end;

                        If valt.LU[i,i]=0.0 Then
                                valt.LU[i,i]:=Tiny;

                        for j:= i+1 to m do
                                begin
                                        if valt.LU[i,i]<>0 then
                                                valt.LU[j,i]:=valt.LU[j,i]/valt.LU[i,i];
                                        for k:=i+1 to m do
                                                valt.LU[j,k]:=valt.LU[j,k]-(valt.LU[j,i]*valt.LU[i,k]);
                                end;
                end;

// LU felbontas kesz---------------------------B. verzio

        End;   // LU_Init vege---------------------------


procedure LU_Frissit(mode:integer;var be:array of integer;var ki:array of integer;n:integer); stdcall;
        var i,j,k,l,PR:integer;
        seged:RealType;
        U:array of array of RealType;
        Count:array of RealType;
        f:array of RealType;

        Begin;
        SetLength(U,size);
        SetLength(Count,size);
        SetLength(f,size);
        PR:=size;

        for i:=1 to size do
                SetLength(U[i],size);

        for i:=1 to size do      // az U matrix masolata
                for j:= 1 to size do
                Begin
                        if i<=j then
                                U[i,j]:=valt.LU[i,j]
                        else  U[i,j]:=0;
                End;

if mode=1 then  //ha a mode 1, akkor Bartels-Golub frissites lesz
        BEGIN
        for i:=1 to n do   //a frissites ciklusa
                Begin

                for j:=1 to size do
                        f[j]:=OriginalProblem.Matrix[j,be[i]];

                for k:=1 to size-1 do  // az f uj oszlopvektor szorzasa L-lel
                       begin
                        if valt.LUP[k]<>-1 then
                                Begin
                                seged:=f[k];
                                f[k]:=f[valt.LUP[k]];
                                f[valt.LUP[k]]:=seged;
                                End;
                        for l:=1 to size do
                                begin
                                count[l]:=0;
                                for j:= 1 to size do
                                        begin
                                        if l=j then
                                                count[l]:=count[l] + f[l];
                                        if l>j then
                                                begin
                                                if j=k then
                                                        count[l]:=(-1*valt.LU[l,j])*f[j];
                                                end;
                                        end;
                                end;
                        for l:=1 to size do
                                begin
                                        f[l]:= count[l];
                                end;
                        end;

                for j:=1 to size do     // az uj bazisoszlop szorzasa a Bartels-Golub frissitesekkel
                        Begin
                        Count[j]:=0;
                        for k:= 1 to size do
                                Begin
                                if valt.BG[j,k]<>0 then
                                        Count[j]:=Count[j]+f[k]*valt.BG[j,k];
                                End;
                        End;

                for j:=1 to size do     // az eredeti sorrendu bazisleiras frissitese, es a kimeno
                        Begin           //oszlop bazisbeli indexenek megkeresese (PR)
                        if valt.OBL[j]=ki[i] then
                                valt.OBL[j]:=be[i];
                        if valt.MBL[j]=ki[i] then
                                PR:=j;

                        End;

                if PR<>size then        //az uj bazisoszlop beszurasa az uj helyre
                        Begin           //modositott bazisleiras frissitese (bazisoszlopok indexenek valos sorrendje)
                        for j:=PR+1 to size do
                                for k:=1 to size do
                                        Begin
                                        U[k,j-1]:=U[k,j];
                                        valt.MBL[j-1]:=valt.MBL[j];
                                        End;
                        End;
                for j:=1 to size do
                        Begin
                        U[j,size]:=Count[j];
                        valt.MBL[size]:=be[i];
                        End;

                for j:=1 to size-1 do //az U matrix atalakitasa, es a Bartels-Golub frissitese
                        Begin
                        if U[j,j]=0 then    //sorcsere
                                Begin
                                //seged:=0.0;
                                for k:=j to size do
                                        Begin
                                        seged:=U[j,k];
                                        U[j,k]:=U[j+1,k];
                                        U[j+1,k]:=seged;
                                        //seged:=0.0;
                                        seged:=valt.BG[j,k];
                                        valt.BG[j,k]:=valt.BG[j+1,k];
                                        valt.BG[j+1,k]:=seged;
                                        End;
                                End;
                        if U[j+1,j]<>0 then    //atlo alatti elem eltuntetese
                                Begin
                               // seged:=0.0;
                                seged:=-1*U[j+1,j]/U[j,j];
                                for k:=1 to size do
                                        Begin
                                        U[j+1,k]:=seged*U[j,k]+U[j+1,k];
                                        valt.BG[j+1,k]:=seged*valt.BG[j,k]+valt.BG[j+1,k];
                                        End;
                                End;
                        End;
                End;  // a frissites ciklusanak vege---- i
        END; //mode vege

         for i:=1 to size do      // az U matrix frissitese az LU-ban
                for j:=1 to size do
                Begin
                        if i<=j then
                                valt.LU[i,j]:=U[i,j];
                End;

        End;

procedure SimTabColJav(var tomb:array of RealType;k:Integer); stdcall;
        var i,j:integer;
        seged:array of RealType;

        Begin  //-------------------------------------

        SetLength(seged,size);

        for i:=1 to size do
                Begin
                        tomb[i]:=OriginalProblem.Matrix[i,k] //a tomb-be teszem a k. oszlopvektort
                End;
// L * y = A k. oszlopa
        for i:=2 to size do
                for j:=1 to i-1 do
                Begin
                        tomb[i]:=tomb[i]-tomb[j]*valt.LU[i,j];
                End;

//az y frissitese a Bartels-Golub szerint
        for i:=1 to size do
                Begin
                seged[i]:=0;
                for j:=1 to size do
                        Begin
                        if valt.BG[i,j]<>0 then
                                seged[i]:=tomb[j]*valt.BG[i,j]+seged[i];
                        End;
                End;

// U * x = y kiszámolása
        seged[size]:=seged[size]/valt.LU[size,size];
        for i:=size-1 downto 1 do
                Begin
                for j:=size downto i+1 do
                        seged[i]:=seged[i]-seged[j]*valt.LU[i,j];
                seged[i]:=seged[i]/valt.LU[i,i];
                End;

//a tomb elemeit helyes sorrendbe kell tenni
        for i:=1 to size do
                for j:=1 to size do
                        if valt.MBL[j]=valt.OBL[i] then
                                tomb[i]:=seged[j];

        End;  //  SimTabColJav vege--------------------------

exports
        LU_Init,
        LU_Frissit,
        SimTabColJav;

begin
end.