c   Метод Годунова-Булгакова исследования устойчивости решения
c
      SUBROUTINE STABIL(N,A,STAB,QAPPA)
      IMPLICIT REAL*8(A-H,O-Z)
      implicit integer*4 (i-n) 
      parameter (N1=20, N2=20, N3=40)
c
      DIMENSION A(N,N)
c
      allocatable Q (:,:)
      allocatable E (:,:)
      allocatable S (:,:)
      allocatable T (:,:)
      allocatable R (:,:)

      allocate (Q(N,N))
      allocate (E(N,N))
      allocate (S(N,N))
      allocate (T(N,N))
      allocate (R(N,N))
c
      DEPS=1.D-7
      ERSTAB=1.d-10
      L=1
      ANORM=0d0
c
        DO I=1,N
        SUMMA=0d0
         DO J=1,N
         SUMMA=SUMMA+DABS(A(I,J))
         T(I,J)=0.d0
         E(I,J)=0.d0
         END DO
        IF(ANORM.LT.SUMMA) ANORM=SUMMA
        T(I,I)=1d0
        E(I,I)=1d0
        END DO
      A1=1d0/ANORM
      S(1:N,1:N)=A(1:N,1:N)*A1
   11 CONTINUE
      A2=1d0/L
      ANORM=0.d0
       DO I=1,N
       SUMMA=0d0
        DO J=1,N
      SMULT=0d0
       DO M=1,N
       SMULT=SMULT+S(I,M)*T(M,J)
       END DO
      SUMMA=SUMMA+DABS(SMULT)
        R(I,J)=A2*SMULT
        E(I,J)=E(I,J)+R(I,J)
        END DO
       IF(ANORM.LT.SUMMA) ANORM=SUMMA
       END DO
      T(1:N,1:N)=R(1:N,1:N)      
      IF(ANORM.LT.DEPS.OR.L.EQ.N1) GOTO 10
      L=L+1
      GOTO 11
   10   L=1
       DO I=1,N
       Q(I,1:N)=0d0
       R(I,1:N)=0d0
       Q(I,I)=1d0
       R(I,I)=1d0
       END DO
C
   20 CONTINUE
       ANORM=0d0
       A2=1d0/(L+1)
       DO I=1,N
       SUMMA=0d0
        DO J=1,N
      SM1=0d0
      SM2=0d0
         DO M=1,N
          SM1=SM1+Q(I,M)*S(M,J)
          SM2=SM2+S(M,I)*Q(M,J)
        END DO
      SUMMA=SUMMA+DABS(SM2)
      T(I,J)=(SM1+SM2)*A2
      R(I,J)=R(I,J)+T(I,J)
      END DO
       IF(ANORM.LT.SUMMA) ANORM=SUMMA
       END DO
c
      Q(1:N,1:N)=T(1:N,1:N)
c
      IF(ANORM.LT.DEPS.OR.L.EQ.N2) GOTO 25
      L=L+1
      GOTO 20
   25 N0=DLOG(STAB*DLOG(2d0*DSQRT(STAB)))/DLOG(2d0)
      N4=N3+N0
      L=1
      QAPP0=0d0
      Q(1:N,1:N)=E(1:N,1:N)
   21    CONTINUE
         DO I=1,N
           DO J=1,N
         SMULT=0d0
          DO M=1,N
            SMULT=SMULT+R(I,M)*Q(M,J)
          END DO
          T(I,J)=SMULT
          END DO
      END DO
      ANORM=0d0
       DO I=1,N
       SUMMA=0d0
        DO J=1,N
      SMULT=0d0
       DO M=1,N
       SMULT=SMULT+Q(M,I)*T(M,J)
       END DO
      R(I,J)=R(I,J)+SMULT
      SUMMA=SUMMA+DABS(R(I,J))
      END DO
       IF(ANORM.LT.SUMMA) ANORM=SUMMA
       END DO
      QAPPA=ANORM+ANORM
      DELTA=DABS(QAPPA-QAPP0)/QAPPA
      IF(QAPPA.GT.STAB.OR.L.EQ.N4.OR.DELTA.LT.ERSTAB)GOTO 35
      ANORM=0d0
      DO I=1,N
        SUMMA=0d0
        DO J=1,N
      SMULT=0d0
       DO M=1,N
       SMULT=SMULT+Q(I,M)*Q(M,J)
       END DO
      SUMMA=SUMMA+DABS(SMULT)
      E(I,J)=SMULT
      END DO
       IF(ANORM.LT.SUMMA) ANORM=SUMMA
       END DO
      Q(1:N,1:N)=E(1:N,1:N)
      L=L+1
      QAPP0=QAPPA
      GOTO 21
   35 CONTINUE
      deallocate(Q,E,S,T,R)      
c
c Ошибки
      CALL error('STABIL')
      RETURN
      END