【求助】求解带形线性方程组的高斯消元法源代码

kuangtq

我采用二维数组保存整个刚度矩阵元素，进行高斯消元求解，能够得到正确解，但耗时较长（400X400的刚阵）。因此改为半带宽格式存储，高斯方法也进行了相应的改正，计算过程出现问题。不知哪位大大有成熟的以半带宽格式存储的高斯消元法源代码？或大家能不能帮忙看看这个程序何处有问题。谢谢！
  
!***************************************!
!       形成P总体刚度矩阵      !
!***************************************!
  
  SUBROUTINE GLOBAL_P_STIFFNESS() !  
  
    USE DATA_MODULE
  
    INTEGER EI
  
    GLB_P_STIF=0
    GLB_P_R=0
  
!    ! 采用等带宽格式存储  
    DO EI=1, NELEM
      CALL ELEM_P_STIFFNESS(EI)            ! 单元有限元特征式分析
      DO I=1, 3
        II=EL_NODE(I, EI)
        GLB_P_R(II)=GLB_P_R(II)+ELE_R(I)      ! 右端项的集成
        DO J=1, 3
          JJ=EL_NODE(J, EI)
          IF(II<=JJ)THEN
            JJ=JJ-II+1              ! 新列号
            GLB_P_STIF(II, JJ)=GLB_P_STIF(II, JJ)+ELE_STIF(I, J)  ! 形成总刚
          ENDIF
        ENDDO
      ENDDO
    ENDDO
  
  END
  
!***************************************!
!        施加压力边界条件       !
!***************************************!
  
  SUBROUTINE IMPLEMENT_P_BOUNDARY()    
    USE DATA_MODULE
                      
    ! 等带宽存储格式
  
    GLB_P_STIF(NODE_ENTRY, 1)=GLB_P_STIF(NODE_ENTRY, 1)*1E20    ! 入口节点 给定非零压力，对角元素乘大数
    GLB_P_R(NODE_ENTRY)=GLB_P_STIF(NODE_ENTRY, 1)*P_ENTRY      ! 右端项
  
    DO I=1, NNODE
      IF(F(1, I)<1) THEN            ! P(I)=0    熔体前沿 气体区压力为0      
        DO J=1, HB_P
          GLB_P_STIF(I, J)=0.0      ! 修改行
        ENDDO
        DO K=1, HB_P         ! HB_P为半带宽
          M=I-(K-1)            ! K=I-M+1  K为新列号, I为原列号，M为行号    
          IF(M>=1) THEN          
            GLB_P_STIF(M, K)=0.0    ! 修改列  (全局矩阵中I列全改为0，实际在带宽内I列只有HB_P个元素)
          ENDIF
        ENDDO
        GLB_P_STIF(I, 1)=1.0        ! 对角元素改1
        GLB_P_R(I)=0            ! 熔体前沿气体区压力为0
      ENDIF
    ENDDO
  
    DO I=1, NNODE                ! 检查总刚是否正定
      IF(GLB_P_STIF(I, 1)<=0)THEN
        WRITE(6, "(/10X, 'MAIN ELEMENT IS NONPOSITIVE NNODE=', I4, 5X, E12.6)") I, GLB_P_STIF(I, 1)   
        STOP
      ENDIF
    ENDDO
  
  END
  
!*******************************************!
!       带宽存储高斯消去法        !
!*******************************************!
  
  SUBROUTINE GAUS_BAND(A,B,X,N,D)    ! 参见《有限单元法基本原理和数值方法》p201
  
    DIMENSION A(N,N),B(N),X(N)
    DOUBLE PRECISION A,B,X
    INTEGER N, D, D0
  
    DO M=1, N-1                        ! 消去  
      M1=MIN(D-1, N-M)
      DO L=1, M1
        C=A(M, L+1)/A(M, 1)
        IF(ABS(C)>=1E-18) THEN
          DO J=1, D-L
            A(M+L, J)=A(M+L, J)-C*A(M, J+L)      ! 消去A
          ENDDO
          B(M+L)=B(M+L)-C*B(M)            ! 消去B
        ENDIF
      ENDDO
    ENDDO
  
    B(N)=B(N)/A(N, 1)                    ! 回代
  
    DO K=1, N-1
      M=N-K
      D0=MIN(D, K+1)                    ! K+1=N-M+1
      DO J=2, D0
        B(M)=B(M)-A(M, J)*B(M+J-1)
      ENDDO
      B(M)=B(M)/A(M, 1)
    ENDDO
  
    X=B
  
  END

kuangtq

问题已解决！

yexm

为什么不试试用变带宽的方法呢？