C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C fff.f     
C
C Copyright (C) 1981,1992.
C
C A. Ronald Gallant
C P.O. Box 5513 
C Raleigh NC 27650-5513 
C USA   
C
C Permission to use, copy, modify, and distribute this software and 
C its documentation for any purpose and without fee is hereby granted, 
C provided that the above copyright notice appear in all copies and 
C that both that copyright notice and this permission notice appear 
C in supporting documentation.  This software is provided "as is" 
C without any expressed or implied warranty.
C
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     FFFCGH      4/21/81
C
C     PURPOSE
C     GENERATE MULTI-INDEXES AND THE JACOBIAN OF A VECTOR THAT CONSISTS
C     OF THE FOURIER FLEXIBLE FORM, ITS GRADIENT, AND ITS HESSIAN.
C
C     USAGE
C     CALL FFFCGH(N,KC,IS,KA,IAA,JJA,M,DL,X,CGH,LT,IW)
C
C     ARGUMENTS
C     N   - LENGTH OF X, INPUT.  LENGTH MUST BE LESS THAN 11.
C           INTEGER*4
C     KC  - LARGEST NORM OF GENERATED MULTI-INDEXES, INPUT.
C           KC MUST BE LESS THAN 30.
C           INTEGER*4
C     IS  - SWITCH, INPUT.
C           IF IS = 0  THEN NO MULTI-INDEXES ARE DELETED.
C           IF IS = L  THEN MULTI-INDEXES WITH
C                      SUM(KA(I,J): I=1,2,...,L).NE.0
C                      ARE DELETED.
C           IF IS = -1 THEN KA AND IAA ARE NOT COMPUTED AND ARE PRESUMED
C                      TO BE AVAILABLE AS INPUT.
C           INTEGER*4
C     KA  - N BY IAA MATRIX CONTAINING COMPUTED MULTI-INDEXES AS
C           COLUMNS, OUTPUT.  INPUT IF IS=-1.  STORED COLUMNWISE
C           (STORAGE MODE 0)
C           INTEGER*4
C     IAA - NUMBER OF MULTI-INDEXES GENERATED, OUTPUT.  INPUT IF
C           IS=-1.
C           INTEGER*4
C     JJA - VECTOR OF LENGTH IAA, INPUT.  JJA(IA) DEFINES THE ORDER
C           OF THE FOURIER EXPANSION OF THE TERM ASSOCIATED WITH
C           MULTI-INDEX NUMBER IA.  ZERO VALUES ARE PERMITTED.  THEY
C           HAVE THE EFFECT OF DELETING MULTI-INDEX NUMBER IA WHEN
C           COMPUTING CGH.  A VECTOR OF ALL ZEROES WILL CAUSE THE
C           IDENTITY BORDERED ABOVE BY A COLUMN OF ONES AND BELOW
C           BY A BLOCK OF ZEROES TO BE RETURNED FOR CGH.
C           INTEGER*4
C     M   - NUMBER OF LEADING ROWS OF THE GRADIENT TO BE RETAINED,
C           AND NUMBER OF ROWS AND COLUMNS OF THE HESSIAN RETAINED,
C           INPUT.  THE FIRST ROW OF THE JACOBIAN IS THE JACOBIAN
C           OF THE FOURIER FLEXIBLE FORM.  THE NEXT M ROWS ARE
C           ROWS FROM THE JACOBIAN OF THE GRADIENT.  THE NEXT M*M
C           ROWS ARE FROM THE PRINCIPAL SUBMATRIX OF THE HESSIAN
C           THE MULTIPLICATION CGH*THETA WILL RESULT IN A 1+M+M*M
C           VECTOR: FIRST ELEMENT IS THE FFF, NEXT M ELEMENTS ARE
C           THE M LEADING ROWS OF THE GRADIENT OF THE FFF, NEXT
C           M*M ELEMENTS IS THE PRINCIPAL SUBMATRIX OF ORDER M OF
C           THE HESSIAN STORED COLUMNWISE (STORAGE MODE 0).
C           INTEGER* 4
C     DL  - SCALING FACTOR, INPUT.
C           REAL*8
C     X   - VECTOR OF LENGTH N, INPUT.
C           REAL*8
C     CGH - MATRIX OF ORDER 1+M+M*M BY LT, OUTPUT.
C           REAL*8
C     LT  - COMPUTED NUMBER OF COLUMNS OF CGH.
C           INTEGER*4
C     IW  - WORK VECTOR OF LENGTH IAA+3*N.
C           INTEGER*4
C
C
      SUBROUTINE FFFCGH(N,KC,IS,KA,IAA,JJA,M,DL,X,CGH,LT,IW)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
      INTEGER*4 KA(N,1),JJA(1),IW(1)
      REAL*8 CGH(1),X(N)
      INDEX=1
      IMN=INDEX+N
      IMX=IMN+N
      NORM=IMX+N
      NIND=NORM+N
      IF(IS.LT.0) GO TO 100
      CALL Z1FFF(N,KC,IS,KA,IAA,IW(INDEX),IW(IMN),IW(IMX),IW(NORM),
     &IW(NIND))
      LIMIT=NIND+IAA
