C polyfac.f
C The authors of this software are Jih-Jie Chang and J D Carroll.

C Copyright (c) 1993 by AT&T.
C Permission to use, copy, modify, and distribute this software for any
C purpose without fee is hereby granted, provided that this entire
C notice is included in all copies of any software which is or includes
C a copy or modification of this software and in all copies of the
C supporting documentation for such software.
C THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMP-
C LIED WARRANTY.  IN PARTICULAR, NEITHER THE AUTHORS NOR AT&T MAKE ANY
C REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY
C OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE.

C This software comes from the FIRST MDS Package of AT&T Bell Laboratories

C For explanation of the method this software implements, see
C Carroll,J.D. (1969), " Polynomial factor analysis" in Proceedings of
C the 77th Annual Convention of the American Psychological Association,
C 4,103-104.
C----+----@----+----@----+----@----+----@----+----@----+----@----+----@

*PFAC                   POLYNOMIAL FACTOR ANALYSIS
*
*     INPUT
*
*     1) N,NF,M,MAXIT,NOPTY,NV,KPUNX,ID,KAL,NRAND,CRIT,CALF(9I4,I6,2F4.0)
*          N = NUMBER OF OBSERVATIONS.
*          NF = NUMBER OF FACTORS.
*          M = TOTAL NUMBER OF POLYNOMIAL COMPONENTS.
*          MAXIT = MAXIMUM NUMBER OF ITERATIONS.
*          NOPTY = CONTROL PARAMETER DEFINING INPUT MATRIX.
*                  IF=0, SQUARE SYMMETRIC MATRIX C IS INPUT.
*                    =1, NV X N MATRIX Y IS INPUT.
*                    C=Y(TRANSPOSED)*Y, IS COMPUTED.
*                    =2, Y IS INPUT, ROW MEANS ARE SUBTRACTED,
*                    C=Y(TRANSPOSED)*Y.
*                    =3, Y IS INPUT, COLUMN MEANS ARE SUBTRACTED AND
*                    COLUMN VARIANCES STANDARDIZED, C=Y(TRANSPOSED)*Y.
*                    =4, Y IS INPUT, COLUMN MEANS ARE SUBTRACTED,
*                    C=Y(TRANSPOSED)*Y.
*                    =5, Y IS INPUT, ROW MEANS ARE SUBTRACTED AND ROW
*                    VARIANCES STANDARDIZED, C=Y(TRANSPOSED)*Y.
*          NV = NUMBER OF VARIABLES (IF NOPTY IS GREATER THAN 0,
*          OTHERWISE NOT RELEVANT).  CANNOT BE GREATER THAN N.
*          KPUNX = PUNCH OPTION FOR MATRIX X.
*                  IF=0, MATRIX X IS PUNCHED.
*                    =1, MATRIX X IS NOT PUNCHED.
*          KAL = OPTION FOR ALPHA
*                IF=0, ALPHA=CALF*(4.0**(COS(IT)**3))/SQRT(GRL(IT)).
*                IF=1, ALPHA=(A-B)/C.
*          ID = IDENTIFICATION OPTION.
*                  IF=0, NUMBERS ARE PLOTTED.
*                  IF=1, READ IN IDENTIFICATION.
*          NRAND = RANDOM NUMBER OPTION.
*                  IF=0, READ IN X MATRIX OF COORDINATES.
*          OTHERWISE = SIX DIGIT ODD INTEGER, USED TO GENERATE
*                     RANDOM INITIAL X MATRIX.
*          CRIT = CRITERION FOR STOPPING.
*                  IF=0, PROGRAM ASSIGNS .01 TO CRIT.
*          CALF = FLOATING POINT CONSTANT, USED TO CALCULATE ALPHA.
*                 PRESENT ONLY IF KAL = 0 OR BLANK.
*
*     2) TITLE          (12A6)
*
*     3) FORMAT FOR MATRIX C.          (12A6)
*
*     4) MATRIX C
*          N X N MATRIX OF CROSS PRODUCTS.
*          Y(TRANSPOSED)*Y (INPUT ONLY IF NOPTY=0, COMPUTED OTHERWISE).
*
*     5) FORMAT FOR MATRIX Y.          (12A6)
*
*     6) MATRIX Y.
*          NV X N MATRIX OF "SCORES" (RELEVANT ONLY IF NOPTY IS GREATER
*          THAN 0).
*
*     7) FORMAT FOR MATRIX X.          (12A6)
*
*     8) MATRIX X.
*          NF X N MATRIX OF "FACTOR LOADINGS".
*         (INPUT ONLY IF NRAND = 0, OR BLANK, OTHERWISE RANDOM STARTING
*         CONFIGURATION IS GENERATED.)
*
*     9) FORMAT FOR MATRIX NE.          (12A6)
*
*     10) MATRIX NE
*          M X NF MATRIX OF EXPONENTS FOR M POLYNOMIAL COMPONENTS WITH
*          RESPECT TO NF FACTORS.
*
*     11) FORMAT FOR IDENTIFICATION.
*
*     12) IDENTIFICATION ARRAY.
*         (INPUT ONLY IF ID = 1, OTHERWISE, NUMBERS ARE PLOTTED).
*
*     REPEAT NUMBERS 1 THRU 12 FOR EACH ADDITIONAL SET OF DATA.
*
*     FINAL CARD MUST BE BLANK.
*
*
      DIMENSION TITLE(13),FMTC(12),FMTX(12),FMTNE(12)
      DIMENSION VAF(101),GRL(101),ALF(101),COS(101),NBACK(101),X(5,100),
     1DX(5,100),Z(20,100),DZ(20,100),NE(20,5),C(100,100),R(20,20),U(20,1
     200),Y(50,100)
      DATA  TITLE (1)/6H75H   /
       CALL  DIMENS (R,U)
 111  READ 501,N,NF,M,MAXIT,NOPTY,NV,KPUNX,ID,KAL,NRAND,CRIT,CALF
 501  FORMAT (9I4,I6,2F4.0)
      IF(N.EQ.0) GO TO 999
      IF(NRAND.NE.0) PRINT 600, NRAND
 600  FORMAT (9H1NRAND = ,I6)
      IF(CRIT.EQ.0.) CRIT=.01
      READ 504,(TITLE(K),K=2,13)
 504  FORMAT (12A6)
      READ 504,(FMTC(I),I=1,12)
      IF(NOPTY.GT.0) GO TO 503
