C @(#)pca.for	19.1 (ES0-DMD) 02/25/03 13:27:15
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 COMPONENTS ANALYSIS                              C
C              (KARHUNEN-LOEVE EXPANSION).                                C
C                                                                         C
C  To call:   CALL PCA(N,M,DATA,METHOD,IPRINT,A1,W1,W2,A2,IERR)    where  C
C                                                                         C
C                                                                         C
C  N, M  : integer dimensions of ...                                      C
C  DATA  : input data.                                                    C
C          On output, DATA contains in first 7 columns the projections of C
C          the row-points on the first 7 principal components.            C
C  METHOD: analysis option.                                               C
C          = 1: on sums of squares & cross products matrix.               C
C          = 2: on covariance matrix.                                     C
C          = 3: on correlation matrix.                                    C
C  IPRINT: print options.                                                 C
C          = 0: no printed output - arrays/vectors, only, contain items   C
C              calculated.                                                C
C          = 1: eigenvalues, only, output.                                C
C          = 2: printed output, in addition, of correlation (or other)    C
C               matrix, eigenvalues and eigenvectors.                     C
C          = 3: full printing of items calculated.                        C
C  A1    : correlation, covariance or sums of squares & cross-products    C
C          matrix, dimensions M * M.                                      C
C          On output, A1 contains in the first 7 columns the projections  C
C          of the column-points on the first 7 principal components.      C
C  W1,W2 : real vectors of dimension M (see called routines for use).     C
C          On output, W1 contains the cumulative percentage variances     C
C          associated with the principal components.                      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, M, DATA, METHOD, IPRINT (and IERR).                 C
C  Output information is contained in DATA, A1, 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  Regarding precision and tolerances, we believe the eigen-extraction    C
C  routine is reasonably precise; however, in subroutine PROJY we have    C
C  taken a zero eigenvalue as being less than or equal to 0.00005 in      C
C  value.  This is adequate for most practical situations.                C
C                                                                         C
C                                                                         C
C  F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.       C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE PCA(N,M,DATA,METHOD,IPRINT,A,W,FV1,Z,IERR)
        INTEGER N,M,METHOD,IPRINT,IERR,M2
        REAL    DATA(N,M), A(M,M), W(M), FV1(M), Z(M,M)
C
        IF (METHOD.EQ.1) GOTO 100
        IF (METHOD.EQ.2) GOTO 400
C       If method.eq.3 or otherwise ...
        GOTO 700
C
C          Form sums of squares and cross-products matrix.
C
  100   CONTINUE
        CALL SCPCLP(N,M,DATA,A,IERR)
        IF (IERR.NE.0) GOTO 9000
C
        IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A)
C
C          Now do the PCA.
C
        GOTO 1000
C
C          Form covariance matrix.
C
  400   CONTINUE
        CALL COVCLP(N,M,DATA,W,A,IERR)
        IF (IERR.NE.0) GOTO 9000
C
        IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A)

C
C          Now do the PCA.
C
        GOTO 1000
C
C          Construct correlation matrix.
C
  700   CONTINUE
        CALL CORCLP(N,M,DATA,W,FV1,A,IERR)
        IF (IERR.NE.0) GOTO 9000
C
        IF (IPRINT.GT.1) CALL OUTHMP(METHOD,M,A)

C
C          Now do the PCA.
C
        GOTO 1000
C
C          Carry out eigenreduction.
C
 1000   M2 = M
        CALL TRED2P(M,M2,A,W,FV1,Z)
        CALL TQL2P(M,M2,W,FV1,Z,IERR)
        IF (IERR.NE.0) GOTO 9000
C
C          Output eigenvalues and eigenvectors.
C
        IF (IPRINT.GT.0) CALL OUTVLP(N,M,W)
        IF (IPRINT.GT.1) CALL OUTVCP(N,M,Z)
C
C          Determine projections and output them.
C
        CALL PROJXP(N,M,DATA,Z,FV1)
        IF (IPRINT.EQ.3) CALL OUTPXP(N,M,DATA)
        CALL PROJYP(M,W,A,Z,FV1)
        IF (IPRINT.EQ.3) CALL OUTPYP(M,A)
C
 9000   RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Determine correlations of columns.                                     C
C  First determine the means of columns, storing in WORK1.                C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE CORCLP(N,M,DATA,WORK1,WORK2,OUT,ISTAT)
        INTEGER N,M,ISTAT,I,J,J1,J2
        REAL   DATA(N,M), OUT(M,M), WORK1(M), WORK2(M)
        REAL   EPS
        DATA            EPS/1.E-10/
C
        DO 30 J = 1, M
           WORK1(J) = 0.0
           DO 20 I = 1, N
              WORK1(J) = WORK1(J) + DATA(I,J)
   20      CONTINUE
           WORK1(J) = WORK1(J)/FLOAT(N)
   30   CONTINUE