100   CONTINUE
      NCOL=1+M+M*M
      LD=N+1
      MPLUS1=M+1
      MPLUS2=M+2
      DO 116 J=1,LD
      DO 110 I=1,MPLUS1
      CGH(NCOL*(J-1)+I)=0.D0
      IF((I.EQ.1).AND.(J.GT.1)) CGH(NCOL*(J-1)+I)=X(J-1)
110   IF(I.EQ.J) CGH(NCOL*(J-1)+I)=1.D0
      DO 115 I=MPLUS2,NCOL
115   CGH(NCOL*(J-1)+I)=0.D0
116   CONTINUE
      K=LD
      DO 170 IA=1,IAA
      IF(JJA(IA).LE.0) GO TO 170
      DLKX=0.D0
      DO 120 I=1,N
120   DLKX=DLKX+DL*DFLOAT(KA(I,IA))*X(I)
      K=K+1
      CGH(NCOL*(K-1)+1)=1.D0-.5D0*DLKX*DLKX
      DO 130 I=1,M
130   CGH(NCOL*(K-1)+I+1)=-DL*DLKX*DFLOAT(KA(I,IA))
      DO 135 I=1,M
      DO 135 J=1,M
135   CGH(NCOL*(K-1)+MPLUS1+M*(J-1)+I)
     &=-DL*DL*DFLOAT(KA(I,IA)*KA(J,IA))
      JJ=JJA(IA)
      DO 160 J0=1,JJ
      K=K+2
      DJ=DFLOAT(J0)
      SIN=DSIN(DJ*DLKX)
      COS=DCOS(DJ*DLKX)
      CGH(NCOL*(K-1-1)+1)=2.D0*COS
      CGH(NCOL*(K-1)+1)=-2.D0*SIN
      DO 140 I=1,M
      DKA=DFLOAT(KA(I,IA))
      CGH(NCOL*(K-1-1)+I+1)=-2.D0*DJ*DL*SIN*DKA
140   CGH(NCOL*(K-1)+I+1)=-2.D0*DJ*DL*COS*DKA
      DO 150 I=1,M
      DKAI=DFLOAT(KA(I,IA))
      DO 150 J=1,M
      DKAJ=DFLOAT(KA(J,IA))
      CGH(NCOL*(K-1-1)+MPLUS1+M*(J-1)+I)
     &=-2.D0*DL*DL*DJ*DJ*COS*DKAI*DKAJ
150   CGH(NCOL*(K-1)+MPLUS1+M*(J-1)+I)
     &=2.D0*DL*DL*DJ*DJ*SIN*DKAI*DKAJ
