unit Linalg;
{=================================================================}
{  Версия файла 3.0.000                                           }
{  от 4 декабря 2002 г.                                           }
{==========================} interface {==========================}

uses
  SysUtils,WinTypes, WinProcs,
  gltype,glproc,MathObj;

const
  cwr  = 4;  cwi  = 4;  cwra = 3;  cwia = 4;

type
  TWR  = array[1..cwr]  of TExtArray;
  TWI  = array[1..cwi]  of TIntArray;
  TWRA = array[1..cwra] of TPtrExt;
  TWIA = array[1..cwia] of TPtrInt;

  TTridag = class (TObject)
   N            : Integer;        {число уравнений }
   WA           : array[1..6] of PExtArr;
   constructor Create(an:Integer);
   destructor  Destroy;override;
   procedure   ChangeNOfE(an:Integer);
   procedure   Progonka;
  end;

  TSparce = class (TObject)
   Matrix_flag  : Boolean;
   Iter_flag    : Boolean;
   NOfE         : Integer;        {число уравнений }
   WR           : TWR;
   WI           : TWI;
   WRA          : TWRA;
   WIA          : TWIA;
   TB           : RealType; {барьер}
   TU           : RealType; {коэффициент численной устойчивости}
   Rows         : Integer; {число просматриваемых строк}

   constructor Create(n:Integer;fm,fi:Boolean);
   destructor  Destroy;override;
   procedure   ChangeNOfE(N : Integer);

   procedure   Init(n:Integer);
   function    Run(n:Integer):Integer;
   function    RunIter(n:Integer;hb,relerr:RealType):Integer;
   function    LUDecomp(n:Integer;flag:Boolean):Integer;
   procedure   SetMatrixMemo(Row,cRow : Integer);
   procedure   ChangeMemoryAlloc(n:Integer);
   procedure   GetRightPart(n:Integer);
   procedure   FillRow(n,ARow:Integer;var X:PExtArr;var ACol:PIntArr);
   procedure   FillItem(ARow,ACol:Integer;X: RealType);
   procedure   FillRightItem(ACol:Integer;X: RealType);
   procedure   FillRightCol(n:Integer;var X: PExtArr);
   procedure   GetRightItem(ACol:Integer;var X: RealType);
   procedure   GetSolutionItem(ACol:Integer;var X: RealType);
   function    GetSolution:PExtArr;
   procedure   GetSolutionCol(n:Integer;var X: PExtArr);
   procedure   GetSolutionError(n:Integer;var dX,dF: PExtArr);
   procedure   GetSolutionErrorItem(ACol:Integer;var X: RealType);
   function    GetDet:RealType;
  end;

function SparceLUDecomp(NOfE,Rows:Integer;var WR:TWR;var Wi:TWI;
          var U,L:TPtrExt;var CU,RU,CL:TPtrInt;TB,TU:RealType;first:Boolean):Integer;
function  SparceSolve(NOfE:Integer;var X,F:array of RealType;var WI:TWI;
          var U,L:TPtrExt;var CU,RU,CL:TPtrInt):Integer;
function  SparceIter(NOfE:Integer;var X,Xr,Fr,F:array of RealType;var WI:TWI;
          var U,US,L:TPtrExt;var CU,RU,CUS,CL:TPtrInt;eps:RealType):Integer;

FUNCTION  GaussLUfact(n: integer;VAR a: PExtPtrArr; VAR indx: PIntArr;
                 VAR vv: PExtArr;VAR ner:Integer):RealType;
PROCEDURE GaussSolve(n: integer;a: PExtPtrArr;indx: PIntArr; VAR b,x: PExtArr);

function  popolam(x0:RealType;const Pt:array of Pointer;fun:TFun):RealType;


var  SAE    : TSparce;

{----------------------------------------------------------------------}
                       IMPLEMENTATION
{--------------------------------------------------------------------------}
uses math;
{const Counter : Integer = 0;}

