regrid_solvers.F90

1! This file is part of MOM6, the Modular Ocean Model version 6.
2! See the LICENSE file for licensing information.
3! SPDX-License-Identifier: Apache-2.0
4
5!> Solvers of linear systems.
7
8use mom_error_handler, only : mom_error, fatal
9
10implicit none ; private
11
13
14contains
15
16!> Solve the linear system AX = R by Gaussian elimination
17!!
18!! This routine uses Gauss's algorithm to transform the system's original
19!! matrix into an upper triangular matrix. Back substitution yields the answer.
20!! The matrix A must be square, with the first index varing down the column.
21subroutine solve_linear_system( A, R, X, N, answer_date )
22 integer, intent(in) :: n !< The size of the system
23 real, dimension(N,N), intent(inout) :: a !< The matrix being inverted in arbitrary units [A] on
24 !! input, but internally modified to become nondimensional
25 !! during the solver.
26 real, dimension(N), intent(inout) :: r !< system right-hand side in arbitrary units [A B] on
27 !! input, but internally modified to have units of [B]
28 !! during the solver
29 real, dimension(N), intent(inout) :: x !< solution vector in arbitrary units [B]
30 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
31 ! Local variables
32 real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed [A]
33 real :: factor ! The factor that eliminates the leading nonzero element in a row [A-1]
34 real :: pivot, i_pivot ! The pivot value and its reciprocal, in [A] and [A-1]
35 real :: swap_a, swap_b ! Swap space in various units [various]
36 logical :: found_pivot ! If true, a pivot has been found
37 logical :: old_answers ! If true, use expressions that give the original (2008 through 2018) MOM6 answers
38 integer :: i, j, k
39
40 old_answers = .true. ; if (present(answer_date)) old_answers = (answer_date < 20190101)
41
42 ! Loop on rows to transform the problem into multiplication by an upper-right matrix.
43 do i = 1,n-1
44
45
46 ! Start to look for a pivot in the current row, i. If the pivot in row i is not valid,
47 ! keep looking for a valid pivot by searching the entries of column i in rows below row i.
48 ! Once a valid pivot is found (say in row k), rows i and k are swaped.
49 found_pivot = .false.
50 k = i
51 do while ( ( .NOT. found_pivot ) .AND. ( k <= n ) )
52 if ( abs( a(k,i) ) > eps ) then ! A valid pivot has been found
53 found_pivot = .true.
54 else ! Seek a valid pivot in the next row
55 k = k + 1
56 endif
57 enddo ! end loop to find pivot
58
59 ! If no pivot could be found, the system is singular.
60 if ( .NOT. found_pivot ) then
61 write(0,*) ' A=',a
62 call mom_error( fatal, 'The linear system is singular !' )
63 endif
64
65 ! If the pivot is in a row that is different than row i, that is if
66 ! k is different than i, we need to swap those two rows
67 if ( k /= i ) then
68 do j = 1,n
69 swap_a = a(i,j) ; a(i,j) = a(k,j) ; a(k,j) = swap_a
70 enddo
71 swap_b = r(i) ; r(i) = r(k) ; r(k) = swap_b
72 endif
73
74 ! Transform pivot to 1 by dividing the entire row (right-hand side included) by the pivot
75 if (old_answers) then
76 pivot = a(i,i)
77 do j = i,n ; a(i,j) = a(i,j) / pivot ; enddo
78 r(i) = r(i) / pivot
79 else
80 i_pivot = 1.0 / a(i,i)
81 a(i,i) = 1.0
82 do j = i+1,n ; a(i,j) = a(i,j) * i_pivot ; enddo
83 r(i) = r(i) * i_pivot
84 endif
85
86 ! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1
87
88 ! Put zeros in column for all rows below that contain the pivot (which is row i)
89 do k = i+1,n ! k is the row index
90 factor = a(k,i)
91 ! A(k,i) = 0.0 ! These elements are not used again, so this line can be skipped for speed.
92 do j = i+1,n ! j is the column index
93 a(k,j) = a(k,j) - factor * a(i,j)
94 enddo
95 r(k) = r(k) - factor * r(i)
96 enddo
97
98 enddo ! end loop on i
99
100 ! Solve system by back substituting in what is now an upper-right matrix.
101 x(n) = r(n) / a(n,n) ! The last row is now trivially solved.
102 do i = n-1,1,-1 ! loop on rows, starting from second to last row
103 x(i) = r(i)
104 do j = i+1,n
105 x(i) = x(i) - a(i,j) * x(j)
106 enddo
107 if (old_answers) x(i) = x(i) / a(i,i)
108 enddo
109
110end subroutine solve_linear_system
111
112!> Solve the linear system AX = R by Gaussian elimination
113!!
114!! This routine uses Gauss's algorithm to transform the system's original
115!! matrix into an upper triangular matrix. Back substitution then yields the answer.
116!! The matrix A must be square, with the first index varing along the row.
117subroutine linear_solver( N, A, R, X )
118 integer, intent(in) :: n !< The size of the system
119 real, dimension(N,N), intent(inout) :: a !< The matrix being inverted in arbitrary units [A] on
120 !! input, but internally modified to become nondimensional
121 !! during the solver.
122 real, dimension(N), intent(inout) :: r !< system right-hand side in [A B] on input, but internally
123 !! modified to have units of [B] during the solver
124 real, dimension(N), intent(inout) :: x !< solution vector [B]
125
126 ! Local variables
127 real, parameter :: eps = 0.0 ! Minimum pivot magnitude allowed [A]
128 real :: factor ! The factor that eliminates the leading nonzero element in a row [A-1].
129 real :: i_pivot ! The reciprocal of the pivot value [A-1]
130 real :: swap ! Swap space used in various units [various]
131 integer :: i, j, k
132
133 ! Loop on rows to transform the problem into multiplication by an upper-right matrix.
134 do i=1,n-1
135 ! Seek a pivot for column i starting in row i, and continuing into the remaining rows. If the
136 ! pivot is in a row other than i, swap them. If no valid pivot is found, i = N+1 after this loop.
137 do k=i,n ; if ( abs( a(i,k) ) > eps ) exit ; enddo ! end loop to find pivot
138 if ( k > n ) then ! No pivot could be found and the system is singular.
139 write(0,*) ' A=',a
140 call mom_error( fatal, 'The linear system is singular !' )
141 endif
142
143 ! If the pivot is in a row that is different than row i, swap those two rows, noting that both
144 ! rows start with i-1 zero values.
145 if ( k /= i ) then
146 do j=i,n ; swap = a(j,i) ; a(j,i) = a(j,k) ; a(j,k) = swap ; enddo
147 swap = r(i) ; r(i) = r(k) ; r(k) = swap
148 endif
149
150 ! Transform the pivot to 1 by dividing the entire row (right-hand side included) by the pivot
151 i_pivot = 1.0 / a(i,i)
152 a(i,i) = 1.0
153 do j=i+1,n ; a(j,i) = a(j,i) * i_pivot ; enddo
154 r(i) = r(i) * i_pivot
155
156 ! Put zeros in column for all rows below that contain the pivot (which is row i)
157 do k=i+1,n ! k is the row index
158 factor = a(i,k)
159 ! A(i,k) = 0.0 ! These elements are not used again, so this line can be skipped for speed.
160 do j=i+1,n ; a(j,k) = a(j,k) - factor * a(j,i) ; enddo
161 r(k) = r(k) - factor * r(i)
162 enddo
163
164 enddo ! end loop on i
165
166 if (a(n,n) == 0.0) then
167 ! no pivot could be found, and the sytem is singular
168 call mom_error(fatal, 'The final pivot in linear_solver is zero.')
169 endif
170
171 ! Solve the system by back substituting into what is now an upper-right matrix.
172 x(n) = r(n) / a(n,n) ! The last row is now trivially solved.
173 do i=n-1,1,-1 ! loop on rows, starting from second to last row
174 x(i) = r(i)
175 do j=i+1,n ; x(i) = x(i) - a(j,i) * x(j) ; enddo
176 enddo
177
178end subroutine linear_solver
179
180
181!> Solve the tridiagonal system AX = R
182!!
183!! This routine uses Thomas's algorithm to solve the tridiagonal system AX = R.
184!! (A is made up of lower, middle and upper diagonals)
185subroutine solve_tridiagonal_system( Al, Ad, Au, R, X, N, answer_date )
186 integer, intent(in) :: n !< The size of the system
187 real, dimension(N), intent(in) :: ad !< Matrix center diagonal in arbitrary units [A]
188 real, dimension(N), intent(in) :: al !< Matrix lower diagonal [A]
189 real, dimension(N), intent(in) :: au !< Matrix upper diagonal [A]
190 real, dimension(N), intent(in) :: r !< system right-hand side in arbitrary units [A B]
191 real, dimension(N), intent(out) :: x !< solution vector in arbitrary units [B]
192 integer, optional, intent(in) :: answer_date !< The vintage of the expressions to use
193 ! Local variables
194 real, dimension(N) :: pivot ! The pivot value [A]
195 real, dimension(N) :: al_piv ! The lower diagonal divided by the pivot value [nondim]
196 real, dimension(N) :: c1 ! Au / pivot for the backward sweep [nondim]
197 real :: i_pivot ! The inverse of the most recent pivot [A-1]
198 integer :: k ! Loop index
199 logical :: old_answers ! If true, use expressions that give the original (2008 through 2018) MOM6 answers
200
201 old_answers = .true. ; if (present(answer_date)) old_answers = (answer_date < 20190101)
202
203 if (old_answers) then
204 ! This version gives the same answers as the original (2008 through 2018) MOM6 code
205 ! Factorization and forward sweep
206 pivot(1) = ad(1)
207 x(1) = r(1)
208 do k = 2,n
209 al_piv(k) = al(k) / pivot(k-1)
210 pivot(k) = ad(k) - al_piv(k) * au(k-1)
211 x(k) = r(k) - al_piv(k) * x(k-1)
212 enddo
213
214 ! Backward sweep
215 x(n) = r(n) / pivot(n) ! This should be X(N) / pivot(N), but is OK if Al(N) = 0.
216 do k = n-1,1,-1
217 x(k) = ( x(k) - au(k)*x(k+1) ) / pivot(k)
218 enddo
219 else
220 ! This is a more typical implementation of a tridiagonal solver than the one above.
221 ! It is mathematically equivalent but differs at roundoff, which can cascade up to larger values.
222
223 ! Factorization and forward sweep
224 i_pivot = 1.0 / ad(1)
225 x(1) = r(1) * i_pivot
226 do k = 2,n
227 c1(k-1) = au(k-1) * i_pivot
228 i_pivot = 1.0 / (ad(k) - al(k) * c1(k-1))
229 x(k) = (r(k) - al(k) * x(k-1)) * i_pivot
230 enddo
231 ! Backward sweep
232 do k = n-1,1,-1
233 x(k) = x(k) - c1(k) * x(k+1)
234 enddo
235
236 endif
237
238end subroutine solve_tridiagonal_system
239
240
241!> Solve the tridiagonal system AX = R
242!!
243!! This routine uses a variant of Thomas's algorithm to solve the tridiagonal system AX = R, in
244!! a form that is guaranteed to avoid dividing by a zero pivot. The matrix A is made up of
245!! lower (Al) and upper diagonals (Au) and a central diagonal Ad = Ac+Al+Au, where
246!! Al, Au, and Ac are all positive (or negative) definite. However when Ac is smaller than
247!! roundoff compared with (Al+Au), the answers are prone to inaccuracy.
248subroutine solve_diag_dominant_tridiag( Al, Ac, Au, R, X, N )
249 integer, intent(in) :: n !< The size of the system
250 real, dimension(N), intent(in) :: ac !< Matrix center diagonal offset from Al + Au in arbitrary units [A]
251 real, dimension(N), intent(in) :: al !< Matrix lower diagonal [A]
252 real, dimension(N), intent(in) :: au !< Matrix upper diagonal [A]
253 real, dimension(N), intent(in) :: r !< system right-hand side in arbitrary units [A B]
254 real, dimension(N), intent(out) :: x !< solution vector in arbitrary units [B]
255 ! Local variables
256 real, dimension(N) :: c1 ! Au / pivot for the backward sweep [nondim]
257 real :: d1 ! The next value of 1.0 - c1 [nondim]
258 real :: i_pivot ! The inverse of the most recent pivot [A-1]
259 real :: denom_t1 ! The first term in the denominator of the inverse of the pivot [A]
260 integer :: k ! Loop index
261
262 ! Factorization and forward sweep, in a form that will never give a division by a
263 ! zero pivot for positive definite Ac, Al, and Au.
264 i_pivot = 1.0 / (ac(1) + au(1))
265 d1 = ac(1) * i_pivot
266 c1(1) = au(1) * i_pivot
267 x(1) = r(1) * i_pivot
268 do k=2,n-1
269 denom_t1 = ac(k) + d1 * al(k)
270 i_pivot = 1.0 / (denom_t1 + au(k))
271 d1 = denom_t1 * i_pivot
272 c1(k) = au(k) * i_pivot
273 x(k) = (r(k) - al(k) * x(k-1)) * i_pivot
274 enddo
275 i_pivot = 1.0 / (ac(n) + d1 * al(n))
276 x(n) = (r(n) - al(n) * x(n-1)) * i_pivot
277 ! Backward sweep
278 do k=n-1,1,-1
279 x(k) = x(k) - c1(k) * x(k+1)
280 enddo
281
282end subroutine solve_diag_dominant_tridiag
283
284
285!> \namespace regrid_solvers
286!!
287!! Date of creation: 2008.06.12
288!! L. White
289!!
290!! This module contains solvers of linear systems.
291!! These routines have now been updated for greater efficiency, especially in special cases.
292
293end module regrid_solvers