I received a code to work on and it has a QR decomposition subroutine as shown below. When checking this QR decomposition using the error calculation I added at the end of it, I found the error increases as the simulation continues to very high magnitudes after hundreds of simulation steps.

The QR decomposition is used to solve the coefficients used in triquadratic interpolation, adapted from trilinear interpolation as found on Wikipedia: https://en.wikipedia.org/wiki/Trilinear_interpolation. In the code, we generate a 5x5x5 grid of values, and uses triquadratic interpolation to derive a spatial function form of the values to obtain f(x,y,z) = a_0 + a_1x + a_2y + a_3*z + …

The function f(x,y,z) can be written as Ax = b for the pre-existing values on each grid point, such that b is the 125×1 vector of existing values, A is a 125×20 matrix calculated using the coordinates, and x is the 20×1 coefficient vector composed of (a_0, a_1, a_2, …).

Then, x is calculated using the QR decomposition algorithm shown below. When I looked up QR decomposition, I found some discussions here that mentioned certain approaches can become increasingly unstable, but approaches such as householder reflection should be very stable.

My first question: I am having trouble identifying the type of approach used in my code, it doesn’t normalize the variables and I can’t match the process to other implementations found online. Can anyone tell me if this implementation is already a householder reflection?

My second question: If the implementation is correct and should be stable, can anyone give me some ideas to check why the solution would become increasingly unstable? Or if it is not householder or simply a bad implementation, would an example householder implementation that I found online help improve the accuracy?

  ! QR decomposition solving for vector B in AX = B
  subroutine qr_eqsystem(A, B, X, M, N, error)
  
  implicit none
  
  integer, intent(in) :: m, n
  real, intent(in) :: A(m, n), B(m)
  
  real, intent(out) :: x(n)
  real, intent(out) :: error
  
  !locals
  real :: q(m, n), r(n, n), aux
  integer :: i, j, k
  
  do i = 1,n
    do k = 1,m
      q(k, i) = A(k, i)
    enddo
  enddo
  
  do i = 1,n
    do j = 1,n
      r(j, i) = 0.0
    enddo
  enddo
  
  do i = 1,n
    do j = 1,i-1
      r(j, i) = 0.0
      do k = 1,m
        r(j, i) = r(j, i) + q(k, j)*A(k, i)
      enddo
      
      do k = 1,m
        q(k, i) = q(k, i) - r(j, i)*q(k,j)
      enddo
    enddo
    
    r(i, i) = 0.0
    do k = 1,m
      r(i, i) = r(i, i) + q(k, i)*q(k, i)
    enddo
    r(i, i) = sqrt(r(i, i))
    
    if(r(i,i).ne.0.0)then
      do k = 1,m
        q(k, i) = q(k, i) / r(i, i)
      enddo
    endif
    
  enddo
  
  do i = 1,n
    x(i) = 0.0
    do j = 1,m
      x(i) = x(i) + q(j, i)*B(j)
    enddo
  enddo
  
  if(r(n, n).ne.0.0)then
    x(n) = x(n) / r(n,n)
  endif
  
  do i = n-1,1,-1
    do j = i+1,n
      x(i) = x(i) - r(i, j) * x(j)
    enddo
    
    if(r(i, i).ne.0.0)then
       x(i) = x(i) / r(i, i)
    endif
  enddo
  
  error = 0
  
  ! calculate AX and compare with B as an error, report total raw error
  do i = 1,m
    aux = 0.0
    do j = 1,n
      aux = aux + A(i,j)*X(j)
    enddo
    error = error + abs(aux - B(i))
  enddo
  
  endsubroutine qr_eqsystem