!  Subroutine My_Discr
  !------------------------------------------------------------------------------
INTEGER :: NGROUP, NVAR, NMAX, N  
          ! So nhom, So bien, Dung luong MAX moi nhom
INTEGER, ALLOCATABLE :: NG(:)    ! Dung luong mau thuc te

REAL, ALLOCATABLE :: X(:,:,:), XX(:,:)     ! So lieu (OBS, VAR, GROUP)
REAL, ALLOCATABLE :: MEAN(:,:), CMEAN(:) 
        ! Trung binh tung nhom va Trung binh chung
REAL, ALLOCATABLE :: CG(:,:,:), C(:,:), C1(:,:)
        ! Covariance Matrixes
REAL, ALLOCATABLE :: P(:), COEF(:,:), Tmp(:)
REAL, ALLOCATABLE :: F(:,:), X1(:,:)
INTEGER, ALLOCATABLE :: LOBS(:), LFSC(:)
! Doc SL
OPEN (1, FILE='my_discr.txt')
READ (1,*) NGROUP,NVAR

ALLOCATE( NG(NGROUP), MEAN(NGROUP, NVAR), CMEAN(NVAR) ) 
ALLOCATE( CG(NVAR, NVAR, NGROUP), C(NVAR, NVAR),C1(NVAR,NVAR))
ALLOCATE( P(NGROUP), COEF(NGROUP,0:NVAR), Tmp(NVAR) )

READ (1,*) (NG(K),K=1,NGROUP)
   N=0         ! Tong so mau = Sum (NG(i)|i=1,NGROUP)
   Do i=1,NGROUP
      N=N+NG(i)
   Enddo
!NMAX  = NG(1)
!Do K=2,NGROUP    ! Dung luong mau CUC DAI cua cac nhom
!   IF (NG(K) > NMAX) NMAX = NG(K)
!Enddo
!ALLOCATE( X(NMAX, NVAR, NGROUP) )
   Allocate (XX(N,NVAR),X1(N, NVAR),LOBS(N),F(N, NGROUP),LFSC(N))
   i=0
   Do K=1,NGROUP
      Do j=1,NG(K)
         !READ (1,*) (X(I,J,K),J=1,NVAR)
         i=i+1; LOBS(i)=K                  ! Nhan cua nhom cua cac mau
         Read (1,*) (XX(i,j1), j1=1,Nvar)
      Enddo
   Enddo
!N = 0   ! Tong dung luong mau
!Do i=1,Ngroup
!   N = N + NG(i)
!Enddo
!ALLOCATE( F(N, NGROUP), X1(N, NVAR), LOBS(N), LFSC(N) )
! Xac dinh cac nhom ban dau
!L = 0
!Do K=1,NGROUP
!   Do I=1,NG(K)
!          L = L + 1
!          LOBS(L) = K
!   Enddo
!Enddo
   Do J = 1, NVAR
      CMEAN(J) = 0.       ! Trung binh chung
      n1=1; n2=NG(1)
      Do k = 1, Ngroup
          MEAN(K,J) = 0.  ! Trung binh tung nhom
          Do i=n1,n2 !1, NG(K)
             MEAN(K,J) = MEAN(K,J) + XX(i,j)  !X(I,J,K)
          Enddo
          CMEAN(J) = CMEAN(J) + MEAN(K,J)
          MEAN(K,J) = MEAN(K,J)/REAL(NG(K))
          if (k<Ngroup) then; n1=n2+1; n2=n2+NG(k+1); else; exit; endif
      Enddo
print*,MEAN(:,J)
      CMEAN(J) = CMEAN(J) /REAL(N)
   Enddo
! Tinh do lech
   Do J=1, NVAR ! K=1,Ngroup
      n1=1; n2=NG(1)
      Do K=1,Ngroup ! J=1, NVAR
         Do I=n1,n2 !1,NG(K)
            X1(I,J) = XX(I,J) - MEAN(K,J) ! CMEAN(J)
print*,XX(I,J),MEAN(K,J),X1(I,J)
         Enddo
         if (k<Ngroup) then; n1=n2+1; n2=n2+NG(k+1); else; exit; endif
      Enddo
   Enddo
