C @(#)lda.for	19.1 (ES0-DMD) 02/25/03 13:27:14
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                       C
C  Carry out a LINEAR DISCRIMINANT ANALYSIS, assigning ungrouped items  C
C  to the closest group centre, using the Mahalanobis distance.         C
C                                                                       C
C                                                                       C
C  To call:   CALL LDA(N,M,DATA,GP,IPRINT,MEAN,MGP,TOTAL,DIFFVEC,       C
C                                W1,W2,IERR)     where                  C
C                                                                       C
C                                                                       C
C  N, M  : integer dimensions of ...                                    C
C  DATA  : input data (real).                                           C
C          On output the first column of DATA contains projections of   C
C          all N items on the line connecting the two group means.      C
C          Zero is the boundary point.                                  C
C  GP    : Integer vector of length N giving group assignments.  An     C
C          unassigned item has group 0.  Otherwise groups 1 and 2 will  C
C          be looked for, and other values here are not acceptable.     C
C  IPRINT: integer; print options (= 3: full; otherwise none).          C
C  MEAN  : real vector of length M (number of attributes or variables). C
C  MGP   : real array of dimensions 2 by M.                             C
C  TOTAL : real array of dimensions M by M; on output contains inverse  C
C          of total variance/covariance matrix.                         C
C  DIFFVEC: real vector of length M.                                    C
C  W1, W2: real vectors of length M.                                    C
C                                                                       C
C                                                                       C
C  Inputs here are N, M, DATA, GP, IPRINT (and IERR).                   C
C  The principle output information is contained in DATA.               C
C  IERR = 1 means that more than two groups have been specified; IERR   C
C  = 2 means that the total variance-covariance matrx was singular.     C
C                                                                       C
C  Note: we require N > M > 2, to prevent the seeking of the inverse    C
C        of a singular matrix.                                          C
C                                                                       C
C                                                                       C
C  F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.     C
C                                                                       C
C-----------------------------------------------------------------------C
        SUBROUTINE LDA(N,M,DATA,GP,IPRINT,MEAN,MGP,TOTAL,DIFFVEC,W1,W2,
     X                                                          IERR)
        REAL    DATA(N,M), TOTAL(M,M), MEAN(M), MGP(2,M), DIFFVEC(M)
        REAL    W1(M), W2(M),D
        INTEGER GP(N), NOG(2), G,IMAT
        INTEGER N,M,IPRINT,IERR,NEFF,I,J
        INTEGER K,J1,J2
C
C          Form global mean.
C
        IERR = 0
        NEFF = 0
        DO 50 I = 1, N
           IF (GP(I).NE.0) NEFF = NEFF + 1
           IF (GP(I).LE.2) GOTO 40
              IERR = 1
              GOTO 9000
   40   CONTINUE
   50   CONTINUE
C
        DO 200 J = 1, M
           MEAN(J) = 0.0
           DO 100 I = 1, N
              IF (GP(I).NE.0) MEAN(J) = MEAN(J) + DATA(I,J)
  100      CONTINUE
           MEAN(J) = MEAN(J)/FLOAT(NEFF)
  200   CONTINUE
C
C          Form (total) variance-covariance matrix.
C
        DO 500 J1 = 1, M
           DO 400 J2 = 1, M
              TOTAL(J1,J2) = 0.0
              DO 300 I = 1, N
                 IF (GP(I).NE.0) TOTAL(J1,J2) = TOTAL(J1,J2) +
     X              (DATA(I,J1)-MEAN(J1))*(DATA(I,J2)-MEAN(J2))
  300         CONTINUE
              TOTAL(J1,J2) = TOTAL(J1,J2)/FLOAT(NEFF)
  400      CONTINUE
  500   CONTINUE
C
        IMAT = 1
C        IF (IPRINT.EQ.3) CALL OUTMTL(IMAT,M,TOTAL)
C
C          Form group means.
C
        DO 700 J = 1, M
           DO 600 K = 1, 2
              MGP(K,J) = 0.0
  600      CONTINUE
  700   CONTINUE
C
        DO 900 I = 1, N
           G = GP(I)
           IF (G.EQ.0) GOTO 900
           NOG(G) = NOG(G) + 1
           DO 800 J = 1, M
              MGP(G,J) = MGP(G,J) + DATA(I,J)
  800      CONTINUE
  900   CONTINUE
C
        DO 1100 K = 1, 2
           DO 1000 J = 1, M
              MGP(K,J) = MGP(K,J)/NOG(K)
 1000      CONTINUE
 1100   CONTINUE
C
C          Invert variance-covariance matrix.
C
        CALL MINVL(M,TOTAL,D,W1,W2)
        IF (ABS(D).GT.0.00001) GOTO 1150
           IERR = 2
           GOTO 9000
 1150   CONTINUE
C
C          Form difference vector of group mean vectors.
C
        DO 1200 J = 1, M
           DIFFVEC(J) = MGP(1,J) - MGP(2,J)
           MEAN(J) = (MGP(1,J) + MGP(2,J))/2.
 1200   CONTINUE
C
C          Determine projections and output them.
C
        CALL PROJXL(N,M,DATA,MEAN,W1,TOTAL,DIFFVEC)
        IF (IPRINT.EQ.3) CALL OUTPXL(N,M,DATA)
C
C
C
 9000   CONTINUE
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                       C
C  Output a matrix.                                                     C
C                                                                       C
C-----------------------------------------------------------------------C
        SUBROUTINE OUTMTL(IMAT,M,ARRAY)
C
C          Output array.
C
        INTEGER IMAT,M,K1,K2
        REAL ARRAY(M,M)
C
C        IF (IMAT.EQ.1)WRITE (6,900)
        DO 100 K1 = 1, M
           WRITE (6,1000) (ARRAY(K1,K2),K2=1,M)
  100   CONTINUE
C
C   900   FORMAT(' VARIANCE/COVARIANCE MATRIX FOLLOWS.',/)
 1000   FORMAT(5(2X,F8.4))
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                      C
C  Invert a symmetric matrix and calculate its determinant.            C
C                                                                      C
C                                                                      C
C  To call:      CALL MATINV(M,ARRAY,DET,W1,W2)   where                C
C                                                                      C
C                                                                      C
C  M       : dimension of ...                                          C
C  ARRAY   : input matrix which is replaced by its inverse.            C
C  NORDER  : degree of matrix (order of determinant)                   C
C  DET     : determinant of input matrix.                              C
C  W1, W2  : work vectors of dimension M.                              C
C                                                                      C
C                                                                      C
C  Reference: Philip B Bevington, "Data Reduction and Error Analysis   C
C             for the Physical Sciences", McGraw-Hill, New York, 1969, C
C             pp. 300-303.                                             C
C                                                                      C
C                                                                      C
C  F. Murtagh, ST-ECF, Garching-bei-Muenchen, January 1986.            C
C                                                                      C
C----------------------------------------------------------------------C
        SUBROUTINE MINVL(M,ARRAY,DET,IK,JK)
        INTEGER M,K,I,J,L
        REAL    ARRAY(M,M), IK(M), JK(M)
        REAL    DET,SAVE,AMAX
C
        DET = 1.0
        DO 100 K = 1, M
C       Find largest element ARRAY(I,J) in rest of matrix.
        AMAX = 0.0
   21      DO 30 I = K, M
              DO 30 J = K, M
                 IF (ABS(AMAX)-ABS(ARRAY(I,J))) 24,24,30
   24            AMAX = ARRAY(I,J)
                 IK(K) = I
                 JK(K) = J
   30      CONTINUE
C          Interchange rows and columns to put AMAX in ARRAY(K,K).
           IF (AMAX) 41,32,41
   32      DET = 0.0
           GOTO 140
   41      I = IK(K)
           IF (I-K) 21,51,43
   43      DO 50 J = 1, M
              SAVE = ARRAY(K,J)
              ARRAY(K,J) = ARRAY(I,J)
   50      ARRAY(I,J) = -SAVE
   51      J = JK(K)
           IF (J-K) 21,61,53
   53      DO 60 I = 1, M
              SAVE = ARRAY(I,K)
              ARRAY(I,K) = ARRAY(I,J)
   60      ARRAY(I,J) = -SAVE
C          Accumulate elements of inverse matrix.
   61      DO 70 I = 1, M
              IF (I-K) 63,70,63
   63         ARRAY(I,K) = -ARRAY(I,K)/AMAX
   70      CONTINUE
           DO 80 I = 1, M
              DO 80 J = 1, M
                 IF (I-K) 74,80,74
   74            IF (J-K) 75,80,75
   75            ARRAY(I,J) = ARRAY(I,J) + ARRAY(I,K)*ARRAY(K,J)
   80      CONTINUE
           DO 90 J = 1, M
              IF (J-K) 83,90,83
   83         ARRAY(K,J) = ARRAY(K,J)/AMAX
   90      CONTINUE
           ARRAY(K,K) = 1.0/AMAX
  100   DET = DET * AMAX
C       Restore ordering of matrix.
        DO 130 L = 1, M
           K = M - L + 1
           J = IK(K)
           IF (J-K) 111,111,105
  105      DO 110 I = 1, M
              SAVE = ARRAY(I,K)
              ARRAY(I,K) = -ARRAY(I,J)
  110      ARRAY(I,J) = SAVE
  111      I = JK(K)
           IF (I-K) 130,130,113
  113      DO 120 J = 1, M
              SAVE = ARRAY(K,J)
              ARRAY(K,J) = -ARRAY(I,J)
  120      ARRAY(I,J) = SAVE
  130   CONTINUE
  140   RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                       C
C  Output projections of row points.                                    C
C                                                                       C
C-----------------------------------------------------------------------C
        SUBROUTINE OUTPXL(N,M,PRJN)
        INTEGER N,M,NUM,K,J
        REAL    PRJN(N,M)
C
C
        NUM = 1
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, N
           WRITE (6,1030) K,(PRJN(K,J),J=1,NUM)
  100   CONTINUE
C
 1000   FORMAT(1X,'PROJECTIONS OF ROW-POINTS FOLLOW.',/)
 1010   FORMAT(' OBJECT   PROJN')
 1020   FORMAT(' ------  ------')
 1030   FORMAT(I5,2X,F8.4)
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                       C
C  Form projections of row-points on factors.                           C
C                                                                       C
C-----------------------------------------------------------------------C
        SUBROUTINE PROJXL(N,M,DATA,MEAN,VEC,TOTINV,DIFF)
        INTEGER N,M,K,L,I,J1,NUM,J2
        REAL DATA(N,M),MEAN(M),VEC(M),TOTINV(M,M),DIFF(M)
C
        NUM = 1
        DO 300 K = 1, N
           DO 50 L = 1, M
              VEC(L) = DATA(K,L)
   50      CONTINUE
           DO 200 I = 1, NUM
              DATA(K,I) = 0.0
              DO 100 J1 = 1, M
                 DO 75 J2 = 1, M
                    DATA(K,I) = DATA(K,I) + (VEC(J1)-MEAN(J1))*
     X                          TOTINV(J1,J2)*DIFF(J2)
   75            CONTINUE
  100         CONTINUE
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END