{процедура выполнения k-го шага исключения Гаусса}
function   ElimStep(var U,L:TPtrExt;var CU,RU,CL:TPtrInt;var Wi:TWi;ind:Integer;TB:RealType):Integer;
var
  akk          : RealType;
  aii,aik,aki  : RealType;
  aa,akkmax    : RealType;
  i,i1,i2,i3   : Integer;
  ik,ki,ii     : Integer;
  j,j1,jk      : Integer;
  countk,counti: Integer;
  countr       : Integer;
  iflag        : Boolean;
  k,m,n        : Integer;
begin
   Result:=0;

   k:=Wi[1].arr^[ind]-1;//номер ведущей строки
   m:=Wi[2].arr^[ind]-1;//номер ведущего столбца
   n:=Wi[3].arr^[ind];  //индекс ведущего столбца в массиве U[k]}

   if (n < 0) then exit;

   akk      :=U.val(k,n);{главный член к-го шага }

   if akk = 0.0 then begin
    Result:=er_NullDet;
    exit
   end;

   countk   :=U.Ptr(k).count;{число ненулевых членов в к-ом уравнении}
   countr   :=RU.Ptr(m).count;{число ненулевых членов в m-ом столбце матрицы}
   if (countk = 0) or (countr = 0) then exit;

   L.Ptr(k).changecount(countr-1);
   CL.Ptr(k).changecount(countr-1);

   ki:=0;
   for i:=0 to countr-1 do begin
    iflag    :=false;
    ii       :=RU.val(m,i)-1;{номер i-го уравнения}
    if ii = k then continue;
    counti   :=U.Ptr(ii).count;    {число ненулевых членов в i-ом уравнении}

{Найти элемент в ii-ом уравнении , находящийся в столбце m}
    for i1:=0 to counti-1 do
    if CU.val(ii,i1) = m+1 then begin
     ik :=i1;
     aik:=U.val(ii,ik);
     break
    end;
    aa  :=-aik/akk;        {множитель для ii-го уравнения}
    L.z(k,ki,aa);          {элемент нижней треугольной матрицы (k,ii)}
    CL.z(k,ki,ii+1);       {индекс элемента L(k,ii)}
    inc(ki);
   {цикл по всем ненулевым элементам ведущей (k-ой) строки}
    for i1:=0 to countk-1 do begin
      if i1 = n then continue;
      i3:=0;
      {найти элемент в ii-ой строке, стоящий в том же столбце,
       что и текущий элемент k-ой строки; если найден, то заполнения
       нет, изменить значение элемента ii-ой строки,
       если не найден, то есть заполнение}
      for i2:=0 to counti-1 do
      if CU.val(k,i1) = CU.val(ii,i2) then begin
        U.z(ii,i2,U.val(ii,i2) + U.val(k,i1)*aa);
        i3:=0;break
      end else i3:=1;

      {есть заполнение, если это первое заполнение
      для ii-ой строки, оно помещается в исключаемый элемент ii-ой
      строки, если не первое, происходит выделение дополнительной памяти }
      if i3 = 1 then
      {если заполнение меньше барьера, заполнения нет}
       if ABS(U.val(k,i1)*aa) >= TB then
        if iflag then begin
         j1:=U.Ptr(ii).count;
         U.Ptr(ii).changecount(j1+1);
         CU.Ptr(ii).changecount(j1+1);
         RU.Ptr(CU.val(k,i1)-1).changecount(RU.Ptr(CU.val(k,i1)-1).count+1);
         U.z(ii,j1,U.val(k,i1)*aa);
         CU.z(ii,j1,CU.val(k,i1));
         RU.z(CU.val(k,i1)-1,RU.Ptr(CU.val(k,i1)-1).count-1,ii+1);
        end
        else begin
         iflag:=true;
         RU.Ptr(CU.val(k,i1)-1).changecount(RU.Ptr(CU.val(k,i1)-1).count+1);
         U.z(ii,ik,U.val(k,i1)*aa);
         CU.z(ii,ik,CU.val(k,i1));
         RU.z(CU.val(k,i1)-1,RU.Ptr(CU.val(k,i1)-1).count-1,ii+1);
        end;
    end;

    if not iflag then begin
     U.Ptr(ii).AtDelete(ik,ik); //удалить элемент с индексом ik
     CU.Ptr(ii).AtDelete(ik,ik);//удалить элемент с индексом ik
    end;
   end;

   for i:=0 to countk-1 do begin
    ki:=CU.val(k,i)-1;
    for i1:=0 to RU.Ptr(ki).count-1 do
    if RU.val(ki,i1) = k+1 then begin
     RU.Ptr(ki).AtDelete(i1,i1);
     break;
    end;
   end;
   RU.Ptr(m).changecount(0);
