求助三对角矩阵的追赶法

楚香饭

三对角矩阵，用追赶法。（取自徐世良的算法集）

其中，数据文件内容为：

32 -8
12 -22 14
14 -22 16
16 -22 18
18 -22 20
20 -22 22
22 -22 24
24 -22 26
26 -22 28
28 -22
20 0 0 0 0 0 0 0 0 0

解为

解为：  0.490707967461251      -0.537168130154994       -1.26472817663892
  -1.26897912899290      -0.426771667312277       0.672632382050108
   1.06060662506127       0.355643056093564      -0.678092760285062
-0.863027149453715

以下是代码：

[Fortran] 纯文本查看 复制代码

Program www_fcode_cn
  Implicit None
  Integer , parameter :: N = 10 !// 10*10 矩阵
  Integer , parameter :: M = 3*N-2
  Integer , parameter :: DP = Kind(0.0D0)
  Real( Kind = DP ) :: B(M) , D(N)
  integer :: L
  Open( 12 , File = "B.txt" )
  read( 12 , * ) B(:) !// 三对角系数
  read( 12 , * ) D(:) !// 右端向量
  call ATRDE( B , N , M , D , L )
  write(*,*) '解为：',D(:)
End Program www_fcode_cn

SUBROUTINE ATRDE(B,N,M,D,L)
  DIMENSION B(M),D(N)
  DOUBLE PRECISION B,D
  L=1
  IF (M.NE.(3*N-2)) THEN
    L=-1
    WRITE(*,10)
    RETURN
  END IF
10  FORMAT(1X,'  ERR  ')
  DO 20 K=1,N-1
    J=3*K-2
    IF (ABS(B(J))+1.0.EQ.1.0) THEN
      L=0
      WRITE(*,10)
      RETURN
    END IF
    B(J+1)=B(J+1)/B(J)
    D(K)=D(K)/B(J)
    B(J+3)=B(J+3)-B(J+2)*B(J+1)
    D(K+1)=D(K+1)-B(J+2)*D(K)
20  CONTINUE
  IF (ABS(B(3*N-2))+1.0.EQ.1.0) THEN
    L=0
    WRITE(*,10)
    RETURN
  END IF
  D(N)=D(N)/B(3*N-2)
  DO 30 K=N-1,1,-1
    D(K)=D(K)-B(3*K-1)*D(K+1)
30  CONTINUE
  RETURN
  END