File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / generic_integer_linear_solvers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:23 2000 UTC (23 years, 10 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
package body Generic_Integer_Linear_Solvers is
use Common_Divisors;
-- SCALERS :
function Divisors ( a : Matrix ) return Vectors.Vector is
-- DESCRIPTION :
-- Returns a vector containing the gcd of the elements of each row.
v : Vectors.Vector(a'range(1));
begin
for i in a'range(1) loop
v(i) := a(i,a'first(2));
if not Equal(v(i),one)
then for j in (a'first(2)+1)..a'last(2) loop
v(i) := gcd(v(i),a(i,j));
exit when Equal(v(i),one);
end loop;
end if;
end loop;
return v;
end Divisors;
function Scale ( a : Matrix ) return Matrix is
v : Vectors.Vector(a'range(1)) := Divisors(a);
b : Matrix(a'range(1),a'range(2));
begin
for i in b'range(1) loop
if (Equal(v(i),one) or Equal(v(i),zero))
then for j in b'range(2) loop
b(i,j) := a(i,j);
end loop;
else for j in b'range(2) loop
b(i,j) := a(i,j) / v(i);
end loop;
end if;
end loop;
return b;
end Scale;
procedure Scale ( a : in out Matrix; v : out Vectors.Vector ) is
dv : Vectors.Vector(a'range(1)) := Divisors(a);
begin
for i in a'range(1) loop
if (Equal(dv(i),one) or Equal(dv(i),zero))
then null;
else for j in a'range(2) loop
a(i,j) := a(i,j) / dv(i);
end loop;
end if;
end loop;
v := dv;
end Scale;
procedure Scale ( a : in out Matrix ) is
v : Vectors.Vector(a'range(1)) := Divisors(a);
begin
for i in a'range(1) loop
if (Equal(v(i),one) or Equal(v(i),zero))
then null;
else for j in a'range(2) loop
a(i,j) := a(i,j) / v(i);
end loop;
end if;
end loop;
end Scale;
procedure Scale ( a : in out Matrix; row,col : in integer ) is
g : number := a(row,col);
begin
if not Equal(g,one)
then for l in (col+1)..a'last(2) loop
g := gcd(g,a(row,l));
exit when Equal(g,one);
end loop;
end if;
if (not Equal(g,zero)) and (not Equal(g,one))
then for l in row..a'last(2) loop
a(row,l) := a(row,l)/g;
end loop;
end if;
end Scale;
-- STATIC TRIANGULATORS :
procedure Upper_Triangulate ( l : out Matrix; a : in out Matrix ) is
row,pivot : integer;
temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
ll : Matrix(a'range(1),a'range(1));
begin
for i in ll'range(1) loop
for j in ll'range(2) loop
Copy(zero,ll(i,j));
end loop;
Copy(one,ll(i,i));
end loop;
row := a'first(1);
for j in a'first(2)..a'last(2) loop
pivot := row-1; -- find pivot
for i in row..a'last(1) loop
if not Equal(a(i,j),zero)
then pivot := i;
exit;
end if;
end loop;
if pivot >= row
then if pivot /= row -- interchange if necessary
then for k in a'range(2) loop
Copy(a(row,k),temp);
Copy(a(pivot,k),a(row,k));
Copy(temp,a(pivot,k));
end loop;
for k in ll'range(2) loop
Copy(ll(row,k),temp);
Copy(ll(pivot,k),ll(row,k));
Copy(temp,ll(pivot,k));
end loop;
end if;
for i in (row+1)..a'last(1) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(row,j),a(i,j),ka,lb,d);
aa := a(row,j)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in j..a'last(2) loop
Copy(a(row,k),a_rowk); Clear(a(row,k));
-- a_rowk := a(row,k);
Copy(a(i,k),a_ik); Clear(a(i,k)); -- a_ik := a(i,k);
a(row,k) := ka*a_rowk + lb*a_ik;
a(i,k) := -bb*a_rowk + aa*a_ik;
end loop;
for k in ll'range(2) loop
a_rowk := ll(row,k);
a_ik := ll(i,k);
ll(row,k) := ka*a_rowk + lb*a_ik;
ll(i,k) := -bb*a_rowk + aa*a_ik;
end loop;
end if;
end loop;
row := row + 1;
end if;
exit when row > a'last(1);
end loop;
l := ll;
end Upper_Triangulate;
procedure Upper_Triangulate ( a : in out Matrix ) is
row,pivot : integer;
temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
begin
row := a'first(1);
for j in a'first(2)..a'last(2) loop
pivot := row-1; -- find pivot
for i in row..a'last(1) loop
if not Equal(a(i,j),zero)
then pivot := i;
exit;
end if;
end loop;
if pivot >= row
then if pivot /= row -- interchange if necessary
then for k in a'range(2) loop
Copy(a(row,k),temp);
Copy(a(pivot,k),a(row,k));
Copy(temp,a(pivot,k));
end loop;
end if;
for i in (row+1)..a'last(1) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(row,j),a(i,j),ka,lb,d);
aa := a(row,j)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in j..a'last(2) loop
Copy(a(row,k),a_rowk); Clear(a(row,k));
--a_rowk := a(row,k);
Copy(a(i,k),a_ik); Clear(a(i,k)); --a_ik := a(i,k);
a(row,k) := ka*a_rowk + lb*a_ik;
a(i,k) := -bb*a_rowk + aa*a_ik;
end loop;
end if;
end loop;
row := row + 1;
end if;
exit when row > a'last(1);
end loop;
end Upper_Triangulate;
procedure Upper_Triangulate
( a : in out Matrix;
ipvt : in out Standard_Integer_Vectors.Vector ) is
row,pivot,tmppiv : integer;
temp,aa,bb,ka,lb,d,a_rowk,a_ik : number;
begin
row := a'first(1);
for j in a'first(2)..a'last(2) loop
pivot := row-1; -- find pivot
for i in row..a'last(1) loop
if not Equal(a(i,j),zero)
then pivot := i;
exit;
end if;
end loop;
if pivot >= row
then if pivot /= row -- interchange if necessary
then for k in a'range(2) loop
Copy(a(row,k),temp);
Copy(a(pivot,k),a(row,k));
Copy(temp,a(pivot,k));
end loop;
tmppiv := ipvt(row);
ipvt(row) := ipvt(pivot);
ipvt(pivot) := tmppiv;
end if;
for i in (row+1)..a'last(1) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(row,j),a(i,j),ka,lb,d);
aa := a(row,j)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in j..a'last(2) loop
Copy(a(row,k),a_rowk); Clear(a(row,k));
-- a_rowk := a(row,k);
Copy(a(i,k),a_ik); Clear(a(i,k)); -- a_ik := a(i,k);
a(row,k) := ka*a_rowk + lb*a_ik;
a(i,k) := -bb*a_rowk + aa*a_ik;
end loop;
end if;
end loop;
row := row + 1;
end if;
exit when row > a'last(1);
end loop;
end Upper_Triangulate;
procedure Lower_Triangulate ( a : in out Matrix; u : out Matrix ) is
column,pivot : integer;
temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
uu : Matrix(a'range(2),a'range(2));
begin
for i in uu'range(1) loop
for j in uu'range(2) loop
Copy(zero,uu(i,j));
end loop;
Copy(one,uu(i,i));
end loop;
column := a'first(2);
for i in a'first(1)..a'last(1) loop
pivot := column-1; -- find pivot
for j in column..a'last(2) loop
if not Equal(a(i,j),zero)
then pivot := j;
exit;
end if;
end loop;
if pivot >= column
then if pivot /= column -- interchange if necessary
then for k in a'range(1) loop
Copy(a(k,column),temp);
Copy(a(k,pivot),a(k,column));
Copy(temp,a(k,pivot));
end loop;
for k in uu'range(1) loop
Copy(uu(k,column),temp);
Copy(uu(k,pivot),uu(k,column));
Copy(temp,uu(k,pivot));
end loop;
end if;
for j in (column+1)..a'last(2) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(i,column),a(i,j),ka,lb,d);
aa := a(i,column)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in i..a'last(1) loop
Copy(a(k,column),a_kcolumn); Clear(a(k,column));
-- a_kcolumn := a(k,column);
Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
a(k,column) := a_kcolumn*ka + a_kj*lb;
a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
end loop;
for k in uu'range(1) loop
Copy(uu(k,column),a_kcolumn); Clear(uu(k,column));
-- a_kcolumn := uu(k,column);
Copy(uu(k,j),a_kj); Clear(u(k,j)); -- a_kj := uu(k,j);
uu(k,column) := a_kcolumn*ka + a_kj*lb;
uu(k,j) := a_kcolumn*(-bb) + a_kj*aa;
end loop;
end if;
end loop;
column := column + 1;
end if;
exit when column > a'last(2);
end loop;
u := uu;
end Lower_Triangulate;
procedure Lower_Triangulate ( a : in out Matrix ) is
column,pivot : integer;
temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
begin
column := a'first(2);
for i in a'first(1)..a'last(1) loop
pivot := column-1; -- find pivot
for j in column..a'last(2) loop
if not Equal(a(i,j),zero)
then pivot := j;
exit;
end if;
end loop;
if pivot >= column
then if pivot /= column -- interchange if necessary
then for k in a'range(1) loop
Copy(a(k,column),temp);
Copy(a(k,pivot),a(k,column));
Copy(temp,a(k,pivot));
end loop;
end if;
for j in (column+1)..a'last(2) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(i,column),a(i,j),ka,lb,d);
aa := a(i,column)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in i..a'last(1) loop
Copy(a(k,column),a_kcolumn); Clear(a(k,column));
-- a_kcolumn := a(k,column);
Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
a(k,column) := a_kcolumn*ka + a_kj*lb;
a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
end loop;
end if;
end loop;
column := column + 1;
end if;
exit when column > a'last(2);
end loop;
end Lower_Triangulate;
procedure Lower_Triangulate
( a : in out Matrix;
ipvt : in out Standard_Integer_Vectors.Vector ) is
column,pivot,tmppiv : integer;
temp,aa,bb,ka,lb,d,a_kcolumn,a_kj : number;
begin
column := a'first(2);
for i in a'first(1)..a'last(1) loop
pivot := column-1; -- find pivot
for j in column..a'last(2) loop
if not Equal(a(i,j),zero)
then pivot := j;
exit;
end if;
end loop;
if pivot >= column
then if pivot /= column -- interchange if necessary
then for k in a'range(1) loop
Copy(a(k,column),temp);
Copy(a(k,pivot),a(k,column));
Copy(temp,a(k,pivot));
end loop;
tmppiv := ipvt(column);
ipvt(column) := ipvt(pivot);
ipvt(pivot) := tmppiv;
end if;
for j in (column+1)..a'last(2) loop -- make zeroes
if not Equal(a(i,j),zero)
then gcd(a(i,column),a(i,j),ka,lb,d);
aa := a(i,column)/d; bb := a(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in i..a'last(1) loop
Copy(a(k,column),a_kcolumn); Clear(a(k,column));
-- a_kcolumn := a(k,column);
Copy(a(k,j),a_kj); Clear(a(k,j)); -- a_kj := a(k,j);
a(k,column) := a_kcolumn*ka + a_kj*lb;
a(k,j) := a_kcolumn*(-bb) + a_kj*aa;
end loop;
end if;
end loop;
column := column + 1;
end if;
exit when column > a'last(2);
end loop;
end Lower_Triangulate;
-- SELECTORS :
function Det ( a : Matrix ) return number is
-- NOTE :
-- The triangulation is implemented independently to keep track
-- of row interchanges.
res : number := one;
m : matrix(a'range(1),a'range(2));
row,pivot : integer;
temp,aa,bb,ka,lb,d,m_rowk,m_ik : number;
begin
Copy(a,m); -- triangulate m
row := m'first(1);
for j in m'first(2)..m'last(2) loop
pivot := row-1; -- find pivot
for i in row..m'last(1) loop
if not Equal(m(i,j),zero)
then pivot := i;
exit;
end if;
end loop;
if pivot >= row
then if pivot /= row -- interchange if necessary
then for k in m'range(2) loop
Copy(m(row,k),temp);
Copy(m(pivot,k),m(row,k));
Copy(temp,m(pivot,k));
end loop;
Min(res);
end if;
for i in (row+1)..m'last(1) loop -- make zeroes
if not Equal(m(i,j),zero)
then gcd(m(row,j),m(i,j),ka,lb,d);
aa := m(row,j)/d; bb := m(i,j)/d;
if Equal(aa,bb) and then Equal(ka,zero)
then Copy(lb,ka); Copy(zero,lb);
end if;
if Equal(aa,-bb) and then Equal(ka,zero)
then ka := -lb; Copy(zero,lb);
end if;
for k in j..m'last(2) loop
Copy(m(row,k),m_rowk); Clear(m(row,k));
-- m_rowk := m(row,k);
Copy(m(i,k),m_ik); Clear(m(i,k)); -- m_ik := m(i,k);
m(row,k) := ka*m_rowk + lb*m_ik;
m(i,k) := -bb*m_rowk + aa*m_ik;
end loop;
end if;
end loop;
row := row + 1;
end if;
exit when row > m'last(1);
end loop;
for k in m'range(1) loop
Mul(res,m(k,k));
end loop;
Clear(m);
return res;
end Det;
function Per ( i,n : natural; a : matrix;
kk : Standard_Integer_Vectors.Vector ) return number is
begin
if i = n+1
then return one;
else declare
res,acc : number;
kkk : Standard_Integer_Vectors.Vector(kk'range) := kk;
begin
Copy(zero,res);
for j in kk'range loop
if not Equal(a(i,j),zero) and then (kk(j) /= 0)
then kkk(j) := kkk(j) - 1;
acc := Per(i+1,n,a,kkk);
Add(res,acc);
Clear(acc);
kkk(j) := kkk(j) + 1;
end if;
end loop;
return res;
end;
end if;
end Per;
function Per ( i,n : natural; a : matrix;
kk : Standard_Integer_Vectors.Vector; max : number )
return number is
begin
if i = n+1
then return one;
else declare
res,acc : number;
kkk : Standard_Integer_Vectors.Vector(kk'range) := kk;
begin
Copy(zero,res);
for j in kk'range loop
if not Equal(a(i,j),zero) and then (kk(j) /= 0)
then kkk(j) := kkk(j) - 1;
acc := a(i,j)*Per(i+1,n,a,kkk,max);
Add(res,acc);
Clear(acc);
kkk(j) := kkk(j) + 1;
end if;
exit when (res > max) or Equal(res,max);
end loop;
return res;
end;
end if;
end Per;
function Per ( a : matrix; k : Standard_Integer_Vectors.Vector )
return number is
-- ALGORITHM :
-- Row expansion without memory, as developed by C.W. Wampler,
-- see `Bezout Number Calculations for Multi-Homogeneous Polynomial
-- Systems', Appl. Math. Comput. 51:(2-3), 143-157, 1992.
begin
return Per(1,a'last(1),a,k);
end Per;
function Per ( a : matrix; k : Standard_Integer_Vectors.Vector;
max : number ) return number is
-- ALGORITHM :
-- Row expansion without memory, as developed by C.W. Wampler,
-- see `Bezout Number Calculations for Multi-Homogeneous Polynomial
-- Systems', Appl. Math. Comput. 51:(2-3), 143-157, 1992.
begin
return Per(1,a'last(1),a,k,max);
end Per;
function Rank ( a : Matrix ) return natural is
res : natural := 0;
m : Matrix(a'range(1),a'range(2));
column : integer;
begin
Copy(a,m);
Upper_Triangulate(m);
-- compute the length of chain of nonzero elements in m :
-- search first nonzero element in first row of m :
column := m'first(2)-1;
for k in m'range(2) loop
if not Equal(m(m'first(1),k),zero)
then column := k;
end if;
exit when (column = k);
end loop;
if column < m'first(2)
then return 0; -- all elements of m are zero
else for k in m'range(1) loop
exit when column > m'last(2);
if not Equal(m(k,column),zero)
then res := res + 1;
else -- search for next nonzero element on row k :
for l in column+1..m'last(2) loop
if not Equal(m(k,l),zero)
then column := l;
res := res + 1;
end if;
exit when (column = l);
end loop;
end if;
column := column + 1;
end loop;
end if;
Clear(m);
return res;
end Rank;
-- DYNAMIC TRIANGULATOR :
procedure Triangulate ( l : in integer; m : in out matrix;
ipvt : in out Standard_Integer_Vectors.vector;
piv : out integer ) is
-- DESCRIPTION :
-- Updates lth row of m such that m remains upper triangular.
tmp,a,b,lcmab,faca,facb : number;
pivot,index,tmppiv : integer;
begin
Switch(ipvt,l,m); -- pivoting for previous unknowns
index := 1; -- update : make l-1 zeroes in row l
for k in 1..(l-1) loop
if (not Equal(m(l,index),zero))
and then (not Equal(m(k,index),zero)) -- make m(l,index) zero
then a := m(k,index); b := m(l,index);
lcmab := lcm(a,b);
if lcmab < zero then lcmab := -lcmab; end if;
facb := lcmab/b; faca := lcmab/a;
if facb > zero
then for i in index..m'last(2) loop
m(l,i) := facb*m(l,i) - faca*m(k,i);
end loop;
else for i in index..m'last(2) loop
m(l,i) := -facb*m(l,i) + faca*m(k,i);
end loop;
end if;
end if;
if not Equal(m(k,index),zero)
then index := index + 1;
end if;
end loop;
pivot := 0; -- search pivot
for k in l..m'last(2)-1 loop
if not Equal(m(l,k),zero)
then pivot := k;
end if;
exit when pivot /= 0;
end loop;
if pivot > l
then for k in 1..l loop -- interchange columns l and pivot
Copy(m(k,l),tmp);
Copy(m(k,pivot),m(k,l));
Copy(tmp,m(k,pivot));
end loop;
tmppiv := ipvt(l);
ipvt(l) := ipvt(pivot);
ipvt(pivot) := tmppiv;
end if;
piv := pivot;
end Triangulate;
procedure Switch ( ipvt : in Standard_Integer_Vectors.vector;
index : in integer; m : in out matrix ) is
tmpv : Vectors.Vector(m'range(2));
begin
for k in tmpv'range loop
Copy(m(index,k),tmpv(k));
end loop;
for k in tmpv'range loop
Copy(tmpv(ipvt(k)),m(index,k));
end loop;
Vectors.Clear(tmpv);
end Switch;
procedure Switch ( ipvt : in Standard_Integer_Vectors.vector;
first,last : in integer; m : in out matrix) is
tmpv : Vectors.Vector(m'range(2));
begin
for index in first..last loop
for k in tmpv'range loop
Copy(m(index,k),tmpv(k));
end loop;
for k in tmpv'range loop
Copy(tmpv(ipvt(k)),m(index,k));
end loop;
Vectors.Clear(tmpv);
end loop;
end Switch;
procedure Switch ( l,pivot,index : in integer; m : in out matrix ) is
tmp : number;
begin
if l /= pivot
then Copy(m(index,l),tmp);
Copy(m(index,pivot),m(index,l));
Copy(tmp,m(index,pivot));
end if;
end Switch;
procedure Switch ( l,pivot : in integer;
first,last : in integer; m : in out matrix ) is
tmp : number;
begin
if l /= pivot
then for index in first..last loop
Copy(m(index,l),tmp);
Copy(m(index,pivot),m(index,l));
Copy(tmp,m(index,pivot));
end loop;
end if;
end Switch;
-- SOLVERS :
function Check0 ( a : Matrix; x : Vectors.Vector ) return boolean is
-- DESCRIPTION :
-- Returns true if x is a solution of the system a*x = 0.
tmp : Vectors.Vector(a'range(1));
begin
tmp := a*x;
for i in tmp'range loop
if not Equal(tmp(i),zero)
then Vectors.Clear(tmp); return false;
end if;
end loop;
Vectors.Clear(tmp);
return true;
end Check0;
procedure Solve0 ( a : in Matrix; x : in out Vectors.Vector ) is
-- ALGORITHM :
-- An intermediate, generating matrix tmp will be constructed,
-- such that
-- 1) the solution x to tmp*x = 0 is the same of a*x = 0;
-- 2) tmp(i,i) and tmp(i,ind) are the only nonzero entries.
-- Before this construction, it will be checked whether there
-- exists a zero column and the index ind must be determined.
-- After the definition of tmp, the back substitution process
-- yields a solution.
piv,ind : integer;
tmp : Matrix(a'first(1)..(a'last(1)+1),a'range(2));
res : Vectors.Vector(x'range);
pivots : Standard_Integer_Vectors.Vector(x'range);
zero_column : boolean;
begin
-- initialization of tmp,ind and pivots :
for i in tmp'range(1) loop
for j in tmp'range(2) loop
Copy(zero,tmp(i,j));
end loop;
end loop;
for i in pivots'range loop
pivots(i) := i;
end loop;
ind := x'first(1)-1;
for i in a'range(1) loop
piv := pivots'first-1;
for j in a'range(2) loop
if not Equal(a(i,j),zero)
then piv := pivots(j);
pivots(j) := pivots(i);
pivots(i) := piv;
exit;
end if;
end loop;
zero_column := true;
for j in a'first(1)..i loop
Copy(a(j,pivots(i)),tmp(j,i));
if zero_column and then not Equal(tmp(j,i),zero)
then zero_column := false;
end if;
end loop;
if piv < pivots'first or else zero_column or else Equal(tmp(i,i),zero)
then ind := i; exit;
end if;
end loop;
if zero_column
then for i in x'range loop
Copy(zero,x(i));
end loop;
Copy(one,x(ind));
elsif (ind < x'first(1)) and (a'last(1) >= a'last(2))
then for i in x'range loop
Copy(zero,x(i));
end loop;
else
if ind < x'first(1)
then ind := a'last(1)+1;
for j in tmp'range(2) loop
Copy(zero,tmp(ind,j));
end loop;
zero_column := true;
for j in a'range(1) loop
Copy(a(j,pivots(ind)),tmp(j,ind));
if zero_column and then not Equal(tmp(j,ind),zero)
then zero_column := false;
end if;
end loop;
end if;
if zero_column
then for i in x'range loop
Copy(zero,x(i));
end loop;
Copy(one,x(ind));
else
-- construct generating matrix :
for i in reverse (tmp'first(2)+1)..(ind-1) loop -- i = column
for k in tmp'first(1)..(i-1) loop
if not Equal(tmp(k,i),zero) -- make tmp(k,i) zero
then declare
aa,bb,d : number;
begin
d := gcd(tmp(i,i),tmp(k,i));
aa := tmp(i,i)/d; bb := tmp(k,i)/d;
for l in k..(i-1) loop
Mul(tmp(k,l),aa);
end loop;
Copy(zero,tmp(k,i)); --aa*tmp(k,i) - bb*tmp(i,i);
tmp(k,ind) := aa*tmp(k,ind) - bb*tmp(i,ind);
end;
Scale(tmp,k,k); -- to avoid numeric_error
end if;
end loop; -- upper half of ith colum consists of zero entries
end loop;
-- generate x by back substitution :
x(ind) := tmp(x'first,x'first);
for i in (x'first+1)..(ind-1) loop
if not Equal(tmp(i,i),zero)
then x(ind) := lcm(tmp(i,i),x(ind));
end if;
end loop;
for i in x'first..(ind-1) loop
if Equal(tmp(i,i),zero)
then Copy(zero,x(i));
else x(i) := -(tmp(i,ind)*x(ind))/tmp(i,i);
end if;
end loop;
end if;
end if;
for i in res'range loop -- take pivots into account
Copy(zero,res(i));
end loop;
for i in x'first..ind loop
Copy(x(i),res(pivots(i)));
end loop;
Vectors.Copy(res,x);
end Solve0;
procedure Solve1 ( a : in Matrix; x : in out Vectors.Vector;
b : in Vectors.Vector; fail : out boolean ) is
acc : number;
begin
fail := false;
for i in reverse x'range loop
Copy(b(i),x(i));
for j in (i+1)..x'last loop
acc := a(i,j)*x(j);
Sub(x(i),acc);
Clear(acc);
end loop;
if not Equal(x(i),zero) and then not Equal(a(i,i),zero)
then acc := Rmd(x(i),a(i,i));
if Equal(acc,zero)
then Div(x(i),a(i,i));
else fail := true;
end if;
Clear(acc);
if fail
then Vectors.Clear(x); return;
end if;
end if;
end loop;
end Solve1;
procedure Solve1 ( a : in Matrix; b : in out Vectors.Vector;
fail : out boolean ) is
acc : number;
begin
fail := false;
for i in reverse b'range loop
for j in (i+1)..b'last loop
acc := a(i,j)*b(j);
Sub(b(i),acc);
Clear(acc);
end loop;
if not Equal(b(i),zero) and then not Equal(a(i,i),zero)
then acc := Rmd(b(i),a(i,i));
if Equal(acc,zero)
then Div(b(i),a(i,i));
else fail := true;
end if;
Clear(acc);
if fail
then Vectors.Clear(b); return;
end if;
end if;
end loop;
end Solve1;
function Solve ( m : Matrix; ipvt : Standard_Integer_Vectors.Vector )
return Vectors.Vector is
x,res : Vectors.Vector(ipvt'range);
a : matrix(m'first(1)..m'last(1)-1,m'range(2));
ind : integer := a'first(1); -- index for the current row number
cnt0 : natural := 0; -- counts the number of zero rows
begin
for k in a'range(1) loop
if not Equal(m(k,k),zero) -- otherwise : skip zero row !
then for l in a'range(2) loop
Copy(m(k,l),a(ind,l));
end loop;
ind := ind + 1;
else for l in a'range(2) loop
Copy(m(k,l),a(a'last(1) - cnt0,l));
end loop;
cnt0 := cnt0 + 1;
end if;
end loop;
for i in x'range loop
Copy(zero,x(i));
end loop;
Solve0(a,x);
for k in res'range loop
res(ipvt(k)) := x(k);
end loop;
if res(res'last) < zero
then Vectors.Min(res);
end if;
Clear(a);
return res;
end Solve;
end Generic_Integer_Linear_Solvers;