!     Chuong trinh xac dinh ham phan lop
!     Kem theo 3 subroutines: DMATX, DISCR, MINV
!
!      Program phanlop
!      PARAMETER (KMAX=5,MMAX=20,NMAX=200)
!      DIMENSION N(KMAX),CMEAN(MMAX),XBAR(KMAX*MMAX),C((MMAX+1)*KMAX),D(MMAX*MMAX)
!      DIMENSION P(NMAX*KMAX),LG(NMAX*KMAX),X(NMAX*KMAX*MMAX)
!      character *12 fsl,fra

!   Kich thuoc duoc mo ta nhu sau
!   Cua N : >= so nhom K
!   Cua CMEAN : >= so bien M
!   Cua XBAR  : >= tich so M x K
!   Cua C : >= tich (M+1) x K
!   Cua D : >= tich M x M
!   cua P,LG : >= kich thuoc chung cua tap mau K nhom, T
!                (T=N(1)+N(2)+...+N(k))
!   Cua X : >= tich T x M
!   Chu y: X la mang 1 chieu chua SL cua cac nhom. Thu tu cua cac phan tu cua X la:
!       +  N(1) x M phan tu dau tien chua SL nhom 1:
!          N(1) phan tu bien 1, N(1) phan tu bien 2,..., N(1) phan tu bien M,
!       +  N(2) x M phan tu tiep theo chua SL nhom 2:
!          N(2) phan tu bien 1, N(2) phan tu bien 2,..., N(2) phan tu bien M,
!       + ...
!       +  N(k) x M phan tu cuoi cung chua chua SL nhom k:
!          N(k) phan tu bien 1, N(k) phan tu bien 2,..., N(k) phan tu bien M
!=====================================================================================

  Subroutine My_Discr (X,M,K,N,Nspl,F,LG)
    Integer :: M,K,N(K) , LG(Nspl) !,Lgr(Nspl)
    Real    :: X(M*Nspl) , F(K,M+1)
    ! Input: X(M*Nspl):          So lieu dau vao chua tam mau cua cac bien cua tung nhom nhu mo ta tren day 
    !        Nspl=N(1)+...+N(K): Tong dung luong mau
    !        M:                  So bien; 
    !        K:                  So nhom; 
    !        N:                  Mang chua dung luong mau tung nhom
    ! Output: F(K,M+1) chua he so cac phuong trinh tung nhom 
    !       + F(k,1) hang so nhom k
    !       + F(k,i), i=2..M+1 he so cac bien 1..M
    !         LG(Nspl)=LG[1:N(1),N(1)+1:N(2),...,N(k-1)+1:N(k)]: Chua so thu tu nhom xuat hien trong tung mau
    ! Working variables:
    Integer :: i,j,L,n1,n2,i1,i2,j1,j2
    !Real :: CMEAN(M),XBAR(K*M),C((M+1)*K),D(M*M)
    !Real :: P(Nspl) 
    Real, Allocatable    :: CMEAN(:),XBAR(:),C(:),D(:),P(:)
    !Real, Allocatable    :: F(:,:)
    !Integer, Allocatable :: LG(:)
    !--------------------------------------------------------------------

    !if (Allocated (LG))    Deallocate (LG);    if (Allocated (F))    Deallocate (F)
    if (Allocated (CMEAN)) Deallocate (CMEAN); if (Allocated (XBAR)) Deallocate (XBAR)
    if (Allocated (C))     Deallocate (C);     if (Allocated (D))    Deallocate (D)
    if (Allocated (P))     Deallocate (P)

    !Allocate (LG(Nspl),F(K,M+1))
    Allocate (CMEAN(M),XBAR(K*M),C((M+1)*K),D(M*M),P(Nspl))

!print"(20F6.3)",X
      call DMATX(k,m,n,x,xbar,d,cmean)   ! Return Mean of variables for each group (XBAR) and Covariance matrix (D)

      !L=0
      !do  i=1,k
      !   do  j=1,m
      !      L=L+1
      !      cmean(j)=xbar(L)
      !   enddo
      !enddo
      !do  i=1,m
      !   L=i-m
      !   do  j=1,m
      !      L=L+m
      !      cmean(j)=D(L)
      !   enddo
      !enddo
      Call MINV(D,M,DET,CMEAN,C)        ! Return Invert of D-matrix and Determinant of D
      Call DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG) ! Return coeficiences of Discriminant functions (C) and Indics of groups (LG)
      n1=1
      n2=M+1
      do  i=1,k
         i1=0
         Do J=N1,N2
            i1=i1+1; F(i,i1)=C(J)