end;
{----------------------------------------------------------------------
 процедура выбора главного элемента на k-м шаге исключения
-----------------------------------------------------------------------}
function SparceSort(Ind,NOfE,p:Integer;var Wi:TWI;var U:TPtrExt;
          var CU,RU:TPtrInt;TU:RealType):Integer;
var
   ir,ic,r   : Integer;
   i,j,k,m   : Integer;
   amax,xmax : RealType;
   n,ii,marc : Integer;
   row       : Integer;
begin
 Result:=0;
 {число просматриваемых строк}
 n:=min(p,NofE-Ind);
 {отсортировать оставшиеся после k-го шага уравнения так, чтобы n первых
  уравнений имели минимальное число ненулевых элементов}
 for j:=0 to n-1 do begin
  r:=MaxInt;
  row:=Ind+j;
  for i:=Ind+j to NofE-1 do begin
   k:=CU.Ptr(Wi[1].arr^[i]-1).Count;//число Н.Н.Э. в строке row
   if k < r then begin
    r:=k;
    row:=i
   end
  end;
  m:=Wi[1].arr^[Ind+j];
  Wi[1].arr^[Ind+j]:=Wi[1].arr^[row];
  Wi[1].arr^[row]:=m;
 end;

 r:=MaxInt;
 for j:=Ind to Ind+n-1 do begin
  row:=Wi[1].arr^[j]-1;

  m:=CU.Ptr(row).Count-1;
  if m < 0 then begin {число ннэ = 0}
    Result:=er_NullDet;
    ii:=-1;
    ic:=Wi[2].arr^[ind-1]-1;  //номер ведущего столбца;
    ir:=j;
    break
  end;

  {найти для каждой из n просматриваемых строк максимальный по значению элемент}
  amax:=0;
  for i:=0 to m do amax:=max(amax,ABS(U.val(row,i)));
  {если ненулевой член в просматриваемой строке удовлетворяет
   условию устойчивости, то выбрать или не выбрать его в качестве главного
   по критерию Марковица}
  for i:=0 to m do begin
   k:=CU.val(row,i)-1;//номер неизвестной

   marc:=m*(RU.Ptr(k).Count-1);//цена Марковица
   if ABS(TU*U.val(row,i)) >= amax then {удовлетворяет устойчивости}
   if marc < r then begin
   {если цена Марковица элемента меньше текущей минимальной,
    то выбрать данный элемент в качестве ведущего}
    r:=marc;{минимальная найденная цена Марковица}
    ii:=i;{индекс ведущего элемента в строке}
    ic:=k;{номер ведущего столбца в строке}
    ir:=j;{номер ведущей строки}
    xmax:=ABS(U.val(row,i));{максимальный элемент среди просматриваемых строк}
   end
   else if marc = r then begin
   {если цена Марковица равна текущей минимальной цене, но
    элемент больше по значению, чем текущий ведущий, то
    выбрать данный элемент в качестве ведущего}
    if ABS(U.val(row,i)) > xmax then begin
     r:=marc;{минимальная найденная цена Марковица}
     ii:=i;{индекс ведущего элемента в строке}
     ic:=k;{номер ведущего столбца в строке}
     ir:=j;{номер ведущей строки}
     xmax:=ABS(U.val(row,i));{максимальный элемент среди просматриваемых строк}
    end
   end
  end
 end;

 row:=Wi[1].arr^[ir];
 Wi[1].arr^[ir]:=Wi[1].arr^[ind];
 Wi[1].arr^[ind]:=row;
 Wi[2].arr^[ind]:=ic+1;
 Wi[3].arr^[ind]:=ii
