$large
$nofloatcalls
      SUBROUTINE LU(A,RXCHG,M,N,NOTOK)
      LOGICAL NOTOK
      REAL*8 A(M,N),ACC,S
      PARAMETER(ACC=1.D-35)
      INTEGER RXCHG(M), P
C-------------------------------------------------------------- 
C     Overwrites A with L and U.  Assume M<=N.  Returns with M=Order
C     if NOTOK true.
C-------------------------------------------------------------- 
      NOTOK=.FALSE.
      DO 10 K = 1,M
         P=1
         S = 0.D0
         DO 11 J=K,N
            IF(ABS(A(K,J)).GT.ABS(S)) THEN
               P=J
               S = A(K,J)
            ENDIF
11       CONTINUE
         IF(ABS(S) .LT. ACC) THEN
            NOTOK=.TRUE.
            M = K-1
            RETURN
         ENDIF
         RXCHG(K) = P
         CALL SWAP(A(1,K),A(1,P),M)
         CALL SCALAR(A(K+1,K),1/A(K,K),M-K)
         DO 20 J = K+1, N
            DO 21 I = K+1,M
21             A(I,J) = A(I,J) - A(K,J)*A(I,K)
20       CONTINUE
10    CONTINUE
      RETURN
      END
C
      SUBROUTINE RLNDET(LU,RXCHG,N,LND,PLUMIN)
      REAL*8 LU(N,*),LND
      INTEGER PLUMIN,RXCHG(N)
C-------------------------------------------------------------- 
C     note that there is no check for nonsingularity.
C-------------------------------------------------------------- 
      LND = 0.0
      PLUMIN = 1
      DO 10 I=1,N
         PLUMIN = PLUMIN*INT(SIGN(1.1D0,LU(I,I)))
         LND = LND + LOG(ABS(LU(I,I)))
10    CONTINUE
      DO 11 I=1,N-1
11       IF(RXCHG(I).NE.I) PLUMIN = -PLUMIN
      RETURN
      END
C
      SUBROUTINE LINV(L,N)
      REAL*8 L(N,*),S
C-------------------------------------------------------------- 
C     We assume L has ones down the diagonal, so we leave these and 
C     the upper triangle alone.
C-------------------------------------------------------------- 
      DO 10 I = N,1,-1
         DO 11 J = I-1,1,-1
            S=L(I,J)
            DO 12 K = I-1,J+1,-1
12             S = S + L(I,K)*L(K,J)
            L(I,J) = -S        
11       CONTINUE
10    CONTINUE
      END
C
      SUBROUTINE UINV(U,N)
C-------------------------------------------------------------- 
C     We leave lower triangle alone.
C-------------------------------------------------------------- 
      REAL*8 U(N,*),S
      DO 10 I = 1, N
         U(I,I) = 1.D0/U(I,I)
         DO 11 J = 1, I-1
            S=U(J,J)*U(J,I)
            DO 12 K = J+1,I-1
12             S = S+ U(J,K)*U(K,I)
            U(J,I) = -S*U(I,I)
11       CONTINUE
10    CONTINUE
      RETURN
      END
C
      SUBROUTINE RECINV(A,LND,SGNLND,N,NOTOK)
      REAL*8 A(N,*),LND
      LOGICAL NOTOK
      INTEGER SGNLND,RXCHG(100)
      CALL LU(A, RXCHG, N, N, NOTOK)
      IF(.NOT.NOTOK) THEN
         CALL RLNDET(A,RXCHG,N,LND,SGNLND)
         CALL LINV(A,N)
         CALL UINV(A,N)
         CALL UXL(A,N)
         CALL RSHUFL(A,RXCHG,N)
      ENDIF
      RETURN
      END
C
      SUBROUTINE UXL(A,N)
      REAL*8 A(N,*), S
C----------------------------------------------------------------
C     Upper and lower triangular matrices are packed in same A.
C     We form U*L and store it in A, using no extra storage space.
C----------------------------------------------------------------
      DO 10 I=1,N
         DO 11 J = 1, I-1
            S=0.0D0
            DO 12 K = I, N
12             S = S+ A(I,K)*A(K,J)
            A(I,J) = S
11       CONTINUE
         DO 13 J=I,N
            S=A(I,J)
            DO 14 K = J+1, N
14             S = S + A(I,K)*A(K,J)
            A(I,J) = S
13       CONTINUE
10    CONTINUE
      RETURN
      END
C
      SUBROUTINE RSHUFL(A,RXCHG,N)
      REAL*8 A(N,*)
      INTEGER RXCHG(N)
      DO 10 ICOL = 1,N
10       CALL VSHUFL(A(1,ICOL),RXCHG,N)
      RETURN
      END
C
      SUBROUTINE VSHUFL(V,RXCHG,N)
      REAL*8 V(N), S
      INTEGER RXCHG(N)
      DO 10 I = N-1,1,-1
         S = V(I)
         V(I) = V(RXCHG(I))
         V(RXCHG(I)) = S
10    CONTINUE
      END
C
      SUBROUTINE SWAP(A,B,N)
      REAL*8 A(*), B(*), WX
      DO 10 I = 1, N
         WX = A(I)
         A(I) = B(I)
         B(I) = WX
10    CONTINUE
      RETURN
      END
C
      REAL*8 FUNCTION TRAATB(A,B,M,N)
      REAL*8 A(M,N),B(M,N),DOT,S
      S=DOT(A(1,1),B(1,1),M)
      DO 10 I=2, N
10       S=S + DOT(A(1,I),B(1,I),M)
      TRAATB=S
      RETURN
      END
C
      SUBROUTINE ATB(A,B,C,M,N,L)
C-------------------------------------------------------------- 
C     A is MxN, B is MxL, thus C is NxL
C-------------------------------------------------------------- 
      REAL*8 A(M,N), B(M,L), C(N,L), DOT
      DO 10 I=1,N
         DO 11 J=1,L
11          C(I,J) = DOT(A(1,I),B(1,J),M)
10    CONTINUE
      RETURN
      END
C
      SUBROUTINE ABT(A,B,C,M,N,L)
      REAL*8 A(M,N), B(L,N), C(M,L), S
      S=0.D0
      CALL MATFIL(C,S,M*L)
      DO 10 JACOL=1,N
         DO 11 JBCOL = 1,N
            DO 12 J = 1, L
               S=B(J,JBCOL)
               DO 13 I= 1, M
13                C(I,J) = C(I,J) + S*A(I,JACOL)
12          CONTINUE
11       CONTINUE
10    CONTINUE
      RETURN
      END
