!C   ARGUMENTS    X      - INPUT/OUTPUT MATRIX OF DIMENSION NV BY THE            ^M
!C                           SUM OF THE ND(I). ON INPUT, X CONTAINS THE          ^M
!C                           GROUP DATA.                                         ^M
!C                         ON OUTPUT, EACH COLUMN DATA VECTOR IS REPLACED        ^M
!C                           BY THE TRANSFORMED NNV BY 1 DATA VECTOR.            ^M
!C                           (SEE REMARK 1).                                     ^M
!C   REMARKS  1.  ON INPUT, THE MATRIX X CONTAINS THE ND(1)+ND(2)+...            ^M
!C                ...+ND(NG) SAMPLE (COLUMN) VECTORS. THE FIRST ND(1)            ^M
!C                COLUMNS OF MATRIX X SHOULD CONTAIN THE DATA FOR GROUP          ^M
!C                ONE. THE NEXT ND(2) COLUMNS OF MATRIX X SHOULD CONTAIN         ^M
!C                THE DATA FOR GROUP TWO, AND SO ON.                             ^M
!C                ON OUTPUT, EACH ORIGINAL OBSERVATION X(I) (WHERE X(I)          ^M
!C                REPRESENTS A COLUMN VECTOR OF INPUT MATRIX X) IS               ^M
!C                REPLACED BY ITS TRANSFORMED OBSERVATION ( A NNV BY 1           ^M
!C                DATA VECTOR). THE TRANSFORMED OBSERVATION IS GIVEN BY          ^M
!C                THE FORMULA (C TRANSPOSE) * X(I) WHERE C IS THE OUTPUT         ^M
!C                MATRIX C.                                                      ^M

      REAL, ALLOCATABLE :: X(:,:)
!C                IX     - INPUT ROW DIMENSION OF MATRIX X EXACTLY AS            ^M
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE         ^M
!C                           CALLING PROGRAM.                                    ^M
!C                NG     - INPUT NUMBER OF DISTINCT GROUPS.                      ^M
!C                NV     - INPUT NUMBER OF DATA VARIABLES.                       ^M
!C                ND     - INPUT VECTOR OF LENGTH NG.  ND(I) CONTAINS            ^M
!C                           THE NUMBER OF DATA VECTORS FROM THE I-TH            ^M
!C                           GROUP IN THE DATA MATRIX X.                         ^M
      Integer :: IX, NG, NV
      Integer, ALLOCATABLE :: ND(:)

!C                XM     - OUTPUT MATRIX OF DIMENSION NV BY NG+1.                ^M
!C                           THE I-TH COLUMN OF XM CONTAINS THE                  ^M
!C                           MEAN VECTOR FOR THE I-TH GROUP.  COLUMN             ^M
!C                           NG+1 CONTAINS THE OVERALL MEAN VECTOR.              ^M
!C                IXM    - INPUT ROW DIMENSION OF MATRIX XM EXACTLY AS           ^M
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE         ^M
!C                           CALLING PROGRAM.                                    ^M
!C                NNV    - OUTPUT NUMBER OF NEW VARIABLES COMPUTED.              ^M
!C                           NNV IS GIVEN BY THE MINIMUM OF NV AND NG-1.         ^M
      REAL, ALLOCATABLE :: XM(:,:)
      Integer :: IXM, NNV

