结构有限元一维变带宽存储LD分解及求解代码
.考虑到编写结构有限元程序的参考，给出结构有限元一维变带宽存储LD分解及求解代码：
C
C This subroutine is used to solve finite element static equations
Cin core , using compacted storage and column reduction scheme .
C
CINPUT VARIABLES
C A => Stiffness matrix stored in compacted form
C V => Right-hand-side load vecter
C MAXA => Vector dontaining addresses of diagonal elements
C of stiffness matrix inA
C NN => Number of equations in the system
C NWK => Number of elements below skyline of matrix
C KK => EQ.1 triangularization of matrix
C EQ.2 reduction and back-substitution of load vector
C
C OUTPUT VARIABLES
C A => D and L - factors of stiffness matrix
C V => Displacement vector
C
C
SUBROUTINE COLSOL(A,V,MAXA,NN,KK)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(1),V(1),MAXA(1)
C----------------------------------------------------------------------
C PERFORM L*D*L(T) FACTORIZATION OF STIFFNESS MATRIX
C----------------------------------------------------------------------
IF(KK-2)40,150,150
40 DO 140 N=1,NN
KN = MAXA(N)
KL = KN + 1
KU = MAXA(N+1) - 1
KH = KU - KL
IF(KH) 110,90,50
50 K= N - KH
IO = 0
KLT= KU
DO 80 J=1,KH
IO = IO + 1
KLT= KLT - 1
KI = MAXA(K)
ND = MAXA(K+1) - KI - 1
IF(ND) 80,80,60
60 KK = MIN(IO,ND)
C= 0.
DO 70 L=1,KK
70 C = C + A(KI+L)*A(KLT+L)
A(KLT)=A(KLT) - C
80 K = K + 1
90 K = N
B = 0.
DO 100 KK=KL,KU
K = K - 1
KI= MAXA(K)
C = A(KK)/A(KI)
B = B + C*A(KK)
100 A(KK) = C
A(KN) = A(KN) - B
110 IF(A(KN)) 120,120,140
120 WRITE(*,2000) N
STOP
140 CONTINUE
RETURN
C----------------------------------------------------------------------
C REDUCE RIGHT-HAND-SIDE LOAD VECTOR
C----------------------------------------------------------------------
150 DO 180 N=1,NN
KL = MAXA(N) + 1
KU = MAXA(N+1) - 1
IF(KU-KL) 180,160,160
160 K = N
C = 0.
DO 170 KK=KL,KU
K = K - 1
170 C = C + A(KK)*V(K)
V(N) = V(N) - C
180 CONTINUE
C----------------------------------------------------------------------
C BACK-SUBSTITUTE
C----------------------------------------------------------------------
DO 200 N=1,NN
K = MAXA(N)
200 V(N) = V(N)/A(K)
IF(NN.EQ.1) RETURN
N = NN
DO 230 L=2,NN
KL = MAXA(N) + 1
KU = MAXA(N+1) - 1
IF(KU-KL) 230,210,210
210 K = N
DO 220 KK=KL,KU
K = K - 1
220 V(K) = V(K) - A(KK)*V(N)
230 N = N - 1
RETURN
C
2000 FORMAT(//5x,'! Error in *COLSOL* ',
&5x,' Stiffness matrix not positive definite'//
&5x,'Nonpositive pivot for equation = ',i4//)
END
C 1998.06.28 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 1998.06.28
C |
C |
C |
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
SUBROUTINE GASJ(A,F,IK,JK,N)
IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION A(N,N),F(N),IK(N),JK(N)
EPS = 1.0E-9
DO 200 K=1,N
C = 0.
DO 130 I=1,N
IF(K.EQ.1) GOTO 110
DO 100 L=1,K-1
IF(I.EQ.IK(L)) GOTO 130
100 CONTINUE
110 DO 120 J=1,N
IF(ABS(A(I,J)).LE.ABS(C)) GOTO 120
C = A(I,J)
I0= I
J0= J
120 CONTINUE
130 CONTINUE
IF(ABS(C).GT.EPS) GOTO 140
WRITE(*,'(47H***ERROR IN *GASJ*THIS MATRAX IS SINGULAR !)')
STOP
140 IK(K) = I0
JK(K) = J0
DO 170 J=1,N
DO 150 L=1,K
IF(J.EQ.JK(L)) GOTO 170
150 CONTINUE
IF(ABS(A(I0,J)).LT.1.E-5) GOTO 170
D = A(I0,J)/C
A(I0,J) = D
DO 160 I=1,N
IF(I.EQ.I0) GOTO 160
A(I,J) = A(I,J) - A(I,J0)*D
160 CONTINUE
170 CONTINUE
D = F(I0)/C
F(I0) = D
DO 180 I=1,N
IF(I.EQ.I0) GOTO 180
F(I) = F(I) - A(I,J0)*D
180 CONTINUE
DO 190 I=1,N
190 A(I,J0) = 0.
200 CONTINUE
DO 210 K=1,N
210 A(JK(K),1) = F(IK(K))
DO 220 K=1,N
220 F(K) = A(K,1)
RETURN
END
[ 本帖最后由 欧阳中华 于 2010-7-31 14:33 编辑 ]
页:
[1]