*
*     READ MATRIX C.  (IF NOPTY=0).
*
      DO 510 I=1,N
 510  READ FMTC,(C(I,J),J=1,N)
      GO TO 502
*
*     READ MATRIX Y.  (IF NOPTY IS GREATER THAN 0).
*
 503  DO 520 I=1,NV
 520  READ FMTC,(C(I,J),J=1,N)
      PRINT 705
 705  FORMAT (17H0MATRIX Y (INPUT))
      DO 710 I=1,NV
 710  PRINT 521,I,(C(I,J),J=1,N)
 502  IF(NRAND.EQ.0) GO TO 505
*
*     GENERATE RANDOM INITIAL X MATRIX.
*
       CALL RANUM (X,N,NF,NRAND)
      DO 95 J=1,N
      SSQ=0.0
      DO 90 I=1,NF
 90   SSQ = SSQ+X(I,J)**2
 95   X(I,J) = ABS(X(I,J)/SQRT(SSQ))
      XN = N
      DO 80 I=1,NF
      SUM = 0.
      DO 85 J=1,N
 85   SUM = SUM + X(I,J)
      XMEAN = SUM/XN
      DO 80 J=1,N
 80   X(I,J) = X(I,J) - XMEAN
      GO TO 506
*
*     READ MATRIX X.
*
 505  READ 504,(FMTX(I),I=1,12)
      DO 530 I=1,NF
 530  READ FMTX,(X(I,J),J=1,N)
*
*     READ MATRIX NE.
*
 506  READ 504,(FMTNE(K),K=1,12)
      DO 540 I=1,M
 540  READ FMTNE,(NE(I,J),J=1,NF)
      PRINT 543
 543  FORMAT (19H0MATRIX NE (M X NF))
      DO 545 I=1,M
 545  PRINT 544,I,(NE(I,J),J=1,NF)
 544  FORMAT (1H0,I2,3X,5I10)