160   CONTINUE
170   CONTINUE
      LT=K
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     FFFCG        4/21/81
C
C     PURPOSE
C     GENERATE MULTI-INDEXES AND THE JACOBIAN OF A VECTOR THAT CONSISTS
C     OF THE FOURIER FLEXIBLE FORM AND ITS GRADIENT.
C
C     USAGE
C     CALL FFFCG(N,KC,IS,KA,IAA,JJA,M,DL,X,CG,LT,IW)
C
C     ARGUMENTS
C     N   - LENGTH OF X, INPUT.  LENGTH MUST BE LESS THAN 11.
C           INTEGER*4
C     KC  - LARGEST NORM OF GENERATED MULTI-INDEXES, INPUT.
C           KC MUST BE LESS THAN 30.
C           INTEGER*4
C     IS  - SWITCH, INPUT.
C           IF IS = 0  THEN NO MULTI-INDEXES ARE DELETED.
C           IF IS = L  THEN MULTI-INDEXES WITH
C                      SUM(KA(I,J): I=1,2,...,L).NE.0
C                      ARE DELETED.
C           IF IS = -1 THEN KA AND IAA ARE NOT COMPUTED AND ARE PRESUMED
C                      TO BE AVAILABLE AS INPUT.
C           INTEGER*4
C     KA  - N BY IAA MATRIX CONTAINING COMPUTED MULTI-INDEXES AS
C           COLUMNS, OUTPUT.  INPUT IF IS=-1.  STORED COLUMNWISE
C           (STORAGE MODE 0)
C           INTEGER*4
C     IAA - NUMBER OF MULTI-INDEXES GENERATED, OUTPUT.  INPUT IF
C           IS=-1.
C           INTEGER*4
C     JJA - VECTOR OF LENGTH IAA, INPUT.  JJA(IA) DEFINES THE ORDER
C           OF THE FOURIER EXPANSION OF THE TERM ASSOCIATED WITH
C           MULTI-INDEX NUMBER IA.  ZERO VALUES ARE PERMITTED.  THEY
C           HAVE THE EFFECT OF DELETING MULTI-INDEX NUMBER IA WHEN
C           COMPUTING CG.  A VECTOR OF ALL ZEROES WILL CAUSE THE
C           IDENTITY TO BE RETURNED FOR CG.
C           INTEGER*4
C     M   - NUMBER OF LEADING ROWS OF THE GRADIENT TO BE RETAINED,
C           INPUT.  THE FIRST ROW OF THE JACOBIAN IS THE JACOBIAN
C           OF THE FOURIER FLEXIBLE FORM.  THE REMANINING ROWS ARE
C           ROWS FROM THE JACOBIAN OF THE GRADIENT.
C           INTEGER* 4
C     DL  - SCALING FACTOR, INPUT.
C           REAL*8
C     X   - VECTOR OF LENGTH N, INPUT.
C           REAL*8
C     CG -  MATRIX OF ORDER M+1 BY LT, OUTPUT.
C           REAL*8
C     LT  - COMPUTED NUMBER OF COLUMNS OF CG.
C           INTEGER*4
C     IW  - WORK VECTOR OF LENGTH IAA+3*N.
C           INTEGER*4
C
C
      SUBROUTINE FFFCG(N,KC,IS,KA,IAA,JJA,M,DL,X,CG,LT,IW)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
      INTEGER*4 KA(N,1),JJA(1),IW(1)
      REAL*8 CG(1),X(N)
      INDEX=1
      IMN=INDEX+N
      IMX=IMN+N
      NORM=IMX+N
      NIND=NORM+N
      IF(IS.LT.0) GO TO 100
      CALL Z1FFF(N,KC,IS,KA,IAA,IW(INDEX),IW(IMN),IW(IMX),IW(NORM),
     &IW(NIND))
      LIMIT=NIND+IAA
100   CONTINUE
      MPLUS1=M+1
      LD=N+1
      DO 110 I=1,MPLUS1
      DO 110 J=1,LD
      CG(MPLUS1*(J-1)+I)=0.D0
      IF((I.EQ.1).AND.(J.GT.1)) CG(MPLUS1*(J-1)+I)=X(J-1)
110   IF(I.EQ.J) CG(MPLUS1*(J-1)+I)=1.D0
      K=LD
      DO 160 IA=1,IAA
      IF(JJA(IA).LE.0) GO TO 160
      DLKX=0.D0
      DO 120 I=1,N
120   DLKX=DLKX+DL*DFLOAT(KA(I,IA))*X(I)
      K=K+1
      CG(MPLUS1*(K-1)+1)=1.D0-.5D0*DLKX*DLKX
      DO 130 I=1,M
130   CG(MPLUS1*(K-1)+I+1)=-DL*DLKX*DFLOAT(KA(I,IA))
      JJ=JJA(IA)
      DO 150 J=1,JJ
      K=K+2
      DJ=DFLOAT(J)
      SIN=DSIN(DJ*DLKX)
      COS=DCOS(DJ*DLKX)
      CG(MPLUS1*(K-1-1)+1)=2.D0*COS
      CG(MPLUS1*(K-1)+1)=-2.D0*SIN
      DO 140 I=1,M
      DKA=DFLOAT(KA(I,IA))
      CG(MPLUS1*(K-1-1)+I+1)=-DL*2.D0*DJ*SIN*DKA