!print*,"F=",i,i1,F(i,i1)
         Enddo
         n1=n1+(M+1)
         n2=n2+(M+1)
      enddo
      !n1=1
      !n2=N(1)
      !do  i=1,k
         !write(3,'(A,I3)') ' Nhom ',I
      !   L=0
         !do  j=n1,n2
         !   L=L+1
         !   write(3,'(1x,i4,2x,i6)') L,LG(j)
         !enddo
      !   if (i<k) then
      !      n1=n1+N(i)
      !      n2=n2+N(i+1)
      !   else
      !      exit
      !   endif
      !enddo
  End Subroutine My_Discr

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C ---------    CAC CHUONG TRINH CON    -------------

!     call DMATX(k,m,n,x,xbar,d,cmean)
!     CMEAN(M),XBAR(K*M),D(M*M)

      SUBROUTINE DMATX(K,M,N,X,XBAR,D,CMEAN)
      dimension N(1),X(1),XBAR(1),D(1),CMEAN(1)
!     Input:  K,M,N,X
!     Output: XBAR,D

      MM=M*M
      do 100 i=1,MM
  100 D(i)=0
      N4=0
      L=0
      LM=0
      do 160 NG=1,K
         N1=N(NG)
         FN=N1
         do 130 j=1,M
            LM=LM+1
            XBAR(LM)=0
            do 120 i=1,N1
               L=L+1
  120       XBAR(LM)=XBAR(LM)+X(L)

  130    XBAR(LM)=XBAR(LM)/FN    ! Trung binh cac bien trong tung nhom
         LMEAN=LM-M
         do 150 i=1,N1
            LL=N4+i-N1
            do 140 j=1,M
               LL=LL+N1
               N2=LMEAN+j
  140       CMEAN(j)=X(LL)-XBAR(N2) ! Do lech (Anomaly) so voi trung binh cua nhom
            LL=0
            do 150 j=1,M
            do 150 jj=1,M
               LL=LL+1
  150       D(LL)=D(LL)+CMEAN(J)*CMEAN(JJ) ! Tinh Covariances = Transpose(X-Xtb nhom)*(X-Xtb nhom)
  160 N4=N4+N1*M

      LL=-K              ! Tru di so nhom
      DO 170 I=1,K
  170 LL=LL+N(I)
      FN=LL
      do 180 I=1,MM
  180 D(I)=D(I)/FN       ! Covariances (FN = Sum (N(i)|i=1,K) - K
      RETURN
      END
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

! In main Prog: Call DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)
!               CMEAN(M),XBAR(K*M),C((M+1)*K),D(M*M), P(Nspl), LG(Nspl)

      SUBROUTINE DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)
      dimension N(1),X(1),XBAR(1),D(1),CMEAN(1),C(1),P(1),LG(1)
!     Input:  K,M,N(K),X,XBAR,D (nghich dao cua Covariance - Output cua Minv)
!     Output: C(M+1), LG(Sum(N(i)|i=1,k)), P

      N1=N(1)
      do 100 i=2,k
  100 N1=N1+N(I)
      FN1=N1                             ! Tong dung luong mau
      do 110 i=1,k
  110 P(i)=N(I)                          ! Dung luong mau tung nhom
      do 130 i=1,M
      CMEAN(I)=0
      N1=I-M
      do 120 j=1,k
      N1=N1+M
  120 CMEAN(I)=CMEAN(I)+P(J)*XBAR(N1)    ! XBAR: Trung binh cac bien trong tung nhom
  130 CMEAN(I)=CMEAN(I)/FN1              ! Trung binh chung cua cac bien
      L=0
      do 140 i=1,k
      do 140 j=1,M
      L=L+1
  140 C(L)=XBAR(L)-CMEAN(J)              ! Do lech trung binh cac nhom so voi Trung binh chung
      V=0
      L=0
      do 160 j=1,M
      do 160 i=1,M
         N1=I-M
         N2=J-M
         Sum=0
         do 150 ij=1,k
            N1=N1+M
            N2=N2+M
  150    Sum=Sum+P(IJ)*C(N1)*C(N2)
         L=L+1
  160 V=V+D(L)*Sum
      N2=0
      do 190 KA=1,k
         do 170 i=1,M
            N2=N2+1
  170    P(I)=XBAR(N2)
         IQ=(M+1)*(KA-1)+1
         Sum=0
         do 180 j=1,M
            N1=J-M
            do 180 L=1,M
               N1=N1+M
  180    Sum=Sum+D(N1)*P(J)*P(L)
         C(IQ)=-(Sum/2)           ! Output He so C0 cua ham phan lop
         do 190 i=1,M
            N1=I-M
            IQ=IQ+1
            C(IQ)=0
            do 190 j=1,M
               N1=N1+M
  190 C(IQ)=C(IQ)+D(N1)*P(J)      ! Output Cac he so Cj cua ham phan lop
      LBASE=0
      N1=0
      do 270 KG=1,K
         NN=N(KG)
         do 260 i=1,NN
            L=I-NN+LBASE
            do 200 j=1,M
               L=L+NN
  200       D(J)=X(L)             ! SL ban dau
            N2=0
            do 220 KA=1,K
               N2=N2+1
               Sum=C(N2)          ! He so C0
               do 210 j=1,M
                  N2=N2+1
  210          Sum=Sum+C(N2)*D(J) ! Tinh Gia tri ham PL cua tung nhom (Tong cua tich cac he so voi cac gia tri cua bien dau vao cong voi hang so)
                                  ! = C0 + Sum (C(j)*X(j)|j=1,M)
  220       XBAR(KA)=Sum
            L=1
            Sum=XBAR(1)           ! Gia su L = Nhom 1
            do 240 j=2,k
               IF(Sum-XBAR(j)) 230,240,240   ! Neu XBAR(j)>Sum gan L = Nhom J
  230          L=J
               Sum=XBAR(J)
  240       CONTINUE
            PL=0
            do 250 j=1,k
  250       PL=PL+EXP(XBAR(J)-Sum)
            N1=N1+1
            LG(N1)=L              ! Output LG(:)=1,2,..,k
  260    P(N1)=1/PL               ! Output Xac suat ???
  270 LBASE=LBASE+NN*M
      RETURN
      END
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!     Call MINV(D,M,DET,CMEAN,C)

      SUBROUTINE MINV(A,N,D,L,M)
      DIMENSION A(1),L(1),M(1)