*
*     NORMALIZE MATRIX Y.
*
      IF(NOPTY.EQ.0) GO TO 1100
      IF(NOPTY.EQ.1) GO TO 1000
      IF(NOPTY.EQ.3) GO TO 900
      DO 820 J=1,NV
      YMEAN=0.0
      DO 810 K=1,N
 810  YMEAN=YMEAN+C(J,K)
      YMEAN=YMEAN/FLOAT(N)
      DO 820 K=1,N
 820  C(J,K)=C(J,K)-YMEAN
      IF(NOPTY.EQ.2) GO TO 1000
      DO 830 J=1,NV
      CNORM=0.0
      DO 825 K=1,N
 825  CNORM=CNORM+C(J,K)**2
      CNORM=1./SQRT(CNORM)
      DO 830 K=1,N
 830  C(J,K)=CNORM*C(J,K)
      GO TO 1000
 900  DO 920 K=1,N
      YMEAN=0.0
      DO 910 J=1,NV
 910  YMEAN=YMEAN+C(J,K)
      YMEAN=YMEAN/FLOAT(NV)
      DO 920 J=1,NV
 920  C(J,K)=C(J,K)-YMEAN
      IF(NOPTY.EQ.4) GO TO 1000
      DO 940 K=1,N
      CNORM=0.0
      DO 930 J=1,NV
 930  CNORM=CNORM+C(J,K)**2
      CNORM=1./SQRT(CNORM)
      DO 940 J=1,NV
 940  C(J,K)=CNORM*C(J,K)
 1000 PRINT 951
 951  FORMAT (1H0,19HNORMALIZED Y MATRIX)
      DO 952I=1,NV
 952  PRINT 521,I,(C(I,J),J=1,N)
      IF(NV.LE.50) NNV=NV
      IF(NV.GT.50) NNV=50
      DO 950 K=1,N
      DO 950 J=1,NNV
 950  Y(J,K)=C(J,K)
*
*     COMPUTE MATRIX C.
*
      DO 1002 K=1,N
      DO 1001 KK=K,N
      VAF(KK)=0.0
      DO 1001 J=1,NV
 1001 VAF(KK)=VAF(KK)+C(J,K)*C(J,KK)
      DO 1002 KK=K,N
 1002 C(KK,K)=VAF(KK)
      DO 1003 K=1,N
      DO 1003 KK=K,N
 1003 C(K,KK)=C(KK,K)
 1100 CDSUM=0.0
      DO 1101 K=1,N
 1101 CDSUM=CDSUM+C(K,K)
      DO 1111 K=1,N
      DO 1111 KK=1,N
 1111 C(K,KK)=C(K,KK)/CDSUM
      PRINT 715
 715  FORMAT (9H0MATRIX C)
      DO 720 I=1,N
 720  PRINT 521,I,(C(I,J),J=1,N)
      ENF=NF
      ALPHA=0.0
      MAXIT1=MAXIT+1
      DO 500 IT=1,MAXIT1
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9999,TIME                                                   A  *****
 9999 FORMAT(14H0A     TIME = ,F8.5)
      VAF(IT)=0.0
      GRL(IT)=0.0
      ALF(IT)=0.0
      COS(IT)=0.0
      NBACK(IT)=0
      ITM1=IT-1
*
*     TAKE STEP OF SIZE ALPHA, AND NORMALIZE X MATRIX.
*
 1    DO 6 J=1,NF
      CNORM=0.0
      DO 5 K=1,N
      X(J,K) = X(J,K)+ALPHA*DX(J,K)
 5    CNORM = CNORM+X(J,K)**2
      CNORM = 1./SQRT(CNORM*ENF)
      DO 6 K=1,N
 6    X(J,K) = CNORM*X(J,K)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9998,TIME                                                   B  *****
 9998 FORMAT(14H0B     TIME = ,F8.5)
      PRINT 430
 430  FORMAT (9H0X MATRIX)
      DO 435 J=1,NF
 435  PRINT 521,J,(X(J,K),K=1,N)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9997,TIME                                                   C  *****
 9997 FORMAT(14H0C     TIME = ,F8.5)
*
*     DEFINE Z
*
      DO 10 L=1,M
      DO 10 K=1,N
      Z(L,K) = 1.0
      DO 10 J=1,NF
      IF(NE(L,J).EQ.0) GO TO 10
      Z(L,K) = Z(L,K)*X(J,K)**NE(L,J)
 10   CONTINUE
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9996,TIME                                                   D  *****
 9996 FORMAT(14H0D     TIME = ,F8.5)
*
*     COMPUTE Z*Z(TRANSPOSED)=R
*
      DO 20 L=1,M
      DO 20 LL=1,L
      R(L,LL)=0.0
      DO 15 K=1,N
 15   R(L,LL) = R(L,LL)+Z(L,K)*Z(LL,K)
 20   R(LL,L)=R(L,LL)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9995,TIME                                                   E  *****
 9995 FORMAT(14H0E     TIME = ,F8.5)
*
*     INVERT R=Z*Z(TRANSPOSED)
*
      MM=1
       CALL  MATINV (R,M,MM,U)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9994,TIME                                                   F  *****
 9994 FORMAT(14H0F     TIME = ,F8.5)
      IF(MM.GT.1) GO TO 888
*
*     COMPUTE U=(Z*Z(TRANSPOSED)) INVERSED*Z
*
      DO 30 L=1,M
      DO 30 K=1,N
      U(L,K)=0.0
      DO 30 LL=1,M
 30   U(L,K) =  U(L,K) + R(L,LL)*Z(LL,K)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9993,TIME                                                   G  *****
 9993 FORMAT(14H0G     TIME = ,F8.5)
