! Chuong trinh tinh ham phan lop:
   ! Cho vao mang X(N,NVAR) chua N=NG(1)+NG(2)+... gia tri cua NVAR bien cua NGROUP nhom
   ! Chuong trinh tra ve cac he so COEF cua NGROUP ham phan lop va chi so cac lop LG cua cac quan trac 
   ! Chuong trinh da duoc test voi SL va so sanh voi bo chuong trinh cua IMSL

   Subroutine Discr_2023(NGROUP,N,NVAR,NG,X,COEF,LG)
      ! Input:  NGROUP: So nhom; NVAR: So bien; NG(NGROUP): Mang chua so mau tuong ung cua tung nhom
      !         NG(1): Nhom 1, NG(2): Nhom 2,..., NG(NGROUP): Nhom NGROUP
      !         N: Tong so mau cua tat ca cac nhom: N=NG(1)+NG(2)+...+NG(NGROUP)
      !         X(N,NVAR): SL quan trac cua cac nhom tuong ung:
      !         X(1:NG(1),NVAR): SL nhom 1; X(NG(1)+1:NG(1)+NG(2),NVAR): SL nhom 2,...
      ! Output: COEF(NGROUP,NVAR+1): Cac he so ham phan lop C0, C1, C2,...,Cnvar cua tung nhom:
      !         COEF(1,1:NVAR+1): Nhom 1, COEF(2,1:NVAR+1): Nhom 2,...
      !         Chu y: He so tu do C0 ung voi COEF(:,1)
      !
      ! Khi su dung cho du bao ham phan lop cua cac nhom k doi voi mot quan trac i nao do se duoc xac dinh boi:
      !         F(i,k) = COEF(k,1) + SUM [COEF(k,j)*X(i,j)|j=1,Nvar]
      ! Chi so nhom cua quan trac i duoc gan boi k0 voi F(i,k0)=MAX{F(i,k)|k=1,NGROUP}
      !------------------------------------------------------------------------------
      INTEGER :: NGROUP, NVAR, N          ! So nhom, So bien, Tong Dung luong
      INTEGER :: NG(NGROUP)    ! Dung luong mau thuc te
      REAL    :: X(N,NVAR)    ! So lieu ban dau
      ! Output:
      Real    :: COEF(NGROUP,NVAR+1)
      INTEGER :: LG(N) 

      REAL, ALLOCATABLE :: MEAN(:,:)       ! Trung binh tung nhom va Trung binh chung
      REAL, ALLOCATABLE :: C1(:,:) !, C1(:,:) ! Covariance Matrixes
      REAL, ALLOCATABLE :: Tmp(:), F(:,:) ! P(:), Tmp(:)
      REAL, ALLOCATABLE :: X1(:,:)
      Integer :: i,j,k,i1,j1,j2,n1,n2

      !print*
      !print*, NGROUP,NVAR,N, NG(:)
      !i=0
      !Do K=1,NGROUP
      !   Do j=1,NG(K)
      !      i=i+1; print*,i,(X(i,j1), j1=1,Nvar)
      !   Enddo
      !Enddo


      if (Allocated (MEAN)) Deallocate (MEAN); !if (Allocated (C)) Deallocate (C)
      if (Allocated (C1)) Deallocate (C1); if (Allocated (X1)) Deallocate (X1)
      if (Allocated (F))  Deallocate (F);  if (Allocated (Tmp)) Deallocate (Tmp)

      Allocate (MEAN(NGROUP,NVAR),C1(NVAR,NVAR),X1(N,NVAR)) ! C(NVAR,NVAR)
      Allocate (F(N,NGROUP),Tmp(NVAR))
   
      Do J = 1, NVAR
         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) + X(i,j) 
            Enddo
            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)
      Enddo
      ! Tinh do lech
      Do J=1, NVAR
         n1=1; n2=NG(1)
         Do K=1,Ngroup 
            Do I=n1,n2
               X1(I,J) = X(I,J) - MEAN(K,J) ! CMEAN(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
            C1(J1,J2) = 0.
            Do i=1,N
               C1(J1,J2) = C1(J1,J2) + X1(I,J1)*X1(I,J2)
            Enddo
            C1(J1,J2) = C1(J1,J2)/REAL(N-Ngroup)
         Enddo
      Enddo

      ! Tinh ma tran nghich dao cua C
      CALL MINVERT(C1,NVAR)

      ! Tinh cac he so cua ham phan biet
      Do K=1, NGroup
         Do J1=1,NVAR
            COEF(K,J1+1) = 0.
            Do J2=1,NVAR
               COEF(K,J1+1) = COEF(K,J1+1) + MEAN(K,J2)*C1(J2,J1)
            Enddo
         Enddo
         !
         DO J1=1,NVAR
            Tmp(J1) = 0.
            Do J2=1,NVAR
               Tmp(J1) = Tmp(J1) + MEAN(K,J2)*C1(J2,J1)
            Enddo
         Enddo
         COEF(K,1) = 0.
         Do J2=1,NVAR
            COEF(K,1) = COEF(K,1) + Tmp(J2)*MEAN(K,J2)
         Enddo
         COEF(K,1) = -0.5*COEF(K,1) !+ LOG(P(K))
         print*,'Cac he so ham phan biet ', K
         print '(8F9.4)',COEF(K,1:)
      Enddo
      ! Tinh ham phan biet (F) va chi so cac nhom (LG) cho tap mau phu thuoc
      Do I=1,N ! 1
         Do K=1,NGROUP
            F(I,K) = COEF(K,1)
            Do J=1,NVAR
               F(I,K) = F(I,K) + COEF(K,J+1)*X(I,J)
            Enddo
         Enddo
         LG(I) = 1
         Temp = F(I,1)
         Do K=2,NGROUP
            if (F(I,K)>Temp) THEN
               Temp = F(I,K)
               LG(I) = K
            endif
         Enddo
      Enddo
      !Do I=1,N ! 1
      !   print '(2F10.4,2I4)', (F(I,K),K=1,NGROUP), LG(I)
      !Enddo
    END Subroutine Discr_2023 
    !------------------------------------------------------------------------------------
    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
       !
       Integer :: N, M,i,j,k
       REAL :: A(N,N), R
       
       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 MINVERT
    !------------------------------------------------------------------------------------