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