C C PROGRAM EMCEE C C Eigenvector Method with Constraint Condition C for C COD and Dispersion Correction C C C First Coded By N. Nakamura (04/11/1993) C Second Coded By Masanori SATOH (25/5/1998) C Third Coded By K.Harada (24,Jan,2000) C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER COMMENT*35, OUTFILE*12, RESPNCFL*12 REAL*8 TPNB, TPNA REAL*8 E(300),V(300,300),W(300,300) INTEGER IC, IE, IM, DRCTN REAL*8 M(300,300), MTRN(300,300), MM(300,300), EINV(300) REAL*8 R(300) REAL*8 Z(300), RESRMS, ETARMS REAL*8 ETA(300),PI REAL*8 C(300,300),CC(300,300) REAL*8 P(300,300),WP(300,300) REAL*8 QQ(300,300),Q(300,300),QCC(300,300),QCCW(300,300) REAL*8 DM1(300,300),DM2(300,300),DM(300,300) REAL*8 DV1(300),DV2(300),ZL(300) INTEGER LCST, LP(300) REAL*8 PP(300,300),EPINV(300),CCTRN(300,300) REAL*8 CX(300),OK2(300),MX(300) C PI=ACOS(-1.0) C C Reading command and optics data file from SAD C File name may be "fort.89", "fort.90" please C OPEN(UNIT=89,FILE='fort.89',STATUS='OLD',FORM='FORMATTED') 10 FORMAT(A35,1F15.9) 20 FORMAT(A35,A12) 30 FORMAT(1F15.9,1F15.9,1F15.9, 1F15.9) 40 FORMAT(A35,1F15.9) 50 FORMAT(1F15.9) 60 FORMAT(A35) READ(89,10) COMMENT, TPNB Write(6,10) COMMENT, TPNB IE = INT(TPNB) READ(89,10) COMMENT, TPNB Write(6,10) COMMENT, TPNB IC = INT(TPNB) READ(89,10) COMMENT, TPNB Write(6,10) COMMENT, TPNB IM = INT(TPNB) IF(IE.GT.IC)IE=IC READ(89,10) COMMENT, TPNB Write(6,10) COMMENT, TPNB DRCTN = INT(TPNB) READ(89,40) COMMENT, CSTNB Write(6,40) COMMENT, CSTNB LCST = INT(CSTNB) READ(89,60) COMMENT Write(6,60) COMMENT DO 260 i=1, LCST READ(89,50) CSTPT LP(i) = INT(CSTPT) Write(6,*) LP(i) 260 CONTINUE READ(89,20) COMMENT, RESPNCFL Write(6,20) COMMENT, RESPNCFL READ(89,20) COMMENT, OUTFILE Write(6,20) COMMENT, OUTFILE DO 270 i=1, IM READ(89,50) ETA(I) C Write(6,*) i," ",ETA(I) 270 CONTINUE CLOSE(89) OPEN(UNIT=90,FILE=RESPNCFL,STATUS='OLD',FORM='FORMATTED') Write(6,*) 'Reading' DO 250 K = 1, IC DO 70 L = 1, IM READ(90,30) A, B, TPNB, TPNA C Write(6,30) A, B, TPNB, TPNA C Write(6,*) K, L IF (DRCTN. EQ. 1) THEN M(L,K) = TPNB ELSE M(L,K) = TPNA ENDIF 70 CONTINUE 250 CONTINUE CLOSE(90) DO 80 I=1,IC DO 90 J=1,IM MTRN(I,J)=M(J,I) 90 CONTINUE 80 CONTINUE CALL MATMUL(MTRN,M,MM,IC,IM,IC) CCCCCCCCC Constructiong CSTRNT CCCCCCCCCCCCC Write(6,*) "Restraint Point Number=",LCST C Write(6,*) "Restraint Point =:" C Do 130 J = 1, LCST C Write(6,*) " Position Monitor No.",LP(J) C 130 CONTINUE DO 370 I=1,LCST DO 380 J=1,IC CC(I,J)=M(LP(I),J) 380 CONTINUE 370 CONTINUE DO 100 I=1,IC DO 110 J=1,LCST CCTRN(I,J)=CC(J,I) 110 CONTINUE 100 CONTINUE DO 360 I=1,LCST ZL(I)=-ETA(LP(I)) 360 CONTINUE CCCCCCCCC Givens-Householder CCCCCCCCCCCC CALL GIVHO(MM,E,V,IC,IE,IE) DO 140 I=1,IE EINV(I)=1.0D00/E(I) C Write(6,*) I, EINV(I),E(I) 140 CONTINUE DO 150 I=1,IC DO 160 J=1,IC W(I,J)=0.0 DO 170 K=1,IE W(I,J)=W(I,J)-EINV(K)*V(I,K)*V(J,K) C Write(6,*) W(I,J) 170 CONTINUE 160 CONTINUE 150 CONTINUE CALL MATMUL(CC,W,PP,LCST,IC,IC) CALL MATMUL(PP,CCTRN,P,LCST,IC,LCST) CALL GIVHO(P,E,V,LCST,LCST,LCST) C Write(6,*) 'For C-matrix' DO 310 I=1,LCST EINV(I)=1.0D00/E(I) C Write(6,*) I, EINV(I), E(I) 310 CONTINUE DO 280 I=1,LCST DO 290 J=1,LCST WP(I,J)=0.0 DO 300 K=1,LCST WP(I,J)=WP(I,J)+EINV(K)*V(I,K)*V(J,K) 300 CONTINUE 290 CONTINUE 280 CONTINUE CALL MATMUL(W,CCTRN,QQ,IC,IC,LCST) CALL MATMUL(QQ,WP,Q,IC,LCST,LCST) CALL MATMUL(Q,CC,QCC,IC,LCST,IC) CALL MATMUL(QCC,W,QCCW,IC,IC,IC) CALL MATMUL(QCCW,MTRN,DM2,IC,IC,IM) DO 180 I=1,IC DO 190 J=1,IC W(I,J)=-W(I,J) 190 CONTINUE 180 CONTINUE CALL MATMUL(W,MTRN,DM1,IC,IC,IM) CALL MATADD(DM1,DM2,DM,IC,IM) CALL MATMUL(DM,ETA,DV1,IC,IM,1) DO 200 I=1,IC DO 210 J=1,LCST Q(I,J)=-Q(I,J) 210 CONTINUE 200 CONTINUE CALL MATMUL(Q,ZL,DV2,IC,LCST,1) CALL MATADD(DV1,DV2,R,IC,1) Write(6,*) 'CSTRNT point : LP : CX + Z = 0' CALL MATMUL(CC,R,CX,L,IC,1) CALL MATADD(CX,ZL,OK2,L,1) DO 320 I=1, LCST Write(6,390) LP(I),CX(I),'+',ZL(I),'=',OK2(I) 320 CONTINUE 390 FORMAT(I4,2X,1F12.9,1X,1A,1X,1F12.9,1X,1A,1X,1D15.7) Write(6,*) 'Is correction well done?' CALL MATMUL(M,R,MX,IM,IC,1) RESRMS = 0 ETARMS = 0 Do 350 I=1, IM ETARMS = ETARMS + ETA(I)*ETA(I) RESRMS = RESRMS + (MX(I)-ETA(I))*(MX(I)-ETA(I)) C Write(6,*) MX(I), ETA(I), MX(I) - ETA(I) 350 CONTINUE Write(6,400) ' RMS of Arg before correction : ',Sqrt(ETARMS)/IM Write(6,400) ' RMS of Arg after correction : ',Sqrt(RESRMS)/IM 340 CONTINUE 330 CONTINUE 400 FORMAT(1A32,1D15.7) OPEN(UNIT=20,FILE=OUTFILE,STATUS='UNKNOWN',FORM='FORMATTED') Write(20,*) IC DO 230 A=1, IC I = INT(A) C WRITE(20,30) A, B, R(I,J) Write(20,*) R(I) C WRITE(6,*) I,"th steer", R(I) 230 CONTINUE CLOSE(20) 999 CONTINUE STOP END C ************************************************************* C ************************************************************* SUBROUTINE GIVHO(A,E,V,N,NEV,NVEC) C ************************************************************* C ************************************************************* C PURPOSE C COMPUTE EIGENVALUES AND EIGENVECTORS OF A REAL NxN SYMMETRIC C MATRIX BY THE GIVENS-HOUSEHOLDER METHOD C DESCRIPTION OF PARAMETERS C A - NxN INPUT MATRIX C NEV - NUMBER OF LARGEST EIGEN VALUES COMPUTED AND STORED IN E. C NVEC - NUMBER OF EIGEN VECTORS COMPUTED AND STORED IN V. NVEC<=NEV C OVINV - INVERSE OF OVERFLOW LIMIT FOR VAX COMPUTER C IY - INITIAL VALUE GENERATING RANDOM VECTOR BY FUNC. RAND C COMMENT C DIMENSION STATEMENT VALID FOR NORDER UP TO 100 C THIS PROGRAM IS BASICALLY VALID FOR NON-ZERO B(J). C BUT EIGENVALUES AND EIGENVECTORS MAY BE OBTAINED APPROXIMATELY, C IF ZERO IS REPLACED WITH A SMALL VALUE(1.0E-5 IN THIS CASE) C IF YOU HAVE A TROUBLE, YOU SHOULD USE AN ALTERNATE WAY C (SEE J.M. ORTEGA, "ON STURM SEQUENCES FOR TRIDIAGONAL MATRICES," C J. ASSOC. COMP. MACH., 7, 260-263(1960).) C C CODED BY N. NAKAMURA(25/10/1993) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 A(300,300),E(300),V(300,300) REAL*8 B(300),C(300),P(300),Q(300),R(300),W(300),Y(300) DIMENSION IN(300) INTEGER AG,IY, N, NEW, NVEC REAL*8 KAP,NORM,LAMBDA,L,MULT LOGICAL FIRST,IN CNN DATA EPS,EPS2,OVINV,IY/1.D-6,1.D-2,1.D-30,31415967D0/ CNN NM1=N-1 NM2=N-2 C BEGIN THE REDUCTION TO TRIDIAGONAL FORM. C THE DIAGONAL AND OFF-DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX ARE C STORED IN C AND B. IF(N.LE.2)GO TO 9 DO 8 I=1,NM2 IP1=I+1 SS=0. DO 1 J=IP1,N SS=SS+A(J,I)*A(J,I) C Write(6,*) N, SS 1 CONTINUE S=SQRT(SS) IF(A(IP1,I).LT.0.)S=-S C(I)=A(I,I) B(I)=-S C IF S=0 THEN ALPHA MUST BE ZERO. ALPHA=0. IF(S.EQ.0.0)GO TO 8 ALPHA=1./(SS+A(IP1,I)*S) T=A(IP1,I)+S A(IP1,I)=T W(I+1)=T IP2=I+2 DO 2 J=IP2,N 2 W(J)=A(J,I) DO 4 J=IP1,N T=0.0 DO 3 K=IP1,N 3 T=T+A(J,K)*W(K) 4 P(J)=T*ALPHA KAP=0 DO 5 K=IP1,N 5 KAP=KAP+W(K)*P(K) KAP=.5*KAP*ALPHA DO 6 K=IP1,N 6 Q(K)=P(K)-KAP*W(K) DO 7 J=IP1,N C COMPUTE THE NEW A. BY SYMMETRY ONLY HALF THE ELEMENTS NEED BE DONE. DO 7 K=J,N A(J,K)=A(J,K)-(Q(J)*W(K)+Q(K)*W(J)) 7 A(K,J)=A(J,K) 8 A(I,I)=ALPHA 9 C(N-1)=A(N-1,N-1) C(N)=A(N,N) B(N-1)=A(N-1,N) C WRITE(6,*)'C=',(C(I),I=1,N) C WRITE(6,*)'B=',(B(I),I=1,N-1) C WRITE(6,*)'A=',(A(I,I),I=1,N) DO 100 I=1,N-1 IF(B(I).EQ.0.0)THEN B(I)=1.E-5 C WRITE(6,*)'B=0' ENDIF 100 CONTINUE C THIS COMPLETES THE REDUCTION TO TRIDIAGONAL FORM. C BEGIN THE EIGENVALUE COMPUTATION. C COMPUTE NORM OF TRIDIAGONAL MATRIX, SQUARES OF THE B(I). NORM=ABS(C(1))+ABS(B(1)) C write(6,*) NORM DO 10 I=2,N T=ABS(C(I))+ABS(B(I))+ABS(B(I-1)) 10 NORM = MAX(NORM,T) C 10 NORM=AMAX1(NORM,T) C WRITE(6,*)'NORM=',NORM DO 11 I=1,N 11 W(I)=B(I)**2 C SET K AND UPPER AND LOWER ESTIMATES. K=1 U=NORM DO 12 I=1,NEV E(I)=-NORM 12 CONTINUE C write(6,*) NORM C BEGIN MAIN LOOP. 13 L=E(K) 14 LAMBDA=.5*(L+U) C write(6,*) LAMBDA C THE CONVERGENCE TEST AS IMPLEMENTED HERE ALLOWS THE COMPUTATION TO C PROCEED UNTIL THE INTERVAL(L,U) CAN BE MADE NO SMALLER. IF((LAMBDA.EQ.L).OR.(LAMBDA.EQ.U))GO TO 30 C BEGIN COMPUTATION OF NUMBER OF SIGN AGREEMENTS,AG AG=0 I=1 16 S=C(I)-LAMBDA 18 IF(S.GE.0.0)AG=AG+1 IF(S.EQ.0.0)GO TO 20 I=I+1 IF(I.GT.N)GO TO 22 CNN SNUM=S*(C(I)-LAMBDA)-W(I-1) IF(SNUM.EQ.0.)THEN S=0.0 GO TO 18 ENDIF SINV=S/SNUM IF(ABS(SINV).GE.OVINV)THEN S=C(I)-LAMBDA-W(I-1)/S GO TO 18 ENDIF CNN C IF OVERFLOW OCURRS, PROCEED AS IF S HAD BEEN ZERO. AG=AG+1 I=I-1 20 I=I+2 IF(I.LE.N)GO TO 16 C THE COMPUTATION OF AG IS COMPLETE. ADJUST INTERVAL. 22 IF(AG.GE.K)GO TO 24 U=LAMBDA GO TO 14 24 L=LAMBDA M=MIN0(AG,NEV) DO 26 I=K,M 26 E(I)=LAMBDA GO TO 14 C THE K-TH EIGENVALUE IS COMPUTED. STORE IN E(K) AND PROCEED. 30 E(K)=LAMBDA K=K+1 IF(K.LE.NEV)GO TO 13 C THIS COMPLETES THE EIGENVALUE COMPUTATION. C IF NVEC NOT ZERO, BEGIN THE EIGENVECTOR COMPUTATION. IF(NVEC.NE.0)GO TO 40 RETURN 40 DO 90 I=1,NVEC C INITIALIZATION FOR THIS EIGENVECTOR. DO 44 J=1,N P(J)=0. Q(J)=B(J) R(J)=C(J)-E(I) C WRITE(6,*)'J,P,Q,R=',J,P(J),Q(J),R(J) 44 Y(J)=1. Y(N+1)=0. Y(N+2)=0. FIRST=.TRUE. C BEGIN THE GAUSSIAN ELIMINATION REDUCTION TO TRIANGULAR FORM. DO 50 J=1,NM1 IF(ABS(R(J)).LT.ABS(B(J)))GO TO 46 C WRITE(6,*)'J,R(J)=',J,R(J) C IF(R(J).EQ.0.)R(J)=OVINV MULT=B(J)/R(J) IN(J)=.FALSE. GO TO 48 46 MULT=R(J)/B(J) IN(J)=.TRUE. R(J)=B(J) T=R(J+1) R(J+1)=Q(J) Q(J)=T P(J)=Q(J+1) Q(J+1)=0. C STORE THE MULTIPLIERS FOR LATTER USE. 48 W(J)=MULT Q(J+1)=Q(J+1)-MULT*P(J) R(J+1)=R(J+1)-MULT*Q(J) IF(ABS(R(J)).EQ.0.0)R(J)=OVINV 50 CONTINUE IF(ABS(R(N)).EQ.0.0)R(N)=OVINV C DO 110 J=1,N C 110 WRITE(6,*)'J,P,Q,R=',J,P(J),Q(J),R(J) C THE COMPLETES THE REDUCTION. BEGIN THE BACK SUBSTITUTION. C IF E(I) IS A CLUSTER, GENERATE A RANDOM VECTOR. C WRITE(6,*)'NORM=',NORM C IF(I.NE.1.AND.ABS(E(I)-E(I-1)).LT.NORM*EPS)E(I)=E(I-1)-NORM*EPS*2 IF(I.EQ.1.OR.ABS(E(I)-E(I-1)).GE.NORM*EPS)GO TO 54 51 CONTINUE C WRITE(6,*)'RAN,I=',I DO 52 J=1,N 52 Y(J)=RAND(IY) C WRITE(6,*)'Y(RAN)=',(Y(J),J=1,N) 54 CONTINUE DO 66 JI=1,N K=N-JI+1 C SAVE Y(K) TO SCALE IF OVERFLOW. T=Y(K) C Y(N+1) AND Y(N+2) ARE ZERO. 62 CONTINUE YNUM=T-Y(K+1)*Q(K)-Y(K+2)*P(K) C WRITE(6,*)'YNUM=',YNUM IF(ABS(YNUM).LE.OVINV)THEN Y(K)=0.0 GO TO 66 ENDIF YINV=R(K)/YNUM C WRITE(6,*)'K,R,YINV=',K,R(K),YINV IF(ABS(YINV).GE.OVINV)THEN Y(K)=(T-Y(K+1)*Q(K)-Y(K+2)*P(K))/R(K) GO TO 66 ENDIF C WRITE(6,*)'Y62(OVFL)=',(Y(J),J=1,N) DO 64 J=1,N 64 Y(J)=Y(J)*1.E-5 T=T*1.E-5 C WRITE(6,*)'Y62*1E-5=',(Y(J),J=1,N) GO TO 62 66 CONTINUE C THIS COMPLETES BACK SUBSTITUTION. IF(.NOT.FIRST)GO TO 74 FIRST=.FALSE. C DO THE REDUCTION ON THE NEW RIGHT HAND SIDE. DO 70 J=1,NM1 IF(IN(J))GO TO 68 Y(J+1)=Y(J+1)-W(J)*Y(J) GO TO 70 68 T=Y(J) Y(J)=Y(J+1) Y(J+1)=T-W(J)*Y(J+1) 70 CONTINUE C THIS COMPLETES REDUCTION OF RIGHT SIDE. DO ONE MORE ITERATION. GO TO 54 C Y IS NOW AN EIGENVECTOR OF TRIDIAGONAL MATRIX. BEGIN BACK C TRANSFORMATION. 74 DO 78 J=1,NM2 T=0.0 K=N-J-1 M=K+1 DO 76 KK=M,N T=T+A(KK,K)*Y(KK) 76 Continue T=A(K,K)*T DO 78 KK=M,N Y(KK)=Y(KK)-T*A(KK,K) 78 Continue C BACK TRANSFORMATION COMPLETE. Y NOW EIGENVECTOR OF A. NORMALIZE. T=ABS(Y(1)) K=1 DO 80 J=2,N S=ABS(Y(J)) IF(S.LE.T)GO TO 80 T=S K=J 80 CONTINUE T=1.0/Y(K) DO 82 J=1,N V(J,I)=Y(J)*T 82 CONTINUE C WRITE(6,*)'Y=',(Y(J),J=1,N) C I-TH COLUMN OF V IS I-TH EIGENVECTOR NORMALIZED WITH MAX ELEMENT =1 C NORMALIZE EIGENVECTORS WITH NORM=1 TV=0.0 DO 84 J=1,N TV=TV+V(J,I)*V(J,I) 84 CONTINUE TV=1./SQRT(TV) DO 86 J=1,N V(J,I)=V(J,I)*TV 86 CONTINUE IF(I.EQ.1)GO TO 90 IF(ABS(E(I)-E(I-1)).LT.NORM*EPS)THEN PROD=0.0 DO 88 J=1,N C WRITE(6,*)'V(I-1),V(I)=',V(J,I-1),V(J,I) PROD=PROD+V(J,I-1)*V(J,I) 88 CONTINUE C WRITE(6,*)'PROD=',PROD IF(ABS(PROD).GE.EPS2)THEN FIRST=.TRUE. GO TO 51 ENDIF ENDIF C 90 CONTINUE C RETURN END C ************************************************************* C ************************************************************* C Matrix Calculation Series C ************************************************************* C ************************************************************* SUBROUTINE MATADD(A,B,C,M,N) C PURPOSE C COMPUTE SUM OF TWO MxN MATRICE. C DESCRIPTION OF PARAMETER C A,B - MxN MATRIX C C - C=A+B C COMMENT C DIMENSION STATEMENT VALID FOR N AND M UP TO 100 C CODED BY N. NAKAMURA (26/10/1993) REAL*8 A(300,300),B(300,300),C(300,300) DO 10 I=1,M DO 10 J=1,N C(I,J)=A(I,J)+B(I,J) 10 CONTINUE RETURN END C SUBROUTINE MATRED(A,B,C,M,N) C PURPOSE C COMPUTE DIFFERENCE OF TWO MxN MATRICE. C DESCRIPTION OF PARAMETER C A,B - MxN MATRIX C C - C=A-B C COMMENT C DIMENSION STATEMENT VALID FOR N AND M UP TO 100 C CODED BY N. NAKAMURA (26/10/1993) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 A(300,300),B(300,300),C(300,300) DO 10 I=1,M DO 10 J=1,N C(I,J)=A(I,J)-B(I,J) 10 CONTINUE RETURN END C SUBROUTINE MATMUL(A,B,C,M,L,N) C PURPOSE C COMPUTE SUM OF TWO MxN MATRICE. C DESCRIPTION OF PARAMETER C A - MxL MATRIX C B - LxN MATRIX C C - C=A*B,MxN MATRIX C COMMENT C DIMENSION STATEMENT VALID FOR N,L AND M UP TO 100 C CODED BY N. NAKAMURA (26/10/1993) IMPLICIT REAL*8 (A-H,O-Z) REAL*8 A(300,300),B(300,300),C(300,300) C REAL I,J,K DO 10 I=1,M DO 10 J=1,N C(I,J)=0.D0 DO 10 K=1,L C(I,J)=C(I,J)+A(I,K)*B(K,J) 10 CONTINUE RETURN END