c_invert_tridiag_matrix Subroutine

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

Arguments

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

Source Code

subroutine c_invert_tridiag_matrix(N,sub_,diag_,over_,Inv)
  integer                            :: i,j,N
  complex(8),dimension(N)            :: diag_
  complex(8),dimension(N-1)          :: sub_
  complex(8),dimension(N-1)          :: over_
  complex(8),dimension(N)            :: Inv
  complex(8)                         :: Foo,Cleft,Cright
  complex(8),dimension(N)            :: dleft
  complex(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) !over_(i-1)/dleft(i-1)*sub_(i)
     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) !sub_(i+1)/dright(i+1)*over_(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 c_invert_tridiag_matrix