subroutine d_backsubplus(matrix, invector, outvector, n, ier) implicit none integer n, ier double precision matrix(n,n), invector(n), outvector(n) integer i, k do i=1, n outvector(i)=invector(i) do k=1, i-1 outvector(i) = outvector(i) - matrix(i,k) * outvector(k) end do if (matrix(i,i) .eq. 0.0) then write (0,*) "d_backsubplus: matrix passed not invertible" ier=1 return end if outvector(i) = outvector(i) / matrix(i,i) end do ier=0 end ! subroutine d_backsubplus c ****************************************************************** subroutine d_backsubminus(matrix, invector, outvector, n, ier) implicit none integer n, ier double precision matrix(n,n), invector(n), outvector(n) integer i, k do i=n, 1, -1 outvector(i)=invector(i) do k=i+1, n outvector(i) = outvector(i) - matrix(i,k) * outvector(k) end do if (matrix(i,i) .eq. 0.0) then write (0,*) "d_backsubplus: matrix passed not invertible" ier=1 return end if outvector(i) = outvector(i) / matrix(i,i) end do ier=0 end ! subroutine d_backsubminus c ******************************************************************