LinPack.mod.m2pp
631 lines (580 with data), 19.8 kB
IMPLEMENTATION MODULE LinPack;
(*========================================================================*)
(* HINWEIS : Bitte nur die Datei LinPack.mod.m2pp editieren *)
(*========================================================================*)
(* Es sind 3 Versionen enthalten die mit *)
(* *)
(* m2pp -D __{Parameter}__ < LinPack.mod.m2pp > LibPack.mod *)
(* *)
(* mit Parameter = {INLINE|BLAS|BLASF77} erzeugt werden koennen. *)
(* *)
(* INLINE : Alle Schleifen werden expliziet ausgefuehrt *)
(* BLAS : Schleifen werden durch Routinen aus LibDBlas ausgefuehrt *)
(* BLASF77: Schleifen werden durch Routinen aus der Fortran 77 *)
(* Bibiliothek BLAS ausgefuehrt - diese muessen mit gelinkt *)
(* werden - Schnittstelle ist in LibDBlasF77.def beschreiben. *)
(* Zusaetzliche Fortran-Bibiliotheken die gelinkt werden *)
(* muessen sind compilerspezifisch. Level 1 DBLAS mit *)
(* GNU Fortran und ATLAS sind mit XDS M2 getestet. *)
(*========================================================================*)
(* Routines dgefa,dgesl,dgedi and dppfa,dppfi *)
(* from LinPack Fortran 77 libary. *)
(*------------------------------------------------------------------------*)
(* Last Change: *)
(* *)
(* 30.11.15, MRi: Adotion of code form linpack benchmark found in *)
(* Stony Brook Modula-2 examples back to double indexes *)
(* by inspecting provided code and Fortran original source *)
(* 15.04.16, MRi: Eliminating all address arithmetic which is now done in *)
(* LibDBlas. Change from row to colume order. *)
(* 17.04.16, MRi: Eliminated final errors ;-) *)
(* 18.04.16, MRi: Introduced preprocessor derectives for inline-code, *)
(* Modula-2 BLAS and Fortran BLAS. Tests look good ... *)
(* 21.04.16, MRi: In dgesl FOR k - Schleifen mit BY -1 eingefuegt, *)
(* FOR kb ... DO k=n-nk ... eleminiert *)
(* 22.04.16, MRi: Added procedure dgedi (untested) *)
(* 23.04.16, MRi: inline code in dgedi, error corrected - looks good ;-) *)
(* 24.04.16, MRi: Added procedures dppfa and dppdi *)
(* 01.11.17, MRi: Adoped to new underscore-free interface to Fortran BLAS *)
(*------------------------------------------------------------------------*)
(* Offene Punkte: *)
(* *)
(* - Weiter testen ;-) *)
(*------------------------------------------------------------------------*)
(* Testcodes in TestLinSol.mod and TestPP.mod *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LinPack.mod.m2pp,v 1.3 2018/01/16 09:19:51 mriedl Exp mriedl $ *)
(*------------------------------------------------------------------------*)
<* IF (__INLINE__) THEN *>
(* HINT: This version uses inlined code. *)
(* This is most probably slower than using BLAS level 1 routines. *)
<* END *>
<* IF (__BLAS__) THEN *>
(* HINT: This version uses Modula-2 BLAS level 1 code. *)
<* END *>
<* IF (__BLASF77__) THEN *>
(* HINT: This version uses Fortran DBLAS level 1 code. *)
(* FORTRAN 77 DOUBLE PRECISION | Fortran 90 REAL (kind=8) must be *)
(* the same type as Modula-2 LONGREAL *)
(* It can be linked against e.g. ATALS Fortran libralies. *)
<* END *>
(*------------------------------------------------------------------------*)
FROM LongMath IMPORT sqrt;
<* IF (__BLAS__) THEN *>
FROM LibDBlas IMPORT dscal,daxpy,ddot,idamax,dswap;
<* END *>
<* IF (__BLASF77__) THEN *>
FROM LibDBlasF77 IMPORT dscal,daxpy,ddot,idamax,dswap;
<* END *>
PROCEDURE dgefa(VAR A : ARRAY OF ARRAY OF LONGREAL;
lda : INTEGER;
n : INTEGER;
VAR ipvt : ARRAY OF INTEGER;
VAR info : INTEGER);
VAR t : LONGREAL;
j,k,kp1,l,nm1 : INTEGER;
<* IF (__INLINE__) THEN *>
dmax : LONGREAL;
ii : INTEGER;
<* END *>
<* IF (__BLASF77__) THEN *>
ione,num : INTEGER;
<* END *>
BEGIN
(* gaussian elimination with partial pivoting *)
IF (n > lda) THEN
(* ... *)
END;
<* IF (__BLASF77__) THEN *>
ione:=1; (* for F77 BLAS routines *)
<* END *>
info := 0;
nm1 := n - 1;
IF (nm1 >= 0) THEN
FOR k:=0 TO nm1-1 DO
kp1 := k + 1;
(* find l := pivot index *)
<* IF (__INLINE__) THEN *>
dmax:=ABS(A[k,k]);
l:=k;
FOR ii:=k TO n-1 DO
IF (ABS(A[k,ii]) > dmax) THEN
dmax:= ABS(A[k,ii]); l:=ii;
END;
END;
<* END *>
<* IF (__BLAS__) THEN *>
l := idamax(n-k,A[k,k],1) + k - 1;
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-k;
l := idamax(num,A[k,k],ione) + k - 1; (* -1 wg. F77 *)
<* END *>
ipvt[k] := l;
(* zero pivot implies this column already triangularized *)
IF (A[k,l] = 0.0) THEN
info := k;
ELSE
(* interchange if necessary *)
IF (l # k) THEN
t := A[k,l];
A[k,l] := A[k,k];
A[k,k] := t;
END;
(* compute multipliers *)
t := -1.0 / A[k,k];
<* IF (__INLINE__) THEN *>
FOR ii:=k+1 TO n-1 DO A[k,ii]:=A[k,ii]*t; END;
<* END *>
<* IF (__BLAS__) THEN *>
dscal(n-k-1,t,A[k,k+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-k-1;
dscal(num,t,A[k,k+1],ione);
<* END *>
(* row elimination with column indexing *)
FOR j:=kp1 TO n-1 DO
t := A[j,l];
IF (l # k) THEN
A[j,l] := A[j,k];
A[j,k] := t;
END;
<* IF (__INLINE__) THEN *>
FOR ii:=k+1 TO n-1 DO A[j,ii]:=A[j,ii] + A[k,ii]*t ; END;
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(n-k-1,t,A[k,k+1],1,A[j,k+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-k-1;
daxpy(num,t,A[k,k+1],ione,A[j,k+1],ione);
<* END *>
END;
END; (* IF (A[k,l] *)
END; (* FOR k *)
END; (* IF nm1 *)
ipvt[n-1] := n-1;
IF (A[n-1,n-1] = 0.0) THEN
info := n-1;
END;
END dgefa;
PROCEDURE dgesl(VAR A : ARRAY OF ARRAY OF LONGREAL;
lda : INTEGER;
n : INTEGER;
VAR ipvt : ARRAY OF INTEGER;
VAR B : ARRAY OF LONGREAL;
job : INTEGER);
(*----------------------------------------------------------------*)
(* dgesl solves the double precision system *)
(* a * x := b or trans(a) * x := b *)
(* *)
(* uses daxpy,ddot level 1 blas routine *)
(*----------------------------------------------------------------*)
VAR t : LONGREAL;
k,l,nm1 : INTEGER;
<* IF (__INLINE__) THEN *>
dot : LONGREAL;
ii : INTEGER;
<* END *>
<* IF (__BLASF77__) THEN *>
ione,num : INTEGER;
<* END *>
BEGIN
IF (n > lda) THEN
(* ... *)
END;
<* IF (__BLASF77__) THEN *>
ione:=1; (* for F77 BLAS routines *)
<* END *>
nm1 := n - 1;
IF (job # 0) THEN
(* job = nonzero, solve trans(a) * x = b *)
(* first solve trans(u)*y = b *)
FOR k:=0 TO n-1 DO
<* IF (__INLINE__) THEN *>
(* t := LibDBlasM2.ddot(k,A[k],1,B,1); *)
t:=0.0; FOR ii:=0 TO k-1 DO t:=t + A[k,ii]*B[ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
t := ddot(k,A[k,0],1,B[0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=k;
t := ddot(num,A[k,0],ione,B[0],ione);
<* END *>
B[k] := (B[k] - t) / A[k,k];
END;
(* now solve trans(l)*x = y *)
IF (nm1 >= 1) THEN
(*
FOR kb:=1 TO nm1-1 DO
k := n - kb - 1;
*)
FOR k:=nm1-1 TO 0 BY -1 DO (* TO 0 statt TO 1, MRi 26.12.16 *)
<* IF (__INLINE__) THEN *>
dot:=0.0; FOR ii:=k+1 TO n-1 DO dot:=dot + A[k,ii]*B[ii]; END;
B[k]:=B[k] + dot;
<* END *>
<* IF (__BLAS__) THEN *>
B[k]:=B[k] + ddot(n-k-1,A[k,k+1],1,B[k+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-k-1;
B[k]:=B[k] + ddot(num,A[k,k+1],ione,B[k+1],ione);
<* END *>
l := ipvt[k];
IF (l # k) THEN
t := B[l];
B[l] := B[k];
B[k] := t;
END;
END;
END;
ELSE (* job = 0 , solve a * x = b *)
IF (nm1 >= 1) THEN (* first solve l*y = b *)
FOR k:=0 TO nm1-1 DO
l := ipvt[k];
t := B[l];
IF (l # k) THEN
B[l] := B[k];
B[k] := t;
END;
<* IF (__INLINE__) THEN *>
FOR ii:=k+1 TO n-1 DO B[ii]:=B[ii] + t*A[k,ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(n-k-1,t,A[k,k+1],1,B[k+1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n-k-1;
daxpy(num,t,A[k,k+1],ione,B[k+1],ione);
<* END *>
END;
END;
(* now solve u*x = y *)
(*
FOR kb:=0 TO n-1 DO
k := n - 1 - kb;
*)
FOR k:=n-1 TO 0 BY -1 DO
B[k] := B[k] / A[k,k];
t := -B[k];
<* IF (__INLINE__) THEN *>
(* LibDBlasM2.daxpy(k,t,A[k],1,B,1); *)
FOR ii:=0 TO k-1 DO B[ii]:=B[ii] + t*A[k,ii]; END;
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(k,t,A[k,0],1,B[0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=k;
daxpy(num,t,A[k,0],ione,B[0],ione);
<* END *>
END;
END;
END dgesl;
PROCEDURE dgedi(VAR A : ARRAY OF ARRAY OF LONGREAL;
lda : INTEGER;
n : INTEGER;
VAR ipvt : ARRAY OF INTEGER;
VAR det : ARRAY OF LONGREAL; (* [0..2] *)
VAR work : ARRAY OF LONGREAL;
job : INTEGER);
VAR i,j,k,l : INTEGER;
t : LONGREAL;
exit : BOOLEAN;
<* IF (__BLASF77__) THEN *>
ione,num : INTEGER;
<* END *>
<* IF (__INLINE__) THEN *>
ii : INTEGER;
<* END *>
BEGIN
IF (lda < n) THEN
(* Fehlerbehandlung *)
END;
IF (job DIV 10 # 0) THEN (* Compute the determinant *)
det[0] := 1.0;
det[1] := 0.0;
det[2] := 1.0;
i:=0; exit:= FALSE;
WHILE (i <= n-1) AND NOT exit DO (* FOR i:=1 TO n DO *)
IF (ipvt[i] # i) THEN
det[0] := -det[0];
det[2] := -det[2];
END;
det[0] := det[0] * A[i,i];
det[2] := det[2] * A[i,i];
IF (det[0] = 0.0) THEN
exit:=TRUE; (* EXIT *)
ELSE
WHILE (ABS(det[0]) < 1.0) DO
det[0] := det[0] * 10.0;
det[1] := det[1] - 1.0;
END;
WHILE (ABS(det[0]) > 10.0) DO
det[0] := det[0] / 10.0;
det[1] := det[1] + 1.0;
END;
INC(i);
END;
END; (* WHILE *)
END;
(* Compute inverse(U) *)
<* IF (__BLASF77__) THEN *>
ione:=1; (* for F77 BLAS routines *)
<* END *>
IF ((job MOD 10) # 0) THEN
FOR k:=0 TO n-1 DO
A[k,k] := 1.0 / A[k,k];
t := -A[k,k];
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO k-1 DO A[k,ii]:=A[k,ii]*t; END;
(* LibDBlasM2.dscal(k,t,A[k],1); *)
(* LibDBlas.dscal(k,t,A[k,0],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
dscal(k,t,A[k,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=k;
dscal(num,t,A[k,0],ione);
<* END *>
FOR j:=k+1 TO n-1 DO
t := A[j,k];
A[j,k] := 0.0;
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO k DO A[j,ii]:=A[j,ii] + t*A[k,ii]; END;
(* LibDBlasM2.daxpy(k+1,t,A[k],1,A[j],1); *)
(* LibDBlas.daxpy(k+1,t,A[k,0],1,A[j,0],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(k+1,t,A[k,0],1,A[j,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=k+1;
daxpy(num,t,A[k,0],ione,A[j,0],ione);
<* END *>
END;
END;
(* Form inverse(U) * inverse(L) *)
FOR k:=n-2 TO 0 BY -1 DO
FOR i:=k+1 TO n-1 DO
work[i] := A[k,i];
A[k,i] := 0.0;
END;
FOR j:=k+1 TO n-1 DO
t := work[j];
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO A[k,ii]:=A[k,ii] + t*A[j,ii]; END;
(* LibDBlasM2.daxpy(n,t,A[j],1,A[k],1); *)
(* LibDBlas.daxpy(n,t,A[j,0],1,A[k,0],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(n,t,A[j,0],1,A[k,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n;
daxpy(num,t,A[j,0],ione,A[k,0],ione);
<* END *>
END;
l := ipvt[k];
IF (l # k) THEN
<* IF (__INLINE__) THEN *>
FOR ii:=0 TO n-1 DO
t := A[k,ii];
A[k,ii] := A[l,ii];
A[l,ii] := t;
END;
(* LibDBlasM2.dswap(n,A[k],1,A[l],1); *)
(* LibDBlas.dswap(n,A[k,0],1,A[l,0],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
dswap(n,A[k,0],1,A[l,0],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=n;
dswap(num,A[k,0],ione,A[l,0],ione);
<* END *>
END;
END;
END;
END dgedi;
PROCEDURE dppfa(VAR AP : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
n : INTEGER;
VAR Info : INTEGER);
(* factor a symmetric positive definite matrix *)
VAR j,jj,k,kj,kk : INTEGER;
s,t : LONGREAL;
<* IF (__INLINE__) THEN *>
dot : LONGREAL;
ij,jjk : INTEGER;
<* END *>
<* IF (__BLASF77__) THEN *>
num,ione : INTEGER;
<* END *>
BEGIN
Info:=0;
<* IF (__BLASF77__) THEN *>
ione:=1; (* Wg. F77 Blas *)
<* END *>
jj:=0;
FOR j:=1 TO n DO
s := 0.0;
kj := jj;
kk := 0;
FOR k:=1 TO j-1 DO
<* IF (__INLINE__) THEN *>
dot:=0.0; jjk:=jj;
FOR ij:=kk TO kk+k-2 DO dot:=dot + AP[ij]*AP[jjk]; INC(jjk); END;
t := AP[kj] - dot;
(* t := AP[kj] - ddot(k-1,AP[kk],1,AP[jj],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
t := AP[kj] - ddot(k-1,AP[kk],1,AP[jj],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=k-1;
t := AP[kj] - ddot(num,AP[kk],ione,AP[jj],ione);
<* END *>
INC(kk,k);
t := t / AP[kk-1];
AP[kj] := t;
s:=s + t*t;
INC(kj);
END;
INC(jj,j);
s := AP[jj-1] - s;
IF (s <= 0.0) THEN
Info:=j;
RETURN;
END;
AP[jj-1] := sqrt(s);
END;
END dppfa;
PROCEDURE dppdi(VAR AP : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
n : INTEGER;
VAR det : ARRAY OF LONGREAL; (* 2 Elemente *)
job : INTEGER);
VAR i,ii,j,j1,jj,k,k1,kj,kk : INTEGER;
t : LONGREAL;
exit : BOOLEAN;
<* IF (__INLINE__) THEN *>
ij,j1j,k1j : INTEGER;
<* END *>
<* IF (__BLASF77__) THEN *>
num,ione : INTEGER;
<* END *>
BEGIN
IF (job DIV 10 # 0) THEN (* Compute the determinant *)
det[0] := 1.0;
det[1] := 0.0;
det[2] := 1.0;
ii := 0;
i:=1; exit:= FALSE;
WHILE (i <= n) AND NOT exit DO (* FOR i:=1 TO n DO *)
INC(ii,i);
det[0] := det[0] * AP[ii-1]*AP[ii-1];
det[2] := det[2] * AP[ii-1]*AP[ii-1];
IF (det[0] = 0.0) THEN
exit:=TRUE; (* EXIT *)
ELSE
WHILE (det[0] < 1.0) DO
det[0] := det[0] * 10.0;
det[1] := det[1] - 1.0;
END;
WHILE (det[0] > 10.0) DO
det[0] := det[0] / 10.0;
det[1] := det[1] + 1.0;
END;
INC(i);
END;
END; (* WHILE *)
END;
(* Compute inverse(R) *)
IF ((job MOD 10) # 0) THEN
<* IF (__BLASF77__) THEN *>
ione:=1; (* Wg. F77 Blas *)
<* END *>
kk := 0;
FOR k:=1 TO n DO
k1 := kk;
INC(kk,k);
AP[kk-1] := 1.0 / AP[kk-1];
t := -AP[kk-1];
<* IF (__INLINE__) THEN *>
FOR ij:=k1 TO k1+k-2 DO AP[ij]:=AP[ij]*t; END;
(* dscal(k-1,t,AP[k1],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
dscal(k-1,t,AP[k1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num := k-1;
dscal(num,t,AP[k1],ione);
<* END *>
j1 := kk;
kj := kk + k - 1;
FOR j:=k+1 TO n DO
t := AP[kj];
AP[kj] := 0.0;
<* IF (__INLINE__) THEN *>
k1j:=k1;
FOR ij:=j1 TO j1+k-1 DO AP[ij]:=AP[ij] + t*AP[k1j]; INC(k1j); END;
(* daxpy(k,t,AP[k1],1,AP[j1],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(k,t,AP[k1],1,AP[j1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num := k;
daxpy(num,t,AP[k1],ione,AP[j1],ione);
<* END *>
INC(j1,j);
INC(kj,j);
END;
END;
(* Form inverse(R) * (inverse(R))' *)
jj := 0;
FOR j:=1 TO n DO
j1 := jj;
kj := jj;
k1 := 0;
FOR k:=1 TO j-1 DO
t := AP[kj];
<* IF (__INLINE__) THEN *>
j1j:=j1;
FOR ij:=k1 TO k1+k-1 DO AP[ij]:=AP[ij] + t*AP[j1j]; INC(j1j); END;
(* daxpy(k,t,AP[j1],1,AP[k1],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
daxpy(k,t,AP[j1],1,AP[k1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num := k;
daxpy(num,t,AP[j1],ione,AP[k1],ione);
<* END *>
INC(k1,k);
INC(kj);
END;
INC(jj,j);
t := AP[jj-1];
<* IF (__INLINE__) THEN *>
FOR ij:=j1 TO jj-1 DO AP[ij]:=AP[ij]*t; END;
(* dscal(j,t,AP[j1],1); *)
<* END *>
<* IF (__BLAS__) THEN *>
dscal(j,t,AP[j1],1);
<* END *>
<* IF (__BLASF77__) THEN *>
num:=j;
dscal(num,t,AP[j1],ione);
<* END *>
END;
END;
END dppdi;
END LinPack.