tri_sle.F90

Go to the documentation of this file.

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!  Subroutine :  t r i _ s l e
!
!> @file
!!
!! Solution of a system of linear equations Ax=b with tridiagonal matrix A.
!!
!! @section Copyright
!!
!! Copyright 2009-2013 Ralf Greve
!!
!! @section License
!!
!! This file is part of SICOPOLIS.
!!
!! SICOPOLIS is free software: you can redistribute it and/or modify
!! it under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! SICOPOLIS is distributed in the hope that it will be useful,
!! but WITHOUT ANY WARRANTY; without even the implied warranty of
!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!! GNU General Public License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with SICOPOLIS.  If not, see <http://www.gnu.org/licenses/>.
!<
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!-------------------------------------------------------------------------------
!> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
!! @param[in]  a0       a0(j) is element A_(j,j-1) of Matrix A
!! @param[in]  a1       a1(j) is element A_(j,j)   of Matrix A
!! @param[in]  a2       a2(j) is element A_(j,j+1) of Matrix A
!! @param[in]  b        inhomogeneity vector
!! @param[in]  nzeilen  size of matrix A (indices run from 0 (!!!) to nzeilen)
!! @param[out] x        Solution vector.
!<------------------------------------------------------------------------------
subroutine tri_sle(a0, a1, a2, x, b, nzeilen)

use sico_types

implicit none

real(dp), dimension(0:*) :: a0, a1, a2, x, b
real(dp), allocatable, dimension(:) :: help_x
integer(i4b) :: nzeilen
integer(i4b) :: n

!--------  Generate an upper triangular matrix
!                      ('obere Dreiecksmatrix') --------

do n=1, nzeilen
   a1(n)   = a1(n) - a0(n)/a1(n-1)*a2(n-1)
end do

do n=1, nzeilen
   b(n)    = b(n) - a0(n)/a1(n-1)*b(n-1)
!         a0(n)  = 0.0_dp , not needed in the following, therefore
!                           not set
end do

!-------- Iterative solution of the new system --------

!      x(nzeilen) = b(nzeilen)/a1(nzeilen)

!      do n=nzeilen-1, 0, -1
!         x(n) = (b(n)-a2(n)*x(n+1))/a1(n)
!      end do

allocate(help_x(0:nzeilen))

help_x(0) = b(nzeilen)/a1(nzeilen)

do n=1, nzeilen
   help_x(n) = b(nzeilen-n)/a1(nzeilen-n) &
              -a2(nzeilen-n)/a1(nzeilen-n)*help_x(n-1)
end do

do n=0, nzeilen
   x(n) = help_x(nzeilen-n)
end do

!       (The trick with the help_x was introduced in order to avoid
!        the negative step in the original, blanked-out loop.)

deallocate(help_x)

!  WARNING: Subroutine does not check for elements of the main
!           diagonal becoming zero. In this case it crashes even
!           though the system may be solvable. Otherwise ok.

end subroutine tri_sle
!