program endgueltig;

{ Attention: This program only runs for }
{ interval H-matrices }
{ with 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, sum1,sum2,klammer : 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;


{calculate a crude enclosure of [x*]}

{first 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]);


{then, 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);

{ jetzt berechnen wir |[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;

writeln;

zaehl:=0;
repeat

     zaehl:= zaehl+1;
     y_alt:=y_neu;
   for i:= 1 to n do
   begin

      if i=1 then begin
	 sum1:=0;
	 for k:=2 to n do sum1:=sum1 + M[1,k]*y_alt[k];
	 klammer:= -sum1 - q[1];
      end
      else if i=n then begin
	 sum2:=0;
	 for k:=1 to n-1 do sum2:=sum2+M[n,k]*y_neu[k];
	 klammer:= -sum2-q[n];
      end
      else begin
	 sum1:=0;
	 sum2:=0;
	 for k:=i+1 to n do sum1:=sum1+M[i,k]*y_alt[k];
	 for k:= 1 to i-1 do sum2:= sum2+M[i,k]*y_neu[k];
	 klammer:=-sum1-sum2-q[i];
      end;
      y_neu[i]:=imax(klammer/M[i,i])**y_alt[i];
   end;




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.