!C     Chuong trinh xac dinh ham phan lop
!C     Day la chuong trinh chinh, can phai
!C     co 3 chuong trinh con la
!C     DMATX, DISCR, MINV
!C
      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)
!      dimension N(5),CMEAN(10),XBAR(50),C(55),D(100)
!      DIMENSION P(250),LG(250),X(2500)
      character *120 fsl,fra
!C   Kich thuoc duoc mo ta nhu sau
!C   Cua N : >= so nhom K
!C   Cua CMEAN : >= so bien M
!C   Cua XBAR  : >= tich so M x K
!C   Cua C : >= tich (M+1) x K
!C   Cua D : >= tich M x M
!C   cua P,LG : >= kich thuoc chung cua tap mau K nhom, T
!C                (T=N(1)+N(2)+...+N(k))
!C   Cua X : >= tich T x M
!C==================================================================
!C...CAU TRUC FILE SO LIEU NHU SAU:
!C...DONG 1: GOM K+1 SO NGUYEN DUONG, XEP THEO THU TU LA
!C           K (SO NHOM), M (SO BIEN), N1,.., NK (DUNG LUONG 
!C                                         MAU CUA TUNG NHOM)
!C...DONG 2 TRO DI GOM N1+...+NK DONG, MOI DONG GOM M SO LA
!C          SO LIEU QUAN TRAC CUA M BIEN. CAC DONG DUOC XEP THEO THU 
!C          TU TU TREN XUONG DUOI TUONG UNG TU LOP 1 (N1) DEN LOP K (NK)
!C
!      write(*,'(a\)')' Cho ten File SL : '
!      read(*,'(a)')fsl
!      write(*,'(a\)')' Cho ten File Ra : '
!      read(*,'(a)') fra
      fsl="Disc_Inp/Test.txt"
      fra="DISC_test.txt"
      open(1,File=fsl,status='old')
      open(3,File=fra,status='Unknown')
      do 400 i=1,25
  400 write(*,*)
  100 read(1,*) K,M,(N(i),i=1,K)
      write(3,2) K,M
    2 format(' Chuong trinh Ham phan lop '/'So nhom',I5/'So bien',I5)
!C      do 110 i=1,k
!C  110 write(3,'(2i6)') I,N(i)
      L=0
      do 130 i=1,k
      n1=n(i)
      do 120 j=1,n1
      read(1,*)(cmean(ij),ij=1,M)
      L=L+1
      n2=L-n1
      do 120 ij=1,M
      n2=n2+n1
  120 x(n2)=cmean(ij)
  130 L=n2
      call DMATX(k,m,n,x,xbar,d,cmean)
      L=0
      do 150 i=1,k
      do 140 j=1,m
      L=L+1
  140 cmean(j)=xbar(L)
!C  150 write(3,'('' Nhom'',i4,'' Tbinh ''/8f15.5)') I,(cmean(j),j=1,M)
!C      write(3,'('' Ma tran hiep phuong sai '')')
  150 continue
      do 170 i=1,m
      L=i-m
      do 160 j=1,m
      L=L+m
  160 cmean(j)=D(L)
!C  170 write(3,8)i,(CMEAN(j),j=1,M)
  170 continue
    8 format(' Hang ',i3/(8f15.5))
      Call MINV(D,M,DET,CMEAN,C)
      Call DISCR(K,M,N,X,XBAR,D,CMEAN,V,C,P,LG)
!C      write(3,9)(CMEAN(I),I=1,M)
    9 format(' Trung binh chung '/(8f15.5))
!C      write(3,10) V
   10 format(' Hang so D-Mahalanov ',f15.5)
      n1=1
      n2=M+1
      do 180 i=1,k
      write(3,11) I,(C(J),J=N1,N2)
   11 FORMAT(' Ham phan lop ',i3/6x,'Hang so    Cac he so'/9f12.3)
      n1=n1+(M+1)
      n2=n2+(M+1)
  180 continue
!C      write(3,12)
   12 format('Evoluation of classification function'
     \' for each observation')
      n1=1
      n2=N(1)
      do 210 i=1,k
      write(3,13) I
   13 format(' Nhom ',i3)
!C   13 format(' Nhom ',i3/14x,'  Xac suat ung voi ham     So thu tu'/
!C     \' Quan trac        phan lop lon nhat       ham lon nhat')
      L=0
      do 190 j=n1,n2
      L=L+1
  190 write(3,14) L,LG(j)
   14 format(1x,i4,2x,i6)
      if (i-k) 200,300,300
  200 n1=n1+N(i)
      n2=n2+N(i+1)
  210 continue
  300 continue
      end

!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C ---------    CAC CHUONG TRINH CON    -------------
!C
      SUBROUTINE DMATX(K,M,N,X,XBAR,D,CMEAN)
      dimension N(1),X(1),XBAR(1),D(1),CMEAN(1)
      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
      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)
      LL=0
      do 150 j=1,M
      do 150 jj=1,M
      LL=LL+1
  150 D(LL)=D(LL)+CMEAN(J)*CMEAN(JJ)
  160 N4=N4+N1*M
      LL=-K
      DO 170 I=1,K
  170 LL=LL+N(I)
      FN=LL
      do 180 I=1,MM
  180 D(I)=D(I)/FN
      RETURN
      END
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C
      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)
      N1=N(1)
      do 100 i=2,k
  100 N1=N1+N(I)
      FN1=N1
      do 110 i=1,k
  110 P(i)=N(I)
      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)
  130 CMEAN(I)=CMEAN(I)/FN1
      L=0
      do 140 i=1,k
      do 140 j=1,M
      L=L+1
  140 C(L)=XBAR(L)-CMEAN(J)
      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)
      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)
      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)
      N2=0
      do 220 KA=1,K
      N2=N2+1
      Sum=C(N2)
      do 210 j=1,M
      N2=N2+1
  210 Sum=Sum+C(N2)*D(J)
  220 XBAR(KA)=Sum
      L=1
      Sum=XBAR(1)
      do 240 j=2,k
      IF(Sum-XBAR(j)) 230,240,240
  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
  260 P(N1)=1/PL
  270 LBASE=LBASE+NN*M
      RETURN
      END
!CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
!C
      SUBROUTINE MINV(A,N,D,L,M)
      DIMENSION A(1),L(1),M(1)
      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 of file
!--------------------------------------------------

