      SUBROUTINE CHOLSKY (IDA, NMAT, M, N, A, NRHS, IDB, B)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
C   CHOLESKY DECOMPOSITION/SUBSTITUTION SUBROUTINE.
C
C   11/28/84  D H BAILEY  MODIFIED FOR NAS KERNEL TEST
C
      DIMENSION A(0:IDA, -M:0, 0:N), 
     $     B(0:NRHS, 0:IDB, 0:N), EPSS(0:256)
      DATA EPS/1E-13/
C
C  CHOLESKY DECOMPOSITION
C
      DO 1 J = 0, N
        I0 = MAX ( -M, -J )
C
C  OFF DIAGONAL ELEMENTS
C
        DO 2 I = I0, -1
          DO 3 JJ = I0 - I, -1
            DO 3 L = 0, NMAT
3             A(L,I,J) = A(L,I,J) - A(L,JJ,I+J) * A(L,I+JJ,J)
          DO 2 L = 0, NMAT
2           A(L,I,J) = A(L,I,J) * A(L,0,I+J)
C
C  STORE INVERSE OF DIAGONAL ELEMENTS
C
        DO 4 L = 0, NMAT
4         EPSS(L) = EPS * A(L,0,J)
        DO 5 JJ = I0, -1
          DO 5 L = 0, NMAT
5           A(L,0,J) = A(L,0,J) - A(L,JJ,J) ** 2
        DO 1 L = 0, NMAT
1         A(L,0,J) = 1. / SQRT ( ABS (EPSS(L) + A(L,0,J)) )
C
C  SOLUTION
C
      DO 6 I = 0, NRHS
        DO 7 K = 0, N
          DO 8 L = 0, NMAT
8           B(I,L,K) = B(I,L,K) * A(L,0,K)
          DO 7 JJ = 1, MIN (M, N-K)
            DO 7 L = 0, NMAT
7             B(I,L,K+JJ) = B(I,L,K+JJ) -  A(L,-JJ,K+JJ) * B(I,L,K)
C
        DO 6 K = N, 0, -1
          DO 9 L = 0, NMAT
9           B(I,L,K) = B(I,L,K) * A(L,0,K)
          DO 6 JJ = 1, MIN (M, K)
            DO 6 L = 0, NMAT
6             B(I,L,K-JJ) = B(I,L,K-JJ) - A(L,-JJ,K) * B(I,L,K)
C
      RETURN
      END