end;
{-------------------------------------------------------------------------}
function SparceLUDecomp;
var
  i,j             : Integer;
  Messg           : TMSG;
begin

 if first then
  for i:=0 to NOfE-1 do
   for j:=1 to 2 do Wi[j].arr^[i]:=i+1;

 for i:=0 to NOfE-2 do begin
 
  if first then begin
   Result:=SparceSort(i,NOfE,Rows,Wi,U,CU,RU,TU);
   if Result > 0 then exit;
  end;

  {выполнить i-ый шаг гауссова исключения}
  Result:=ElimStep(U,L,CU,RU,CL,Wi,i,TB);
  if Result > 0 then exit;

  if PeekMessage(Messg,0,0,0,PM_REMOVE) then
  SendMessage(Messg.hwnd,Messg.message,Messg.wParam,Messg.lParam);
 end;

 if first then begin
  Result:=SparceSort(NofE-1,NOfE,Rows,Wi,U,CU,RU,TU);
  if Result > 0 then exit;
 end

end;
{------------------------------------------------------------------------}
function SparceSolve;
var
 m,n,k,j,i : integer;
 sum   :Realtype;
begin
 Result:=0;
{определение новых правых частей при гауссовом исключении}
   for i:=1 to NOfE-1 do begin
    k:=Wi[1].arr^[i-1]-1;
    for j:=0 to L.Ptr(k).count-1 do begin
     m:=CL.val(k,j)-1;
     if m >= 0 then F[m]:=F[m]+F[k]*L.val(k,j);
    end;
   end;
{выполнить обратную прогонку,найти начальное приближение вектора
решений}
 for i:=NOfE-1 Downto 0 do begin
  k:=Wi[1].arr^[i]-1;
  m:=Wi[2].arr^[i]-1;
  n:=Wi[3].arr^[i];
  Sum:=F[k];

  for j:=0 to U.Ptr(k).count-1 do
   if j <> n then Sum:=Sum-U.val(k,j)*X[CU.val(k,j)-1];

  if (U.val(k,n) = 0.0)or(n < 0) then begin
   X[m]:=1;
   Result:=er_NullDet
  end else X[m]:=Sum/U.val(k,n)
 end
end;
{------------------------------------------------------------------------
------------------------------------------------------------------------}
function SparceIter;
var
 errdx         : Realtype;
 errdx1        : Realtype;
 i,j,k,iter    : Integer;
 flag          : Boolean;
const
 maxiter = 20;abserr = 1.e-15;
begin
{Итерационный процесс состоит из трех этапов}
{ Этап 2 : d(i) = U^(-1)*L^(-1)*r(i) }
{ Этап 1 : r(i) = f(i) - A*x(i)}
{ Этап 3 : x(i+1) = x(i) + d(i)}
 iter:=0;errdx1:=1.e30;flag:=false;


 REPEAT
  for i:=0 to NOfE-1 do begin
   F[i]:=Fr[i];
   for j:=0 to Us.Ptr(i).count-1 do F[i]:=F[i]-Us.val(i,j)*X[CUs.val(i,j)-1];
  end;

  Result:=SparceSolve(NOfE,Xr,F,Wi,U,L,CU,RU,CL);
  if Result > 0 then exit;

  errdx:=0;
  for k:=0 to NOfE-1 do begin
{    errx:=errx+ABS(X[k]);
    errdx:=errdx+ABS(Xr[k]);}
    errdx:=max(errdx,ABS(Xr[k])/(ABS(X[k])+abserr));
    X[k]:=X[k]+Xr[k];
  end;
  if errdx <= eps then flag:=true;
 { if errdx <= eps*amax1(errx,abserr) then flag:=True;}
  if iter > 0 then
   if errdx >= 2*errdx1 then begin
    {tells(0,'Итерационное уточнение системы ЛАУ расходится');}
    flag:=True
   end;

