MODULE lss_aprx;

USE mv_ari;

CONST RCONST = 1E-15;
      REINS   = 1.0;
      RNULL   = 0.0;


GLOBAL PROCEDURE MINV( VAR A : RMATRIX; VAR err : INTEGER );
{  MINV berechnet Naeherungsinverse der Matrix A }

VAR v     : DYNAMIC ARRAY [LB(A,1)..UB(A,1)] OF INTEGER;
    r     : RVECTOR[LB(A,1)..UB(A,1)];
    h     : REAL;
    ok    : BOOLEAN;
    i,j,k,l,
    au1,ao1,au2,ao2 : INTEGER;
BEGIN
   au1:= LB(A,1);  ao1:= UB(A,1);
   au2:= LB(A,2);  ao2:= UB(A,2);

   ok:= ao1-au1 = ao2-au2;   { quadratische Matrix ? }

   FOR i:= au1 TO ao1 DO v[i]:= i;

   i := au1-1; { Zeilenindex  }
   j := au2-1; { Spaltenindex }

   WHILE (j<ao2) AND ok DO
   BEGIN
      i:= succ(i);  j:= succ(j);  { i und j laufen synchron }
      IF j<ao2 THEN
      BEGIN
         l:= i;
         FOR k:= i+1 TO ao1 DO  IF ABS(a[k,j]) > ABS(a[l,j]) THEN l:= k;
         IF i<>l THEN
         BEGIN
            r:= a[l];  a[l]:= a[i];  a[i]:= r;
            k:= v[l];  v[l]:= v[i];  v[i]:= k;
         END;
      END;
      IF a[i,j]=RNULL THEN ok:= FALSE ELSE
      BEGIN
         h     := REINS / a[i,j];
         a[*,j]:= h*RVECTOR(a[*,j]);
         a[i,j]:= h;
         FOR k:= au2 TO ao2 DO
             IF k<>j THEN
             BEGIN
                FOR l:= au1 TO ao1 DO
                    IF l<>i THEN
                    BEGIN
                       h := a[l,k] - a[l,j]*a[i,k];
                       IF h=RNULL THEN a[l,k]:= RCONST*a[l,k]
                                ELSE a[l,k]:= h;
                    END;
                a[i,k]:= -a[i,k]*a[i,j];
             END;
      END;
   END;
   IF ok THEN
   FOR k:= au1 TO ao1 DO
   BEGIN
      l:= au2+k-au1;
      WHILE v[k]<>k DO
      BEGIN
         i:= v[k];            v[k]  := v[i];   v[i]  := i;
         j:= au2+i-au1;
         r:= RVECTOR(a[*,j]); a[*,j]:= a[*,l]; a[*,l]:= r;
      END;
   END;
   err:= ORD( NOT ok );
END {*** MINV ***};

GLOBAL PROCEDURE MINV1( VAR W : RMATRIX; VAR err : INTEGER );
{ MINV1 berechnet Naeherungsinverse fuer Matrix der Form I+W,
  Ausgabematrix ist ebenfalls von der Form I+W', wobei W E/A-Parameter }
VAR i,j,k,l,
    wu1,wo1,wu2,wo2 : INTEGER;
    h,wjj           : REAL;
BEGIN
   wu1:= LB(W,1);  wo1:= UB(W,1);
   wu2:= LB(W,2);  wo2:= UB(W,2);

   IF wo1-wu1 <> wo2-wu2 THEN err:= 1   { Keine quadratische Matrix !!! }
                         ELSE err:= 0;   

   i:= wu1-1;
   FOR j:= wu2 TO wo2 DO
   BEGIN
      i:= succ(i);
      IF (W[i,j]<>-1) AND (err=0) THEN
      BEGIN
         wjj   := W[i,j];
         h     := REINS /(1+wjj);
         W[*,j]:= h*RVECTOR(W[*,j]);
         W[i,j]:= -wjj*h;
         FOR k:= wu2 TO wo2 DO
         IF k<>j THEN
         BEGIN
            FOR l:= wu1 TO wo1 DO
            IF l<>i THEN  W[l,k]:= W[l,k] - W[l,j]*W[i,k];
            W[i,k]:= - W[i,k]*h;
         END;
      END ELSE err:= 1;
   END;
END;  {*** MINV1 ***}

END.
