$large
$nofloatcalls
C This version plunges ahead with decomposition and inversion despite
C finding small diagonal elements, giving out error messages as it goes.
      SUBROUTINE SYMINV(X,N,NOTPSD)
C  LOWER TRIANGLE OF X IS STORED BY ROWS, N ROWS.
C  INVERSION IS DONE IN PLACE.
      LOGICAL NOTPSD
      DOUBLE PRECISION X(1)
      CALL CHOLES(X,N,NOTPSD)
      CALL TRIINV(X,N)
      CALL TRISQ(X,N)
      RETURN
      END
      
      INTEGER FUNCTION INDX(I,J)
      IF(I.GT.J)THEN
          IU=I
          IL=J
      ELSE
          IU=J
          IL=I
      ENDIF
      INDX=IU*(IU-1)/2+IL
      RETURN
      END
      
      SUBROUTINE CHOLES(X,N,NOTPSD) 
C---------------------------------------------------------------------
C  THIS FINDS TRIANGULAR XT S.T. XT TIMES XT-TRANSPOSE  IS   X. 
C      If X has near-zeroes for eigenvalues, the routine assumes these 
C      are exactly zero.  It leaves the diagonal values causing the 
C      difficulties unaltered, so diagonal elements less than PRECIS 
C      label the problem points, and it sets elements below this 
C      diagonal to zero.  It is set up to handle near-singularity
C      properly on the assumption that diagonal elements are about 1.0.
C--------------------------------------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z) 
      PARAMETER(PRECIS=1.D-15) 
      LOGICAL NOTPSD 
      DOUBLE PRECISION X(1),SUM 
      II=0 
      NOTPSD=.FALSE.
      DO 10 I=1,N 
          II=II+I 
          VALUE=X(II) 
          IF (VALUE.GT.PRECIS) THEN 
            VALUE=SQRT(VALUE) 
            X(II)=VALUE 
            JJ=II+I 
            DO 11 J=I+1,N 
                X(JJ)= X(JJ)/VALUE 
                JJK=JJ 
                JJI=II+I 
                DO 12 K=I+1,J 
                  JJK=JJK+1 
                  X(JJK)=X(JJK)-X(JJ)*X(JJI) 
                  JJI=JJI+K 
12              CONTINUE 
                JJ=JJ+J 
11          CONTINUE 
          ELSE 
            NOTPSD=.TRUE. 
            WRITE(*,*)'Accuracy problem at diagonal element ',I,'.' 
            JJ=II+I 
            DO 21 J=I+1,N 
                X(JJ)=0.D0 
                JJ=JJ+J 
21          CONTINUE 
          ENDIF 
10    CONTINUE 
      RETURN 
      END

      SUBROUTINE SYMVEC(X,Y,Z,N)
C  SYMMETRIC, PACKED X TIMES VECTOR Y, BOTH OF DIMENSION N. RESULT IN Z.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION X(1),Y(N),SUM,Z(N)
      DO 10 I=1,N
      SUM=0.
      IROW=I*(I-1)/2
      DO 12 J=1,I
      IROW=IROW+1
      SUM=X(IROW)*Y(J)+SUM
   12 CONTINUE
      J1=I+1
        DO 13 J=J1,N
          IROW=IROW+I+J-J1
          SUM=X(IROW)*Y(J)+SUM
   13   CONTINUE
   10 Z(I)=SUM
      RETURN
      END

      SUBROUTINE TRIINV(X,N)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION X(1),SUM,PRECIS
      PARAMETER(PRECIS=1.D-15) 
      II=INDX(N,N)
      III=II
C-----------------------------------------------
C      II indexes the end of the current row in the I loop
C      III indexes the element of the inverse which is currently 
C      being calculated.
C---------------------------------------------------
      DO 10 I=N,1,-1
C        HERE'S THE DIAGONAL ELEMENT CASE
      IF(X(II).GT.PRECIS) THEN
          X(II)=1./X(II)
      ELSE
          X(II)=0.0D0
      ENDIF
      III=III-1
      DO 11 J=I-1,1,-1
C     FORM CROSS PRODUCT
          SUM=0.D0
          JJI=II
          JJJ=II-I+J
          DO 12 K=I,J+1,-1
            SUM=SUM+X(JJI)*X(JJJ)
            JJI=JJI-1
            JJJ=JJJ-K+1
12        CONTINUE
          IF(X(JJJ).GT.PRECIS) THEN
            X(III)=-SUM/X(JJJ)
          ELSE
            X(III)=0.D0
          ENDIF
          III=III-1
11    CONTINUE
      II=II-I
10    CONTINUE
      RETURN
      END

      SUBROUTINE TRISQ(X,N)
C  THIS FINDS XTRANSPOSE TIMES X.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION X(1),SUM
      DO 10 I=1,N
      DO 11 J=I,N
          SUM=0.
          DO 12 K=J,N
            SUM=SUM+X(INDX(K,I))*X(INDX(K,J))
12           CONTINUE
          X(INDX(I,J))=SUM
11    CONTINUE
10    CONTINUE
      RETURN
      END

      FUNCTION DOT(A,B,N)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 A(*),B(*)
      DOT=0.
      DO 10 I=1,N
10       DOT=DOT+A(I)*B(I)
      RETURN
      END

      SUBROUTINE SCALAR(MAT,REALNO,N)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 MAT(*)
      DO 10 I=1,N
