条件付固有値法のプログラム解説
プログラム自体はサブルーチンも含めて加速器概要のトップページからダウンロード可能です。
C
C PROGRAM EMCEE
C
C
Eigenvector Method with Constraint Condition
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
COD、分散関数、ベータなど（解説で"y"としたベクトルデータ）の読み込み
ＳＡＤからの出力ファイルの例(fort.89）
number of eigen value to use
: .0000000000
number of corrector to use
: 111.0000000000
number of monitor to use
: 131.0000000000
direction of calculation
: 1.0000000000
number of constraint condition
: .0000000000
number of constraint BMPs below :
responce matrix data file name
: fort.92
output file name
: fort.893
.3703329221
.2011214551
：
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)
応答行列（解説で"M"とした行列データ）の読み込み
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
束縛条件（解説で"C"とした行列）の作成
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) CC …… CT
380 CONTINUE
370 CONTINUE
DO 100 I=1,IC
DO 110 J=1,LCST
CCTRN(I,J)=CC(J,I) CCTRN
…… C
110
CONTINUE
100 CONTINUE
DO 360 I=1,LCST
ZL(I)=-ETA(LP(I)) ZL …… Z
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) 1/ki
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) W
…… A-1
C
Write(6,*) W(I,J)
170 CONTINUE
160 CONTINUE
150 CONTINUE
CALL MATMUL(CC,W,PP,LCST,IC,IC) PP=CC*W
…… CTA-1
CALL MATMUL(PP,CCTRN,P,LCST,IC,LCST) P=PP*CCTRN
…… CTA-1C
CALL GIVHO(P,E,V,LCST,LCST,LCST)
P=CTA-1Cの逆行列を求める
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) QQ
= A-1C
CALL MATMUL(QQ,WP,Q,IC,LCST,LCST) Q
= A-1CP-1
CALL MATMUL(Q,CC,QCC,IC,LCST,IC) QCC
= A-1CP-1CT
CALL MATMUL(QCC,W,QCCW,IC,IC,IC) QCCW
= A-1CP-1CTA-1
CALL MATMUL(QCCW,MTRN,DM2,IC,IC,IM) DM2
= A-1CP-1CTA-1MT
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) DM1
= -A-1MT
CALL MATADD(DM1,DM2,DM,IC,IM) DM
= (-A-1MT+A-1CP-1CTA-1MT)
CALL MATMUL(DM,ETA,DV1,IC,IM,1) DV1
= (-A-1MT+A-1CP-1CTA-1MT)
y
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) DV2 = -(A-1CP-1) z
求まった蹴り角（解説の"x"ベクトル）返し
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 *************************************************************
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
ＳＯＲ施設加速器部門