Fortran Coder

标题: 求助三对角矩阵的追赶法 [打印本页]

作者: csycsy 时间: 2014-6-3 14:09
标题: 求助三对角矩阵的追赶法
哪位大神帮帮忙啊！！急啊

ac345982b2b7d0a2427f1160c9ef76094a369a9c.jpg (33.06 KB, 下载次数: 473)

作者: 楚香饭 时间: 2014-6-3 14:37
三对角矩阵，用追赶法。（取自徐世良的算法集）

其中，数据文件内容为：

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

作者: csycsy 时间: 2014-6-3 14:59

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

其中，数据文件内容为：

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

作者: 楚香饭 时间: 2014-6-3 15:17

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

是的

作者: csycsy 时间: 2014-6-3 16:33

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

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

作者: 楚香饭 时间: 2014-6-3 17:07

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

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

作者: csycsy 时间: 2014-6-3 17:11

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

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

作者: 楚香饭 时间: 2014-6-3 17:18

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

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

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

作者: 楚香饭 时间: 2014-6-3 17:21
解为：  0.753906250000000    0.515625000000000    0.382812500000000
  0.300781250000000    0.246093750000000    0.207495629370629
  0.178485576923077    0.159555288461538    0.128906250000000
  0.164062500000000

作者: csycsy 时间: 2014-6-3 17:22

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

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

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

作者: 楚香饭 时间: 2014-6-3 17:23

csycsy 发表于 2014-6-3 17:22
万分感谢！！不好意思！！！！来重庆请你吃饭！

我还真就在重庆，渝北区，两路街道。

作者: csycsy 时间: 2014-6-4 10:17

chuxf 发表于 2014-6-3 17:23
我还真就在重庆，渝北区，两路街道。

我在沙坪坝，你什么时候有空过来嘛！

作者: csycsy 时间: 2014-6-4 12:06

chuxf 发表于 2014-6-3 17:21
解为： 0.753906250000000 0.515625000000000 0.382812500000000
0.300781250000000 ...

你什么时候有空？

作者: 楚香饭 时间: 2014-6-4 12:18
别太热情了，我开玩笑而已。我家是重庆的，但是每年在重庆呆的时间也不多。

举手之劳，不必请客吃饭那么认真。常来论坛看看就是了~

作者: csycsy 时间: 2014-6-4 12:31

chuxf 发表于 2014-6-4 12:18
别太热情了，我开玩笑而已。我家是重庆的，但是每年在重庆呆的时间也不多。

举手之劳，不必请客吃饭那么认 ...

同为重庆人，。再次万分感谢！！！！

作者: aliouying 时间: 2014-6-6 16:25

csycsy 发表于 2014-6-3 17:22
万分感谢！！不好意思！！！！来重庆请你吃饭！

他就在重庆~

作者: 岸边的鱼 时间: 2014-6-9 22:09
为何设置为此帖子仅作者可见？？？？？？

作者: qqqqq 时间: 2018-11-4 17:18

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

其中，数据文件内容为：

我完全复制你的代码之后，运行时，无法获得计算结果，线程退出了

作者: 楚香饭 时间: 2018-11-4 21:12

qqqqq 发表于 2018-11-4 17:18
我完全复制你的代码之后，运行时，无法获得计算结果，线程退出了

你用的什么操作系统？什么编译器？
是否提供了输入文件？
在命令行下运行，看看错误提示。

作者: qqqqq 时间: 2018-11-6 00:35

楚香饭发表于 2018-11-4 21:12
你用的什么操作系统？什么编译器？
是否提供了输入文件？
在命令行下运行，看看错误提示。 ...

我用的是visual studio 2010的fortran90，编译器我也不太清楚
提供了输入文件，但不知道是不是因为文件里数字提供的计算机算不出来
提示是文件.exe触发断点

作者: fcode 时间: 2018-11-6 08:34
请看黑框错误提示。
（黑框，就是Hello，world输出的地方）

欢迎光临 Fortran Coder (http://bbs.fcode.cn/)