{   tells(0,'RelErr='+FloatToStrF(errdx/amax1(errx,abserr),ffGeneral,6,2)+' d(i)='+FloatToStrF(errdx,ffGeneral,6,2)+
 ' d(i-1)='+FloatToStrF(errdx1,ffGeneral,6,2)+' iter='+IntToStr(iter));}

  errdx1:=errdx;
  inc(iter);

 UNTIL (flag) or (iter >= maxiter);

end;
{---------------------------------------------------------------------}
{                      TSPARCE METHODS                               }
{---------------------------------------------------------------------}
function TSparce.GetDet:RealType;
var i,j:Integer;
    h  : RealType;
begin
 h:=1;
 for i:=0 to NOfE-1 do
  for j:=0 to WRA[1].Ptr(i).Count-1 do
   if WIA[1].val(i,j) = i+1 then begin
     h:=h*abs(WRA[1].val(i,j));
     break
    end;

 Result:=h;
end;
{-------------------------------------------------------------------------}
procedure TSparce.FillRightItem;
begin
 WR[1].arr^[ACol]:=X
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetRightItem;
begin
 X:=WR[1].arr^[ACol]
end;
{-------------------------------------------------------------------------}
procedure TSparce.FillRightCol;
begin
 Move(X^,WR[1].arr^,n*SOfR);
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetSolutionItem;
begin
 X:=WR[3].arr^[ACol]
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetSolutionError;
begin
 Move(WR[4].arr^,dX^,n*SOfR);
 Move(WR[2].arr^,dF^,n*SOfR)
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetSolutionErrorItem;
begin
 X:=WR[4].arr^[ACol];
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetSolutionCol;
begin
 Move(WR[3].arr^,X^,n*SOfR)
end;
{-------------------------------------------------------------------------}
function TSparce.GetSolution;
begin
 Result:=WR[3].arr
end;
{-------------------------------------------------------------------------}
procedure TSparce.FillRow(n,ARow:Integer;var X:PExtArr;var ACol:PIntArr);
var i : Integer;
begin
  for i:=0 to n-1 do begin
   if X^[i] = 0.0 then continue;
   WIA[2].Ptr(ARow).Add(ACol^[i]+1);
   WRA[2].Ptr(ARow).Add(X^[i])
  end
end;
{-------------------------------------------------------------------------}
procedure TSparce.FillItem(ARow,ACol:Integer;X: RealType);
begin
  if X = 0.0 then exit;
  WIA[2].Ptr(ARow).Add(ACol+1);
  WRA[2].Ptr(ARow).Add(X)
end;
{-------------------------------------------------------------------------}
constructor TSparce.Create;
var i : Integer;
begin
  inherited Create;
  NOfE:=N;
  Matrix_flag:=fm;Iter_flag:=fi;
  TU:=1;TB:=0;Rows:=3;

  for i:=1 to cwr do WR[i]:=TExtArray.Create(NOfE);
  for i:=1 to cwi do WI[i]:=TIntArray.Create(NOfE);
  for i:=1 to cwra do WRA[i]:=TPtrExt.Create(NOfE,1);
  for i:=1 to cwia do WIA[i]:=TPtrInt.Create(NOfE,1);
