MAIC-2  Revision 19
 All Classes Files Functions Variables
tri_sle.F90
Go to the documentation of this file.
1 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 !
3 ! Subroutine : t r i _ s l e
4 !
5 !> @file
6 !!
7 !! Solution of a system of linear equations Ax=b with tridiagonal matrix A.
8 !!
9 !! @section Copyright
10 !!
11 !! Copyright 2010-2013 Ralf Greve, Bjoern Grieger, Oliver J. Stenzel
12 !!
13 !! @section License
14 !!
15 !! This file is part of MAIC-2.
16 !!
17 !! MAIC-2 is free software: you can redistribute it and/or modify
18 !! it under the terms of the GNU General Public License as published by
19 !! the Free Software Foundation, either version 3 of the License, or
20 !! (at your option) any later version.
21 !!
22 !! MAIC-2 is distributed in the hope that it will be useful,
23 !! but WITHOUT ANY WARRANTY; without even the implied warranty of
24 !! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 !! GNU General Public License for more details.
26 !!
27 !! You should have received a copy of the GNU General Public License
28 !! along with MAIC-2. If not, see <http://www.gnu.org/licenses/>.
29 !<
30 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 
32 !-------------------------------------------------------------------------------
33 !> Solution of a system of linear equations Ax=b with tridiagonal matrix A.
34 !! @param[in] a0 a0(j) is element A_(j,j-1) of matrix A
35 !! @param[in] a1 a1(j) is element A_(j,j) of matrix A
36 !! @param[in] a2 a2(j) is element A_(j,j+1) of matrix A
37 !! @param[in] b inhomogeneity vector
38 !! @param[in] n_rows size of matrix A (indices run from 0 (!!!) to n_rows)
39 !! @param[out] x solution vector.
40 !<------------------------------------------------------------------------------
41 subroutine tri_sle(a0, a1, a2, x, b, n_rows)
42 
43 use maic2_types
44 
45 implicit none
46 
47 integer(i4b), intent(in) :: n_rows
48 real(dp), dimension(0:*), intent(inout) :: a0, a1, a2, b
49 real(dp), dimension(0:*), intent(out) :: x
50 
51 integer(i4b) :: n
52 real(dp), allocatable, dimension(:) :: help_x
53 
54 !-------- Generation of an upper triangular matrix --------
55 
56 do n=1, n_rows
57  a1(n) = a1(n) - a0(n)/a1(n-1)*a2(n-1)
58 end do
59 
60 do n=1, n_rows
61  b(n) = b(n) - a0(n)/a1(n-1)*b(n-1)
62 ! a0(n) = 0.0_dp , not needed in the following, therefore
63 ! not set
64 end do
65 
66 !-------- Iterative solution of the new system --------
67 
68 ! x(n_rows) = b(n_rows)/a1(n_rows)
69 
70 ! do n=n_rows-1, 0, -1
71 ! x(n) = (b(n)-a2(n)*x(n+1))/a1(n)
72 ! end do
73 
74 allocate(help_x(0:n_rows))
75 
76 help_x(0) = b(n_rows)/a1(n_rows)
77 
78 do n=1, n_rows
79  help_x(n) = b(n_rows-n)/a1(n_rows-n) &
80  -a2(n_rows-n)/a1(n_rows-n)*help_x(n-1)
81 end do
82 
83 do n=0, n_rows
84  x(n) = help_x(n_rows-n)
85 end do
86 
87 ! (The trick with the help_x was introduced in order to avoid
88 ! the negative step in the original, blanked-out loop.)
89 
90 deallocate(help_x)
91 
92 ! WARNING: Subroutine does not check for elements of the main
93 ! diagonal becoming zero. In this case it crashes even
94 ! though the system may be solvable. Otherwise ok.
95 
96 end subroutine tri_sle
97 !