140   CG(MPLUS1*(K-1)+I+1)=-DL*2.D0*DJ*COS*DKA
150   CONTINUE
160   CONTINUE
      LT=K
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
C     FFFG      4/21/81
C
C     PURPOSE
C     GENERATE MULTI-INDEXES AND THE JACOBIAN OF THE GRADIENT OF THE
C     FOURIER FLEXIBLE FORM.
C
C     USAGE
C     CALL FFFG(N,KC,IS,KA,IAA,JJA,M,DL,X,DG,LT,IW)
C
C     ARGUMENTS
C     N   - LENGTH OF X, INPUT.  LENGTH MUST BE LESS THAN 11.
C           INTEGER*4
C     KC  - LARGEST NORM OF GENERATED MULTI-INDEXES, INPUT.
C           KC MUST BE LESS THAN 30.
C           INTEGER*4
C     IS  - SWITCH, INPUT.
C           IF IS = 0  THEN NO MULTI-INDEXES ARE DELETED.
C           IF IS = L  THEN MULTI-INDEXES WITH
C                      SUM(KA(I,J): I=1,2,...,L).NE.0
C                      ARE DELETED.
C           IF IS = -1 THEN KA AND IAA ARE NOT COMPUTED AND ARE PRESUMED
C                      TO BE AVAILABLE AS INPUT.
C           INTEGER*4
C     KA  - N BY IAA MATRIX CONTAINING COMPUTED MULTI-INDEXES AS
C           COLUMNS, OUTPUT.  INPUT IF IS=-1.  STORED COLUMNWISE
C           (STORAGE MODE 0)
C           INTEGER*4
C     IAA - NUMBER OF MULTI-INDEXES GENERATED, OUTPUT.  INPUT IF
C           IS=-1.
C           INTEGER*4
C     JJA - VECTOR OF LENGTH IAA, INPUT.  JJA(IA) DEFINES THE ORDER
C           OF THE FOURIER EXPANSION OF THE TERM ASSOCIATED WITH
C           MULTI-INDEX NUMBER IA.  ZERO VALUES ARE PERMITTED.  THEY
C           HAVE THE EFFECT OF DELETING MULTI-INDEX NUMBER IA WHEN
C           COMPUTING DG.  A VECTOR OF ALL ZEROES WILL CAUSE THE
C           IDENTITY WITH A LEADING COLUMN OF ZEROES TO BE RETURNED
C           FOR DG.
C           INTEGER*4
C     M   - NUMBER OF LEADING ROWS OF THE JACOBIAN TO BE RETAINED.
C           INPUT.
C           INTEGER*4
C     DL  - SCALING FACTOR, INPUT.
C           REAL*8
C     X   - VECTOR OF LENGTH N, INPUT.
C           REAL*8
C     DG  - MATRIX OF ORDER M BY LT, OUTPUT.
C           REAL*8
C     LT  - COMPUTED NUMBER OF COLUMNS OF DG.
C           INTEGER*4
C     IW  - WORK VECTOR OF LENGTH IAA+3*N.
C           INTEGER*4
C
C
      SUBROUTINE FFFG(N,KC,IS,KA,IAA,JJA,M,DL,X,DG,LT,IW)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
      INTEGER*4 KA(N,1),JJA(1),IW(1)
      REAL*8 DG(M,1),X(N)
      INDEX=1
      IMN=INDEX+N
      IMX=IMN+N
      NORM=IMX+N
      NIND=NORM+N
      IF(IS.LT.0) GO TO 100
      CALL Z1FFF(N,KC,IS,KA,IAA,IW(INDEX),IW(IMN),IW(IMX),IW(NORM),
     &IW(NIND))
      LIMIT=NIND+IAA
100   CONTINUE
      K=N+1
      DO 110 I=1,M
      DO 110 J=1,K
110   DG(I,J)=0.D0
      DO 115 I=1,M
115   DG(I,I+1)=1.D0
      DO 160 IA=1,IAA
      IF(JJA(IA).LE.0) GO TO 160
      DLKX=0.D0
      DO 120 I=1,N
120   DLKX=DLKX+DL*DFLOAT(KA(I,IA))*X(I)
      K=K+1
      DO 130 I=1,M
130   DG(I,K)=-DL*DLKX*DFLOAT(KA(I,IA))
      JJ=JJA(IA)
      DO 150 J=1,JJ
      K=K+2
      DJ=DFLOAT(J)
      ASIN=-DL*2.D0*DJ*DSIN(DJ*DLKX)
      ACOS=-DL*2.D0*DJ*DCOS(DJ*DLKX)
      DO 140 I=1,M
      DKA=DFLOAT(KA(I,IA))
      DG(I,K-1)=ASIN*DKA
