消主元解方程的程序问题

fcode

你的代码风格太古老了，应该更新一下了。

我对代码做了现代化风格的改造，主要错误是应该先对 u 操作，然后对 rm 操作。

给你的建议：
1.要学习 Fortran95 或以上的新语法
2.放弃 Compaq Fortran，放弃固定格式，放弃 Common！！

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

Program Main
  Implicit None
  Integer , Parameter :: nub0 = 5
  Real :: rm , x , u , rma(nub0,nub0)
  Common /a4/rm(nub0, nub0), x(nub0), u(nub0)
  call RANDOM_SEED()
  call RANDOM_NUMBER(rm)
  call RANDOM_NUMBER(u)
  rma = rm
  write(*,'("rm0:",/,5(5es14.6,/))') rma!//原始系数矩阵
  write(*,'("U0:",/,5es14.6)') u!//原始右端向量
  call ch( nub0 )  
  write(*,'("X:",/,5es14.6)') x!//解
  write(*,'("rm_up:",/,5(5es14.6,/))') rm!//上三角化
  write(*,'("U_Check:",/,5es14.6)') matmul( rma , reshape(x,[nub0,1]) )!//验证
End Program Main

Subroutine ch(z)
  Implicit None
  Integer , Parameter :: nub0 = 5
  Real :: rm , x , u
  integer :: i , z , j
  Common /a4/rm(nub0, nub0), x(nub0), u(nub0)
  Do i = 1, z
    u(i)=u(i) / rm(i,i)
    rm(i,i:z) = rm(i,i:z) / rm(i,i)    
    Do j = i + 1, z
      u(j)=u(j)-u(i)*rm(j,i)
      rm(j,i:z) = rm(j,i:z) - rm(i,i:z) * rm(j,i)     
    End Do
  End Do
  x(z) = u(z)/rm(z, z)
  Do i = z - 1, 1, -1    
    x(i) = u(i) - sum( x(i+1:z) * rm(i,i+1:z) )
  End Do
End Subroutine ch