end;
{-------------------------------------------------------------------------}
procedure TSparce.GetRightPart;
begin
 Move(WR[1].arr^,WR[2].arr^,n*SOfR)
end;
{--------------------------------------------------------------------------}
destructor TSparce.Destroy;
var
   i,j : integer;
begin
 for i:=1 to cwr do WR[i].Free;
 for i:=1 to cwi do WI[i].Free;

 for j:=1 to cwra do begin
  for i:=0 to NOfE-1 do WRA[j].Ptr(i).Free;
  WRA[j].Free
 end;

 for j:=1 to cwia do begin
  for i:=0 to NOfE-1 do WIA[j].Ptr(i).Free;
  WIA[j].Free
 end;
 inherited Destroy
end;
{-------------------------------------------------------------------------}
procedure TSparce.ChangeNOfE;
var i,j : Integer;
begin
  for i:=1 to cwr do WR[i].ChangeCount(N);
  for i:=1 to cwi do WI[i].ChangeCount(N);
  for i:=1 to cwra do begin
   for j:=N to NofE-1 do WRA[i].Ptr(j).Free;
   for j:=NofE to N-1 do WRA[i].Add(1)
  end;
  for i:=1 to cwia do begin
   for j:=N to NofE-1 do WIA[i].Ptr(j).Free;
   for j:=NofE to N-1 do WIA[i].Add(1)
  end;
  NofE:=N
end;
{--------------------------------------------------------------------------}
procedure TSparce.SetMatrixMemo;
begin
  WRA[1].Ptr(Row).ChangeCount(cRow);
  WRA[2].Ptr(Row).ChangeCount(cRow);
  WIA[1].Ptr(Row).ChangeCount(cRow);
  WIA[2].Ptr(Row).ChangeCount(cRow);
  WRA[1].Ptr(Row).Count:=0;
  WRA[2].Ptr(Row).Count:=0;
  WIA[1].Ptr(Row).Count:=0;
  WIA[2].Ptr(Row).Count:=0;
end;
{--------------------------------------------------------------------------}
procedure TSparce.Init;
var i,j,k : Integer;
begin
  ChangeMemoryAlloc(n);

  for j:=0 to N-1 do WI[4].arr^[j]:=0;
  for j:=0 to N-1 do
   for i:=0 to WIA[1].Ptr(j).count-1 do Inc(WI[4].arr^[WIA[1].val(j,i)-1]);

  for j:=0 to N-1 do begin
   WIA[3].Ptr(j).changecount(WI[4].arr^[j]);
   WI[4].arr^[j]:=0
  end;

  for j:=0 to N-1 do
   for i:=0 to WIA[1].Ptr(j).count-1 do begin
    k:=WIA[1].val(j,i)-1;
    WIA[3].z(k,WI[4].arr^[k],j+1);
    Inc(WI[4].arr^[k])
   end
end;
{--------------------------------------------------------------------------}
function TSparce.LUDecomp;
begin
  Result:=SparceLUDecomp(N,Rows,WR,WI,WRA[1],WRA[3],WIA[1],WIA[3],WIA[4],TB,TU,flag);
end;
{--------------------------------------------------------------------------}
procedure TSparce.ChangeMemoryAlloc;
var j:integer;
begin
 for j:=0 to N-1 do begin
   WRA[1].Ptr(j).changecount(WRA[2].Ptr(j).count);
   WIA[1].Ptr(j).changecount(WRA[2].Ptr(j).count);
 end;
 for j:=0 to N-1 do begin
   Move(WRA[2].Ptr(j).arr^,WRA[1].Ptr(j).arr^,WRA[1].Ptr(j).count*SOfR);
   Move(WIA[2].Ptr(j).arr^,WIA[1].Ptr(j).arr^,WRA[1].Ptr(j).count*SOfI);
 end
