C @(#)ca.for	17.1.1.1 (ES0-DMD) 01/25/02 17:16:37
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 CORRESPONDENCE ANALYSIS.                                      C
C                                                                            C
C                                                                            C
C  To call:   CALL CA(N,M,DATA,IPRINT,A1,W1,W2,A2,R,C,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 factors.                            C
C  IPRINT: print options.                                                    C
C          = 1: no printed output - arrays/vectors, only, contain items      C
C              calculated.                                                   C
C          = 2: printed output, in addition, of eigenvalues and eigenvectors.C
C          = 3: full printing of items calculated.                           C
C  A1    : matrix to be diagonalized, of dimensions M * M.                   C
C          On output, A1 contains in the first 7 columns the projections of  C
C          the column-points on the first 7 factors.                         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  R, C  : vectors of length N and M (respectively) containing, on output,   C
C          the masses of the row and column points.                          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, W1, R and C.                 C
C  All printed outputting is carried out in easily recognizable subroutines  C
C  called from the first subroutine which follows.                           C
C                                                                            C
C  Note the presence of the trivial first eigenvalue, and accompanying       C
C  eigenvector (all values 1, ignoring the sign).                            C
C                                                                            C
C  Errors: IERR = 1, no convergence in eigenextraction routine; IERR = 2,    C
C  values of a row or a column (hence, mass of the row or column) are all    C
C  zero.
C                                                                            C
C                                                                            C
C  F. Murtagh, ST-ECF/ESA/ESO, Garching-bei-Muenchen, January 1986.          C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE CA(N,M,DATA,IPRINT,A,W,FV1,Z,R,C,IERR)
        REAL    DATA(N,M), A(M,M), W(M), FV1(M), Z(M,M), R(N), C(M)
        REAL    TOT
        INTEGER N,M,IPRINT,IERR,I,J,J1,J2,M2
C
C          Form row sums and total
C
        TOT = 0.0
        DO 200 I = 1, N
           R(I) = 0.0
           DO 100 J = 1, M
              TOT = TOT + DATA(I,J)
              R(I) = R(I) + DATA(I,J)
  100      CONTINUE
  200   CONTINUE
C
C          Form column sums and then means
C
        DO 400 J = 1, M
           C(J) = 0.0
           DO 300 I = 1, N
              C(J) = C(J) + DATA(I,J)
  300      CONTINUE
           IF (C(J).GT.0.0) GOTO 350
C             WRITE (6,*) ' ALL VALUES IN COL.',J,' ARE 0: PLS. REMOVE'
              IERR = 2
              RETURN
  350      CONTINUE
           C(J) = C(J)/TOT
  400   CONTINUE
C
C          Form row means and make data into frequencies
C
        DO 600 I = 1, N
           IF (R(I).GT.0.0) GOTO 450
C             WRITE (6,*) ' ALL VALUES IN ROW',I,' ARE 0: PLS. REMOVE'
              IERR = 2
              RETURN
  450      CONTINUE
           R(I) = R(I)/TOT
           DO 500 J = 1, M
              DATA(I,J) = DATA(I,J)/TOT
  500      CONTINUE
  600   CONTINUE
C
C          Form matrix to be diagonalized
C
        DO 900 J1 = 1, M
           DO 800 J2 = 1, M
              A(J1,J2) = 0.0
              DO 700 I = 1, N
                 A(J1,J2) = A(J1,J2) + DATA(I,J1)*DATA(I,J2)/
     X                                (R(I)*SQRT(C(J1)*C(J2)))
  700         CONTINUE
  800      CONTINUE
  900   CONTINUE
C
C          Carry out eigenreduction.
C
        M2 = M
        CALL TRED2C(M,M2,A,W,FV1,Z)
        CALL TQL2C(M,M2,W,FV1,Z,IERR)
        IF (IERR.NE.0) GOTO 9000
C
C          Output eigenvalues and eigenvectors.
C
        IF (IPRINT.GT.1) CALL OUTVC(M,W)
C
C          Determine factors from the eigenvectors of the symmetric
C          matrix which has been diagonalozed.
C
        DO 1100 J1 = 1, M
           DO 1000 J2 = 1, M
              Z(J1,J2) = Z(J1,J2)/SQRT(C(J1))
 1000      CONTINUE
 1100   CONTINUE
C
        IF (IPRINT.GT.1) CALL OUTVC(M,Z)
C
C          Determine projections and contributions and output them.
C
        CALL PROJXC(N,M,DATA,Z,FV1,R)
        IF (IPRINT.EQ.3) CALL OUTPXC(N,M,DATA)
        CALL PROJYC(M,W,A,Z,FV1,C)
        IF (IPRINT.EQ.3) CALL OUTPYC(M,A)
        IF (IPRINT.EQ.3) CALL OUTCTX(N,M,DATA,R)
        IF (IPRINT.EQ.3) CALL OUTCTY(M,A,C)
C
 9000   RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                            C
C Reduce a real, symmetric matrix to a symmetric, tridiagonal                C
C 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. 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 producting 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 TRED2C(NM,N,A,D,E,Z)
C
        REAL A(NM,N),D(N),E(N),Z(NM,N)
        REAL SCALE,H,F,G,HH
        INTEGER NM,N,I,J,II,L,K,JP1
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-ECFESA/ESO, Garching-bei-Muenchen, January 1986.            C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE TQL2C(NM,N,D,E,Z,IERR)
C
        REAL    D(N), E(N), Z(NM,N)
        REAL    EPS,F,H,B,G,R,P,C,S
        INTEGER NM,N,IERR,I,J,L,M,L1,MML,II
        INTEGER K
        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 bottom
C             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 OUTMTC(N,M,ARRAY)
C
        INTEGER   N,M,K1,K2
        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 eigenvalues in order of decreasing value.                          C
C  The first eigenvalue in Correspondence Analysis is trivially              C
C  of value 1, and so is ignored.                                            C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTVLC(NVALS,VALS)
        INTEGER   NVALS,K
        REAL       VALS(NVALS),TOT,CUM,VPC,VCPC
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                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTVC(NDIM,VECS)
        INTEGER         NDIM,K1,K2,NUM
        REAL       VECS(NDIM,NDIM)
C
        NUM = MIN0(NDIM,7)
        WRITE (6,1000)
        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
 1000   FORMAT(1H0,'FACTORS 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 (up to) 7 factors.              C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTPXC(N,M,PRJN)
        INTEGER N,M,K,NUM,J
        REAL    PRJN(N,M)
C
        NUM = MIN0(M,7)
        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(1H0,'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 column points on first 7 factors.                  C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTPYC(M,PRJNS)
        INTEGER K,M,J,NUM
        REAL    PRJNS(M,M)
C
        NUM = MIN0(M,7)
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, M
           WRITE (6,1030) K,(PRJNS(K,J),J=1,NUM)
  100   CONTINUE
C
 1000   FORMAT(1H0,'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 factors.                       C
C   (The first - trivial - factor is ignored.)                               C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE PROJXC(N,M,DATA,EVEC,VEC,R)
        INTEGER N,M,K,I,J,L,NUM
        REAL    DATA(N,M), EVEC(M,M), VEC(M), R(N)
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) * EVEC(J,M-I)
  100         CONTINUE
              IF (R(K).GT.0.0) DATA(K,I) = DATA(K,I)/R(K)
              IF (R(K).EQ.0.0) DATA(K,I) = 0.0
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                            C
C   Form projections of column points on first 7 factors.                    C
C   (The first - trivial - factor is ignored.)                               C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE PROJYC(M,EVALS,A,Z,VEC,C)
        INTEGER M,L,J1,J2,J3,NUM
        REAL    EVALS(M), A(M,M), Z(M,M), VEC(M), C(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)*
     X                                   SQRT(C(J3))
  100         CONTINUE
              IF ((EVALS(M-J2).GT.0.0).AND.(C(J1).GT.0.0))
     X                  A(J1,J2) = A(J1,J2)/SQRT(EVALS(M-J2)*C(J1))
              IF ((EVALS(M-J2).EQ.0.0).OR.(C(J1).EQ.0.0)) A(J1,J2)=0.0
  200      CONTINUE
  300   CONTINUE
C
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                            C
C   Output contributions of row-points to first 7 factors.                   C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTCTX(N,M,PRJN,R)
        INTEGER N,M,NUM,K,J
        REAL    PRJN(N,M), R(N), XDUM(7)
C
        NUM = MIN0(M,7)
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, N
           DO 50  J = 1, NUM
              XDUM(J) = R(K)*PRJN(K,J)**2
   50      CONTINUE
           WRITE (6,1030) K,(XDUM(J),J=1,NUM)
  100   CONTINUE
C
 1000   FORMAT(1H0,'ABS. CONTRIBUTIONS OF OBJECTS TO FACTORS FOLLOW.',/)
 1010   FORMAT(' OBJECT  CNTR-1  CNTR-2  CNTR-3  CNTR-4  CNTR-5  CNTR-6
     X  CNTR-7')
 1020   FORMAT(' ------  ------  ------  ------  ------  ------  ------
     X  ------')
 1030   FORMAT(I5,2X,7F8.4)
        RETURN
        END
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++C
C                                                                            C
C   Output contributions of column-points to first 7 factors.                C
C                                                                            C
C----------------------------------------------------------------------------C
        SUBROUTINE OUTCTY(M,PRJNS,C)
        INTEGER M,K,J,NUM
        REAL    PRJNS(M,M), C(M), XDUM(7)
C
        NUM = MIN0(M,7)
        WRITE (6,1000)
        WRITE (6,1010)
        WRITE (6,1020)
        DO 100 K = 1, M
           DO 50 J = 1, NUM
              XDUM(J) = C(K)*PRJNS(K,J)**2
   50      CONTINUE
           WRITE (6,1030) K,(XDUM(J),J=1,NUM)
  100   CONTINUE
C
 1000   FORMAT
     X  (1H0,'ABS. CONTRIBUTIONS OF VARIABLES TO FACTORS FOLLOW.',/)
 1010   FORMAT('  VBLE.  CNTR-1  CNTR-2  CNTR-3  CNTR-4  CNTR-5  CNTR-6
     X  CNTR-7')
 1020   FORMAT(' ------  ------  ------  ------  ------  ------  ------
     X  ------')
 1030   FORMAT(I5,2X,7F8.4)
        RETURN
        END