--- a
+++ b/LinPack.mod.m2pp
@@ -0,0 +1,630 @@
+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.