C @(#)cmds.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 PRINCIPAL COORDINATES ANALYSIS                             C
C              (CLASSICAL MULTIDIMENSIONAL SCALING).                      C
C                                                                         C
C  To call:   CALL PCA(N,A,IPRINT,W1,W2,A2,IERR)    where                 C
C                                                                         C
C                                                                         C
C  N     : integer number of objects.                                     C
C  A     : input distances array, real, of dimensions N by N.             C
C          On output, A contains in the first 7 columns the projections   C
C          of the objects on the first 7 principal components.            C
C  IPRINT: print options.                                                 C
C          = 3: full printing of items calculated.                        C
C          = 2: printing of everything except the input distance matrix.  C
C          Otherwise: no printing.                                        C
C  W1,W2 : real vectors of dimension M (see called routines for use).     C
C          On output, W1 contains the eigenvalues (in increasing order of C
C          magnitude).                                                    C
C  A2    : real array of dimensions M * M (see called routines for use).  C
C  IERR  : error indicator (normally zero).                               C
C                                                                         C
C                                                                         C
C  Inputs here are N, A, IPRINT (and IERR).                               C
C  Output information is contained in A, and W1.                          C
C  All printed outputs are carried out in easily recognizable subroutines C
C  called from the first subroutine following.                            C
C                                                                         C
C                                                                         C
C  F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.       C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE CMDS(N,A,IPRINT,W,FV1,Z,IERR)
        INTEGER N,I1,I2,IPRINT,IERR,N2
        REAL    A(N,N), W(N), FV1(N), Z(N,N)
        REAL    TOT
C
        IF (IPRINT.EQ.3) CALL OUTMAT(N,A)
C
        TOT = 0.0
        DO 200 I1 = 1, N
           W(I1) = 0.0
           DO 100 I2 = 1, N
              W(I1) = W(I1) + A(I1,I2)
              TOT = TOT + A(I1,I2)
  100      CONTINUE
           W(I1) = W(I1)/FLOAT(N)
  200   CONTINUE
        TOT = TOT/(FLOAT(N)*FLOAT(N))
C
        DO 300 I1 = 1, N
           DO 300 I2 = 1, N
              A(I1,I2) = -0.5 * (A(I1,I2) - W(I1) - W(I2) - TOT)
  300   CONTINUE
C
C          Carry out the eigenreduction.
C
        N2 = N
        CALL TRED2(N,N2,A,W,FV1,Z)
        CALL TQL2(N,N2,W,FV1,Z,IERR)
        IF (IERR.NE.0) GOTO 9000
C
C          Output eigenvalues and eigenvectors.
C
        IF (IPRINT.GE.2) CALL OUTEVL(N,W)
        IF (IPRINT.GE.2) CALL OUTEVC(N,Z)
C
C          Determine projections and output them.
C
        CALL PROJN(N,W,A,Z,FV1)
        IF (IPRINT.GE.2) CALL OUTPR(N,A)
C
 9000   RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C Reduce a real, symmetric matrix to a symmetric, tridiagonal matrix.     C
C                                                                         C
C To call:    CALL TRED2(NM,N,A,D,E,Z)    where                           C
C                                                                         C
C NM = row dimension of A and Z;                                          C
C N = order of matrix A (will always be <= NM);                           C
C A = symmetric matrix of order N to be reduced to tridiagonal form;      C
C D = vector of dim. N containing, on output, diagonal elts. of trid.     C
C     matrix.                                                             C
C E = working vector of dim. at least N-1 to contain subdiagonal elts.;   C
C Z = matrix of dims. NM by N containing, on output, orthogonal           C
C     transformation matrix producing the reduction.                      C
C                                                                         C
C Normally a call to TQL2 will follow the call to TRED2 in order to       C
C produce all eigenvectors and eigenvalues of matrix A.                   C
C                                                                         C
C Algorithm used: Martin et al., Num. Math. 11, 181-195, 1968.            C
C                                                                         C
C Reference: Smith et al., Matrix Eigensystem Routines - EISPACK          C
C Guide, Lecture Notes in Computer Science 6, Springer-Verlag, 1976,      C
C pp. 489-494.                                                            C
C                                                                         C
C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.        C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE TRED2(NM,N,A,D,E,Z)
        INTEGER NM,N,I,J,L,II,K,JP1
        REAL A(NM,N),D(N),E(N),Z(NM,N)
        REAL SCALE,F,G,H,HH
C
        DO 100 I = 1, N
           DO 100 J = 1, I
              Z(I,J) = A(I,J)
  100   CONTINUE
        IF (N.EQ.1) GOTO 320
        DO 300 II = 2, N
           I = N + 2 - II
           L = I - 1
           H = 0.0
           SCALE = 0.0
           IF (L.LT.2) GOTO 130
           DO 120 K = 1, L
              SCALE = SCALE + ABS(Z(I,K))
  120      CONTINUE
           IF (SCALE.NE.0.0) GOTO 140
  130      E(I) = Z(I,L)
           GOTO 290
  140      DO 150 K = 1, L
              Z(I,K) = Z(I,K)/SCALE
              H = H + Z(I,K)*Z(I,K)
  150      CONTINUE
C
           F = Z(I,L)
           G = -SIGN(SQRT(H),F)
           E(I) = SCALE * G
           H = H - F * G
           Z(I,L) = F - G
           F = 0.0
C
           DO 240 J = 1, L
              Z(J,I) = Z(I,J)/H
              G = 0.0
C             Form element of A*U.
              DO 180 K = 1, J
                 G = G + Z(J,K)*Z(I,K)
  180         CONTINUE
              JP1 = J + 1
              IF (L.LT.JP1) GOTO 220
              DO 200 K = JP1, L
                 G = G + Z(K,J)*Z(I,K)
  200         CONTINUE
C             Form element of P where P = I - U U' / H .
  220         E(J) = G/H
              F = F + E(J) * Z(I,J)
  240      CONTINUE
           HH = F/(H + H)
C          Form reduced A.
           DO 260 J = 1, L
              F = Z(I,J)
              G = E(J) - HH * F
              E(J) = G
              DO 250 K = 1, J
                 Z(J,K) = Z(J,K) - F*E(K) - G*Z(I,K)
  250         CONTINUE
  260      CONTINUE
  290      D(I) = H
  300   CONTINUE
  320   D(1) = 0.0
        E(1) = 0.0
C       Accumulation of transformation matrices.
        DO 500 I = 1, N
           L = I - 1
           IF (D(I).EQ.0.0) GOTO 380
           DO 360 J = 1, L
              G = 0.0
              DO 340 K = 1, L
                 G = G + Z(I,K) * Z(K,J)
  340         CONTINUE
              DO 350 K = 1, L
                 Z(K,J) = Z(K,J) - G * Z(K,I)
  350         CONTINUE
  360      CONTINUE
  380      D(I) = Z(I,I)
           Z(I,I) = 1.0
           IF (L.LT.1) GOTO 500
           DO 400 J = 1, L
              Z(I,J) = 0.0
              Z(J,I) = 0.0
  400      CONTINUE
  500   CONTINUE
C
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C Determine eigenvalues and eigenvectors of a symmetric,                  C
C tridiagonal matrix.                                                     C
C                                                                         C
C To call:    CALL TQL2(NM,N,D,E,Z,IERR)    where                         C
C                                                                         C
C NM = row dimension of Z;                                                C
C N = order of matrix Z;                                                  C
C D = vector of dim. N containing, on output, eigenvalues;                C
C E = working vector of dim. at least N-1;                                C
C Z = matrix of dims. NM by N containing, on output, eigenvectors;        C
C IERR = error, normally 0, but 1 if no convergence.                      C
C                                                                         C
C Normally the call to TQL2 will be preceded by a call to TRED2 in        C
C order to set up the tridiagonal matrix.                                 C
C                                                                         C
C Algorithm used: QL method of Bowdler et al., Num. Math. 11,             C
C 293-306, 1968.                                                          C
C                                                                         C
C Reference: Smith et al., Matrix Eigensystem Routines - EISPACK          C
C Guide, Lecture Notes in Computer Science 6, Springer-Verlag, 1976,      C
C pp. 468-474.                                                            C
C                                                                         C
C F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.        C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE TQL2(NM,N,D,E,Z,IERR)
        INTEGER N,NM,I,J,IERR,M,L1,L
        INTEGER II,K,MML
        REAL    D(N), E(N), Z(NM,N),EPS,F,B,H
        REAL    G,P,R,C,S
        DATA    EPS/1.E-12/
C
        IERR = 0
        IF (N.EQ.1) GOTO 1001
        DO 100 I = 2, N
           E(I-1) = E(I)
  100   CONTINUE
        F = 0.0
        B = 0.0
        E(N) = 0.0
C
        DO 240 L = 1, N
           J = 0
           H = EPS * (ABS(D(L)) + ABS(E(L)))
           IF (B.LT.H) B = H
C          Look for small sub-diagonal element.
           DO 110 M = L, N
              IF (ABS(E(M)).LE.B) GOTO 120
C             E(N) is always 0, so there is no exit through the
C             bottom of the loop.
  110      CONTINUE
  120      IF (M.EQ.L) GOTO 220
  130      IF (J.EQ.30) GOTO 1000
           J = J + 1
C          Form shift.
           L1 = L + 1
           G = D(L)
           P = (D(L1)-G)/(2.0*E(L))
           R = SQRT(P*P+1.0)
           D(L) = E(L)/(P+SIGN(R,P))
           H = G-D(L)
C
           DO 140 I = L1, N
              D(I) = D(I) - H
  140      CONTINUE
C
           F = F + H
C          QL transformation.
           P = D(M)
           C = 1.0
           S = 0.0
           MML = M - L
C
           DO 200 II = 1, MML
              I = M - II
              G = C * E(I)
              H = C * P
              IF (ABS(P).LT.ABS(E(I))) GOTO 150
              C = E(I)/P
              R = SQRT(C*C+1.0)
              E(I+1) = S * P * R
              S = C/R
              C = 1.0/R
              GOTO 160
  150         C = P/E(I)
              R = SQRT(C*C+1.0)
              E(I+1) = S * E(I) * R
              S = 1.0/R
              C = C * S
  160         P = C * D(I) - S * G
              D(I+1) = H + S * (C * G + S * D(I))
C             Form vector.
              DO 180 K = 1, N
                 H = Z(K,I+1)
                 Z(K,I+1) = S * Z(K,I) + C * H
                 Z(K,I) = C * Z(K,I) - S * H
  180         CONTINUE
  200      CONTINUE
           E(L) = S * P
           D(L) = C * P
           IF (ABS(E(L)).GT.B) GOTO 130
  220      D(L) = D(L) + F
  240   CONTINUE
C
C       Order eigenvectors and eigenvalues.
        DO 300 II = 2, N
           I = II - 1
           K = I
           P = D(I)
           DO 260 J = II, N
              IF (D(J).GE.P) GOTO 260
              K = J
              P = D(J)
  260      CONTINUE
           IF (K.EQ.I) GOTO 300
           D(K) = D(I)
           D(I) = P
           DO 280 J = 1, N
              P = Z(J,I)
              Z(J,I) = Z(J,K)
              Z(J,K) = P
  280      CONTINUE
  300   CONTINUE
C
        GOTO 1001
C       Set error - no convergence to an eigenvalue after 30 iterations.
 1000   IERR = L
 1001   RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output array.                                                          C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTMAT(N,ARRAY)
        INTEGER N,K1,K2
        REAL ARRAY(N,N)
C
        DO 100 K1 = 1, N
           WRITE (6,1000) (ARRAY(K1,K2),K2=1,N)
  100   CONTINUE
C
 1000   FORMAT(10(2X,F8.4))
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output eigenvalues in order of decreasing value.                       C
C  Ignore first (trivial) eigenvalue.                                     C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTEVL(NVALS,VALS)
        INTEGER K,NVALS
        REAL VALS(NVALS),CUM,VPC,VCPC,TOT
C
        TOT = 0.0
        DO 100 K = 1, NVALS-1
           TOT = TOT + VALS(K)
  100   CONTINUE
C
        WRITE (6,1000)
        CUM = 0.0
        K = NVALS
        WRITE (6,1010)
        WRITE (6,1020)
  200   CONTINUE
        K = K - 1
        CUM = CUM + VALS(K)
        VPC = VALS(K) * 100.0 / TOT
        VCPC = CUM * 100.0 / TOT
        WRITE (6,1030) VALS(K),VPC,VCPC
        IF (K.GT.1) GOTO 200
C
        RETURN
 1000   FORMAT(1H0,'EIGENVALUES FOLLOW.',/)
 1010   FORMAT
     X  (' Eigenvalues       As Percentages    Cumul. Percentages')
 1020   FORMAT
     X  (' -----------       --------------    ------------------')
 1030   FORMAT(F10.4,9X,F10.4,10X,F10.4)
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C         Output FIRST SEVEN eigenvectors associated with eigenvalues     C
C         in descending order.                                            C
C         Ignore first (trivial) eigenvector.                             C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTEVC(NDIM,VECS)
        INTEGER NDIM,K1,K2,NUM,KDUM
        REAL  VECS(NDIM,NDIM)
C
        NUM = MIN0(NDIM,7)
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K1 = 1, NDIM
           KDUM = NDIM - K2
           WRITE (6,1030) K1,(VECS(K1,KDUM),K2=1,NUM)
  100   CONTINUE
C
        RETURN
 1000   FORMAT(1H0,'EIGENVECTORS FOLLOW.',/)
 1010   FORMAT('  VBLE.   EV-1    EV-2    EV-3    EV-4    EV-5    EV-6
     X   EV-7')
 1020   FORMAT(' ------  ------  ------  ------  ------  ------  ------
     X------')
 1030   FORMAT(I5,2X,7F8.4)
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output projections of objects on first 7 principal components.         C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTPR(N,PRJNS)
        INTEGER N,K,J,NUM
        REAL    PRJNS(N,N)
C
        NUM = MIN0(N,7)
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, N
           WRITE (6,1030) K,(PRJNS(K,J),J=1,NUM)
  100   CONTINUE
C
 1000   FORMAT(1H0,'PROJECTIONS OF OBJECTS FOLLOW.',/)
 1010   FORMAT('  VBLE.  PROJ-1  PROJ-2  PROJ-3  PROJ-4  PROJ-5  PROJ-6
     X  PROJ-7')
 1020   FORMAT(' ------  ------  ------  ------  ------  ------  ------
     X  ------')
 1030   FORMAT(I5,2X,7F8.4)
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Determine projections of objects on 7 principal components.            C
C  Ignore first (trivial) axis.                                           C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE PROJN(N,EVALS,A,Z,VEC)
        INTEGER N,J1,J2,L,NUM,J3,JDUM
        REAL    EVALS(N), A(N,N), Z(N,N), VEC(N)
C
        NUM = MIN0(N,7)
        DO 300 J1 = 1, N
           DO 50 L = 1, N
              VEC(L) = A(J1,L)
   50      CONTINUE
           DO 200 J2 = 1, NUM
              A(J1,J2) = 0.0
              JDUM = N - J2
              DO 100 J3 = 1, N
                 A(J1,J2) = A(J1,J2) + VEC(J3)*Z(J3,JDUM)
  100         CONTINUE
              IF (EVALS(N-J2).GT.0.0) A(J1,J2) =
     X                                  A(J1,J2)/SQRT(EVALS(JDUM))
              IF (EVALS(N-J2).LE.0.0) A(J1,J2) = 0.0
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END
