Fortran Coder

查看: 4605|回复: 1

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

[复制链接]

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

63

帖子

9

主题

0

精华

专家

超凡脱俗

F 币
474 元
贡献
237 点
发表于 2019-10-31 13:32:48 | 显示全部楼层
本帖最后由 Jackdaw 于 2019-10-31 13:34 编辑

[Fortran] 纯文本查看 复制代码
    !> @details 下三角矩阵求逆(稠密存储)
    !> @author Jackdaw毛毛
    !> @date 2019 年 8 月 26 日
    subroutine LowerTriMatrixInv(mat,matInv)
        real(fp),intent(in)    :: mat(:,:)    !< 输入矩阵,下三角,稠密存储
        real(fp),intent(inout) :: matInv(:,:) !< 输出矩阵,下三角,稠密存储,与输入矩阵关系:互为逆矩阵
        integer :: i,j,n 
        n = size(mat,dim=1)  

        matInv = 0.0_fp 
        do i = 1, n 
            matInv(i,i) = 1.0_fp / mat(i,i) 
        end do 
        do i = 2, n 
            do j = 1, i-1
                !do k = j, i-1 
                !   matInv(i,j) = matInv(i,j) + mat(i,k)*matInv(k,j)
                !end do 
                matInv(i,j) = dot_product( mat(i,j:i-1),matInv(j:i-1,j) )
                matInv(i,j) = matInv(i,j) * (-matInv(i,i)) 
            end do 
        end do 
    end subroutine 
    !> @details 高斯赛德尔迭代法计算 \f$ Ax=b \f$,收敛速度比雅克比迭代更快。
    !> @author Jackdaw毛毛
    !> @date 2019 年 8 月 26 日
    subroutine useGaussSeidelIterSolveAx_eq_b(matrixA,vectorB,vectorX)
        real(fp),intent(in)    :: matrixA(:,:) !< 系数矩阵(迭代收敛的充分条件:严格对角占优或对称正定)
        real(fp),intent(in)    :: vectorB(:)   !< 右端项
        real(fp),intent(inout) :: vectorX(:)   !< 解向量 
        
        real(fp),allocatable   :: G(:,:),tmpM(:,:),tmpB(:),tmpX(:)
        integer ::i,n 
        
        n = size(matrixA,dim=1)
        allocate(G(n,n),tmpM(n,n),tmpX(n),tmpB(n))
        ! 迭代格式
        ! x(k+1) = D^-1(b-Ux(k)-Lx(k+1)) 
        ! x(k+1) = (L+D)^-1(-Ux(k) + b )
        ! 构造迭代矩阵
        tmpM = 0.0_fp 
        do i = 1, n 
            tmpM(i:n,i) = matrixA(i:n,i) ! L+D
        end do 
        call LowerTriMatrixInv(tmpM,G)   ! 求逆,存在G里面
        tmpM = 0.0_fp 
        do i = 1, n-1
            tmpM(i,i+1:n) = matrixA(i,i+1:n) ! U
        end do 
        tmpB = matmul(G,vectorB) ! 常向量
        G = matmul(G,tmpM) ! 迭代矩阵
        deallocate(tmpM) 
        do i = 1, 1000
            tmpX = tmpB - matmul(G,vectorX) 
            if( maxval(abs(vectorX-tmpX)) .lt. 1.e-10_fp ) then 
                vectorX = vectorX + tmpX
                vectorX = vectorX / 2.0_fp
                exit
            else
                vectorX = tmpX
            endif 
        end do 
        if(i .ge. 1000 ) write(*,*) "GaussSeidel-iteration failed."
        deallocate( G,tmpX,tmpB )
    end subroutine  
    
    !> @details 高斯赛德尔迭代法计算 \f$ Ax=b \f$,测试算例。
    !> @author Jackdaw毛毛
    !> @date 2019 年 8 月 26 日
    subroutine useGaussSeidelIterSolveAx_eq_b_TEST
        real(fp) :: a(2,2),b(2),x(2)
        a(1,:) = [ 3.0_fp, 1.0_fp ]
        a(2,:) = [ 1.0_fp, 2.0_fp ] 
        b      = [ 5.0_fp, 5.0_fp ]
        x      = 0.0_fp
        call useGaussSeidelIterSolveAx_eq_b(a,b,x) 
        write(*,*) "解析解 x1 = 1.0, x2 = 2.0"
        write(*,*) "数值解",x
    end subroutine useGaussSeidelIterSolveAx_eq_b_TEST 


天下英雄出我辈,一入江湖岁月催。

鸿图霸业谈笑间,不胜人生一场醉。
您需要登录后才可以回帖 登录 | 极速注册

本版积分规则

捐赠本站|Archiver|关于我们 About Us|小黑屋|Fcode ( 京ICP备18005632-2号 )

GMT+8, 2024-3-28 20:42

Powered by Tencent X3.4

© 2013-2024 Tencent

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