我在 Fortran (2008) 中进行了一些修改并编写了以下内容来解决 上三角形 ,。$Rx=b$$R\in\mathbb{R}^{n\times n}$$x,b\in\mathbb{R}^n$

我的代码如下所示：

function backsub_upper(R, b) result(x)

    ! R is a square upper triangular matrix
    ! b is the vector to be solved for
    ! x is the solution to Rx=b

    real(kind=wp), dimension(:,:), intent(in)   :: R
    real(kind=wp), dimension(:) , intent(in)    :: b

    real(kind=wp), dimension(:)                 :: x(size(b)), bc(size(b))
    integer(kind=sp)                            :: i,j

    x = 0.0_wp
    bc(:) = b(:)
    write(*,*)"solving"
    do j=1,size(bc)
        i = size(bc)+1-j
        x(i) = bc(i)/R(i,i)

        if (i/=1) bc(1:i-1) = bc(1:i-1) - x(i)*R(1:i-1,i)
    end do
    write(*,*)"done"

end function backsub_upper

我用随机（上三角）和随机进行了测试。效果很好，但超过 n=~50 时，错误开始变得疯狂。我可以通过强制执行一些对角优势来降低误差，但只是想知道是否有更好的算法可以更好地控制一般矩阵的误差？$R$$b$

我意识到随机矩阵不是“通用的”，但这个功能让我担心