高斯消元法求解线性齐次方程组

兰溪之水

在别人的基础上扩展加强，修正的，高斯消元法求解线性齐次方程组！！！！！

1. module LinearAlgebra
   
2. implicit none
   
3. contains
   
4. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
5. !!!!!!!!                     !!!!!!!!!
   
6. !!!!!!!! Gauss_Jordan法 !!!!!!!!!
   
7. !!!!!!!! !!!!!!!!!
   
8. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
9. subroutine Gauss_Jordan(A,S,ANS)
   
10. implicit none
    
11. real :: A(:,:)
    
12. real :: S(:)
    
13. real :: ANS(:)
    
14. real, allocatable :: B(:,:)
    
15. integer :: i, N
    
16. N = size(A,1)
    
17. allocate(B(N,N))
    
18. ! 保存原先的矩阵A,及数组S
    
19. B=A
    
20. ANS=S
    
21. ! 把B化成对角矩阵
    
22. call Upper(B,ANS,N) ! 把B化成上三角矩阵
    
23. call Lower(B,ANS,N) ! 把B化成下三角矩阵
    
24. 
    
25. ! 求解
    
26. forall(i=1:N)
    
27. ANS(i)=ANS(i)/B(i,i)
    
28. end forall
    
29. return
    
30. end subroutine Gauss_Jordan
    
31. ! 输出等式子程序
    
32. subroutine output(M,S)
    
33. implicit none
    
34. real :: M(:,:), S(:)
    
35. integer :: N,i,j
    
36. N = size(M,1)
    
37. ! write中加上advance="no",可以中止断行发生,使下一次的
    
38. ! write接续在同一行当中.
    
39. do i=1,N
    
40. write(*,"(1x,f5.2,a1)", advance="NO") M(i,1),'A'
    
41. do j=2,N
    
42. if ( M(i,j) < 0 ) then
    
43. write(*,"('-',f5.2,a1)",advance="NO") -M(i,j),char(64+j)
    
44. else
    
45. write(*,"('+',f5.2,a1)",advance="NO") M(i,j),char(64+j)
    
46. end if
    
47. end do
    
48. write(*,"('=',f8.4)") S(i)
    
49. end do
    
50. return
    
51. end subroutine output
    
52. ! 求上三角矩阵的子程序
    
53. subroutine Upper(M,S,N)
    
54. implicit none
    
55. integer :: N
    
56. real :: M(N,N)
    
57. real :: S(N)
    
58. integer :: I,J
    
59. real :: E
    
60. do I=1,N-1
    
61. do J=I+1,N
    
62. E=M(J,I)/M(I,I)
    
63. M(J,I:N)=M(J,I:N)-M(I,I:N)*E
    
64. S(J)=S(J)-S(I)*E
    
65. end do
    
66. end do
    
67. return
    
68. end subroutine Upper
    
69. ! 求下三角矩阵的子程序
    
70. subroutine Lower(M,S,N)
    
71. implicit none
    
72. integer :: N
    
73. real :: M(N,N)
    
74. real :: S(N)
    
75. integer :: I,J
    
76. real :: E
    
77. do I=N,2,-1
    
78. do J=I-1,1,-1
    
79. E=M(J,I)/M(I,I)
    
80. M(J,1:N)=M(J,1:N)-M(I,1:N)*E
    
81. S(J)=S(J)-S(I)*E
    
82. end do
    
83. end do
    
84. return
    
85. end subroutine Lower
    
86. end module
    
87. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
88. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
89. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
90. !!!!!!!!! !!!!!!!!!
    
91. !!!!!!!!! 求解齐次线性方程组 !!!!!!!!!
    
92. !!!!!!!!! !!!!!!!!!
    
93. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
94. program main
    
95. use LinearAlgebra
    
96. implicit none
    
97. integer:: N ! Size of Matrix
    
98. real, allocatable :: A(:,:),S(:),ans(:)
    
99. integer :: i
    
100. write(*,*)"请输入要求方程组的元数:"
     
101. read(*,*) N
     
102. allocate(A(N,N))
     
103. do i=1,N
     
104. write(*,"(A6,I2,A30)")"输入第",i,"条方程变量的系数(用空格隔开)："
     
105. read(*,*) A(i,:)
     
106. end do
     
107. allocate(S(N))
     
108. write(*,"(A6,I2,A36)")"输入第",i,"条方程等号右边常数系数(用空格隔开)："
     
109. read(*,*) S
     
110. allocate(ans(N))
     
111. 
     
112. write(*,*) 'Equation:'
     
113. call output(A,S)
     
114. call Gauss_Jordan(A,S,ANS)
     
115. write(*,*) 'Ans:'
     
116. do i=1,N
     
117. write(*,"(1x,a1,'=',F8.4)") char(64+i),ANS(i)
     
118. end do
     
119. stop
     
120. end program
     
121. 
     

复制代码

mofangbao

哈哈 不会是被我逼迫出来的吧

天目神眉

哈哈哈，这个是不是和那个追赶法有点类似，貌似以前小研究过，但愿不如兰溪版主这么高深

兰溪之水

天目神眉 发表于 2011-10-22 12:19

                               
登录/注册后可看大图

哈哈哈，这个是不是和那个追赶法有点类似，貌似以前小研究过，但愿不如兰溪版主这么高深

好吧，但愿。。。

兰溪之水

mofangbao 发表于 2011-10-22 11:15

                               
登录/注册后可看大图

哈哈 不会是被我逼迫出来的吧

呃。。。算是吧，最近在调整状态啊~

天目神眉

兰溪之水 发表于 2011-10-22 12:36

                               
登录/注册后可看大图

好吧，但愿。。。

，写错字鸟，是“但是”

尽头的尽头

兰溪给力啊

兰溪之水

尽头的尽头 发表于 2011-10-22 12:50

                               
登录/注册后可看大图

兰溪给力啊

因为我站在巨人的头上~

尽头的尽头

兰溪之水 发表于 2011-10-22 13:05

                               
登录/注册后可看大图

因为我站在巨人的头上~

兰溪谦虚了，哈哈

心语