IMPLEMENTATION MODULE NumAlLib1;
(*------------------------------------------------------------------------*)
(* Einige Routinen aus der NUMAL Numerical Algol library (Stichting CWI) *)
(* Teil 1 (numal5p1) *)
(* Some routinen from the NUMAL Numerical Algol library (Stichting CWI) *)
(*------------------------------------------------------------------------*)
(* Letzte Veraenderung *)
(* *)
(* 13.12.16, MRi: Prozedur VecVec und. ChlDec1 eingefuegt *)
(* 16.12.16, MRi: SeqVec,SymMatVec *)
(* 16.12.16, MRi: ChlInv1 eingefuegt *)
(* 29.07.17, MRi: Erstellen der ersten Versionen von ChlSol1 *)
(* 31.07.17, MRi: Korrektur in SeqVec wg. offener Felder *)
(*------------------------------------------------------------------------*)
(* Testprogram in TestChol *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: NumAlLib1.mod,v 1.1 2018/01/16 09:20:42 mriedl Exp $ *)
FROM SYSTEM IMPORT TSIZE;
FROM Storage IMPORT ALLOCATE,DEALLOCATE;
FROM LongMath IMPORT sqrt;
PROCEDURE VecVec( L,U : CARDINAL;
Shift : CARDINAL;
VAR A : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Scalar product of the vector given in array A[L:U] and *)
(* array B[Shift+L : Shift+U]. *)
(* *)
(* The meaning of the formal parameters is: *)
(* *)
(* L,U : lower and upper bound of the running subscript *)
(* Shift : index-shifting parameter of the vector b; *)
(* A,B: : array a[L:U], b[L+Shift : U+Shift]. *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
VAR k : CARDINAL;
s : LONGREAL;
BEGIN
s:= 0.0;
FOR k:=L TO U DO s:=s + A[k-1]*B[Shift+k-1]; END;
RETURN s;
END VecVec;
PROCEDURE SeqVec( L,U : CARDINAL;
IL : CARDINAL;
Shift : CARDINAL;
VAR A : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* Scalar product of the vectors given in array *)
(* A[IL:il + (U+L-1)*(U-L)/2] and array B[Shift+L : Shift+U], *)
(* where the elements of the first vector are *)
(* A[IL+(j+L-1)*(j-L)/2] for j = l, ..., u. *)
(* *)
(* The meaning of the formal parameters is: *)
(* *)
(* L,U : lower and upper bound of the running subscript; *)
(* IL : lower bound of the vector a; *)
(* Shift : index-shifting parameter of the vector b; *)
(* A,B : array A[p : q], B[l + shift, U + Shift]; *)
(* the values of p and q should satisfy p <= IL and *)
(* q >= IL+(U+L-1)*(U-L)/2). *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
VAR s : LONGREAL;
l,il : CARDINAL;
BEGIN
s:= 0.0; il:=IL;
FOR l:=L TO U DO
s:=s + A[il-1]*B[l + Shift - 1]; INC(il,l);
END;
RETURN s;
END SeqVec;
PROCEDURE SymMatVec( L,U : CARDINAL;
I : CARDINAL;
VAR A : ARRAY OF LONGREAL;
VAR B : ARRAY OF LONGREAL) : LONGREAL;
(*----------------------------------------------------------------*)
(* The scalarproduct of (a part of) A vector and *)
(* (a part of) a row of a symmetric matrix, whose uppertriangle *)
(* is given columnwise in an one-dimensional array. *)
(* *)
(* The value of the scalar product of the vectors given *)
(* in array A[p:q] and array B[l:u], where the elements *)
(* of the first vector are: if l<i then A[(i-1)*i//2 + j], *)
(* j=l,..., min(u, i-1) and A[(j-1)*j//2 + i], j=i,..., u, *)
(* respectively , otherwise A[(j-1)*j//2 + i], j=l,..., u. *)
(* *)
(* The meaning of the formal parameters is: *)
(* *)
(* L,U : lower and upper bound of the vector B, respectively *)
(* L >=1; *)
(* I : row index of the matrix A; I >= 1; *)
(* A : a one-dimensional array A[p:q] with *)
(* if i > l then p=(i-1)*i//2 + l else p=(l-1)*l//2 + i *)
(* and if i > u then q=(i-1)*i//2 + u else q=(u-1)*u//2+i *)
(* B : a one-dimensional array B[l:u]; *)
(* *)
(* Source: NUMAL Numerical Algol library (Stichting CWI) *)
(*----------------------------------------------------------------*)
VAR k,m : CARDINAL;
u : CARDINAL;
BEGIN
IF (L > I) THEN m:=L ELSE m:=I; END;
k:= m*(m-1) DIV 2;
IF (I <= U) THEN u:=I-1 ELSE u:=U; END;
RETURN VecVec(L,u,k,B,A) + SeqVec(m,U,k + I,0,A,B);
END SymMatVec;
(*------------------------------------------------------------------------*)
PROCEDURE ChlDec1(VAR A : ARRAY OF LONGREAL;
N : CARDINAL;
VAR iFehl : INTEGER;
eps : LONGREAL);
VAR j,k,kk,kj,low,up : CARDINAL;
r,EpsNorm : LONGREAL;
BEGIN
IF (N = 0) THEN iFehl:= -1; RETURN; END;
r:= 0.0; kk:= 0;
FOR k:=1 TO N DO
INC(kk,k);
IF (A[kk-1] > r) THEN r:= A[kk-1]; END;
END;
EpsNorm := eps*r;
kk:=0; k:=1;
LOOP (* FOR k:=1 TO n DO *)
INC(kk,k);
low := kk - k + 1;
up := kk - 1;
r:= A[kk-1] - VecVec(low,up,0,A,A);
IF (r <= EpsNorm) THEN
iFehl := k - 1;
EXIT;
END;
r:= sqrt(r); A[kk-1]:=r;
kj:= kk + k;
FOR j:=k+1 TO N DO
A[kj-1]:= (A[kj-1] - VecVec(low,up,kj - kk,A,A)) / r;
INC(kj,j);
END;
INC(k);
IF (k > N) THEN iFehl:=0; EXIT; END;
END;
END ChlDec1;
PROCEDURE ChlSol1( N : CARDINAL;
VAR A : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
VAR B : ARRAY OF LONGREAL);
VAR i,ii : CARDINAL;
BEGIN
ii:=0;
FOR i:=1 TO N DO
ii:=ii + i;
B[i-1]:= (B[i-1] - VecVec(1, i - 1, ii - i, B, A)) / A[ii-1];
END;
FOR i:=N TO 1 BY -1 DO
B[i-1]:= (B[i-1] - SeqVec(i+1,N,ii+i,0,A,B)) / A[ii-1];
ii:=ii - i;
END
END ChlSol1;
PROCEDURE ChlInv1(VAR A : ARRAY OF LONGREAL;
N : CARDINAL);
VAR i,ii,i1,j,ij : CARDINAL;
R : LONGREAL;
U : POINTER TO ARRAY [1..MAX(INTEGER)] OF LONGREAL;
BEGIN
ALLOCATE(U,N*TSIZE(LONGREAL));
ii:= (N + 1) * N DIV 2;
FOR i:=N TO 1 BY -1 DO
R:= 1.0 / A[ii-1];
i1:= i + 1;
ij:= ii + i;
FOR j:= i1 TO N DO
U^[j]:= A[ij-1];
ij:=ij + j;
END;
FOR j:= N TO i1 BY -1 DO
ij:=ij - j;
A[ij-1]:= -SymMatVec(i1,N,j,A,U^)*R;
END;
A[ii-1]:= (R - SeqVec(i1,N,ii+i,0,A,U^))*R;
ii:=ii - i;
END;
DEALLOCATE(U,N*TSIZE(LONGREAL));
END ChlInv1;
END NumAlLib1.