*
*     COMPUTE  DZ = 2.0((UC*U(TRANSPOSED))Z-UC)
*
      DO 40 L=1,M
      DO 40 J=1,N
      DZ(L,J)=0.0
      DO 40 I=1,N
 40   DZ(L,J) = DZ(L,J)+U(L,I)*C(I,J)
      DO 50 L=1,M
      DO 50 LL=L,M
      R(L,LL) = 0.0
      DO 45 J=1,N
 45    R(L,LL) = R(L,LL)+DZ(L,J)*U(LL,J)
 50   R(LL,L) = R(L,LL)
      DO 100 K=1,N
      DO 100 L=1,M
      A=0.0
      DO 60 I=1,M
 60   A = A+R(L,I)*Z(I,K)
 100  DZ(L,K) = 2.0*(A-DZ(L,K))
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9992,TIME                                                   H  *****
 9992 FORMAT(14H0H     TIME = ,F8.5)
*
*     COMPUTE TRACE  U*C*Z(TRANSPOSED)=VAF(IT)
*
      VAF(IT)=0.0
      DO 35 L=1,M
      DO 35 K=1,N
      DO 35 KK=1,N
 35   VAF(IT) = VAF(IT)+U(L,K)*C(K,KK)*Z(L,KK)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9991,TIME                                                   I  *****
 9991 FORMAT(14H0I     TIME = ,F8.5)
      PRINT 470,IT,VAF(IT)
 470  FORMAT (5H0VAF(,I3,4H) = ,F10.5)
*
 38   ALF(IT) = ALPHA
      PRINT 475,ITM1,NBACK(ITM1)
 475  FORMAT (7H0NBACK(,I2,4H) = ,I3)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9990,TIME                                                   J  *****
 9990 FORMAT(14H0J     TIME = ,F8.5)
*
*     COMPUTE DX
*
 39   CRPR=0.0
      GR=0.0
      DO 200 J=1,NF
      DO 200 K=1,N
      DXC=0.0
      DO 170 L=1,M
 170  DXC = DXC+DZ(L,K)*FLOAT(NE(L,J))*Z(L,K)/X(J,K)
      CRPR = CRPR+DXC*DX(J,K)
      GR = GR+DXC**2
 200  DX(J,K)=DXC
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9989,TIME                                                   K  *****
 9989 FORMAT(14H0K     TIME = ,F8.5)
      GRL(IT) = SQRT(GR)
      PRINT 485,IT,GRL(IT)
 485  FORMAT (5H0GRL(,I3,4H) = ,F10.5)
      IF(IT.EQ.1) GO TO 299
      COS(IT)=CRPR/(GRL(ITM1)*GRL(IT))
 299  IF(GRL(IT).LE.CRIT) GO TO 777
      PRINT 490
 490  FORMAT (10H0DX MATRIX)
      DO 495 J=1,NF
 495  PRINT 521,J,(DX(J,K),K=1,N)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9988,TIME                                                   L  *****
 9988 FORMAT(14H0L     TIME = ,F8.5)
*
*     BEGIN COMPUTATION OF "STEP-SIZE" ALPHA.
*
*     COMPUTE DELTA Z (CALLED DZ).
*
      IF (KAL.EQ.0) GO TO 265
      IF(IT.EQ.1) GO TO 263
      ALPHA = ALPHA*ABS(GRL(ITM1)/(GRL(ITM1)-GRL(IT)*COS(IT)))
      GO TO 500
 263  ALPHA = CALF
      GO TO 500
 265  DO 300 L=1,M
      DO 300 K=1,N
      DDZ=0.0
      DO 270 J=1,NF
 270  DDZ = DDZ+FLOAT(NE(L,J))*DX(J,K)/X(J,K)
 300  DZ(L,K) = Z(L,K)*DDZ
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9987,TIME                                                   M  *****
 9987 FORMAT(14H0M     TIME = ,F8.5)
