共轭梯度与带状矩阵向量乘法组合模块误差原因排查求助

我实现了两个Fortran模块：

• 共轭梯度法（conjugate gradient method）
• 带状矩阵-向量乘法（band matrix-vector multiplication）

两个模块单独运行时，仅因加法非结合性存在1e-14至1e-15级别的机器误差。但组合运行后，自定义matVecMult模块与常规矩阵向量乘法的结果差异显著，无法定位误差原因，恳请协助排查。

相关数据与模块代码 #

核心数据 #

• 稠密带状矩阵：dnsa
• 稀疏矩阵：a
• 非零对角线索引向量：

cols = [1 2 3 4 5 6 9]

CG求解器模块 #

module LinSolver
    use denseTimesVec
    implicit none
    contains
    subroutine cg(a,dnsa,n,b,cols,x)
        implicit none
        integer*8, intent (in) :: n
        real*8, intent (in) :: a(:,:)
        real*8, intent (in) :: dnsa(:,:)
        real*8, intent (inout):: b(:)
        real*8, intent (out) :: x(:)
        integer*8, intent (in) :: cols(:)


        ! local variables
        integer :: i, j, k
        real*8, dimension(n) :: r, p, Ap, ApDNS
        real*8, parameter :: eps = 1.e-6
        real*8 :: rsold, rsnew, alpha, pAp

        call matVecMult(n,dnsa,x,cols,ApDNS)
        do i = 1,n
            Ap(i) = 0
            do j = 1,n
                Ap(i) = Ap(i) + a(i,j) * x(j)
            enddo
        enddo
        print*,cols

        r = b - Ap
        p = r

        rsold = 0;
        do j = 1,n
            rsold = rsold + r(j)*r(j);
        enddo
        
        
        do i=1,n
            call matVecMult(n,dnsa,p,cols,ApDNS)
            do k = 1,n
                Ap(k) = 0;
                do j = 1,n
                    Ap(k) = Ap(k) + A(k,j) * p(j);
                enddo
            enddo
            
            pAp = 0;
            do j = 1,n
                pAp = pAp + Ap(j) * p(j);
            enddo
            alpha = rsold / pAp;
            x = x + alpha * p;
            r = r - alpha * Ap;
            rsnew = 0;
            do j = 1,n
                rsnew = rsnew + r(j)*r(j);
            enddo
            if (sqrt(rsnew)<eps) exit
            p = r + (rsnew / rsold) * p
            rsold = rsnew
        enddo
    end subroutine cg
end module LinSolver

带状矩阵-向量乘法模块 #

module denseTimesVec
    
    implicit none
    contains
    subroutine matVecMult(n,a,x,cols,b)
    
    implicit none        
    integer*8, intent (in) :: n
    real*8, intent (out) :: b(:)
    real*8, intent (in):: a(:,:), x(:)
    integer*8, intent (in) :: cols(:)
    integer*8   ::  i, j, width
    real*8, dimension (:), allocatable :: xTimesUpper, xTimesLower
    
    allocate ( xTimesUpper(n), xTimesLower(n)) 
    
    
    do j = 1,size(cols)
        width = cols(j)-1;
        do i = 1,n-width
            if (width == 0) then
                b(i) = b(i)+a(i,1)*x(i);
            else
                xTimesUpper = x(1+width:n);
                xTimesLower = x(1:n-width);
                
                ! multiplying upper/lower part
                b(i) = b(i) + a(i,j) * xTimesUpper(i); ! OK if compared to triuA*x-b
                b(i+width) = b(i+width) + a(i,j)*xTimesLower(i);
            endif
        enddo
    enddo
    
    end subroutine matVecMult
    
    end module denseTimesVec

内容的提问来源于stack exchange，提问作者2Napasa