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;
	real procedure DSBAND (A,m,l,b,n); value m,l,b,n; integer m,l,b,n; array A;
	begin	integer i,j,k; real w;
	ELIM:	for k:= 1 step 1 until m do 
		begin	l:= if k + b <= m then l + 1 else m; MAX (i,k,l,abs (A[i,1]));
			if i != k then for j:= 1 step 1 until n do 
			begin w:= - A[i,j]; A[i,j]:= A[k,j]; A[k,j]:= w end;
			for i:= k + 1 step 1 until l do 
			begin	w:= A[i,1] / A[k,1];
				for j:= 2 step 1 until n do
				A[i, if j <= b then j - 1 else j]:= A[i,j] - w * A[k,j];
				A[i,b]:= 0
			end
		end ELIM;
		DSBAND:= PROD (k,1,m,A[k,1]);
	BACK:	for j:= b + 1 step 1 until n do for k:= m step - 1 until 1 do
		A[k,j]:= (A[k,j] - SUM (i,2, if k + b <= m then b else m - k + 1,A[k,i] * A[k + i - 1,j]))/A[k,1]
	end DSBAND;
	procedure PSP1 (A,i,j,n,B,BB,D); value n; integer i,j,n; real A; array B,BB,D;
	begin	integer p,r; real w,x;
	HA:	for r:= 1 step 1 until n do
		begin	j:= i:= r; D[r]:= A; BB[r]:= SUM (j,r+1,n,A ** 2);
			B[r]:= sqrt (BB[r]); if B[r] < &apos;-200 then goto HB;
			j:= r+1; if A > 0 then B[r]:= -B[r]; A:= A-B[r]; w:= A * B[r];
			for j:= r+1 step 1 until n do D[j]:= A;
			for p:= r+1 step 1 until n do
			begin	j:= p;
				B[p]:= (SUM (i,r+1,p-1,A * D[i]) + SUM (j,p,n,A * D[j]))/w
			end;
			x:= SUM (p,r+1,n,D[p] * B[p])/(2 * w);
			for j:= r+1 step 1 until n do B[j]:= D[j] * x + B[j];
			for i:= r+1 step 1 until n do for j:= i step 1 until n do
			A:= D[i] * B[j] + B[i] * D[j] + A;
		HB:
		end
	end PSP1;
	procedure PREP (D,B,n,eps,E); value n; integer n; real eps; array D,B,E;
	begin	integer r; real s,x;
		s:= 0;
		for r:= 1 step 1 until n do
		begin	x:= (if r = 1 then 0 else abs (B[r-1]))
			+ abs (D[r]) + abs (B[r]);
			if x > s then s:= x
		end;
		E[1]:= 0; E[2]:= eps * s; E[3]:= 1.1 * s
	end PREP;
	real procedure SEIVA (D,BB,n,n1,n2,k,E); value n; integer n,n1,n2,k; array D,BB,E;
	begin	integer r,t; real x;
		real procedure SDET (q,q2); value q,q2; integer q,q2;
		begin	integer p; real d0,d1,d2;
			p:= 0; SDET:= E[5];
			d1:= t; d2:= (x-D[q]) * d1; goto DB;
		DA:	q:= q+1; d0:= d1; d1:= d2; d2:= (x-D[q]) * d1 - BB[q-1] * d0;
		DB:	if d2 > 0 == d1 <= 0 then p:= p+1;
			if p <= r then
			begin	if q < q2 then goto DA;
				if x < E[4] then E[4]:= x;
				if abs (d2) > E[5] then E[5]:= abs (d2);
				SDET:= d2
			end
		end SDET;
	GA:	if k = 0 then n2:= 0;
		if k = n2 then
		begin	E[4]:= E[3]; E[5]:= 0; n1:= n2+1;
		GC:	n2:= n2+1; if 1-BB[n2]/E[3] ** 2 != 1 then goto GC
		end;
		k:= k+1; r:= k-n1+1; t:= EVEN (r);
		SEIVA:= ZERO (x,E[4],-E[3],SDET (n1,n2),E)
	end SEIVA;
	procedure SEIVEC (D,B,n,n1,n2,x,V);
		value n,n1,n2,x; integer n,n1,n2; real x; array D,B,V;
	begin	integer i,p,q; real x1;
	WA:	p:= n1-1; q:= n2+1; V[p]:= V[q]:= 1;
	WB:	i:= p:= p + 1; if p = n2 then goto WD;
		V[p]:= (if p = n1 then (x - D[p]) else ((x - D[p]) * V[p - 1] - B[p - 1] * V[p - 2])) / B[p];
		if abs (V[p]) >= abs (V[p - 1]) then goto WB;
		if p >= q then goto WD;
	WC:	i:= q:= q - 1; if q = n1 then goto WD;
		V[q]:= (if q = n2 then (x - D[q]) else ((x - D[q]) * V[q + 1] - B[q] * V[q + 2])) / B[q - 1];
		if abs (V[q]) >= abs (V[q + 1]) then goto WC;
		if p < q then goto WB;
	WD:	V[i]:= 1/sqrt (SUM (p, n1 - 1, i - 2, V[p] ** 2)/V[i - 1] ** 2 + 1
		        + SUM (p, i + 2, n2 + 1, V[p] ** 2)/V[i + 1] ** 2);
		x1:= V[i]/V[i - 1]; for p:= i - 1 step -1 until n1 do V[p]:= V[p - 1] * x1;
		x1:= V[i]/V[i + 1]; for p:= i + 1 step  1 until n2 do V[p]:= V[p + 1] * x1;
		for p:= 1 step 1 until n1 - 1, n2 + 1 step 1 until n do V[p]:= 0
	end SEIVEC;
	procedure TRASF1 (A,i,j,n,B,V); value n; integer i,j,n; real A; array B,V;
	begin	real x1,f1;
		for i:= n-1 step -1 until 1 do
		begin	if abs (B[i]) >= &apos;-200 then
			begin	j:= i+1; x1:= A;
				f1:= SUM (j,i+1,n,A*V[j])/(x1*B[i]);
				for j:= i+1 step 1 until n do V[j]:= A*f1+V[j]
			end
		end
	end TRASF1;
	procedure SYMEVE (A,i,j,n,eps,OVA,OVEC);
		value n; integer i,j,n; real A,eps; procedure OVA,OVEC;
	begin	integer k,n1,n2,p; real x; array B,BB,D[1:n],V[0:n+1],E[1:5];
		PSP1 (A,i,j,n,B,BB,D); PREP (D,B,n,eps,E); k:= 0;
	next:	x:= SEIVA (D,BB,n,n1,n2,k,E); OVA (x);
		SEIVEC (D,B,n,n1,n2,x,V); TRASF1 (A,i,j,n,B,V); OVEC (p,1,n,V[p]);
		if k < n then goto next
	end SYMEVE;
	real procedure CSQRT (a,b,rp,ip); value a,b; real a,b,rp,ip;
	begin	rp:= sqrt ((abs (a) + sqrt (a * a + b * b)) / 2.0);
		ip:= b:= b / (2.0 * rp);
		if a < 0 then begin ip:= if b >= 0 then rp else -rp; rp:= abs (b) end;
		CSQRT:= rp
	end	CSQRT;
	real procedure CZERO (x,y,r,s,e); real x,y,r,s; array e;
	begin	real a,b,c,d,e1,e2,e3,g,h,r0,r1,r2,s0,s1,s2,t,u;
		e1:= e[1] ** 2; e2:= e[2] ** 2; e3:= e[3] ** 2; a:= x;
		g:= sqrt (a * a + y ** 2) * .1 + .1; h:= 0; c:= -.5; d:= 0;
		x:= a + g; r0:= r; s0:= s; x:= a - g; r1:= r; s1:= s;
		r0:= r0 - r1; s0:= s0 - s1; x:= a;
		r2:= r; s2:= s;
	LL:	r1:= r1 - r2; s1:= s1 - s2;
		t:= r0 * c - s0 * d - r1; u:= r0 * d + s0 * c - s1;
		a:= (t - r1) * c - (u - s1) * d - r1; b:= (t - r1) * d + (u - s1) * c - s1;
		r0:= t * c - u * d; s0:= t * d + u * c;
		t:= -2.0 * ((1.0 + c) * r2 - d * s2); u:= -2.0 * ((1.0 + c) * s2 + d * r2);
		CSQRT (a * a - b * b + 2.0 * (t * r0 - u * s0),2.0 * (a * b + t * s0 + u * r0),c,d);
		if a * c + b * d < 0 then begin c:= - c; d:= - d end;
		a:= a + c; b:= b + d;
		c:= (a * t + b * u) / (a * a + b * b); d:= (a * u - b * t) / (a * a + b * b);
		a:= sqrt (c * c + d * d) / 10;
		if a > 1 then begin c:= c/a; d:= d/a end;
		a:= g * c - h * d; h:= g * d + h * c; g:= a; x:= x + g; y:= y + h;
		r0:= r1; s0:= s1; r1:= r2; s1:= s2; r2:= r; s2:= s;
		if g * g + h * h > (x ** 2 + y ** 2) * e1 + e2 | r2 * r2 + s2 * s2 > e3 then goto LL;
		CZERO:= x
	end	CZERO;
	real procedure CPROD (k,a,b,rk,sk,rp,ip); value b; integer k,a,b; real rk,sk,rp,ip;
	begin	real p,q,r,s; p:= 1; ip:= q:= 0;
		for k:= a step 1 until b do
		begin	r:= rk; s:= sk; ip:= p * s + q * r;
			p:= p * r - q * s; q:= ip
		end;
		CPROD:= rp:= p
	end	CPROD;
	real procedure CPOL (A,k,n,x,y,rp,ip); value x,y,n; real A,x,y,rp,ip; integer k,n;
	begin	real p,q,b0,b1,b2;
		p:= 2 * x; q:= x * x + y * y; b2:= b1:= 0;
		for k:= 0 step 1 until n - 1 do
		begin b0:= b1; b1:= b2; b2:= p * b1 - q * b0 + A end;
		ip:= y * b2; k:= n; CPOL:= rp:= x * b2 - q * b1 + A
	end	CPOL;
	real procedure SYMDET (A,i,j,n); value n; integer i,j,n; real A;
	begin	integer k; real d,r; array v[1:n];
		d:= 1;
		for k:= 1 step 1 until n do
		begin	j:= k; for i:= 1 step 1 until k do v[i]:= A;
			r:= v[k] - SUM (i,1,k-1,v[i] ** 2);
			d:= r * d; i:= k; A:= r:= 1 / sqrt (r);
			for j:= k+1 step 1 until n do
			begin	i:= k; A:= (A - SUM (i,1,k-1,A * v[i])) * r end
		end	LU;
		SYMDET:= d
	end	SYMDET;
	procedure SYMSOL (A,i,j,n,b); value n; integer i,j,n; real A; array b;
	begin	for j:= 1 step 1 until n do b[j]:= (b[j] - SUM (i,1,j-1,A * b[i])) * A;
		for i:= n step -1 until 1 do
		begin	j:= i; b[i]:= A * (b[i] - SUM (j,i+1,n,A * b[j])) end
	end	SYMSOL;
	procedure SYMINV (A,i,j,n); value n; integer i,j,n; real A;
	begin	integer k; array v[1:n];
		for k:= 1 step 1 until n do
		begin	i:= j:= k; v[k]:= A;
			for j:= k+1 step 1 until n do
			begin	i:= k; A:= v[j]:= - SUM (i,k,j-1,A * v[i]) * A end;
			for i:= 1 step 1 until k do
			begin	j:= k; A:= SUM (j,k,n,A * v[j]) end
		end
	end	SYMINV;
	procedure syminv (A,i,j,n); value n; integer i,j,n; real A;
	begin	integer k; array v[1:n];
		for k:= 1 step 1 until n do
		begin	i:= j:= k; v[k]:= A;
			for j:= k+1 step 1 until n do v[j]:= - SUM (i,k,j-1,A * v[i]) * A;
			i:= j:= k; A:= SUM (j,k,n,v[j] ** 2)
		end
	end	syminv;
	comment	AP 224
		SYMDET1:=  determinant of the n-th order  symmetric  positive definite matrix M  which is defined as
		follows:  the actual parameter for A  -  being a subscripted real variable whose indices  (or index)
	depend(s) on the actual parameters for i and j  -  is the (i, j)-th element of M for each i and j satisfying
	1 <= i <= j <= n. Thus one needs to give only  the upper triangle of M.  In order to avoid waste of space,  one
	may give this  triangle  in a one-dimensional  array.  E.g.,  if the upper triangle  of M  is given in array
	C[1 : n * (n + 1) % 2]  columnwise,  i.e.  the  columns  one  after  the  other,  and  the successive values
	(j - 1) * j % 2  have been recorded  in an auxiliary  integer  array J[1 : n],  then the appropriate call of
	SYMDET1 reads:
						    SYMDET1 (C[i + J[j]], i, j, n).

	The method used is the square root method  of Cholesky,  yielding an upper triangle which,  premultiplied by
	its transpose,  gives the original matrix.  SYMDET1 replaces the elements of M by the corresponding elements
	of this upper triangle. It uses the non-local real procedure SUM (= AP 119);

	real procedure SYMDET1 (A,i,j,n); value n; integer i,j,n; real A;

	begin	integer k; real d,r; array v[1:n];

		d:= 1;
		for k:= 1 step 1 until n do

		begin	j:= k; for i:= 1 step 1 until k do v[i]:= A;
			i:= k; A:= r:= sqrt (v[k] - SUM (i,1,k-1,v[i] ** 2));
			d:= r * d;
			for j:= k+1 step 1 until n do
			begin i:= k; A:= (A - SUM (i,1,k-1,A * v[i])) / r end
		end LU;

		SYMDET1:= d ** 2

	end SYMDET1;

	comment	AP 225
		SYMSOL1 replaces the vector given in  array b[1 : n], by the solution vector x of the linear system:
		U transpose * U * x = b,  where U is an upper triangle which is defined by the actual parameters for
	A,  i,  j and n  in the same way as the upper triangle of M in SYMDET1 (= AP 224).  Consequently,  a call of
	SYMSOL1, following a call of SYMDET1  with the same actual parameters for A,  i,  j  and  n,  has the effect
	that b is replaced by the solution vector x of the linear system M * x = b.
	SYMSOL1 leaves the elements A unaltered. It uses the non-local real procedure SUM (= AP 119);

	procedure SYMSOL1(A,i,j,n,b); value n; integer i,j,n; real A; array b;
	begin	real r;
		for j:= 1 step 1 until n do
		begin	i:= j; r:= A;
			b[j]:= (b[j] - SUM (i,1,j-1,A * b[i])) / r
		end;
		for i:= n step -1 until 1 do
		begin	j:= i; r:= A;
			b[i]:= (b[i] - SUM (j,i+1,n,A * b[j])) / r
		end
	end SYMSOL1;

	comment	AP 226
		SYMINV1 replaces the matrix elements A by the corresponding upper triangular elements of the inverse
		of  U transpose  * U,  where U is an upper triangle which is defined by the actual parameters in the
	same way as the upper triangle of  M  in SYMDET1 (= AP 224).  Consequently,  a call of SYMINV1,  following a
	call of SYMDET1  with the same actual parameters,  has the effect  that  the upper triangle of the symmetric
	positive definite matrix M is replaced by the upper triangle of the inverse of M.
	SYMINV1 uses the non-local real procedure SUM (= AP 119);

	procedure SYMINV1 (A,i,j,n); value n; integer i,j,n; real A;
	begin	integer k; real r; array v[1:n];
		for k:= 1 step 1 until n do
		begin	i:= j:= k; A:= v[k]:= 1 / A;
			for j:= k+1 step 1 until n do
			begin	i:= j; r:= A; i:= k;
				A:= v[j]:= - SUM (i,k,j-1,A * v[i]) / r
			end;
			for i:= 1 step 1 until k do
			begin	j:= k; A:= SUM (j,k,n,A * v[j]) end
		end
	end SYMINV1;

	comment	AP 227
		syminv1 calculates the main diagonal of the inverse of U transpose * U, where U is an upper triangle
		which is defined by the actual parameters for A, i, j and n in the same way as the upper triangle of
	M in SYMDET1 (= AP 224).  The calculated diagonal elements are delivered in array d[1 : n].  Consequently, a
	call of syminv1,  following a call of SYMDET1 with the same actual parameters for A,  i,  j and n,   has the
	effect that the diagonal elements of the inverse of the symmetric positive definite matrix  M  are delivered
	in array d.
	syminv1 leaves the elements A unaltered. It uses the non-local real procedure SUM (= AP 119);

	procedure syminv1(A,i,j,n,d); value n; integer i,j,n; real A; array d;
	begin	integer k; real r;
		for k:= 1 step 1 until n do
		begin	i:= j:= k; d[k]:= 1 / A;
			for j:= k+1 step 1 until n do
			begin	i:= j; r:= A;
				d[j]:= - SUM (i,k,j-1,A * d[i]) / r
			end;
			d[k]:= SUM (j,k,n,d[j] ** 2)
		end
	end syminv1;

	comment	AP 228
		SYMDET2:= determinant of the n-th order symmetric positive definite matrix,  given in  integer array
		A[1 : n * (n + 1) % 2] in such a way that,  for all i and j satisfying 1 <= i <= j <= n,  the (i, j)-th
	element is A[i + (j - 1) * j % 2].  The method used is the square root method of Cholesky, yielding an upper
	triangle U which, premultiplied by its transpose, gives alfa * matrix A.  The elements of U are written over
	the corresponding elements of A.  The scaling factor alfa must be chosen so that the maximal element of U is
	just within the integer capacity,  in order to obtain a reasonably accurate representation of U.  In view of
	the definiteness of A this means that alfa must be slightly less (but not too critically,  on account of the
	inexactness of the arithmetic)  than the square of the integer capacity divided by the maximal element of A.
	Also,  one may use SYMDET2 with real array A,  in which case 1.0 is the most obvious value of alfa.  If A is
	negative definite, one may use SYMDET2 with alfa negative.
	SYMDET2 uses the non-local real procedure INPROD (= AP 120);

	real procedure SYMDET2 (A,n,alfa);

	value n,alfa; integer n; real alfa; integer array A;

	begin	integer i,j,k,kk,kj; real d;

		d:= 1; kk:= 0;
		for k:= 1 step 1 until n do

		begin	kk:= kk+k; A[kk]:=
			sqrt(A[kk] * alfa - INPROD(i,1-k,-1,A[kk+i],A[kk+i]));
			d:= A[kk] * d; kj:= kk;
			for j:= k+1 step 1 until n do

			begin	kj:= kj+j-1;
				A[kj]:= (A[kj] * alfa
				- INPROD (i,1-k,-1,A[kj+i],A[kk+i])) / A[kk]
			end
		end LU;

		SYMDET2:= d ** 2 / alfa ** n

	end SYMDET2;

	comment	AP 229
		SYMSOL2  replaces the vector given in  real array b[1 : n],  by the solution vector x  of the linear
		system:  U transpose * U * x = alfa * b,  where U is an upper triangle,  given in integer  (or real)
	array A[1 : n * (n + 1) % 2]  in such a way that,  for all i and j satisfying  1 <= i <= j <= n,  the (i, j)-th
	element is A[i + (j - 1) * j % 2].  The scaling factor alfa is chosen in relation to the scaling of U.  Con-
	sequently,  the call SYMSOL2 (A, n, alfa, b) following the call SYMDET2 (A, n, alfa)  (viz.: AP 228) has the
	effect that b is replaced by the solution vector x of the linear system A * x = b.
	SYMSOL2 leaves the elements of A unaltered. It uses the non-local real procedure SUM (= AP 119);

	procedure SYMSOL2 (A,n,alfa,b);

	value n,alfa; integer n; real alfa; integer array A; real array b;

	begin	integer i,j,j0; integer array J[1:n];

		j0:= 0;

		for j:= 1 step 1 until n do

		begin	b[j]:=
			(b[j] * alfa - SUM (i,1,j-1,A[i+j0] * b[i]))/A[j+j0];

			J[j]:= j0; j0:= j0+j
		end;

		for i:= n step -1 until 1 do
		b[i]:= (b[i] - SUM (j,i+1,n,A[i + J[j]] * b[j]))/A[i + J[i]]
	end SYMSOL2;