Power-Programmierung/CD1.mdf/fortran/library/ssp/statcorr/canor.for

Text File | 1985-11-29 | 6KB | 213 lines

C
C     ..................................................................
C
C        SUBROUTINE CANOR
C
C        PURPOSE
C           COMPUTE THE CANONICAL CORRELATIONS BETWEEN TWO SETS OF
C           VARIABLES.  CANOR IS NORMALLY PRECEDED BY A CALL TO SUBROU-
C           TINE CORRE.
C
C        USAGE
C           CALL CANOR (N,MP,MQ,RR,ROOTS,WLAM,CANR,CHISQ,NDF,COEFR,
C                       COEFL,R)
C
C        DESCRIPTION OF PARAMETERS
C           N     - NUMBER OF OBSERVATIONS
C           MP    - NUMBER OF LEFT HAND VARIABLES
C           MQ    - NUMBER OF RIGHT HAND VARIABLES
C           RR    - INPUT MATRIX (ONLY UPPER TRIANGULAR PORTION OF THE
C                   SYMMETRIC MATRIX OF M X M, WHERE M = MP + MQ)
C                   CONTAINING CORRELATION COEFFICIENTS.  (STORAGE MODE
C                   OF 1)
C           ROOTS - OUTPUT VECTOR OF LENGTH MQ CONTAINING EIGENVALUES
C                   COMPUTED IN THE NROOT SUBROUTINE.
C           WLAM  - OUTPUT VECTOR OF LENGTH MQ CONTAINING LAMBDA.
C           CANR  - OUTPUT VECTOR OF LENGTH MQ CONTAINING CANONICAL
C                   CORRELATIONS.
C           CHISQ - OUTPUT VECTOR OF LENGTH MQ CONTAINING THE
C                   VALUES OF CHI-SQUARES.
C           NDF   - OUTPUT VECTOR OF LENGTH MQ CONTAINING THE DEGREES
C                   OF FREEDOM ASSOCIATED WITH CHI-SQUARES.
C           COEFR - OUTPUT MATRIX (MQ X MQ) CONTAINING MQ SETS OF
C                   RIGHT HAND COEFFICIENTS COLUMNWISE.
C           COEFL - OUTPUT MATRIX (MP X MQ) CONTAINING MQ SETS OF
C                   LEFT HAND COEFFICIENTS COLUMNWISE.
C           R     - WORK MATRIX (M X M)
C
C        REMARKS
C           THE NUMBER OF LEFT HAND VARIABLES (MP) SHOULD BE GREATER
C           THAN OR EQUAL TO THE NUMBER OF RIGHT HAND VARIABLES (MQ).
C           THE VALUES OF CANONICAL CORRELATION, LAMBDA, CHI-SQUARE,
C           DEGREES OF FREEDOM, AND CANONICAL COEFFICIENTS ARE COMPUTED
C           ONLY FOR THOSE EIGENVALUES IN ROOTS WHICH ARE GREATER THAN
C           ZERO.
C
C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
C           MINV
C           NROOT  (WHICH, IN TURN, CALLS THE SUBROUTINE EIGEN.)
C
C        METHOD
C           REFER TO W. W. COOLEY AND P. R. LOHNES, 'MULTIVARIATE PRO-
C           CEDURES FOR THE BEHAVIORAL SCIENCES', JOHN WILEY AND SONS,
C           1962, CHAPTER 3.
C
C     ..................................................................
C
      SUBROUTINE CANOR (N,MP,MQ,RR,ROOTS,WLAM,CANR,CHISQ,NDF,COEFR,
     1                  COEFL,R)
      DIMENSION RR(1),ROOTS(1),WLAM(1),CANR(1),CHISQ(1),NDF(1),COEFR(1),
     1          COEFL(1),R(1)