140   DG(I,K)=ACOS*DKA
150   CONTINUE
160   CONTINUE
      LT=K
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
      SUBROUTINE Z1FFF(LENGTH,KC,IS,INDICE,NUM,ID1,ID2,ID3,ID4,NIND)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
C     PROGRAM TO GENERATE MULTI-INDEXES FOR THE FOURIER FLEXIBLE FORM.
C     WRITTTEN BY JOHN MONAHAN.
C     LENGTH IS THE NUMBER OF COMMODIDTIES OR LENGTH OF THE MULTI-INDEX.
C     KC=||MULTI-INDEX||*  OR NORM OF THE MULTI-INDEX.
C     KCMUST BE LESS THAN 30.
      LOGICAL Z2FFF,OK
      INTEGER*4 INDEX(10),PT,IMN(10),IMX(10),NORM(10)
      INTEGER*4 INDICE(LENGTH,1),NIND(1)
      IF(KC.GE.30) RETURN
      DO 10I=1,10
  10  INDEX(I)=0
      NUM=0
      DO 9L=1,LENGTH
      LP1=L+1
      DO 8I=1,10
      IMN(I)=-KC
  8   IMX(I)=KC
      IMN(L)=1
      NORM(LP1)=0
      DO 1I=1,L
      J=LP1-I
      NORM(J)=NORM(J+1)+IABS(IMN(J))
  1   INDEX(I)=IMN(I)
  2   PT=1
      INDEX(PT)=IMN(PT)
      NORM(PT)=NORM(PT+1)+IABS(INDEX(PT))
  3   CONTINUE
C THIS IS THE CENTER OF THE LOOPING
      IF(NORM(1).GT.KC) GO TO 4
      OK=Z2FFF(INDEX,KC,L,IS,LENGTH)
      IF(.NOT.OK) GO TO 4
      NUM=NUM+1
      NIND(NUM)=NORM(1)
      DO 7I=1,LENGTH
  7   INDICE(I,NUM)=INDEX(LENGTH+1-I)
C  THIS ENDS THE CENTER LOOP
  4   INDEX(PT)=INDEX(PT)+1
      NORM(PT)=NORM(PT+1)+IABS(INDEX(PT))
      IF(INDEX(PT).LE.IMX(PT)) GO TO (3,2,2,2,2,2,2,2,2,2),PT
      INDEX(PT)=IMN(PT)
      PT=PT+1
      IF(PT.LE.L) GO TO 4
  9   CONTINUE
      CALL Z3FFF(NIND,INDICE,NUM,LENGTH)
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
      LOGICAL FUNCTION Z2FFF(INDEX,KC,L,IS,LENGTH)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
      INTEGER*4 INDEX(1),P,PRIME(10)
      DATA PRIME/2,3,5,7,11,13,17,19,23,29/
      Z2FFF=.TRUE.
      IF(IS.EQ.0) GO TO 1
      ISUM=0
      DO 5 I=1,IS
  5   ISUM=ISUM+INDEX(LENGTH+1-I)
      IF(ISUM.NE.0) GO TO 4
  1   DO 2 I=1,10
      P=PRIME(I)
      IF(P.GT.KC) RETURN
      DO 3 J=1,L
      IF(IABS(INDEX(J)).EQ.1) RETURN
      IF(MOD(INDEX(J),P).NE.0) GO TO 2
  3   CONTINUE
      Z2FFF=.FALSE.
      RETURN
  2   CONTINUE
      RETURN
  4   CONTINUE
      Z2FFF=.FALSE. 
      RETURN
      END
C....*...1.........2.........3.........4.........5.........6.........7.*.......8
      SUBROUTINE Z3FFF(A,B,N,K)
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4(I-N)
      SAVE
      INTEGER*4 A(N),B(K,N),R
      M=N
20    M=M/2
      IF(M)30,40,30
30    KK=N-M
      J=1
41    I=J
49    L=I+M
      IF(A(I)-A(L))60,60,50
50    R=A(I)
      A(I)=A(L)
      A(L)=R
      DO 55 JJ=1,K
      R=B(JJ,I)
      B(JJ,I)=B(JJ,L)
55    B(JJ,L)=R
      I=I-M
      IF(I-1)60,49,49
60    J=J+1
      IF(J-KK)41,41,20
40    CONTINUE
      RETURN
      END