C
C          Next det. the std. devns. of cols., storing in WORK2.
C
        DO 50 J = 1, M
           WORK2(J) = 0.0
           DO 40 I = 1, N
              WORK2(J) = WORK2(J) + 
     X                   (DATA(I,J)-WORK1(J))*(DATA(I,J)-WORK1(J))
   40      CONTINUE
           WORK2(J) = WORK2(J)/FLOAT(N)
           WORK2(J) = SQRT(WORK2(J))
           IF (WORK2(J).LE.EPS) WORK2(J) = 1.0
   50   CONTINUE
C
C          Now centre and reduce the column points.
C
        DO 70 I = 1, N
           DO 60 J = 1, M
              DATA(I,J) = (DATA(I,J)-WORK1(J))/(SQRT(FLOAT(N))*WORK2(J))
   60      CONTINUE
   70   CONTINUE
C
C          Finally calculate the cross product of the redefined data matrix.
C
        DO 100 J1 = 1, M-1
           OUT(J1,J1) = 1.0
           DO 90 J2 = J1+1, M
              OUT(J1,J2) = 0.0
              DO 80 I = 1, N
                 OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2)
   80         CONTINUE
              OUT(J2,J1) = OUT(J1,J2)
   90      CONTINUE
  100   CONTINUE
        OUT(M,M) = 1.0
C
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Determine covariances of columns.                                      C
C  First determine the means of columns, storing in WORK.                 C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE COVCLP(N,M,DATA,WORK,OUT,ISTAT)
        INTEGER N,M,ISTAT,I,J,J1,J2
        REAL    DATA(N,M), OUT(M,M), WORK(M)
C
        DO 30 J = 1, M
           WORK(J) = 0.0
           DO 20 I = 1, N
              WORK(J) = WORK(J) + DATA(I,J)
   20      CONTINUE
           WORK(J) = WORK(J)/FLOAT(N)
   30   CONTINUE
C
C          Now centre the column points.
C
        DO 50 I = 1, N
           DO 40 J = 1, M
              DATA(I,J) = DATA(I,J)-WORK(J)
   40      CONTINUE
   50   CONTINUE
C
C          Finally calculate the cross product matrix of the redefined
C          data matrix.
C
        DO 80 J1 = 1, M
           DO 70 J2 = J1, M
              OUT(J1,J2) = 0.0
              DO 60 I = 1, N
                 OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2)
   60         CONTINUE
              OUT(J2,J1) = OUT(J1,J2)
   70      CONTINUE
   80   CONTINUE
C
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Determine sums of squares and cross-products of columns.               C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE SCPCLP(N,M,DATA,OUT,ISTAT)
        INTEGER M,N,J1,J2,I,ISTAT
        REAL    DATA(N,M), OUT(M,M)
C
        DO 30 J1 = 1, M
           DO 20 J2 = J1, M
              OUT(J1,J2) = 0.0
              DO 10 I = 1, N
                 OUT(J1,J2) = OUT(J1,J2) + DATA(I,J1)*DATA(I,J2)
   10         CONTINUE
              OUT(J2,J1) = OUT(J1,J2)
   20      CONTINUE
   30   CONTINUE
C
        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 TRED2P(NM,N,A,D,E,Z)
        INTEGER I,J,N,L,K,II,NM,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 TQL2P(NM,N,D,E,Z,IERR)
        INTEGER I,N,IERR,M,NM,L,J,L1,II,MML,K
        REAL    D(N), E(N), Z(NM,N),EPS,F,B,H
        REAL    P,R,G,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 = 1
 1001   RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output array.                                                          C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTMTP(N,M,ARRAY)
        INTEGER K1,K2,N,M
        REAL ARRAY(N,M)
C
        DO 100 K1 = 1, N
           WRITE (6,1000) (ARRAY(K1,K2),K2=1,M)
  100   CONTINUE
C
 1000   FORMAT(10(2X,F8.4))
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output half of (symmetric) array.                                      C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTHMP(ITYPE,NDIM,ARRAY)
        INTEGER ITYPE,NDIM,K1,K2,STATUS
        REAL ARRAY(NDIM,NDIM)
        CHARACTER*80 MYLINE
C
        IF (ITYPE.EQ.1) WRITE (MYLINE,1000)
        IF (ITYPE.EQ.2) WRITE (MYLINE,2000)
        IF (ITYPE.EQ.3) WRITE (MYLINE,3000)
        CALL STTPUT(MYLINE,STATUS)
C
        DO 100 K1 = 1, NDIM
           WRITE (MYLINE,4000) (ARRAY(K1,K2),K2=1,K1)
           CALL STTPUT(MYLINE,STATUS)
  100   CONTINUE
C
 1000   FORMAT(1X ,'SUMS OF SQUARES & CROSS-PRODUCTS MATRIX FOLLOWS.')
 2000   FORMAT(1X ,'COVARIANCE MATRIX FOLLOWS.')
 3000   FORMAT(1X ,'CORRELATION MATRIX FOLLOWS.')
 4000   FORMAT(8(2X,F8.4))
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Output eigenvalues in order of decreasing value.                       C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTVLP(N,NVALS,VALS)
        INTEGER K,NVALS,N,M,ISTAT
        REAL       VALS(NVALS),TOT,CUM
        REAL       VPC,VCPC
        CHARACTER       LINE*80