!     Input:  A, N
!     Output: A^-1, Det(A)

      D=1.
      NK=-N
      DO 80 K=1,N
      NK=NK+N
      L(K)=K
      M(K)=K
      KK=NK+K
      BIGA=A(KK)
      DO 20 J=K,N
      IZ=N*(J-1)
      DO 20 I=K,N
      IJ=IZ+I
   10 IF(ABS(BIGA)-ABS(A(IJ))) 15,20,20   !
   15 BIGA=A(IJ)
      L(K)=I
      M(K)=J
   20 CONTINUE
      J=L(K)
      IF(J-K)35,35,25                     !
   25 KI=K-N
      DO 30 I=1,N
      KI=KI+N
      HOLD=-A(KI)
      JI=KI-K+J
      A(KI)=A(JI)
   30 A(JI)=HOLD
   35 I=M(K)
      IF(I-K)45,45,38                     !
   38 JP=N*(I-1)
      DO 40 J=1,N
      JK=NK+J
      JI=JP+J
      HOLD=-A(JK)
      A(JK)=A(JI)
   40 A(JI)=HOLD
   45 IF(BIGA)48,46,48                    !
   46 D=0.
      RETURN
   48 DO 55 I=1,N
      IF(I-K)50,55,50                     !
   50 IK=NK+I
      A(IK)=A(IK)/(-BIGA)
   55 CONTINUE
      DO 65 I=1,N
      IK=NK+I
      HOLD=A(IK)
      IJ=I-N
      DO 65 J=1,N
      IJ=IJ+N
      IF(I-K)60,65,60                     !
   60 IF(J-K)62,65,62                     !
   62 KJ=IJ-I+K
      A(IJ)=HOLD*A(KJ)+A(IJ)
   65 CONTINUE
      KJ=K-N
      DO 75 J=1,N
      KJ=KJ+N
      IF(J-K)70,75,70                     !
   70 A(KJ)=A(KJ)/BIGA
   75 CONTINUE
      D=D*BIGA
      A(KK)=1./BIGA
   80 CONTINUE
      K=N
  100 K=(K-1)
      IF(K)150,150,105                     !
  105 I=L(K)
      IF(I-K)120,120,108                   !
  108 JQ=N*(K-1)
      JR=N*(I-1)
      DO 110 J=1,N
      JK=JQ+J
      HOLD=A(JK)
      JI=JR+J
      A(JK)=-A(JI)
  110 A(JI)=HOLD
  120 J=M(K)
      IF(J-K)100,100,125                   !
  125 KI=K-N
      DO 130 I=1,N
      KI=KI+N
      HOLD=A(KI)
      JI=KI-K+J
      A(KI)=-A(JI)
  130 A(JI)=HOLD
      GOTO 100
  150 RETURN
      END
!END