!C                E      - OUTPUT VECTOR OF LENGTH NV CONTAINING THE             ^M
!C                           EIGENVALUES IN THE FIRST NNV LOCATIONS IN           ^M
!C                           DESCENDING ORDER. THE REMAINDER OF E IS             ^M
!C                           USED AS WORK AREA. (SEE REMARK 2).                  ^M
!C            2.  THE COLUMN VECTORS IN THE TRANSFORMATION MATRIX C ARE          ^M
!C                IN ORDER OF DECREASING IMPORTANCE, THAT IS, THE EIGEN-         ^M
!C                VECTOR ASSOCIATED WITH THE LARGEST EIGENVALUE IS IN            ^M
!C                THE FIRST COLUMN OF C, AND SO ON.                              ^M
!C                C      - OUTPUT MATRIX OF DIMENSION NV BY NNV.  THE            ^M
!C                           I-TH COLUMN OF C CONTAINS THE VECTOR FOR            ^M
!C                           COMPUTING THE I-TH NEW LINEAR VARIABLE.             ^M
!C                           (SEE REMARK 1).                                     ^M
!C                IC     - INPUT ROW DIMENSION OF MATRIX C EXACTLY AS            ^M
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE         ^M
!C                           CALLING PROGRAM.                                    ^M
      REAL, ALLOCATABLE :: E(:), C(:,:)
      Integer :: IC

!C                SW     - WORK VECTOR OF LENGTH ((NV+1)*NV)/2 OR 4,             ^M
!C                           WHICHEVER IS GREATER. (SEE REMARK 3).               ^M
!C                SB     - WORK MATRIX OF DIMENSION NV BY NV.                    ^M
!C                           (SEE REMARK 3).                                     ^M
!C                ISB    - INPUT ROW DIMENSION OF MATRIX SB EXACTLY AS           ^M
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE         ^M
!C                           CALLING PROGRAM.                                    ^M
!C            3.  IF THE WITHIN AND BETWEEN SUMS OF SQUARES MATRICES, SW         ^M
!C                AND SB, ARE DESIRED, THEY MAY BE PRINTED BY MODIFYING          ^M
!C                THE SOURCE CODE AT THE DESIGNATED COMMENT CARD.                ^M
      REAL, ALLOCATABLE :: SW(:), SB(:,:)
      Integer :: ISB

!C                EX     - COMPLEX WORK VECTOR OF LENGTH NV CONTAINING           ^M
!C                           THE UNSORTED EIGENVALUES.  (SEE REMARK 4).          ^M
!C                         NOTE - THE ROUTINE TREATS EX AS A REAL VECTOR         ^M
!C                           OF LENGTH 2*NV. AN APPROPRIATE EQUIVALENCE          ^M
!C                           STATEMENT MAY BE REQUIRED. SEE DOCUMENT             ^M
!C                           EXAMPLE.                                            ^M
!C                CX     - COMPLEX WORK MATRIX OF DIMENSION NV BY NV             ^M
!C                           CONTAINING THE EIGENVECTORS ASSOCIATED WITH         ^M
!C                           EX.  (SEE REMARK 4).                                ^M
!C                         NOTE - THE ROUTINE TREATS CX AS A REAL VECTOR         ^M
!C                           OF LENGTH 2*NV*NV. AN APPROPRIATE EQUI-             ^M
!C                           VALENCE STATEMENT MAY BE REQUIRED. SEE              ^M
!C                           DOCUMENT EXAMPLE.                                   ^M
!C                ICX    - INPUT ROW DIMENSION OF MATRIX CX EXACTLY AS           ^M
!C                           SPECIFIED IN THE DIMENSION STATEMENT IN THE         ^M
!C                           CALLING PROGRAM.^M
!C            4.  ONE SHOULD PRUDENTLY INSPECT THE COMPLEX EIGENVALUE            ^M
!C                VECTOR AND EIGENVECTOR MATRIX AND CHECK THAT THERE ARE         ^M
!C                ONLY NG-1 OR NV NONZERO EIGENVALUES AND THAT THE IMAG-         ^M
!C                INARY PARTS OF THE EIGENVECTORS ARE TRULY ZERO.                ^M
      REAL, ALLOCATABLE :: EX(:,:), CX(:,:)
      Integer :: ICX, IER
            
!C                IS     - INTEGER WORK VECTOR OF LENGTH NV.                     ^M
      Integer, ALLOCATABLE :: IS(:)

      REAL, ALLOCATABLE :: Y(:), Y1(:), Xtmp(:,:), S(:), Thres(:)   ! Working arrays
      Integer :: NThres
      Real :: Xtc