end;
{-------------------------------------------------------------------------}
function TSparce.Run;
begin {Run}
  if not matrix_flag then GetRightPart(n)
  else begin
   Init(n);
   GetRightPart(n);
   LUDecomp(n,false)
  end;

  Result:=SparceSolve(N,WR[3].arr^,WR[2].arr^,WI,WRA[1],WRA[3],WIA[1],WIA[3],WIA[4]);
  if Result > 0 then exit;
  if iter_flag then
  Result:=SparceIter(N,WR[3].arr^,WR[4].arr^,WR[1].arr^,WR[2].arr^,WI,WRA[1],WRA[2],
                     WRA[3],WIA[1],WIA[3],WIA[2],WIA[4],0.0001)
end;
{-------------------------------------------------------------------------}
function TSparce.RunIter(n:Integer;hb,relerr:RealType):Integer;
var i,j,it     : integer;
    sum,err,f    : RealType;
const Itmax = 50;abserr = 1.e-10;
begin
  Result:=0;
  for i:=0 to N-1 do WR[4].arr^[i]:=0;
  IT:=0;
  repeat
   for i:=0 to N-1 do begin
    sum:=0;
    for j:=0 to WRA[2].Ptr(i).Count-1 do
     sum:=sum+WRA[2].val(i,j)*WR[4].arr^[WIA[2].val(i,j)-1];
    WR[3].arr^[i]:=WR[1].arr^[i]+hb*sum;
   end;

   err:=0;
   for i:=0 to N-1 do begin
    f:=max(abs(WR[4].arr^[i])+abs(WR[3].arr^[i]),abserr);
    err:=max(err,abs(WR[4].arr^[i]-WR[3].arr^[i])/f);
    WR[4].arr^[i]:=WR[3].arr^[i]
   end;
   Inc(IT)
  until (err <= relerr) or (IT > ITMAX);
  if IT > ITMAX then Result:=er_MaxIter
