real procedure DET (A,n,p); value n; integer n; array A; integer array p;
	begin	integer i,j,k; real d,r,s; array v[1: n];
		for i:= 1 step 1 until n do v[i]:= sqrt (INPROD (j,1,n,A[i,j],A[i,j]));
		d:= 1;
		for k:= 1 step 1 until n do
		begin	r:= - 1;
			for i:= k step 1 until n do
			begin	A[i,k]:= A[i,k] - INPROD (j,1,k - 1,A[i,j],A[j,k]);
				s:= abs (A[i,k]) / v[i];
				if s > r then begin r:= s; p[k]:= i end
			end LOWER;
			v[p[k]]:= v[k];
			for j:= 1 step 1 until n do
			begin	r:= A[k,j]; A[k,j]:= if j <= k then A[p[k],j] else
				(A[p[k],j] - INPROD (i,1,k - 1,A[k,i],A[i,j])) / A[k,k];
				if p[k] != k then A[p[k],j]:= - r
			end UPPER;
			d:= A[k,k] * d
		end LU;
		DET:= d
	end DET;
	procedure SOL (LU,b,n,p); value n; integer n; array LU,b; integer array p;
	begin	integer i,k; real r;
		for k:= 1 step 1 until n do
		begin	r:= b[k];
			b[k]:= (b[p[k]] - INPROD (i,1,k - 1,LU[k,i],b[i])) / LU[k,k];
			if p[k] != k then b[p[k]]:= - r
		end;
		for k:= n step - 1 until 1 do
		b[k]:= b[k] - INPROD (i,k + 1,n,LU[k,i],b[i])
	end SOL;
	real procedure DETSOL (A,b,n);
	value n; integer n; array A,b;
	begin	integer array p[1:n];
		DETSOL:= DET (A,n,p); SOL (A,b,n,p)
	end DETSOL;
	procedure INV (LU,n,p); value n; integer n; array LU; integer array p;
	begin	integer i,j,k; real r; array v[1: n];
		for k:= n step - 1 until 1 do
		begin	for j:= k + 1 step 1 until n do
			begin	v[j]:= LU[k,j]; LU[k,j]:= 0 end;
			LU[k,k]:= 1 / LU[k,k];
			for j:= k - 1 step - 1 until 1 do
			LU[k,j]:= - INPROD (i,j + 1,k,LU[k,i],LU[i,j]) / LU[j,j];
			for j:= 1 step 1 until n do
			LU[k,j]:= LU[k,j] - INPROD (i,k + 1,n,v[i],LU[i,j])
		end;
		for k:= n step - 1 until 1 do
		begin	if p[k] != k then for i:= 1 step 1 until n do
			begin	r:= LU[i,k]; LU[i,k]:= - LU[i,p[k]]; LU[i,p[k]]:= r end
		end
	end INV;
	real procedure DETINV (A,n); value n; integer n; array A;
	begin	integer array p[1:n];
		DETINV:= DET (A,n,p); INV (A,n,p)
	end DETINV;