***********************************************************************
* 程序说明：N --自由度 *
* NS --需求的振型数目 *
* RTOL--控制精度 *
* SM --质量矩阵 *
* MS12--质量矩阵控制系数：1--对角阵 *
* 2--满阵 *
* KD12--刚度柔度控制系数：1--柔度矩阵 *
* 2--刚度矩阵 *
* W(N)--频率向量 *
* T(N)--周期向量 *
***********************************************************************
REAL A(1000)
OPEN(1,FILE='INPUT.DAT')
OPEN(2,FILE='OUTPUT.DAT')
READ(1,*)N,NS,MS12,KD12,RTOL
KA=1
KYO=KA+N*N
KSM=KYO+N
KU=KSM+N*N
KYOP=KU+N*NS
KR=KYOP+N*N
KBT=KR+N
KSX=KBT+N
KYK=KSX+N
KW=KYK+N
KT=KW+N
CALL ITER(A(KA),A(KYO),A(KSM),A(KU),A(KYOP),A(KR),
! A(KBT),A(KSX),A(KYK),A(KW),A(KT),N,NS,MS12,KD12,RTOL)
STOP
END
*
*-----------------------------------------------------------------------
SUBROUTINE ITER(A,YO,SM,U,YOP,R,BT,SX,YK,W,T,N,MJ,MS12,KD12,RTOL)
*-----------------------------------------------------------------------
REAL A(N,N),YO(N),SM(N,N),U(N,MJ),YOP(N),R(N),BT(N),SX(N)
! ,YK(N),W(N),T(N)
WRITE(2,*)'THE SYSTEM STIFF MATRIX'
*
DO 5 I=1,N
READ(1,*)(A(I,J),J=1,N)
WRITE(2,*)(A(I,J),J=1,N)
5 CONTINUE
*
WRITE(2,*)'THE SYSTEM MASS MATRIX'
IF (MS12.EQ.2) THEN
DO 6 I=1,N
READ(1,*)(SM(I,J),J=1,N)
WRITE(2,*)(SM(I,J),J=1,N)
6 CONTINUE
*
ELSE
READ(1,*)(SM(I,I),I=1,N)
WRITE(2,*)(SM(I,I),I=1,N)
END IF
*
IF (KD12.EQ.2) THEN
CALL UTDU3(A,SX,N)
END IF
*
DO 300 J=1,MJ
J1=J-1
NITER=0
W1=0
*
DO 10 I=1,N
U(I,J)=1.0
10 CONTINUE
*
DO 30 I=1,N
SUM=0.0
DO 20 K=1,N
SUM=SUM+SM(I,K)*U(K,J)
20 CONTINUE
YO(I)=SUM
30 CONTINUE
*
40 NITER=NITER+1
IF(NITER.GE.40) STOP
*
DO 50 I=1,N
YOP(I)=YO(I)
50 CONTINUE
*
IF(J.EQ.1)GOTO 80
DO 70 L=1,J1
RJ=0.0
DO 60 I=1,N
RJ=RJ+U(I,J)*YO(I)
60 CONTINUE
BT(L)=R(L)*RJ
write(*,*)R(L),rj
70 CONTINUE
*
80 IF (KD12.EQ.2)THEN
CALL DECOMP(A,YO,N)
ELSE
END IF
*
DO 120 I=1,N
U(I,J)=YO(I)
120 CONTINUE
*
IF(J.EQ.1) GOTO 150
DO 140 I=1,N
ST=0.0
DO 130 L=1,J1
ST=ST+BT(L)*U(I,L)
130 CONTINUE
YO(I)=YO(I)-ST
140 CONTINUE
*
150 DO 170 I=1,N
SUM=0.0
DO 160 K=1,N
SUM=SUM+SM(I,K)*YO(K)
160 CONTINUE
YK(I)=SUM
170 CONTINUE
*
ST=0.0
RJ=0.0
DO 180 I=1,N
ST=ST+YO(I)*YOP(I)
RJ=RJ+YO(I)*YK(I)
180 CONTINUE
*
W2=ST/RJ
ST=1.0/SQRT(RJ)
DO 190 I=1,N
YO(I)=YK(I)*ST
U(I,J)=U(I,J)*ST
190 CONTINUE
*
EL=ABS(W2-W1)
W1=W2
EL=EL/W2
IF(EL.GT.RTOL) GOTO 40
*
DO 210 I=1,N
SUM=0.0
DO 200 K=1,N
SUM=SUM+SM(I,K)*U(K,J)
200 CONTINUE
YK(I)=SUM
210 CONTINUE
*
RJ=0.0
DO 220 I=1,N
RJ=RJ+YK(I)*U(I,J)
220 CONTINUE
*
R(J)=1/(W2*RJ)
W(J)=SQRT(W2)
WRITE(2,230)J,W(J)
WRITE(2,240)J,(U(I,J)/U(1,J),I=1,N)
WRITE(*,*)NITER
230 FORMAT('W(',I2,')=',F15.7,5X)
240 FORMAT('A(',I2,')=',5F15.7)
300 CONTINUE
RETURN
END
*
*----------------------------------------------------------------------
SUBROUTINE UTDU3(A,C,N)
*----------------------------------------------------------------------
REAL A(N,N),C(N)
IF (N.EQ.1) GOTO 40
L=N-1
*
DO 20 I=1,L
IF(A(I,I).LE.0.0) GOTO 30
M=I+1
DO 10 J=M,N
C(J)=A(I,J)
A(I,J)=A(I,J)/A(I,I)
DO 10 K=M,J
A(K,J)=A(K,J)-A(I,J)*C(K)
10 CONTINUE
20 CONTINUE
*
IF(A(N,N).GT.0) GOTO 40
I=N
30 WRITE(*,*)I,A(I,J)
STOP
40 RETURN
END
*
*-------------------------------------------------------------------------
SUBROUTINE DECOMP(A,B,N)
*-------------------------------------------------------------------------
REAL A(N,N),B(N)
IF (N.GT.1)GOTO 10
B(1)=B(1)/A(1,1)
GOTO 50
10 L=N-1
*
DO 20 I=2,N
K=I-1
DO 20 J=1,K
B(I)=B(I)-B(J)*A(J,I)
20 CONTINUE
*
DO 30 I=1,N
B(I)=B(I)/A(I,I)
30 CONTINUE
*
DO 40 I=1,L
J1=L-I+1
IP=J1+1
DO 40 J=IP,N
B(J1)=B(J1)-B(J)*A(J1,J)
40 CONTINUE
*
50 RETURN
END |
发表于 2008-3-29 21:33
提升卡
变色卡
显身卡
楼主