结构有限元一维变带宽存储LD分解及求解代码

欧阳中华

.
      考虑到编写结构有限元程序的参考，给出结构有限元一维变带宽存储LD分解及求解代码：

C
C     This subroutine is used to solve finite element static equations  
C  in core , using compacted storage and column reduction scheme .   
C                                                                     
C  INPUT VARIABLES                                                     
C             A[NWK] => Stiffness matrix stored in compacted form         
C                V[NN] => Right-hand-side load vecter                        
C     MAXA[NN+1] => Vector dontaining addresses of diagonal elements   
C                   of stiffness matrix in  A                          
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[NWK] => D and L - factors of stiffness matrix              
C            V[NN] => 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 编辑 ]