d_invert_tridiag_matrix Subroutine

subroutine d_invert_tridiag_matrix(N, sub_, diag_, over_, Inv)

Arguments

Type IntentOptional Attributes Name
integer :: N
real(kind=8), dimension(N-1) :: sub_
real(kind=8), dimension(N) :: diag_
real(kind=8), dimension(N-1) :: over_
real(kind=8), dimension(N) :: Inv

Source Code

subroutine d_invert_tridiag_matrix(N,sub_,diag_,over_,Inv)
  integer                         :: i,j,N
  real(8),dimension(N)            :: diag_
  real(8),dimension(N-1)          :: sub_
  real(8),dimension(N-1)          :: over_
  real(8),dimension(N)            :: Inv
  real(8)                         :: Foo,Cleft,Cright
  real(8),dimension(N)            :: dleft
  real(8),dimension(N)            :: dright
  !
  !DOWNWARD:
  dleft(1) = diag_(1)
  if(dleft(1)==0d0)stop "matrix is ill-conditioned: no inverse exists"
  do i=2,N
     foo      = 1d0/dleft(i-1)
     cleft    = sub_(i-1)*foo
     dleft(i) = diag_(i) - cleft*over_(i-1)
     if(dleft(i)==0d0)stop "matrix is ill-conditioned: no inverse exists"
  enddo
  !
  !UPWARD:
  dright(N) = diag_(N)
  if(dright(N)==0d0)stop "matrix is ill-conditioned: no inverse exists"
  do i=N-1,1,-1
     foo       = 1d0/dright(i+1)
     cright    = over_(i)*foo
     dright(i) = diag_(i) - cright*sub_(i)
     if(dright(i)==0d0)stop "matrix is ill-conditioned: no inverse exists"
  enddo
  !
  do i=1,N
     Foo    =  dleft(i) + dright(i) - diag_(i)
     Inv(i) = 1d0/Foo
  end do
end subroutine d_invert_tridiag_matrix