10       MAT(I)=REALNO*MAT(I)
      END

      SUBROUTINE MMNSCL(A,SCALAR,N)
      IMPLICIT REAL*8(A-H,O-Z)
      REAL*8 A(*),SCALAR
      DO 10 I=1,N
10       A(I)=A(I)-SCALAR
      END

      SUBROUTINE VECTRI(VEC,TRI,RESULT,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 VEC(*),TRI(*),RESULT(*)
      II=1
      DO 2 I=1,N
2        RESULT(I)=0.
      DO 11 IROW=1,N
         DO 10 ICOL=1,IROW
            RESULT(ICOL)=RESULT(ICOL)+VEC(IROW)*
     1            TRI(INDX(IROW,ICOL))
10       CONTINUE
11    CONTINUE
      RETURN
      END

      SUBROUTINE TRIVEC(TRI,VEC,RESULT,N)
      IMPLICIT REAL*8 (A-H,O-Z)
      REAL*8 VEC(*),TRI(*),RESULT(*)
      DO 10 IROW=1,N
      RESULT(IROW)=0.0
      DO 11 K=1,IROW
11          RESULT(IROW)=RESULT(IROW)+VEC(K)*
     1            TRI(INDX(IROW,K))
10    CONTINUE
      RETURN
      END

      SUBROUTINE MATSUB(A,B,C,N)
      REAL*8 A(*),B(*),C(*)
      DO 10 I=1,N
10       C(I) = A(I)-B(I)
      RETURN
      END

      SUBROUTINE MATMOV(A,B,N)
      REAL*8 A(*),B(*)
      DO 10 I=1,N
10       B(I)=A(I)
      END

      SUBROUTINE MATADD(A,B,C,N)
      REAL*8 A(*),B(*),C(*)
      DO 10 I=1,N
10       C(I) = A(I)+B(I)
      RETURN
      END

      SUBROUTINE SYMWRI(A,LUN,N)
      REAL*8 A(*)
      II=0
      DO 10 I=1,N
         WRITE(LUN,*) I,(A(II+J),J=1,I)
         II=II+I
10    CONTINUE
      END
C
      SUBROUTINE MATFIL(A,X,N)
      REAL*8 A(*),X
      DO 10 I=1,N
10       A(I)=X
      END
      SUBROUTINE MMULT(A,B,C,L,M,N)
      REAL*8 A(L,M), B(M,N), C(L,N)
      CALL MATFIL(C,0.0D0,L*N)
      DO 10 K=1,N
         DO 11 J=1,M
            DO 12 I=1,L
               C(I,K) = C(I,K) + A(I,J)*B(J,K)
12          CONTINUE
11       CONTINUE
10    CONTINUE
      END
C
C
      SUBROUTINE RECTRI(REC,TRI,RESULT,N,M)
      REAL*8 REC(N,M),TRI(*),RESULT(N,M), A
      A=0.0
      CALL MATFIL(RESULT,A,N*M)
      II=0
      DO 11 IRUN=1,M
         DO 10 ICOL=1,IRUN
            II=II+1
            DO 12 IROW=1,N
12             RESULT(IROW,ICOL)=RESULT(IROW,ICOL)+REC(IROW,IRUN)*
     1               TRI(II)
10       CONTINUE
11    CONTINUE
      RETURN
      END
      SUBROUTINE IMATFL(IMAT,ITEM,ISIZE)
      INTEGER IMAT(*)
      DO 10 I = 1, ISIZE
10       IMAT(I) = ITEM
      RETURN
      END
      SUBROUTINE XTX(X,XX,NROW,NCOL)
      REAL*8 X(NROW,NCOL),XX(*),DOT
      II=0
      DO 10 I=1,NCOL
         DO 11 J=1,I
            II=II+1
            XX(II) = DOT(X(1,I),X(1,J),NROW)
11       CONTINUE
10    CONTINUE
      RETURN
      END
      SUBROUTINE TRIREC(TRI,REC,RESULT,N,M)
      REAL*8 TRI(*),REC(*),RESULT(*)
      REAL*8 S
C-------------------------------------------------------------------
C     ACTUALLY, REC AND RESULT ARE NXM, BUT WE GET POSSIBLY MORE
C     EFFICIENT INDEXING BY TREATING THEM AS VECTORS AND RUNNING
C     THROUGH THEM WITH A CONSTANT INCREMENT.
C-------------------------------------------------------------------
      II=0
      S=0.0   
      CALL MATFIL(RESULT,S,N*M)
      DO 10 I=1,N
         DO 11 K=1,I
            II=II+1
            IJ = I
            DO 12 KJ = K,M*N,N
               RESULT(IJ) = RESULT(IJ)+TRI(II)*REC(KJ)
               IJ=IJ+N
12          CONTINUE
11       CONTINUE
10    CONTINUE
      END
      FUNCTION LOGDET(A,N)
      REAL*8 A(*),LOGDET
      INTEGER N 
      LOGICAL NOTPSD  
      LOGDET=0.0
      CALL CHOLES(A,N,NOTPSD) 
      IF(NOTPSD) THEN 
         WRITE(*,*) 'DET OF NONPSD MATRIX'
      ELSE  
         II=0 
         DO 10 I=1,N  
            II=II+I 
            LOGDET=LOGDET+LOG(A(II))
10       CONTINUE 
      ENDIF 
      RETURN
      END 
