求助三对角矩阵的追赶法

csycsy

哪位大神帮帮忙啊！！急啊

楚香饭

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

其中，数据文件内容为：

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] syntaxhighlighter_viewsource syntaxhighlighter_copycode

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

csycsy

chuxf 发表于 2014-6-3 14:37
三对角矩阵，用追赶法。（取自徐世良的算法集）

其中，数据文件内容为：

谢谢啊，你这是先把数字文件存为B。txt吗？然后再用这个代码解吗？

楚香饭

csycsy 发表于 2014-6-3 14:59
谢谢啊，你这是先把数字文件存为B。txt吗？然后再用这个代码解吗？

是的

csycsy

chuxf 发表于 2014-6-3 15:17
是的

大神大神，错了错了，，，不好意思啊，那个数据文本为32 -8
12 -22 6
14 -22 4
16 -22 2
18 -22 0
20 -22 -2
22 -22 -4
24 -22 -6
26 -22 -8
28 -22
20 0 0 0 0 0 0 0 0 0


能不能帮我运行下结果？？？？？？？？？？？万分感谢！！！！！！

楚香饭

csycsy 发表于 2014-6-3 16:33
大神大神，错了错了，，，不好意思啊，那个数据文本为32 -8
12 -22 6
14 -22 4

你这要求过分了点。我不是你的保姆，自己运行吧。

csycsy

chuxf 发表于 2014-6-3 17:07
你这要求过分了点。我不是你的保姆，自己运行吧。

大神啊，主要我不是这个专业的，隔行如隔山啊！！万分感谢啊@@！！

楚香饭

csycsy 发表于 2014-6-3 17:11
大神啊，主要我不是这个专业的，隔行如隔山啊！！万分感谢啊@@！！

中国没有大学开设 Fortran 专业。

1.你一开始强调需要源代码。现在给你源代码了，你又不用，改而要结果。这是对回答者的不尊重。
2.这里是讨论技术问题的，不是免费劳动力市场。

楚香饭

解为：  0.753906250000000       0.515625000000000       0.382812500000000
  0.300781250000000       0.246093750000000       0.207495629370629
  0.178485576923077       0.159555288461538       0.128906250000000
  0.164062500000000

csycsy

chuxf 发表于 2014-6-3 17:18
中国没有大学开设 Fortran 专业。

1.你一开始强调需要源代码。现在给你源代码了，你又不用，改而要结果 ...

万分感谢！！不好意思！！！！来重庆请你吃饭！