C
        TOT = 0.0
        DO 100 K = 1, NVALS
           TOT = TOT + VALS(K)
  100   CONTINUE
C
C       WRITE (6,1000)
        WRITE (LINE, 900)
        CALL STTPUT(LINE,ISTAT)
        WRITE(6,900)
        WRITE (LINE,1000)
        CALL STTPUT(LINE,ISTAT)
        WRITE(6,900)
        WRITE (LINE, 900)
        CALL STTPUT(LINE,ISTAT)
        CUM = 0.0
        K = NVALS + 1
C
        M = NVALS
C       (We only want Min(nrows,ncols) eigenvalues output:)
        M = MIN0(N,NVALS)
C
C       WRITE (6,1010)
C       WRITE (6,1020)
        WRITE (LINE,1010)
        CALL STTPUT(LINE,ISTAT)
        WRITE (LINE,1020)
        CALL STTPUT(LINE,ISTAT)
  200   CONTINUE
        K = K - 1
        CUM = CUM + VALS(K)
        VPC = VALS(K) * 100.0 / TOT
        VCPC = CUM * 100.0 / TOT
C       WRITE (6,1030) VALS(K),VPC,VCPC
        WRITE (LINE,1030) VALS(K),VPC,VCPC
        CALL STTPUT(LINE,ISTAT)
C        VALS(K) = VCPC
        IF (K.GT.NVALS-M+1) GOTO 200
C
        RETURN
  900   FORMAT('   ')
 1000   FORMAT(1X,' EIGENVALUES FOLLOW.')
 1010   FORMAT
     X  (' Eigenvalues        As Percentages    Cumul. Percentages')
 1020   FORMAT
     X  (' -----------        --------------    ------------------')
 1030   FORMAT(F13.4,7X,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                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTVCP(N,NDIM,VECS)
        INTEGER NDIM,NUM,N,K1,K2
        REAL    VECS(NDIM,NDIM)
C
        NUM = MIN0(N,NDIM,7)
C
        WRITE(6,985)
        WRITE (6,1000)
        WRITE(6,985)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K1 = 1, NDIM
        WRITE (6,1030) K1,(VECS(K1,NDIM-K2+1),K2=1,NUM)
  100   CONTINUE
C
        RETURN
  985   FORMAT('  ')
 1000   FORMAT(1X ,'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 row-points on first 7 pricipal components.       C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTPXP(N,M,PRJN)
        INTEGER N,M,K,NUM,J
        REAL    PRJN(N,M)
C
        NUM = MIN0(M,7)
        WRITE(6,985)
        WRITE (6,1000)
        WRITE(6,985)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, N
           WRITE (6,1030) K,(PRJN(K,J),J=1,NUM)
  100   CONTINUE
C
  985   FORMAT('  ')
 1000   FORMAT(1X ,'PROJECTIONS OF ROW-POINTS FOLLOW.')
 1010   FORMAT(' OBJECT  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  Output projections of columns on first 7 principal components.         C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE OUTPYP(M,PRJNS)
        INTEGER K,NUM,M,J
        REAL    PRJNS(M,M)
C
        NUM = MIN0(M,7)
        WRITE(6,985)
        WRITE (6,1000)
        WRITE(6,985)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, M
           WRITE (6,1030) K,(PRJNS(K,J),J=1,NUM)
  100   CONTINUE
C
  985   FORMAT('  ') 
 1000   FORMAT(1X ,'PROJECTIONS OF COLUMN-POINTS 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  Form projections of row-points on first 7 principal components.        C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE PROJXP(N,M,DATA,EVEC,VEC)
        INTEGER N,M,K,NUM,L,I,J
        REAL    DATA(N,M), EVEC(M,M), VEC(M)
C
        NUM = MIN0(M,7)
        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 J = 1, M
                 DATA(K,I) = DATA(K,I) + VEC(J) *
     X                                   EVEC(J,M-I+1)
  100         CONTINUE
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                         C
C  Determine projections of column-points on 7 principal components.      C
C                                                                         C
C-------------------------------------------------------------------------C
        SUBROUTINE PROJYP(M,EVALS,A,Z,VEC)
        INTEGER M,NUM,L,J1,J2,J3
        REAL    EVALS(M), A(M,M), Z(M,M), VEC(M)
C
        NUM = MIN0(M,7)
        DO 300 J1 = 1, M
           DO 50 L = 1, M
              VEC(L) = A(J1,L)
   50      CONTINUE
           DO 200 J2 = 1, NUM
              A(J1,J2) = 0.0
              DO 100 J3 = 1, M
                 A(J1,J2) = A(J1,J2) + VEC(J3)*Z(J3,M-J2+1)
  100         CONTINUE
              IF (EVALS(M-J2+1).GT.0.00005) A(J1,J2) =
     X                                  A(J1,J2)/SQRT(EVALS(M-J2+1))
              IF (EVALS(M-J2+1).LE.0.00005) A(J1,J2) = 0.0
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END







