# More stable method of back substitution?

I've been tinkering a little in Fortran (2008) and wrote the following to solve $$Rx=b$$ for $$R\in\mathbb{R}^{n\times n}$$ upper-triangular, $$x,b\in\mathbb{R}^n$$.

My code looks like this:

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


I tested with random $$R$$ (upper triangular) and random $$b$$. Works great but for more than n=~50 the error starts going crazy. I could keep the error down by enforcing some diagonal dominance, but just wondering if there is a better algorithm that control the error better for general matrices?

I realise random matrices are not 'general' but this feature got me worried

• What did you compare your computed solution with to assess its accuracy, if you took random $R$ and $b$? May 21 '20 at 17:15
• I checked the norm2(matmul(R,x)-b) May 23 '20 at 10:08
• This is a bit difficult because actually I am implementing a back substitution of course to do solve upper triangular matrices such as generated from QR decomposition. So @rchilton1980 answer is meaningful. On the other hand Federico Poloni's answer tells me "No. This is as good as it gets". I was actually wondering if I had implemented this in a way that exaggerates round-off errors. Not sure which answer to accept here. May 23 '20 at 10:13
• That is usually called residual, not error. What do you mean exactly by 'going crazy'? Some growth in the residual with the dimension is to be expected, because $R$ and $b$ have more and larger entries. The recommended thing to check is the relative residual norm2(matmul(R,x)-b) / norm2(R)norm2(x), or also norm2(matmul(R,x)-b) / norm2(b) which should be of the same order unless $R$ and $b$ are chosen in a special way. This quantity should always be in the ballpark of machine precision (times a factor $n$ or so), irrespective of condition numbers. May 23 '20 at 11:53
• netlib.org/lapack/lapack-3.1.1/html/dtrsv.f.html
– vibe
May 23 '20 at 12:23