C     Chuong trinh xac dinh ham phan lop
C     Day la chuong trinh chinh, can phai
C     co 3 chuong trinh con tham du la
C     DMATX, DISCR, MINV
C
$debug
      Program phanlop
	PARAMETER (NLOP=10)
      dimension N(10),CMEAN(10),XBAR(100),C(110),D(100)
      DIMENSION P(500),X(5000)
	integer LG(500)
      character *12 fsl,fra
C	character *80 ten1, ten2
	DIMENSION ND(NLOP),ND1(NLOP),ND2(NLOP),ND3(NLOP)
	REAL PP(NLOP,NLOP), P1(NLOP),P2(NLOP)   ! LUU KET QUA DANG GIA
	NAMELIST  /FN_IN_OUT/ FSL,FRA 
C
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	

      CALL TACHSL    !...TACH SO LIEU THANH CAC LOP

	OPEN (21,FILE='INPUT_DISR.TXT',STATUS='OLD')
	READ (21,FN_IN_OUT)
	CLOSE (21)

!      write(*,'(a\)')' Cho ten File SL : '
!      read(*,'(a)')fsl
!      write(*,'(a\)')' Cho ten File Ra : '
!      read(*,'(a)') fra
      open(1,File=fsl,status='unknown')
      open(3,File=fra,status='unknown')
C      read(1,*) 
C	read(1,*) 
      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     \	   'Dung luong mau nhom 1:',I3/'Dung luong mau nhom 2:',I3/
