portedportable, not yet portedexecuted, not portableexecutable, not hit by this run
| 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. | |
| 6 | module regrid_solvers | |
| 7 | ||
| 8 | use MOM_error_handler, only : MOM_error, FATAL | |
| 9 | ||
| 10 | implicit none ; private | |
| 11 | ||
| 12 | public :: solve_linear_system, linear_solver, solve_tridiagonal_system, solve_diag_dominant_tridiag | |
| 13 | ||
| 14 | contains | |
| 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. | |
| 21 | 0 | subroutine 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 | 0 | 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 | 0 | 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 | 0 | found_pivot = .false. |
| 50 | 0 | k = i |
| 51 | 0 | do while ( ( .NOT. found_pivot ) .AND. ( k <= N ) ) |
| 52 | 0 | if ( abs( A(k,i) ) > eps ) then ! A valid pivot has been found |
| 53 | 0 | found_pivot = .true. |
| 54 | else ! Seek a valid pivot in the next row | |
| 55 | 0 | 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 | 0 | if ( .NOT. found_pivot ) then |
| 61 | 0 | write(0,*) ' A=',A |
| 62 | 0 | 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 | 0 | if ( k /= i ) then |
| 68 | 0 | do j = 1,N |
| 69 | 0 | swap_a = A(i,j) ; A(i,j) = A(k,j) ; A(k,j) = swap_a |
| 70 | enddo | |
| 71 | 0 | 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 | 0 | if (old_answers) then |
| 76 | 0 | pivot = A(i,i) |
| 77 | 0 | do j = i,N ; A(i,j) = A(i,j) / pivot ; enddo |
| 78 | 0 | R(i) = R(i) / pivot |
| 79 | else | |
| 80 | 0 | I_pivot = 1.0 / A(i,i) |
| 81 | 0 | A(i,i) = 1.0 |
| 82 | 0 | do j = i+1,N ; A(i,j) = A(i,j) * I_pivot ; enddo |
| 83 | 0 | 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 | 0 | do k = i+1,N ! k is the row index |
| 90 | 0 | 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 | 0 | do j = i+1,N ! j is the column index |
| 93 | 0 | A(k,j) = A(k,j) - factor * A(i,j) |
| 94 | enddo | |
| 95 | 0 | 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 | 0 | X(N) = R(N) / A(N,N) ! The last row is now trivially solved. |
| 102 | 0 | do i = N-1,1,-1 ! loop on rows, starting from second to last row |
| 103 | 0 | X(i) = R(i) |
| 104 | 0 | do j = i+1,N |
| 105 | 0 | X(i) = X(i) - A(i,j) * X(j) |
| 106 | enddo | |
| 107 | 0 | if (old_answers) X(i) = X(i) / A(i,i) |
| 108 | enddo | |
| 109 | ||
| 110 | 0 | end 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. | |
| 117 | 0 | subroutine 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 | 0 | 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 | 0 | do k=i,N ; if ( abs( A(i,k) ) > eps ) exit ; enddo ! end loop to find pivot |
| 138 | 0 | if ( k > N ) then ! No pivot could be found and the system is singular. |
| 139 | 0 | write(0,*) ' A=',A |
| 140 | 0 | 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 | 0 | if ( k /= i ) then |
| 146 | 0 | do j=i,N ; swap = A(j,i) ; A(j,i) = A(j,k) ; A(j,k) = swap ; enddo |
| 147 | 0 | 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 | 0 | I_pivot = 1.0 / A(i,i) |
| 152 | 0 | A(i,i) = 1.0 |
| 153 | 0 | do j=i+1,N ; A(j,i) = A(j,i) * I_pivot ; enddo |
| 154 | 0 | 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 | 0 | do k=i+1,N ! k is the row index |
| 158 | 0 | 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 | 0 | do j=i+1,N ; A(j,k) = A(j,k) - factor * A(j,i) ; enddo |
| 161 | 0 | R(k) = R(k) - factor * R(i) |
| 162 | enddo | |
| 163 | ||
| 164 | enddo ! end loop on i | |
| 165 | ||
| 166 | 0 | if (A(N,N) == 0.0) then |
| 167 | ! no pivot could be found, and the sytem is singular | |
| 168 | 0 | 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 | 0 | X(N) = R(N) / A(N,N) ! The last row is now trivially solved. |
| 173 | 0 | do i=N-1,1,-1 ! loop on rows, starting from second to last row |
| 174 | 0 | X(i) = R(i) |
| 175 | 0 | do j=i+1,N ; X(i) = X(i) - A(j,i) * X(j) ; enddo |
| 176 | enddo | |
| 177 | ||
| 178 | 0 | end 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) | |
| 185 | 0 | subroutine 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 | 0 | real, dimension(N) :: pivot ! The pivot value [A] |
| 195 | 0 | real, dimension(N) :: Al_piv ! The lower diagonal divided by the pivot value [nondim] |
| 196 | 0 | 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 | 0 | old_answers = .true. ; if (present(answer_date)) old_answers = (answer_date < 20190101) |
| 202 | ||
| 203 | 0 | 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 | 0 | pivot(1) = Ad(1) |
| 207 | 0 | X(1) = R(1) |
| 208 | 0 | do k = 2,N |
| 209 | 0 | Al_piv(k) = Al(k) / pivot(k-1) |
| 210 | 0 | pivot(k) = Ad(k) - Al_piv(k) * Au(k-1) |
| 211 | 0 | X(k) = R(k) - Al_piv(k) * X(k-1) |
| 212 | enddo | |
| 213 | ||
| 214 | ! Backward sweep | |
| 215 | 0 | X(N) = R(N) / pivot(N) ! This should be X(N) / pivot(N), but is OK if Al(N) = 0. |
| 216 | 0 | do k = N-1,1,-1 |
| 217 | 0 | 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 | 0 | I_pivot = 1.0 / Ad(1) |
| 225 | 0 | X(1) = R(1) * I_pivot |
| 226 | 0 | do k = 2,N |
| 227 | 0 | c1(K-1) = Au(k-1) * I_pivot |
| 228 | 0 | I_pivot = 1.0 / (Ad(k) - Al(k) * c1(K-1)) |
| 229 | 0 | X(k) = (R(k) - Al(k) * X(k-1)) * I_pivot |
| 230 | enddo | |
| 231 | ! Backward sweep | |
| 232 | 0 | do k = N-1,1,-1 |
| 233 | 0 | X(k) = X(k) - c1(K) * X(k+1) |
| 234 | enddo | |
| 235 | ||
| 236 | endif | |
| 237 | ||
| 238 | 0 | end 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. | |
| 248 | 303504 | subroutine 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 | 607008 | 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 | 303504 | I_pivot = 1.0 / (Ac(1) + Au(1)) |
| 265 | 303504 | d1 = Ac(1) * I_pivot |
| 266 | 303504 | c1(1) = Au(1) * I_pivot |
| 267 | 303504 | X(1) = R(1) * I_pivot |
| 268 | 22762800 | do k=2,N-1 |
| 269 | 22459296 | denom_t1 = Ac(k) + d1 * Al(k) |
| 270 | 22459296 | I_pivot = 1.0 / (denom_t1 + Au(k)) |
| 271 | 22459296 | d1 = denom_t1 * I_pivot |
| 272 | 22459296 | c1(k) = Au(k) * I_pivot |
| 273 | 22762800 | X(k) = (R(k) - Al(k) * X(k-1)) * I_pivot |
| 274 | enddo | |
| 275 | 303504 | I_pivot = 1.0 / (Ac(N) + d1 * Al(N)) |
| 276 | 303504 | X(N) = (R(N) - Al(N) * X(N-1)) * I_pivot |
| 277 | ! Backward sweep | |
| 278 | 23066304 | do k=N-1,1,-1 |
| 279 | 23066304 | X(k) = X(k) - c1(k) * X(k+1) |
| 280 | enddo | |
| 281 | ||
| 282 | 303504 | end 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 | ||
| 293 | end module regrid_solvers |