!Subroutine for calculating the inverse of a matrix.

subroutine QR_INVERSE(Q, R, m, Z)
implicit none

     integer, parameter :: dbl=SELECTED_REAL_KIND(14)
     real(KIND=dbl) :: Q(m,m), Q_T(m,m), R(m,m), Z(m,m), QQ_T(m,m)
     integer :: m, i, j, k

Q_T = TRANSPOSE(Q)

Z = 0.0

do i = m, 1, -1
   do j = 1, m
      do k = i+1, m
         Z(i,j) = Z(i,j) - R(i,k)*Z(k,j)/R(i,i)
      enddo
   Z(i,j) = Z(i,j) + Q_T(i,j)/R(i,i)
   enddo
enddo

end subroutine QR_INVERSE