Fortran Coder

查看: 190|回复: 0

[线性代数] 高斯-赛德尔迭代解线性方程组的问题

[复制链接]

6

帖子

3

主题

0

精华

入门

F 币
41 元
贡献
22 点
发表于 2019-6-27 08:41:01 | 显示全部楼层 |阅读模式
本帖最后由 鹰击长空 于 2019-6-27 08:45 编辑

两个高斯-赛德尔迭代解线性方程组子程序,算法一样,只是单精度和双精度的区别。麻烦帮忙看看是哪里的问题

单精度
[Fortran] 纯文本查看 复制代码
    subroutine GS_solver (A,b,n,x)
        implicit none
        
        integer,parameter :: step = 5000  !最大迭代步数
        real   ,parameter :: err  = 1e-7 !误差
        
        integer :: i,j,k,n
        real    :: A(n,n),b(n),x0(n),x(n)
        real    :: x1(n)
        real    :: s,dx
        
        x0 = 0
        x1 = x0
        x  = x0

        
        do k = 1,step
            dx = 0
            
            do i = 1,n
               x1 = x0
               x1(i) = 0
               s = 0 
               s = dot_product( A(i,:), x1 )
               
               x(i) = ( b(i) - s ) / A(i,i) 
               
               dx = dx + ( x0(i)-x(i) )**2
               x0(i) = x(i) 
            end do
            
            dx = sqrt( dx )
            if ( dx < err ) exit
            
            
            x0 = x

    end do
    
 end subroutine GS_solver


双精度
[Fortran] 纯文本查看 复制代码
    subroutine GSeidelSolver (A,b,n,x)
        implicit none
        
        integer   ,parameter :: MaxStep = 9000   !最大迭代步数
        real(8)   ,parameter :: cov  = 1d-8      !误差
        integer  :: i, j, k, n, step
        real(8)  :: A(n,n), b(n), x(n)
        real(8)  :: x0(n), x1(n)
        real(8)  :: dx, err
        
        x0 = 0d0
        x1 = 0d0
        x  = 0d0
        
        LOOP: do step = 1, MaxStep
            
            dx = 0d0
            do i = 1,n

                x1 = x0
                
                x1(i) = 0
                
                k = 0d0
                k = dot_product( A(i,:), x1 )
                
                x(i) = ( b(i) - k ) / A(i,i)
                dx = dx + ( x(i)-x0(i) )**2
                
                x0(i) = x(i) 
         
           end do
         
        err = dsqrt( dx )
                        
       if ( err < cov )  then
                exit LOOP
        end if

        x0 = x          

  end do LOOP         

 end subroutine GSeidelSolver
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

QQ|捐赠本站|Archiver|关于我们 About Us|群聊|Fcode

GMT+8, 2019-10-17 03:12

Powered by Discuz! X3.2

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表