小弟写了一个求对称正定矩阵的cholesky分解的程序，参考的是徐士良的程序，但是总是出错，不知道是为什么，有时在10*10大小矩阵时能运行，但是对于20*20矩阵就不能用了，求教高手解答。
这是我的程序代码： PROGRAM MAIN
    IMPLICIT NONE
C  creat an original ordering matrix
      INTEGER N,K
      PARAMETER (K=20,n=20)
      INTEGER B,C,I,J,T,W,x
    INTEGER L(K,N),PP(K,N)
    INTEGER P(N)                       定义了好多，有一些都没用到
    DOUBLE PRECISION D(N),A(N)  
    DIMENSION  JS(K),IS(K)
    DOUBLE PRECISION JS,IS
      DOUBLE PRECISION PO(K,K),PPO(K,K)
    DOUBLE PRECISION FI(K,N)
    DOUBLE PRECISION G(N),H(N)
    DOUBLE PRECISION RELATION,DET
    REAL R


    B=1
    C=N
    R=1.0
    DO I=1,K
       CALL NRAB1(R,B,C,N,P)      用随机数程序生成一个K*N的矩阵
       DO J=1,N
         L(I,J)=P(J)
       END DO 
       R=R+2
    END DO 
c     DO I=1,N
c     WRITE(*,30) (L(I,J),J=1,N)
c30     FORMAT(1X,10I6) 
c      WRITE(*,*)
c    END DO
      DO I=1,K
      DO J=1,N
      A(J)=L(I,J)
        END DO
      DO T=1,K
         DO J=1,N
            D(J)=L(T,J)
         END DO
        PO(I,T)=ABS(RELATION(A,D,N))      求矩阵的相关性系数矩阵
      END DO
    END DO

      WRITE(*,30) ((PO(I,J),J=1,K),I=1,K)
30     FORMAT(1X,20F10.6) 
      WRITE(*,*)

      CALL BCHOL(PO,K,DET,x)            求矩阵的cholesky分解，调试过程中该子程序运行到最后几步时就     DO I=1,K                                               出错
    DO J=1,K
       PPO(I,J)=PO(I,J)
    END DO
    END DO
c    WRITE(*,20)((ppo(I,J),J=1,K),I=1,K)
c    WRITE(*,*)
    
c    IF(x.LE.0) THEN
c    WRITE(*,*)
cc    WRITE(*,10)
c    WRITE(*,20)((ppo(I,J),J=1,4),I=1,4)
c    WRITE(*,*)
c    END IF
c10    FORMAT(1X,'MAT L IS:')
c20    FORMAT(1X,10D15.6)
c30    FORMAT(1X,'DET=',D15.6)






    STOP
    END
***********************************************************
      
*************************************************************
      SUBROUTINE NRAB1(R,A,B,N,P)   ! CREAT A PERMUTATION
    INTEGER A,B
    INTEGER P(N)

    S=B-A+1.0
    K=LOG(S-0.5)/LOG(2.0)+1
    L=1
    DO 10 I=1,K
10    L=2*L
      K=1
    S=4.0*L
    I=1
20    IF(I<=L.AND.K<=N) THEN
      R=R+R+R+R+R
    M=R/S
    R=R-M*S
    J=A+R/4.0
    IF(J.LE.B) THEN
    P(K)=J
    K=K+1
    END IF
    I=I+1
    GOTO 20
    END IF 
    RETURN
    END           
******************************************************
      DOUBLE PRECISION FUNCTION RELATION(A,D,N)
    IMPLICIT NONE 
    INTEGER  N
    DOUBLE PRECISION  A(N),D(N)

    INTEGER I,J
    DOUBLE PRECISION SFENZI,SFENMU1,SFENMU2,S
    DOUBLE PRECISION AMEAN,DMEAN
C
      S=0
    DO I=1,N
     S=S+A(I)
    END DO
    AMEAN=S/N
    S=0
    DO I=1,N
     S=S+D(I)
    END DO 
    DMEAN=S/N
C
      SFENZI=0
    SFENMU1=0
    SFENMU2=0
    DO I=1,N
       SFENZI=SFENZI+(A(I)-AMEAN)*(D(I)-DMEAN)
       SFENMU1=SFENMU1+(A(I)-AMEAN)**2
       SFENMU2=SFENMU2+(D(I)-DMEAN)**2
    END DO 
    RELATION=SFENZI/SQRT(SFENMU1*SFENMU2)
    END 
*******************************************************
      SUBROUTINE BCHOL(A,N,DET,L)    
    DIMENSION A(N,N)
    DOUBLE PRECISION A,DET
    INTEGER N
    L=1
    IF(A(1,1)+1.0.EQ.1.0)THEN
    L=0
    WRITE(*,10)
    RETURN
    END IF
10    FORMAT(1X,'FIAL')
      A(1,1)=SQRT(A(1,1))
    DET=A(1,1)
    DO 20 I=2,N
20    A(I,1)=A(I,1)/A(1,1)
      DO 60 J=2,N
      DO 30 K=1,J-1
30       A(J,J)=A(J,J)-A(J,K)*A(J,K)
      IF(A(J,J)+1.0.EQ.1.0)THEN
    L=0
    WRITE(*,10)
    RETURN
    END IF
    A(J,J)=SQRT(A(J,J))
    DET=DET*A(J,J)
    DO 50 I=J+1,N
       DO 40 K=1,J-1
40    A(I,J)=A(I,J)-A(I,K)*A(J,K)
    A(I,J)=A(I,J)/A(J,J)
50    CONTINUE
60    CONTINUE
      DET=DET*DET
       DO 70 I=1,N-1
    DO 70 J=I+1,N
70    A(I,J)=0.0
      RETURN
    END