program vers; use i_ari, mvi_ari, mv_ari; type variablen = record art: (w,z); stelle: integer; end; basis = dynamic array[*] of variablen; procedure GJ(var A:imatrix; var q: ivector; pz,ps:integer); var i,j : integer; begin q[pz]:=q[pz]/A[pz,ps]; for j := lb(A,2) to ub(A,2) do if j<> ps then A[pz,j]:=A[pz,j]/A[pz,ps]; A[pz,ps]:=1; for i:= lb(A,1) to ub(A,1) do if not (i = pz) then begin for j:= lb(A,2) to ub(A,2) do if j<>ps then A[i,j]:=A[i,j]-A[pz,j]*A[i,ps]; q[i]:=q[i] - q[pz]*A[i,ps]; end; for i:= lb(A,1) to ub(A,1) do if i<>pz then A[i,ps]:=0; end; function anf_pivotzeile(q:ivector):integer; { es wird angenommen, dass nicht q>=0 gilt } var stelle, i : integer; vorerst : interval; begin stelle:=lb(q); vorerst:=q[lb(q)]; for i:= lb(q)+1 to ub(q) do if q[i].sup < vorerst.inf then begin vorerst:= q[i]; stelle:=i; end; anf_pivotzeile:=stelle; end; function pivotzeile(A:imatrix;q:ivector;ps:integer):integer; var loesbar : boolean; i,stelle,t : integer; test, vorerst : interval; begin loesbar:=false; for i:= ub(q) downto lb(q) do if A[i,ps].inf > 0 then begin loesbar:=true; stelle:=i; end; if not loesbar then begin writeln('Interval ray termination'); pivotzeile:=0; end else begin vorerst:=q[stelle]/A[stelle,ps]; t:=stelle; for i:= lb(q) to ub(q) do if A[i,ps].sup > 0 then begin test:= q[i].inf/A[i,ps].sup; if test.inf < vorerst.inf then begin vorerst:=test; t:=i; if A[t,ps].inf > 0 then stelle:=t; end; end; if A[t,ps].inf <= 0 then pivotzeile:=stelle else pivotzeile:=t; end; end; procedure einlesen(var A:imatrix; var q : ivector); var i,j : integer; M : imatrix[lb(A,1)..ub(A,1),lb(A,1)..ub(A,1)]; begin A:=0; writeln('Give your interval matrix [M] row-wise'); for i:= lb(M) to ub(M) do for j:= lb(M) to ub(M) do begin read(M[i,j]); A[i,ub(A,1)+j]:=-M[i,j]; end; writeln('Give the interval vector [q]'); for i:= lb(q) to ub(q) do begin read(q[i]); A[i,i]:=1; A[i,ub(A,2)]:=-1 end; writeln; end; procedure haupt(n:integer); var A : imatrix[1..n,1..2*n+1]; q : ivector[1..n]; B : basis[1..n]; i,ps,pz : integer; glob : integer; loesung : ivector[1..2*n]; unten,oben : real; gegangen, jetzt : variablen; vorsicht : boolean; fertig : boolean; label 1; { Zentrale Idee des Algorithmus: } { Hat die Variable z_j die Basis verlassen, } { so ist im naechsten Schritt die Variable w_j } { die Pivotspalte und w_j kommt zur Basis } { Und zwar an die Stelle, die sich als Pivotzeile herausstellt } { In gegangen merken wir uns, die Variable, die gegangen ist } { In jetzt, die die Pivotspalte wird und in die Basis kommt } begin { Initialisierung der Basis } for i:= 1 to n do begin B[i].art:=w; B[i].stelle:=i; end; einlesen(A,q); {schreiben(A,q);} fertig := true; for i:= 1 to n do if q[i].inf <0 then fertig:=false; if fertig then goto 1; { erster pivotschritt } ps:=2*n+1; pz:= anf_pivotzeile(q); glob:=pz; GJ(A,q,pz,ps); {schreiben(A,q);} B[pz].art:= z; B[pz].stelle:=0; gegangen.art:=w; gegangen.stelle:=pz; { z[pz] ist neue Pivotspalte } { w[pz] hat Basis verlassen } while B[glob].stelle=0 do begin if gegangen.art=w then { Pivotspalte } begin ps:=n+gegangen.stelle; jetzt.art:=z end; if gegangen.art=z then { Pivotspalte } begin ps:=gegangen.stelle; jetzt.art:=w; end; jetzt.stelle:=gegangen.stelle; { Pivotspalte merken } pz:=pivotzeile(A,q,ps); GJ(A,q,pz,ps); { Basiswechsel } {schreiben(A,q);} if B[pz].art=z then gegangen.art:=z; { z[pz] hat Basis verlassen } if B[pz].art=w then gegangen.art:=w; { w[pz] hat Basis verlassen } gegangen.stelle:=B[pz].stelle; { Zur Basis kommt die Pivotspalte } B[pz].art:=jetzt.art; B[pz].stelle:=jetzt.stelle; end;{while} { jetzt aufsammeln } 1: loesung:=0; for i := 1 to n do begin if B[i].art=z then loesung[n+B[i].stelle]:=q[i]; if B[i].art=w then loesung[B[i].stelle]:=q[i]; end; for i:= 1 to 2*n do begin unten:= loesung[i].inf; oben:= loesung[i].sup; writeln('[',unten:23:0:-1,',',oben:23:0:1,']'); end; { w - Mz =q koennte zwei Loesungen haben } { wovon eine nichtnegative Eintraege haben } { koennte. Die koennte in dem Aufgeblaehten q } { dann enthalten sein } vorsicht:=false; for i:= 1 to 2*n do if loesung[i].inf<0 then vorsicht:=true; if vorsicht then writeln('ILA fails!'); end; var n : integer; begin writeln('Give the dimension of the LCP'); read(n); haupt(n); end.