DEFINITION MODULE LinPack;
(*------------------------------------------------------------------------*)
(* Routines dgefa,dgesl,dgedi and dppfa,dppfi *)
(* from LinPack Fortran 77 libary. *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id$ *)
PROCEDURE dgefa(VAR A : ARRAY OF ARRAY OF LONGREAL;
lda : INTEGER;
n : INTEGER;
VAR ipvt : ARRAY OF INTEGER;
VAR info : INTEGER);
(*----------------------------------------------------------------*)
(* dgefa factors a double precision matrix by gaussian *)
(* elimination. *)
(* *)
(* dgefa is usually called by dgeco, but it can be called *)
(* directly with a saving in time if rcond is not needed. *)
(* (time for dgeco) := (1 + 9/n)*(time for dgefa) . *)
(* *)
(* on entry *)
(* *)
(* A LONGREAL precision[n][lda] *)
(* the matrix to be factored. *)
(* *)
(* lda integer *)
(* the leading dimension of the array a . *)
(* *)
(* n integer *)
(* the order of the matrix a . *)
(* *)
(* on RETURN *)
(* *)
(* A an upper triangular matrix and the multipliers *)
(* which were used to obtain it. *)
(* the factorization can be written a := l*u where *)
(* l is a product of permutation and unit lower *)
(* triangular matrices and u is upper triangular. *)
(* *)
(* ipvt integer[n] *)
(* an integer vector of pivot indices. *)
(* *)
(* info integer *)
(* := 0 normal value. *)
(* := k If u[k][k] .eq. 0.0 . this is not an error *)
(* condition for this subroutine, but it does *)
(* indicate that dgesl or dgedi will divide by *)
(* zero. *)
(* If called use rcond in dgeco for a reliable *)
(* indication of singularity. *)
(* *)
(* linpack. this version dated 08/14/78 . *)
(* cleve moler, university of New Mexico, argonne national lab. *)
(*----------------------------------------------------------------*)
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 *)
(* using the factors computed by dgeco or dgefa. *)
(* *)
(* on entry *)
(* *)
(* A double precision[n][lda] *)
(* the output from dgeco or dgefa. *)
(* *)
(* lda integer *)
(* the leading dimension of the array a . *)
(* *)
(* n integer *)
(* the order of the matrix a . *)
(* *)
(* ipvt integer[n] *)
(* the pivot vector from dgeco or dgefa. *)
(* *)
(* B double precision[n] *)
(* the right hand side vector. *)
(* *)
(* job integer *)
(* := 0 to solve a*x := b , *)
(* := nonzero to solve trans(a)*x := b where *)
(* trans(a) is the transpose. *)
(* *)
(* n RETURN *)
(* *)
(* B the solution vector x . *)
(* *)
(* error condition *)
(* *)
(* a division by zero will occur if the input factor contains *)
(* a zero on the diagonal. technically this indicates *)
(* singularity but it is often caused by improper arguments or *)
(* improper setting of lda . it will not occur if the *)
(* subroutines are called correctly and if dgeco has set *)
(* rcond .gt. 0.0 or dgefa has set info .eq. 0 . *)
(* *)
(* to compute inverse(a) * c where c is a matrix *)
(* with p columns *)
(* dgeco(a,lda,n,ipvt,rcond,z) *)
(* IF (!rcond is too small){ *)
(* for (j:=0,j<p,j++) *)
(* dgesl(a,lda,n,ipvt,c[j][0],0); *)
(* END; *)
(* cleve moler, university of new mexico, argonne national lab. *)
(*----------------------------------------------------------------*)
PROCEDURE dgedi(VAR A : ARRAY OF ARRAY OF LONGREAL;
lda : INTEGER;
n : INTEGER;
VAR ipvt : ARRAY OF INTEGER;
VAR det : ARRAY OF LONGREAL; (* [0..1] *)
VAR work : ARRAY OF LONGREAL;
job : INTEGER);
(*----------------------------------------------------------------*)
(* dgedi computes the determinant and inverse of a matrix *)
(* using the factors computed by dgeco or dgefa. *)
(* *)
(* on entry *)
(* *)
(* A the output from dgeco or dgefa. *)
(* lda the leading dimension of the array a . *)
(* n the order of the matrix a . *)
(* ipvt the pivot vector from dgeco or dgefa. *)
(* work work vector. contents destroyed. *)
(* job = 11 both determinant and inverse. *)
(* = 01 inverse only. *)
(* = 10 determinant only. *)
(* *)
(* on return *)
(* *)
(* A inverse of original matrix if requested. *)
(* otherwise unchanged. *)
(* det double precision(2) *)
(* determinant of original matrix if requested. *)
(* otherwise not referenced. *)
(* determinant = det(1) * 10.0**det(2) *)
(* with 1.0 .le. dabs(det(1)) .lt. 10.0 *)
(* or det(1) .eq. 0.0 . *)
(* *)
(* error condition *)
(* *)
(* a division by zero will occur if the input factor contains *)
(* a zero on the diagonal and the inverse is requested. *)
(* it will not occur if the subroutines are called correctly *)
(* and if dgeco has set rcond .gt. 0.0 or dgefa has set *)
(* info .eq. 0 . *)
(* *)
(* linpack. this version dated 08/14/78 . *)
(* Cleve Moler, university of new mexico, argonne national lab. *)
(*----------------------------------------------------------------*)
PROCEDURE dppfa(VAR AP : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
n : INTEGER;
VAR Info : INTEGER);
(*----------------------------------------------------------------*)
(* dppfa factors a double precision symmetric positive definite *)
(* matrix stored in packed form. *)
(* *)
(* dppfa is usually called by dppco, but it can be called *)
(* directly with a saving in time if rcond is not needed. *)
(* (time for dppco) = (1 + 18/n)*(time for dppfa) . *)
(* *)
(* on entry *)
(* *)
(* AP The packed form of a symmetric matrix A. The *)
(* columns of the upper triangle are stored *)
(* sequentially in a one-dimensional array of length *)
(* n*(n+1)/2. *)
(* n the order of the matrix a . *)
(* *)
(* on return *)
(* *)
(* AP an upper triangular matrix r , stored in packed *)
(* form, so that a = trans(r)*r . *)
(* info integer *)
(* = 0 for normal return. *)
(* = k if the leading minor of order k is not *)
(* positive definite. *)
(* *)
(* LinPack. This version dated 08/14/78 . *)
(* Cleve Moler, university of new mexico, argonne national lab. *)
(* *)
(* Modula-2 version M.Riedl, 23.04.2016 *)
(*----------------------------------------------------------------*)
PROCEDURE dppdi(VAR AP : ARRAY OF LONGREAL; (* SUPERVEKTOR *)
n : INTEGER;
VAR det : ARRAY OF LONGREAL; (* 2 Elemente *)
job : INTEGER);
(*----------------------------------------------------------------*)
(* dppdi computes the determinant and inverse *)
(* of a double precision symmetric positive definite matrix *)
(* using the factors computed by dppco or dppfa . *)
(* *)
(* on entry *)
(* *)
(* AP the output from dppco or dppfa (n*(n+1)/2) *)
(* n the order of the matrix a . *)
(* job = 11 both determinant and inverse. *)
(* = 01 inverse only. *)
(* = 10 determinant only. *)
(* *)
(* on return *)
(* *)
(* AP the upper triangular half of the inverse . *)
(* the strict lower triangle is unaltered. *)
(* det determinant of original matrix if requested. *)
(* otherwise not referenced. *)
(* determinant = det(1) * 10.0**det(2) *)
(* with 1.0 .le. det(1) .lt. 10.0 *)
(* or det(1) .eq. 0.0 . *)
(* *)
(* error condition *)
(* *)
(* a division by zero will occur if the input factor contains *)
(* a zero on the diagonal and the inverse is requested. *)
(* it will not occur if the subroutines are called correctly *)
(* and if dpoco or dpofa has set info .eq. 0 . *)
(* *)
(* LinPack. This version dated 08/14/78 . *)
(* Cleve Moler, university of new mexico, argonne national lab. *)
(* *)
(* Modula-2 version M.Riedl, 23.04.2016 *)
(*----------------------------------------------------------------*)
END LinPack.