MOM6
regrid_solvers.F90
Go to the documentation of this file.
1 !> Solvers of linear systems.
3 
4 ! This file is part of MOM6. See LICENSE.md for the license.
5 
6 use mom_error_handler, only : mom_error, fatal
7 
8 implicit none ; private
9 
11 
12 contains
13 
14 !> Solve the linear system AX = B by Gaussian elimination
15 !!
16 !! This routine uses Gauss's algorithm to transform the system's original
17 !! matrix into an upper triangular matrix. Back substitution yields the answer.
18 !! The matrix A must be square and its size must be that of the vectors B and X.
19 subroutine solve_linear_system( A, B, X, system_size )
20  real, dimension(:,:), intent(inout) :: a !< The matrix being inverted
21  real, dimension(:), intent(inout) :: b !< system right-hand side
22  real, dimension(:), intent(inout) :: x !< solution vector
23  integer, intent(in) :: system_size !< The size of the system
24  ! Local variables
25  integer :: i, j, k
26  real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed
27  real :: factor
28  real :: pivot
29  real :: swap_a, swap_b
30  logical :: found_pivot ! boolean indicating whether
31  ! a pivot has been found
32  ! Loop on rows
33  do i = 1,system_size-1
34 
35  found_pivot = .false.
36 
37  ! Start to look for a pivot in row i. If the pivot
38  ! in row i -- which is the current row -- is not valid,
39  ! we keep looking for a valid pivot by searching the
40  ! entries of column i in rows below row i. Once a valid
41  ! pivot is found (say in row k), rows i and k are swaped.
42  k = i
43  do while ( ( .NOT. found_pivot ) .AND. ( k <= system_size ) )
44 
45  if ( abs( a(k,i) ) > eps ) then ! a valid pivot is found
46  found_pivot = .true.
47  else ! Go to the next row to see
48  ! if there is a valid pivot there
49  k = k + 1
50  endif
51 
52  enddo ! end loop to find pivot
53 
54  ! If no pivot could be found, the system is singular and we need
55  ! to end the execution
56  if ( .NOT. found_pivot ) then
57  write(0,*) ' A=',a
58  call mom_error( fatal, 'The linear system is singular !' )
59  endif
60 
61  ! If the pivot is in a row that is different than row i, that is if
62  ! k is different than i, we need to swap those two rows
63  if ( k /= i ) then
64  do j = 1,system_size
65  swap_a = a(i,j)
66  a(i,j) = a(k,j)
67  a(k,j) = swap_a
68  enddo
69  swap_b = b(i)
70  b(i) = b(k)
71  b(k) = swap_b
72  endif
73 
74  ! Transform pivot to 1 by dividing the entire row
75  ! (right-hand side included) by the pivot
76  pivot = a(i,i)
77  do j = i,system_size
78  a(i,j) = a(i,j) / pivot
79  enddo
80  b(i) = b(i) / pivot
81 
82  ! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1
83 
84  ! Put zeros in column for all rows below that containing
85  ! pivot (which is row i)
86  do k = (i+1),system_size ! k is the row index
87  factor = a(k,i)
88  do j = (i+1),system_size ! j is the column index
89  a(k,j) = a(k,j) - factor * a(i,j)
90  enddo
91  b(k) = b(k) - factor * b(i)
92  enddo
93 
94  enddo ! end loop on i
95 
96 
97  ! Solve system by back substituting
98  x(system_size) = b(system_size) / a(system_size,system_size)
99  do i = system_size-1,1,-1 ! loop on rows, starting from second to last row
100  x(i) = b(i)
101  do j = (i+1),system_size
102  x(i) = x(i) - a(i,j) * x(j)
103  enddo
104  x(i) = x(i) / a(i,i)
105  enddo
106 
107 end subroutine solve_linear_system
108 
109 !> Solve the tridiagonal system AX = B
110 !!
111 !! This routine uses Thomas's algorithm to solve the tridiagonal system AX = B.
112 !! (A is made up of lower, middle and upper diagonals)
113 subroutine solve_tridiagonal_system( Al, Ad, Au, B, X, system_size )
114  real, dimension(:), intent(inout) :: ad !< Maxtix center diagonal
115  real, dimension(:), intent(inout) :: al !< Matrix lower diagonal
116  real, dimension(:), intent(inout) :: au !< Matrix upper diagonal
117  real, dimension(:), intent(inout) :: b !< system right-hand side
118  real, dimension(:), intent(inout) :: x !< solution vector
119  integer, intent(in) :: system_size !< The size of the system
120  ! Local variables
121  integer :: k ! Loop index
122  integer :: n ! system size
123 
124  n = system_size
125 
126  ! Factorization
127  do k = 1,n-1
128  al(k+1) = al(k+1) / ad(k)
129  ad(k+1) = ad(k+1) - al(k+1) * au(k)
130  enddo
131 
132  ! Forward sweep
133  do k = 2,n
134  b(k) = b(k) - al(k) * b(k-1)
135  enddo
136 
137  ! Backward sweep
138  x(n) = b(n) / ad(n)
139  do k = n-1,1,-1
140  x(k) = ( b(k) - au(k)*x(k+1) ) / ad(k)
141  enddo
142 
143 end subroutine solve_tridiagonal_system
144 
145 !> \namespace regrid_solvers
146 !!
147 !! Date of creation: 2008.06.12
148 !! L. White
149 !!
150 !! This module contains solvers of linear systems.
151 !! These routines could (should ?) be replaced later by more efficient ones.
152 
153 end module regrid_solvers
regrid_solvers
Solvers of linear systems.
Definition: regrid_solvers.F90:2
regrid_solvers::solve_linear_system
subroutine, public solve_linear_system(A, B, X, system_size)
Solve the linear system AX = B by Gaussian elimination.
Definition: regrid_solvers.F90:20
regrid_solvers::solve_tridiagonal_system
subroutine, public solve_tridiagonal_system(Al, Ad, Au, B, X, system_size)
Solve the tridiagonal system AX = B.
Definition: regrid_solvers.F90:114
mom_error_handler::mom_error
subroutine, public mom_error(level, message, all_print)
This provides a convenient interface for writing an mpp_error message with run-time filter based on a...
Definition: MOM_error_handler.F90:72
mom_error_handler
Routines for error handling and I/O management.
Definition: MOM_error_handler.F90:2