C     \	   'Dung luong mau nhom 3:',I3/)
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)
C	write(*,'(10f4.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(6,'('' Nhom'',i4,'' Tbinh ''/8f15.5)') I,(cmean(j),j=1,M)
C      write(6,'('' 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(6,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(6,9)(CMEAN(I),I=1,M)
    9 format(' Trung binh chung '/(8f15.5))
C      write(6,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)
  180 n2=n2+(M+1)
C      write(3,12)
   12 format('Evoluation of classification function'
     \' for each observation')
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
	ND=0
	ND1=0
	ND2=0
	ND3=0
	PP = 0.0 ! GAN BANG 0 CHO MANG KET QUA DANH GIA
	P1 = 0.0
	P2 = 0.0
	TONG = 0.0
	 
      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')
		
!##############################################
		DO II=1,K
             DO J=N1,N2
	          IF (LG(J).EQ.II) then 
                   PP(I,II)=PP(I,II)+1
	          ENDIF
	       ENDDO
	    ENDDO
! ######################################
		L=0
		do 190 j=n1,n2
			L=L+1
		IF (LG(J).EQ.1) then 
		ND1(I)=ND1(I)+1

	    else if (LG(J).EQ.2) then 
		ND2(I)=ND2(I)+1 
	    else 
		ND3(I)=ND3(I)+1
	    endif
  190	    CONTINUE !! 	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
!###########################################
          DCX=0.0
	    DO J=1,K    ! CONG HANG DUOI CUNG (CONG THEO COT)
	       DO II=1,K
	          P1(J)=P1(J)+PP(II,J)
	          TONG=TONG+PP(II,J)
	       ENDDO
	    ENDDO
	    DO II=1,K    ! CONG COT BEN PHAI NHAT (CONG THEO HANG)
	       DO J=1,K
	          P2(II)=P2(II)+PP(II,J)
	       ENDDO
	    ENDDO
	    DO II=1,K
	       P1(II)=P1(II)/TONG*100
	       P2(II)=P2(II)/TONG*100
	       DO J=1,K
	          PP(II,J)=PP(II,J)/TONG*100
	          IF (II.EQ.J) DCX=DCX+PP(II,J)
	       ENDDO
	    ENDDO
!##############################################
      DO I=1,K
!	  WRITE(3,*)'Do chinh xac LOP',I,'='
!	  WRITE(3,*)'ND1(',I,')=',ND1(I),'/',n(i)
!        WRITE(3,*)'ND2(',I,')=',ND2(I),'/',n(i)
!	  WRITE(3,*)'ND3(',I,')=',ND3(I),'/',n(i)
	  WRITE (3,'(10F10.1)') (PP(I,II),II=1,K), P2(I)
	ENDDO
	WRITE (3,'(10F10.1)') (P1(I),I=1,K), TONG
	WRITE (3,'("    DO CHINH XAC CHUNG = ",F10.2)') DCX
      end


C      Chuong trinh con trong bo chuong trinh
C      phan tich phan biet (phan lop)
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


      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

C     Chuong trinh trong bo chuong trinh phan lop

      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
!####################################################
      SUBROUTINE TACHSL
	PARAMETER (NMAX=200, MMAX=12, KMAX=3,MVAR=8, N=40)
	REAL X(NMAX,0:MMAX), XX(NMAX*MMAX,MVAR), YY(NMAX*MMAX)
	REAL XG(NMAX*MMAX,MVAR,KMAX) ! MANG CHUA SO LIEU CAC NHOM
	                             ! (KMAX NHOM)
!	INTEGER YEAR1,YEAR2,YEAR,NAM(NMAX)
	INTEGER IDX_YT(12), IDK_YT(12)       ! 12 MONTH IN YEAR
	INTEGER IDX_NT(12), IDK_NT(12)       ! 12 MONTH IN YEAR
	INTEGER SHIFT
	CHARACTER *20 FNAME_NT(MVAR)

	DATA FNAME_NT/'PREDICTOR\NINO12.TXT','PREDICTOR\NINO3.TXT',
     1              'PREDICTOR\NINO34.TXT','PREDICTOR\NINO4.TXT',
     2              'PREDICTOR\SOI.TXT',   'PREDICTOR\SST1.TXT',
     3              'PREDICTOR\SST2.TXT',  'PREDICTOR\SST3.TXT' /
	DATA SHIFT /1/	! =1 : DICH CHUYEN MOT NAM
	                ! =0 : KHONG DICH CHUYEN
      DATA IDX_YT   / 6, 7, 8, 0, 0, 0, 0, 0, 0, 0, 0, 0/   ! THANG DUOC CHON DB
	DATA IDK_YT   / 0, 0, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1/

	DATA IDX_NT   / 1, 2, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0 /  ! THANG DUOC CHON
	                                                      ! LAM NHAN TO DB
	DATA IDK_NT   / 0, 0, 0,-1,-1,-1,-1,-1,-1,-1,-1,-1 /

                ! IDK(0) = 1 : THERE ARE SEVERAL OF MONTHS MUST BE SHIFT 1 YEAR
	          !        = 0 : NO SHIFT
	          ! 
	          ! IDK(I) = 0 : THIS YEAR (NO SHIFT)
	          ! (I<>0) = 1 : NEXT YEAR (SHIFT)
	          !        =-1 : IGNOR

C...READ PREDICTANT-DATA 

	OPEN (1,FILE='PREDICTANT\CODE_B1T.DAT',STATUS='OLD')
	M=MMAX

	DO I=1,15     ! BO 15 DONG
	   READ(1,*)
	ENDDO

	DO I=1,N
	   READ(1,*) (X(I,J),J=0,M)
	ENDDO
	CLOSE(1)

	IF (SHIFT.EQ.1) THEN ! (1)
      NN=0
      DO K=1,12
         IF (IDK_YT(K).EQ.0) THEN
		  DO I=1,N-1
	         NN=NN+1
	         YY(NN)=X(I,IDX_YT(K))
	      ENDDO
	   ELSE IF (IDK_YT(K).EQ.1) THEN
	      DO I=2,N
	         NN=NN+1
	         YY(NN)=X(I,IDX_YT(K))
            ENDDO
	   ENDIF	      
	ENDDO
	NN_YT=NN
	ELSE   ! (1)
      NN=0
      DO K=1,12
         IF (IDK_YT(K).EQ.0) THEN
	      DO I=1,N
	         NN=NN+1
	         YY(NN)=X(I,IDX_YT(K))
	      ENDDO
	   ELSE IF (IDK_YT(K).EQ.1) THEN
	      DO I=1,N
	         NN=NN+1
	         YY(NN)=X(I,IDX_YT(K))
            ENDDO
	   ENDIF	      
	ENDDO
	NN_YT=NN
	END IF ! (1)
C
C...READ PREDICTOR-DATA 
C
      DO 50 L=1,MVAR     

	OPEN (1,FILE=FNAME_NT(L),STATUS='OLD')
	M=MMAX
	READ(1,*)
	READ(1,*) !N, YEAR1, YEAR2
	READ(1,*)
	DO I=1,N
	   READ(1,*) (X(I,J),J=0,M)
	ENDDO

	IF (SHIFT.EQ.1) THEN ! (2)
      NN=0
      DO K=1,12
         IF (IDK_NT(K).EQ.0) THEN
	      DO I=1,N-1
	         NN=NN+1
	         XX(NN,L)=X(I,IDX_NT(K))
	      ENDDO
	   ELSE IF (IDK_NT(K).EQ.1) THEN
	      DO I=2,N
	         NN=NN+1
	         XX(NN,L)=X(I,IDX_NT(K))
            ENDDO
	   ENDIF	      
	ENDDO
	NN_NT=NN
	ELSE   ! (2)
      NN=0
      DO K=1,12
         IF (IDK_NT(K).EQ.0) THEN
	      DO I=1,N
	         NN=NN+1
	         XX(NN,L)=X(I,IDX_NT(K))
	      ENDDO
	   ELSE IF (IDK_NT(K).EQ.1) THEN
	      DO I=1,N
	         NN=NN+1
	         XX(NN,L)=X(I,IDX_NT(K))
            ENDDO
	   ENDIF	      
	ENDDO
	NN_NT=NN
	END IF ! (2)
	CLOSE(1)
  50  CONTINUE
C
C   TACH NHOM
C
      N1=0
	N2=0
	N3=0
	DO I=1,NN_NT
	   IF (YY(I).EQ.1.0) THEN
	      N1=N1+1
	      DO L=1,MVAR
	         XG(N1,L,1) = XX(I,L)
		  ENDDO
	   ELSE IF (YY(I).EQ.2.0) THEN
	      N2=N2+1
	      DO L=1,MVAR
      	     XG(N2,L,2) = XX(I,L)
	      ENDDO
	   ELSE IF (YY(I).EQ.3.0) THEN
	      N3=N3+1
	      DO L=1,MVAR
	         XG(N3,L,3) = XX(I,L)
	      ENDDO
         ENDIF
  	ENDDO
C
	OPEN (3,FILE='TESTOUT.TXT', STATUS='UNKNOWN')
      WRITE(3,*) KMAX,MVAR,N1,N2,N3
C... LOP 1
!!	DO K=1,KMAX
	   DO I=1,N1
	      WRITE(3,'(8F7.2)') (XG(I,L,1),L=1,MVAR)
	   ENDDO
!!	ENDDO

C... LOP 2
!!	DO K=1,KMAX
	   DO I=1,N2
	      WRITE(3,'(8F7.2)') (XG(I,L,2),L=1,MVAR)
	   ENDDO
!!	ENDDO

C... LOP 3
!!	DO K=1,KMAX
	   DO I=1,N3
	      WRITE(3,'(8F7.2)') (XG(I,L,3),L=1,MVAR)
	   ENDDO
!!	ENDDO
	CLOSE (3)
	END ! TACH SO LIEU