! Tinh ma tran Covariance:
   Do J1=1,NVAR
      Do J2=1,NVAR
         C(J1,J2) = 0.
         Do i=1,N
            C(J1,J2) = C(J1,J2) + X1(I,J1)*X1(I,J2)
         Enddo
         ! Do K=1,Ngroup
         !        CG(J1,J2,K) = 0.
         !        Do I=1,NG(K)
         !           CG(J1,J2,K) = CG(J1,J2,K) + & 
         !               X(I,J1,K)*X(I,J2,K)
         !        Enddo
         !        CG(J1,J2,K) = CG(J1,J2,K)/REAL(NG(K))
         !        C(J1,J2) = C(J1,J2) + CG(J1,J2,K)*NG(K)
         ! Enddo
          C(J1,J2) = C(J1,J2)/REAL(N-Ngroup)
      Enddo
   Enddo

print*,C

! Tinh xac suat tien nghiem
   Do K=1, Ngroup
      P(K) = REAL(NG(K))/REAL(N)
   Enddo
! Tinh ma tran nghich dao cua C
   C1(:,:) = C(:,:)
   CALL MINVERT(C1,NVAR)
! Tinh cac he so cua ham phan biet
   Do K=1, NGroup
      Do J1=1,NVAR
         COEF(K,J1) = 0.
         Do J2=1,NVAR
             COEF(K,J1) = COEF(K,J1) + MEAN(K,J2)*C1(J2,J1)
         Enddo
      Enddo
      COEF(K,0) = 0.
      DO J1=1,NVAR
          Tmp(J1) = 0.
          Do J2=1,NVAR
              Tmp(J1) = Tmp(J1) + MEAN(K,J2)*C1(J2,J1)
          Enddo
      Enddo
      Do J2=1,NVAR
         COEF(K,0) = COEF(K,0) + Tmp(J2)*MEAN(K,J2)
      Enddo
      !COEF(K,0) = -0.5*COEF(K,0) + LOG(P(K))
      COEF(K,0) = -0.5*COEF(K,0) !+ LOG(P(K))
      print*,'Cac he so ham phan biet ', K
      print '(8F9.4)',COEF(K,0:)
   Enddo
! Tinh ham phan biet cho tap mau phu thuoc
!   Do J=1,NVAR
!      N1 = 0
!   Do K=1,Ngroup
!          Do I=1,NG(K)
!             N1 = N1+1
!             X1(N1,J) = X(I,J,K) ! + CMEAN(J)
!          Enddo
!   Enddo
!Enddo
    Do I=1,N ! 1
       Do K=1,NGROUP
          F(I,K) = 0.
          Do J=1,NVAR
             F(I,K) = F(I,K) + COEF(K,J)*XX(I,J)
          Enddo
          F(I,K) = F(I,K) + COEF(K,0)
       Enddo
       LFSC(I) = 1
       Temp = F(I,1)
       Do K=2,NGROUP
          IF (F(I,K)>Temp) THEN
             Temp = F(I,K)
             LFSC(I) = K
          EndIf
       Enddo
    Enddo
    Do I=1,N ! 1
       print '(4F10.4,2I4)', (F(I,K),K=1,NGROUP),LOBS(I), LFSC(I)
    Enddo
END  ! PROGRAM
!-----------------
SUBROUTINE MINVERT(A,N)
! C/TRINH TINH MA TRAN NGHICH DAO CUA MA TRAN A
! INPUT: + A(N,N) MANG So lieu 
!        + N SO HANG/COT CUA MA TRAN A
! OUTPUT: A MA TRAN NGHICH DAO CUA A
!
REAL  A(N,N)
M=N-1
DO I=1,N
   R=A(I,1)
   DO K=1,M
      A(I,K)=A(I,K+1)/R
   ENDDO
   A(I,N)=1.0/R
   DO J=1,N
      IF (I/=J) THEN
         R=A(J,1)
         DO K=1,M
            A(J,K)=A(J,K+1)-R*A(I,K)
         ENDDO
         A(J,N)=-R*A(I,N)
      ENDIF
   ENDDO
ENDDO 
RETURN
END SUBROUTINE