end;
{------------------------------------------------------------------------
                      Объект - прогонка
{------------------------------------------------------------------------}
CONSTRUCTOR TTridag.Create;
var i : Integer;
begin
 inherited Create;
 n:=an;
 for i:=1 to 6 do GetMem(wa[i],n*SOfR);
 wa[1]^[0]:=0; wa[3]^[n-1]:=0;
end;
{------------------------------------------------------------------------}
destructor TTridag.Destroy;
var i : Integer;
begin
 for i:=1 to 6 do FreeMem(wa[i],n*SOfR);
 inherited Destroy
end;
{------------------------------------------------------------------------}
procedure TTridag.ChangeNOfE;
var i : Integer;
begin
 if n <> an then begin
  for i:=1 to 6 do FreeMem(wa[i],n*SOfR);
  n:=an;
  for i:=1 to 6 do GetMem(wa[i],n*SOfR);
 end;
 wa[1]^[0]:=0; wa[3]^[n-1]:=0;
end;
{------------------------------------------------------------------------}
PROCEDURE TTridag.Progonka;
VAR
   j  : integer;
   bet: RealType;
BEGIN
   j:=0;
   bet := wa[2]^[j];
   wa[5]^[j] := wa[4]^[j]/bet;
   FOR j := 1 to n-1 DO BEGIN
      wa[6]^[j] := wa[3]^[j-1]/bet;
      bet := wa[2]^[j]-wa[1]^[j]*wa[6]^[j];
      wa[5]^[j] := (wa[4]^[j]-wa[1]^[j]*wa[5]^[j-1])/bet;
   END;
   FOR j := n-2 DOWNTO 0 DO wa[5]^[j] := wa[5]^[j]-wa[6]^[j+1]*wa[5]^[j+1]
END;

function GaussLUfact(n: integer;VAR a: PExtPtrArr; VAR indx: PIntArr;
                 VAR vv: PExtArr;VAR ner:Integer):RealType;
CONST
   tiny=1.0e-20;
VAR
   k,j,imax,i: integer;
   sum,dum,big,d: realtype;
BEGIN
   ner:=0;
   d := 1.0;
   FOR i := 0 to n-1 DO BEGIN
      big := 0.0;
      FOR j := 0 to n-1 DO IF (abs(a^[i]^[j]) > big) THEN big := abs(a^[i]^[j]);
      IF (big = 0.0) THEN BEGIN
       ner:=3;
       big:=tiny
      END;
      vv^[i] := 1.0/big
   END;
   FOR j := 0 to n-1 DO BEGIN
      IF (j > 0) THEN BEGIN
         FOR i := 0 to j-1 DO BEGIN
            sum := a^[i]^[j];
            IF (i > 0) THEN BEGIN
               FOR k := 0 to i-1 DO sum := sum-a^[i]^[k]*a^[k]^[j];
               a^[i]^[j] := sum
            END
         END
      END;
      big := 0.0;
      FOR i := j to n-1 DO BEGIN
         sum := a^[i]^[j];
         IF (j > 0) THEN BEGIN
            FOR k := 0 to j-1 DO  sum := sum-a^[i]^[k]*a^[k]^[j];
            a^[i]^[j] := sum
         END;
         dum := vv^[i]*abs(sum);
         IF (dum > big) THEN BEGIN
            big := dum;
            imax := i
         END
      END;
      IF (j <> imax) THEN BEGIN
         FOR k := 0 to n-1 DO BEGIN
            dum := a^[imax]^[k];
            a^[imax]^[k] := a^[j]^[k];
            a^[j]^[k] := dum
         END;
         d := -d;
         vv^[imax] := vv^[j]
      END;
      indx^[j] := imax;
      IF (j <> n-1) THEN BEGIN
         IF (a^[j]^[j] = 0.0) THEN a^[j]^[j] := tiny;
         dum := 1.0/a^[j]^[j];
         FOR i := j+1 to n-1 DO a^[i]^[j] := a^[i]^[j]*dum
      END
   END;
   IF (a^[n-1]^[n-1] = 0.0) THEN a^[n-1]^[n-1] := tiny;
   Result:=d
END;

PROCEDURE GaussSolve(n: integer;a: PExtPtrArr;indx: PIntArr; VAR b,x: PExtArr);
VAR
   j,ip,ii,i: integer;
   sum: realtype;
BEGIN
   Move(b^,x^,n*SOfR);
   ii := -1;
   FOR i := 0 to n-1 DO BEGIN
      ip := indx^[i];
      sum := x^[ip];
      x^[ip] := x^[i];
      IF  (ii <> -1) THEN
      FOR j := ii to i-1 DO sum := sum-a^[i]^[j]*x^[j]
      ELSE IF (sum <> 0.0) THEN ii := i;
      x^[i] := sum
   END;
   FOR i := n-1 DOWNTO 0 DO BEGIN
      sum := x^[i];
      IF (i < n-1) THEN
      FOR j := i+1 to n-1 DO sum := sum-a^[i]^[j]*x^[j];
      x^[i] := sum/a^[i]^[i]
   END
END;

function popolam;
var y,y0,y1,dx,x1,x,err : RealType;
const epsx = 0.005;
      epsdx = 0.0001;
begin
     y0:=fun(x0,pt);
     dx:=max(abs(epsx*x0),epsx);
     y1:=fun(x0+dx,pt);

     if (y1*y0 > 0) and (abs(y1) > abs(y0)) then dx:=-dx;
     repeat
      dx:=2*dx;
      x1:=x0+dx;
      y1:=fun(x1,pt);
     until y1*y0 < 0;

     if y0 < 0 then begin
      dx:=x1;
      x1:=x0;
      x0:=dx
     end;

     repeat
      x:=0.5*(x0+x1);
      y:=fun(x,pt);
      if y > 0 then x0:=x else x1:=x;
      err:=abs(x0-x1)/max(abs(x0),abs(x1));
     until err < epsdx;
     result:=x
end;

END.