!      OPEN (1, FILE="data_test_disc.txt") !"Disc_Inp/JMA_VES_TC_All_1-5.txt")
      OPEN (1, FILE="Disc_Inp/JMA_VES_TC_All_1-5.txt")
      READ (1,*) NG, NV, N
      NThres=NG-1

      ALLOCATE (X(NV,N), ND(NG), P(NG), XM(NV,NG), W(NV,NG), D(NG,NG))
      ALLOCATE (Thres(NThres))
      ALLOCATE (Y(NV), Y1(N), Xtmp(N,NV), S(((NV+1)*NV)/2))

      Read (1,*) (Thres(i),i=1,NThres)

      print*,NG, NV, N, Thres(:)


      ND(:)=0
      Do i=1,N
         read (1,*) Xtc, (Xtmp(i,j),j=1,NV)
         Y1(i) = Group (Xtc)
         print "(20F7.3)",Y1(i),(Xtmp(i,j),j=1,NV)
         Do k=1,NG; if (Y1(i)==k) ND(k)=ND(k)+1 ; Enddo
      Enddo
      print*, (ND(k),k=1,NG)
      Close (1)
      ! Sort Y1:
      Do i=1, N-1
         Do j=i+1,N
            if (Y1(j)<Y1(i)) then
               Xtc=Y1(j); Y1(j)=Y1(i); Y1(i)=Xtc
               Y(:)=Xtmp(j,:); Xtmp(j,:)=Xtmp(i,:); Xtmp(i,:)=Y(:)
            endif
         Enddo
      Enddo

      ! Create Data array:
      Do iv=1, NV
         X(iv,:)=Xtmp(:,iv)
      Enddo

      ! PRIOR PROBABILITY:
      Do i=1,NG; P(i)=real(ND(i))/real(N); Enddo
      !Do i=1,NG; P(i)=1.0/Real(NG); Enddo

      ! Parameter:
      IX=NV; IXM=NV; IW=NV; ID=NG;  
      CALL ODNORM (X,IX,NG,NV,ND,P,XM,IXM,S,W,IW,D,ID,IER)

      ! Discriminal Equations:
      Do i=1, NG
         print*,"Group ",i
         print "(20F10.3)", P(i), (W(j,i),j=1,NV)
      Enddo

      ! Probability:
      Do i=1,NG
         print "(50F10.3)", (D(i,j),j=1,NG)
      Enddo

      ! Testing using Discr. functions:
      Do i=1,N
         Do k=1,NG
            Xtmp(i,k)=0
            Do j=1,NV
               Xtmp(i,k)=Xtmp(i,k)+W(j,k)*X(j,i) 
            Enddo
            Xtmp(i,k)=Xtmp(i,k)+P(k)
         Enddo
         Y(i)=GRoup_ID(Xtmp(i,1:NG),NG)
      Enddo
      
      ! Print test:
      i1=0
      Do k=1,NG
         Do j=1,ND(k)
            i1=i1+1
            print*,k, Int(Y(i1))
         Enddo
      Enddo

      CONTAINS
      ! -----------------------------------------------------------------
      Integer Function GRoup_ID (X,N)
         Integer :: N, i, j, ix
         Real :: X(N), Xmax
         Xmax=X(1); ix=1
         Do i=2,N
            if (X(i)>Xmax) ix=i
         Enddo
         GRoup_ID=ix
      End Function GRoup_ID
      ! -----------------------------------------------------------------
      Integer Function Group (x)
       ! Identify which group x is belong to
       Real :: x
       Integer :: k
       Do k=1,NG
          if (k==1) then
             if (x <= Thres(k)) Group=1
          elseif (k==NG) then
             if (x > Thres(k-1)) Group=NG
          else
             if (x > Thres(k-1).AND. x <= Thres(k)) Group=k
          endif
       Enddo
      End Function Group

      END  ! PROGRAM
     ! -----------------------------------------------------------------
     Include "FDA.for"