program endgueltig;

{ Attentiom: This program only runs,}
{ if [M] is an interval H-matrix }
{ with positive diagonal entries }

use i_ari, mvi_ari, mv_ari;


function iga(A:imatrix;b:ivector):ivector[lb(b)..ub(b)];
var k,i,j : integer;
    x,y   : ivector[lb(b)..ub(b)];
    hilf  : interval;
begin
for k:=lb(b) to ub(b)-1 do
begin

  for i:= k+1 to ub(b) do
  begin
      for j:= k+1  to ub(b) do
      begin
           a[i,j]:=a[i,j] - a[k,j] * a[i,k] /a[k,k];
      end;

      {jetzt wird die untere Dreiecksmatrix besetzt}
      a[i,k]:=a[i,k]/a[k,k];

  end;
end;

{Loese nun Ly=b}

for i:= lb(b) to ub(b) do
begin
    hilf:=0;
    for j:= lb(b) to i-1 do
    begin
        hilf:= hilf + a[i,j]*y[j];
    end;
    y[i]:=b[i]- hilf;
end;
{Loese nun Rx=y}

for i:= ub(b) downto lb(b) do
begin
     hilf:=0;
     for j:= i+1 to ub(b) do
     begin
         hilf:= hilf + a[i,j]*x[j];
     end;
     x[i]:= (y[i] - hilf)/a[i,i];
end;

for i:= lb(b) to ub(b) do iga[i]:=x[i];

end;






function imax(a	: interval):interval;
begin
   if (a.sup < 0) or (a.sup=0) then imax:=0
   else if (a.inf > 0) or (a.inf=0) then imax:=a
   else imax:=intval(0,a.sup);
end;

function betrag(a:interval):real;
begin
if a.inf > 0 then betrag:=a.sup
else if a.sup < 0 then betrag:=-a.inf
     else if -a.inf < a.sup then betrag:=a.sup
          else betrag:=-a.inf;
end;


procedure haupt(n:integer;stop:integer);
var
      i,j,k,h,l,zaehl         : integer;
      Hilfmatrix              : imatrix[1..n,1..n];
      ep                      : rvector[1..n];
      x_1,v,u,alpha             : ivector[1..n];
      M,R                     : imatrix[1..n,1..n];
      q                       : ivector[1..n];
      hvector                 : ivector[1..n];
      y_alt,y_neu,ystart      : ivector[1..n];
      y                       : ivector[1..n];
      Loese                   : ivector[1..2*n];
      unten,oben              : real;
      hilf                    : interval;

begin


writeln('Give the interval matrix [M] row-wise');
writeln('[M] has to be an H-matrix with positive diagonal entries');
for i:= 1 to n do
for j:= 1 to n do
read(M[i,j]);
writeln('Give the interval vector [q]');
for i:= 1 to n do
read(q[i]);
writeln;


R:= -M;
for i:= 1 to n do R[i,i]:=0;


{first we calculate a crude enclosure of [x*]}

{calculate P(I-P)-1*alpha}

{with alpha:=q([x_0],[x_1])}

{it is P(I-P)-1*alpha = <[D]>-1*|[R]|*<M>-1*<[D]>*alpha}

{calculate [x_1], if [x_0]=0}

   for i:= 1 to n do x_1[i]:=imax(-q[i]/M[i,i]);


{set alpha := q([x_0],[x_1])=q(0,[x_1])}
{and multiply with <[D]>}

   for j:= 1 to n do alpha[j]:=x_1[j].sup*M[j,j].inf;

{ calculate <[M]>-1* alpha}

for i := 1 to n do
for j:= 1 to n do
if i=j then Hilfmatrix[i,j]:=M[i,j].inf
       else Hilfmatrix[i,j]:=-betrag(M[i,j]);

u:=iga(Hilfmatrix,alpha);

{ calculate |[R]|* <[M]>-1* alpha}

for i := 1 to n do
for j:= 1 to n do Hilfmatrix[i,j]:= betrag(R[i,j]);

hvector:=Hilfmatrix*u;

   for i:= 1 to n do v[i]:=hvector[i]/M[i,i].inf;

for i:= 1 to n do ep[i]:= v[i].sup;


for i:= 1 to n do
ystart[i] := intval(x_1[i].inf -< ep[i], x_1[i].sup +> ep[i] );

   for i:= 1 to n do if ystart[i].inf<0 then ystart[i].inf:=0;

y_neu:=ystart;
for j:= 1 to n do
begin
unten :=y_neu[j].inf;
oben  :=y_neu[j].sup;
writeln('[',unten:23:0:-1,',',oben:23:0:1,']');
end;

writeln;

zaehl:=0;
repeat

     zaehl:= zaehl+1;
     y_alt := y_neu;
     hvector:= R*y_alt-q;
     for i:=1 to n do y_neu[i]:= imax(hvector[i]/M[i,i])**y_alt[i];

until  ((zaehl=stop) or (y_neu=y_alt));

writeln;
write('We had ',zaehl,' iteration steps');
writeln;
writeln('Enclosure of LB : ');
writeln;

for j:= 1 to n do
begin
unten :=y_neu[j].inf;
oben  :=y_neu[j].sup;
writeln('[',unten:23:0:-1,',',oben:23:0:1,']');
end;


end;

var n,stop        : integer;


begin

writeln('Give the dimension of the LCP');
read(n);
writeln('How many iteration steps should be done at most ?');
read(stop);
haupt(n,stop);


end.