c_invert_tridiag_block_matrix Subroutine

subroutine c_invert_tridiag_block_matrix(Nb,N,sub_,diag_,over_,Ainv)
  integer                        :: ib,i,j,Nb,N
  complex(8),dimension(Nb,N,N)   :: diag_
  complex(8),dimension(Nb-1,N,N) :: sub_
  complex(8),dimension(Nb-1,N,N) :: over_
  complex(8),dimension(Nb,N,N)   :: Ainv
  complex(8),dimension(N,N)      :: Foo,Cleft,Cright
  complex(8),dimension(Nb,N,N)   :: Dleft
  complex(8),dimension(Nb,N,N)   :: Dright
  !
  !DOWNWARD:
  dleft(1,:,:) = diag_(1,:,:)
  do i=2,Nb
     foo  = dleft(i-1,:,:) ; call inv(foo)
     cleft= matmul(sub_(i-1,:,:),foo)
     dleft(i,:,:) = diag_(i,:,:) - matmul(cleft,over_(i-1,:,:))
  enddo
  !
  !BACKWARD:
  dright(Nb,:,:) = diag_(Nb,:,:)
  do i=Nb-1,1,-1
     foo   = dright(i+1,:,:) ; call inv(foo)
     cright= matmul(over_(i,:,:),foo)
     dright(i,:,:) = diag_(i,:,:) - matmul(cright,sub_(i,:,:))
  enddo
  !
  do ib=1,Nb
     Ainv(ib,:,:)    =  dleft(ib,:,:) + dright(ib,:,:) - diag_(ib,:,:)
     call inv(Ainv(ib,:,:))
  end do
end subroutine c_invert_tridiag_block_matrix