*
*     COMPUTE
*     A = 1/2 TRACE(Z*DZ(TRANSPOSED)+DZ*Z(TRANSPOSED))
*     *(U*C*U(TRANSPOSED)),
*     B = TRACE DZ*C*U(TRANSPOSED), AND
*     CC = TRACE (DZ*DZ(TRANSPOSED))*(U*C*U(TRANSPOSED)).
*     ALPHA = (A-B)/C
*
      A=0.0
      B=0.0
      CC=0.0
      DO 400 L=1,M
      DO 400 LL=1,L
      C1=0.0
      C2=0.0
      C3=0.0
      C4=0.0
      DO 350 K=1,N
      C1 = C1+Z(L,K)*DZ(LL,K)+DZ(L,K)*Z(LL,K)
      C2 = C2+DZ(L,K)*DZ(LL,K)
      DO 350 KK=1,N
      C3 = C3+U(L,K)*C(K,KK)*U(LL,KK)
 350  C4 = C4+DZ(L,K)*C(K,KK)*U(LL,KK)
      IF(L.EQ.LL) GO TO 380
      A=A+C1*C3
      CC=CC+2.0*C2*C3
      GO TO 400
 380  A=A+C1*C3/2.0
      B=B+C4
      CC=CC+C2*C3
 400  CONTINUE
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9986,TIME                                                   N  *****
 9986 FORMAT(14H0N     TIME = ,F8.5)
      ALPHA=(B-A)/CC
 500  CONTINUE
 5001 IF(NOPTY.EQ.0) GO TO 606
*
*     COMPUTE BHAT TRANSPOSED = U * Y (TRANSPOSED).
*
      DO 601 J=1,M
      DO 601 K=1,NNV
      DZ(J,K)=0.0
      DO 601 L=1,N
 601  DZ(J,K) = DZ(J,K)+U(J,L)*Y(K,L)
      PRINT 603
 603  FORMAT (1H0,58HMATRIX BHAT TRANSPOSED (MATRIX OF REGRESSION COEFFI
     1CIENTS))
      DO 610 I=1,M
 610  PRINT 521,I,(DZ(I,J),J=1,NNV)
*
*     COMPUTE YHAT = BHAT * Z.
*
      DO 602 J=1,NNV
      DO 602 K=1,N
      Y(J,K)=0.0
      DO 602 L=1,M
 602  Y(J,K) = Y(J,K)+DZ(L,J)*Z(L,K)
      PRINT 604
 604  FORMAT (1H0,11HMATRIX YHAT)
      DO 605 I=1,NNV
 605  PRINT 521,I,(Y(I,J),J=1,N)
 606  TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9985,TIME                                                   O  *****
 9985 FORMAT(14H0O     TIME = ,F8.5)
*
*     PLOT "FINAL CONFIGURATION."
*
       CALL PLOT (ID,N,NF,TITLE,X)
*
*     OUTPUT AND FINIS.
*
 777  PRINT 519,(TITLE(K),K=2,13)
 519  FORMAT (1H1,12A6)
      PRINT 511,ITM1
 511  FORMAT (14H0 ITERATION = ,I3)
      PRINT 512
 512  FORMAT (11H0 ITERATION,7X,3HVAF,7X,3HGRL,7X,3HCOS,7X,3HALF)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9984,TIME                                                   P  *****
 9984 FORMAT(14H0P     TIME = ,F8.5)
      DO 550 K=1,IT
      KK=K-1
 550  PRINT 513,KK,VAF(K),GRL(K),COS(K),ALF(K)
 513  FORMAT (1H ,4X,I3,3X,4F10.5)
      TIME = ITIMEZ(0)
      TIME = TIME/10.                                                       0000
      PRINT 9983,TIME                                                   Q  *****
 9983 FORMAT(14H0Q     TIME = ,F8.5)
      PRINT 514
 514  FORMAT (56H0 MATRIX X (NUMBER OF FACTORS BY NUMBER OF OBSERVATIONS
     1))
      DO 560 I=1,NF
 560  PRINT 521,I,(X(I,J),J=1,N)
      PRINT 516
 516  FORMAT (57H0 MATRIX DX (NUMBER OF FACTORS BY NUMBER OF OBSERVATION
     1S))
      DO 570 I=1,NF
 570  PRINT 521,I,(DX(I,J),J=1,N)
*
 521  FORMAT (1H ,I3,10E12.4/(1H ,3X,10E12.4))
*
      PRINT 517
 517  FORMAT (59H0 MATRIX Z (NUMBER OF COMPONENTS BY NUMBER OF OBSERVATI
     1ONS))
      DO 580 I=1,M
 580  PRINT 521,I,(Z(I,J),J=1,N)
      IF(KPUNX.EQ.1) GO TO 999
      PUNCH 504,(TITLE(K),K=2,13)
      PUNCH 522
 522  FORMAT (11H(2X,5E13.5))
      DO 590 I=1,NF
 590  PUNCH 518,I,(X(I,J),J=1,N)
 518  FORMAT (I2,5E13.5/(2X,5E13.5))
      GO TO 111
 888  PRINT 889
 889  FORMAT (29H0 R REMAINS UNINVERTED. STOP.)
 999  STOP
      END