C
C        ...............................................................
C
C        IF A DOUBLE PRECISION VERSION OF THIS ROUTINE IS DESIRED, THE
C        C IN COLUMN 1 SHOULD BE REMOVED FROM THE DOUBLE PRECISION
C        STATEMENT WHICH FOLLOWS.
C
C     DOUBLE PRECISION RR,ROOTS,WLAM,CANR,CHISQ,COEFR,COEFL,R,DET,SUM
C
C        THE C MUST ALSO BE REMOVED FROM DOUBLE PRECISION STATEMENTS
C        APPEARING IN OTHER ROUTINES USED IN CONJUNCTION WITH THIS
C        ROUTINE.
C
C        THE DOUBLE PRECISION VERSION OF THIS SUBROUTINE MUST ALSO
C        CONTAIN DOUBLE PRECISION FORTRAN FUNCTIONS.  SQRT IN STATEMENT
C        165 MUST BE CHANGED TO DSQRT.  ALOG IN STATEMENT 175 MUST BE
C        CHANGED TO DLOG.
C
C        ...............................................................
C
C     PARTITION INTERCORRELATIONS AMONG LEFT HAND VARIABLES, BETWEEN
C     LEFT AND RIGHT HAND VARIABLES, AND AMONG RIGHT HAND VARIABLES.
C
      M=MP+MQ
      N1=0
      DO 105 I=1,M
      DO 105 J=1,M
      IF(I-J) 102, 103, 103
  102 L=I+(J*J-J)/2
      GO TO 104
  103 L=J+(I*I-I)/2
  104 N1=N1+1
  105 R(N1)=RR(L)
      L=MP
      DO 108 J=2,MP
      N1=M*(J-1)
      DO 108 I=1,MP
      L=L+1
      N1=N1+1
  108 R(L)=R(N1)
      N2=MP+1
      L=0
      DO 110 J=N2,M
      N1=M*(J-1)
      DO 110 I=1,MP
      L=L+1
      N1=N1+1
  110 COEFL(L)=R(N1)
      L=0
      DO 120 J=N2,M
      N1=M*(J-1)+MP
      DO 120 I=N2,M
      L=L+1
      N1=N1+1
  120 COEFR(L)=R(N1)
C
C     SOLVE THE CANONICAL EQUATION
C
      L=MP*MP+1
      K=L+MP
      CALL MINV (R,MP,DET,R(L),R(K))
C
C        CALCULATE T = INVERSE OF R11 * R12
C
      DO 140 I=1,MP
      N2=0
      DO 130 J=1,MQ
      N1=I-MP
      ROOTS(J)=0.0
      DO 130 K=1,MP
      N1=N1+MP
      N2=N2+1
  130 ROOTS(J)=ROOTS(J)+R(N1)*COEFL(N2)
      L=I-MP
      DO 140 J=1,MQ
      L=L+MP
  140 R(L)=ROOTS(J)
C
C        CALCULATE A = R21 * T
C
      L=MP*MQ
      N3=L+1
      DO 160 J=1,MQ
      N1=0
      DO 160 I=1,MQ
      N2=MP*(J-1)
      SUM=0.0
      DO 150 K=1,MP
      N1=N1+1
      N2=N2+1
  150 SUM=SUM+COEFL(N1)*R(N2)
      L=L+1
  160 R(L)=SUM
C
C        CALCULATE EIGENVALUES WITH ASSOCIATED EIGENVECTORS OF THE
C        INVERSE OF R22 * A
C
      L=L+1
      CALL NROOT (MQ,R(N3),COEFR,ROOTS,R(L))
C
C     FOR EACH VALUE OF I = 1, 2, ..., MQ, CALCULATE THE FOLLOWING
C     STATISTICS
C
      DO 210 I=1,MQ
C
C        TEST WHETHER EIGENVALUE IS GREATER THAN ZERO
C
      IF(ROOTS(I)) 220, 220, 165
C
C        CANONICAL CORRELATION
C
  165 CANR(I)= SQRT(ROOTS(I))
C
C        CHI-SQUARE
C
      WLAM(I)=1.0
      DO 170 J=I,MQ
  170 WLAM(I)=WLAM(I)*(1.0-ROOTS(J))
      FN=N
      FMP=MP
      FMQ=MQ
  175 CHISQ(I)=-(FN-0.5*(FMP+FMQ+1.0))*ALOG(WLAM(I))
C
C        DEGREES OF FREEDOM FOR CHI-SQUARE
C
      N1=I-1
      NDF(I)=(MP-N1)*(MQ-N1)
C
C        I-TH SET OF RIGHT HAND COEFFICIENTS
C
      N1=MQ*(I-1)
      N2=MQ*(I-1)+L-1
      DO 180 J=1,MQ
      N1=N1+1
      N2=N2+1
  180 COEFR(N1)=R(N2)
C
C        I-TH SET OF LEFT HAND COEFFICIENTS
C
      DO 200 J=1,MP
      N1=J-MP
      N2=MQ*(I-1)
      K=MP*(I-1)+J
      COEFL(K)=0.0
      DO 190 JJ=1,MQ
      N1=N1+MP
      N2=N2+1
  190 COEFL(K)=COEFL(K)+R(N1)*COEFR(N2)
  200 COEFL(K)=COEFL(K)/CANR(I)
  210 CONTINUE
  220 RETURN
      END