\hypertarget{namespaceregrid__solvers}{}\section{regrid\+\_\+solvers Module Reference}
\label{namespaceregrid__solvers}\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}}


\subsection{Detailed Description}
Solvers of linear systems. 

Date of creation\+: 2008.\+06.\+12 L. White.

This module contains solvers of linear systems. These routines have now been updated for greater efficiency, especially in special cases. \subsection*{Functions/\+Subroutines}
\begin{DoxyCompactItemize}
\item 
subroutine, public \hyperlink{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}{solve\+\_\+linear\+\_\+system} (A, R, X, N, answers\+\_\+2018)
\begin{DoxyCompactList}\small\item\em Solve the linear system AX = R by Gaussian elimination. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}{linear\+\_\+solver} (N, A, R, X)
\begin{DoxyCompactList}\small\item\em Solve the linear system AX = R by Gaussian elimination. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}{solve\+\_\+tridiagonal\+\_\+system} (Al, Ad, Au, R, X, N, answers\+\_\+2018)
\begin{DoxyCompactList}\small\item\em Solve the tridiagonal system AX = R. \end{DoxyCompactList}\item 
subroutine, public \hyperlink{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag} (Al, Ac, Au, R, X, N)
\begin{DoxyCompactList}\small\item\em Solve the tridiagonal system AX = R. \end{DoxyCompactList}\end{DoxyCompactItemize}


\subsection{Function/\+Subroutine Documentation}
\mbox{\Hypertarget{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}\label{namespaceregrid__solvers_a2ed09d1e0857610eb3638601e20d2a2c}} 
\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!linear\+\_\+solver@{linear\+\_\+solver}}
\index{linear\+\_\+solver@{linear\+\_\+solver}!regrid\+\_\+solvers@{regrid\+\_\+solvers}}
\subsubsection{\texorpdfstring{linear\+\_\+solver()}{linear\_solver()}}
{\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::linear\+\_\+solver (\begin{DoxyParamCaption}\item[{integer, intent(in)}]{N,  }\item[{real, dimension(n,n), intent(inout)}]{A,  }\item[{real, dimension(n), intent(inout)}]{R,  }\item[{real, dimension(n), intent(inout)}]{X }\end{DoxyParamCaption})}



Solve the linear system AX = R by Gaussian elimination. 

This routine uses Gauss\textquotesingle{}s algorithm to transform the system\textquotesingle{}s original matrix into an upper triangular matrix. Back substitution then yields the answer. The matrix A must be square, with the first index varing along the row.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in,out}  & {\em a} & The matrix being inverted \mbox{[}nondim\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em r} & system right-\/hand side \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em x} & solution vector \mbox{[}A\mbox{]} \\
\hline
\end{DoxyParams}


Definition at line 112 of file regrid\+\_\+solvers.\+F90.


\begin{DoxyCode}
112   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{  !< The size of the system}
113   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: a\textcolor{comment}{  !< The matrix being inverted [nondim]}
114   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: r\textcolor{comment}{  !< system right-hand side [A]}
115   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: x\textcolor{comment}{  !< solution vector [A]}
116 
117   \textcolor{comment}{! Local variables}
118   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter} :: eps = 0.0   \textcolor{comment}{! Minimum pivot magnitude allowed}
119   \textcolor{keywordtype}{real}    :: factor       \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.}
120   \textcolor{keywordtype}{real}    :: i\_pivot      \textcolor{comment}{! The reciprocal of the pivot value [inverse of the input units of a row of A]}
121   \textcolor{keywordtype}{real}    :: swap
122   \textcolor{keywordtype}{logical} :: found\_pivot  \textcolor{comment}{! If true, a pivot has been found}
123   \textcolor{keywordtype}{integer} :: i, j, k
124 
125   \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-right matrix.}
126   \textcolor{keywordflow}{do} i=1,n-1
127     \textcolor{comment}{! Seek a pivot for column i starting in row i, and continuing into the remaining rows.  If the}
128     \textcolor{comment}{! pivot is in a row other than i, swap them.  If no valid pivot is found, i = N+1 after this loop.}
129     \textcolor{keywordflow}{do} k=i,n ; \textcolor{keywordflow}{if} ( abs( a(i,k) ) > eps ) \textcolor{keywordflow}{exit} ;\textcolor{keywordflow}{ enddo} \textcolor{comment}{! end loop to find pivot}
130     \textcolor{keywordflow}{if} ( k > n ) \textcolor{keywordflow}{then}  \textcolor{comment}{! No pivot could be found and the system is singular.}
131       \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a
132       \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system is singular !'} )
133 \textcolor{keywordflow}{    endif}
134 
135     \textcolor{comment}{! If the pivot is in a row that is different than row i, swap those two rows, noting that both}
136     \textcolor{comment}{! rows start with i-1 zero values.}
137     \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then}
138       \textcolor{keywordflow}{do} j=i,n ; swap = a(j,i) ; a(j,i) = a(j,k) ; a(j,k) = swap ;\textcolor{keywordflow}{ enddo}
139       swap = r(i) ; r(i) = r(k) ; r(k) = swap
140 \textcolor{keywordflow}{    endif}
141 
142     \textcolor{comment}{! Transform the pivot to 1 by dividing the entire row (right-hand side included) by the pivot}
143     i\_pivot = 1.0 / a(i,i)
144     a(i,i) = 1.0
145     \textcolor{keywordflow}{do} j=i+1,n ; a(j,i) = a(j,i) * i\_pivot ;\textcolor{keywordflow}{ enddo}
146     r(i) = r(i) * i\_pivot
147 
148     \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)}
149     \textcolor{keywordflow}{do} k=i+1,n    \textcolor{comment}{! k is the row index}
150       factor = a(i,k)
151       \textcolor{comment}{! A(i,k) = 0.0  ! These elements are not used again, so this line can be skipped for speed.}
152       \textcolor{keywordflow}{do} j=i+1,n ; a(j,k) = a(j,k) - factor * a(j,i) ;\textcolor{keywordflow}{ enddo}
153       r(k) = r(k) - factor * r(i)
154 \textcolor{keywordflow}{    enddo}
155 
156 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on i}
157 
158   \textcolor{comment}{! Solve the system by back substituting into what is now an upper-right matrix.}
159   x(n) = r(n) / a(n,n)  \textcolor{comment}{! The last row is now trivially solved.}
160   \textcolor{keywordflow}{do} i=n-1,1,-1 \textcolor{comment}{! loop on rows, starting from second to last row}
161     x(i) = r(i)
162     \textcolor{keywordflow}{do} j=i+1,n ; x(i) = x(i) - a(j,i) * x(j) ;\textcolor{keywordflow}{ enddo}
163 \textcolor{keywordflow}{  enddo}
164 
\end{DoxyCode}
\mbox{\Hypertarget{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}\label{namespaceregrid__solvers_a3f1d27aaab0a19c9abd07c96bb704bb6}} 
\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+diag\+\_\+dominant\+\_\+tridiag@{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag}}
\index{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag@{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag}!regrid\+\_\+solvers@{regrid\+\_\+solvers}}
\subsubsection{\texorpdfstring{solve\+\_\+diag\+\_\+dominant\+\_\+tridiag()}{solve\_diag\_dominant\_tridiag()}}
{\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::solve\+\_\+diag\+\_\+dominant\+\_\+tridiag (\begin{DoxyParamCaption}\item[{real, dimension(n), intent(in)}]{Al,  }\item[{real, dimension(n), intent(in)}]{Ac,  }\item[{real, dimension(n), intent(in)}]{Au,  }\item[{real, dimension(n), intent(in)}]{R,  }\item[{real, dimension(n), intent(out)}]{X,  }\item[{integer, intent(in)}]{N }\end{DoxyParamCaption})}



Solve the tridiagonal system AX = R. 

This routine uses a variant of Thomas\textquotesingle{}s algorithm to solve the tridiagonal system AX = R, in a form that is guaranteed to avoid dividing by a zero pivot. The matrix A is made up of lower (Al) and upper diagonals (Au) and a central diagonal Ad = Ac+\+Al+\+Au, where Al, Au, and Ac are all positive (or negative) definite. However when Ac is smaller than roundoff compared with (Al+\+Au), the answers are prone to inaccuracy.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in}  & {\em ac} & Matrix center diagonal offset from Al + Au\\
\hline
\mbox{\tt in}  & {\em al} & Matrix lower diagonal\\
\hline
\mbox{\tt in}  & {\em au} & Matrix upper diagonal\\
\hline
\mbox{\tt in}  & {\em r} & system right-\/hand side\\
\hline
\mbox{\tt out}  & {\em x} & solution vector \\
\hline
\end{DoxyParams}


Definition at line 235 of file regrid\+\_\+solvers.\+F90.


\begin{DoxyCode}
235   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)}  :: n\textcolor{comment}{   !< The size of the system}
236   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: ac\textcolor{comment}{  !< Matrix center diagonal offset from Al + Au}
237   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: al\textcolor{comment}{  !< Matrix lower diagonal}
238   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: au\textcolor{comment}{  !< Matrix upper diagonal}
239   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: r\textcolor{comment}{   !< system right-hand side}
240   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{   !< solution vector}
241   \textcolor{comment}{! Local variables}
242   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: c1       \textcolor{comment}{! Au / pivot for the backward sweep}
243   \textcolor{keywordtype}{real}               :: d1       \textcolor{comment}{! The next value of 1.0 - c1}
244   \textcolor{keywordtype}{real}               :: i\_pivot  \textcolor{comment}{! The inverse of the most recent pivot}
245   \textcolor{keywordtype}{real}               :: denom\_t1 \textcolor{comment}{! The first term in the denominator of the inverse of the pivot.}
246   \textcolor{keywordtype}{integer}            :: k        \textcolor{comment}{! Loop index}
247 
248   \textcolor{comment}{! Factorization and forward sweep, in a form that will never give a division by a}
249   \textcolor{comment}{! zero pivot for positive definite Ac, Al, and Au.}
250   i\_pivot = 1.0 / (ac(1) + au(1))
251   d1 = ac(1) * i\_pivot
252   c1(1) = au(1) * i\_pivot
253   x(1) = r(1) * i\_pivot
254   \textcolor{keywordflow}{do} k=2,n-1
255     denom\_t1 = ac(k) + d1 * al(k)
256     i\_pivot = 1.0 / (denom\_t1 + au(k))
257     d1 = denom\_t1 * i\_pivot
258     c1(k) = au(k) * i\_pivot
259     x(k) = (r(k) - al(k) * x(k-1)) * i\_pivot
260 \textcolor{keywordflow}{  enddo}
261   i\_pivot = 1.0 / (ac(n) + d1 * al(n))
262   x(n) = (r(n) - al(n) * x(n-1)) * i\_pivot
263   \textcolor{comment}{! Backward sweep}
264   \textcolor{keywordflow}{do} k=n-1,1,-1
265     x(k) = x(k) - c1(k) * x(k+1)
266 \textcolor{keywordflow}{  enddo}
267 
\end{DoxyCode}
\mbox{\Hypertarget{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}\label{namespaceregrid__solvers_a36827238948de49ff0dd0be54cfaf719}} 
\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+linear\+\_\+system@{solve\+\_\+linear\+\_\+system}}
\index{solve\+\_\+linear\+\_\+system@{solve\+\_\+linear\+\_\+system}!regrid\+\_\+solvers@{regrid\+\_\+solvers}}
\subsubsection{\texorpdfstring{solve\+\_\+linear\+\_\+system()}{solve\_linear\_system()}}
{\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::solve\+\_\+linear\+\_\+system (\begin{DoxyParamCaption}\item[{real, dimension(n,n), intent(inout)}]{A,  }\item[{real, dimension(n), intent(inout)}]{R,  }\item[{real, dimension(n), intent(inout)}]{X,  }\item[{integer, intent(in)}]{N,  }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})}



Solve the linear system AX = R by Gaussian elimination. 

This routine uses Gauss\textquotesingle{}s algorithm to transform the system\textquotesingle{}s original matrix into an upper triangular matrix. Back substitution yields the answer. The matrix A must be square, with the first index varing down the column.


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in,out}  & {\em a} & The matrix being inverted \mbox{[}nondim\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em r} & system right-\/hand side \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in,out}  & {\em x} & solution vector \mbox{[}A\mbox{]}\\
\hline
\mbox{\tt in}  & {\em answers\+\_\+2018} & If true or absent use older, less efficient expressions. \\
\hline
\end{DoxyParams}


Definition at line 20 of file regrid\+\_\+solvers.\+F90.


\begin{DoxyCode}
20   \textcolor{keywordtype}{integer},              \textcolor{keywordtype}{intent(in)}    :: n\textcolor{comment}{  !< The size of the system}
21   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N,N)}, \textcolor{keywordtype}{intent(inout)} :: a\textcolor{comment}{  !< The matrix being inverted [nondim]}
22   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: r\textcolor{comment}{  !< system right-hand side [A]}
23   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)},   \textcolor{keywordtype}{intent(inout)} :: x\textcolor{comment}{  !< solution vector [A]}
24   \textcolor{keywordtype}{logical},    \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}    :: answers\_2018\textcolor{comment}{ !< If true or absent use older, less efficient
       expressions.}
25   \textcolor{comment}{! Local variables}
26   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{parameter}       :: eps = 0.0        \textcolor{comment}{! Minimum pivot magnitude allowed}
27   \textcolor{keywordtype}{real}    :: factor       \textcolor{comment}{! The factor that eliminates the leading nonzero element in a row.}
28   \textcolor{keywordtype}{real}    :: pivot, i\_pivot \textcolor{comment}{! The pivot value and its reciprocal [nondim]}
29   \textcolor{keywordtype}{real}    :: swap\_a, swap\_b
30   \textcolor{keywordtype}{logical} :: found\_pivot  \textcolor{comment}{! If true, a pivot has been found}
31   \textcolor{keywordtype}{logical} :: old\_answers  \textcolor{comment}{! If true, use expressions that give the original (2008 through 2018) MOM6
       answers}
32   \textcolor{keywordtype}{integer} :: i, j, k
33 
34   old\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) old\_answers = answers\_2018
35 
36   \textcolor{comment}{! Loop on rows to transform the problem into multiplication by an upper-right matrix.}
37   \textcolor{keywordflow}{do} i = 1,n-1
38 
39 
40     \textcolor{comment}{! Start to look for a pivot in the current row, i.  If the pivot in row i is not valid,}
41     \textcolor{comment}{! keep looking for a valid pivot by searching the entries of column i in rows below row i.}
42     \textcolor{comment}{! Once a valid pivot is found (say in row k), rows i and k are swaped.}
43     found\_pivot = .false.
44     k = i
45     \textcolor{keywordflow}{do} \textcolor{keywordflow}{while} ( ( .NOT. found\_pivot ) .AND. ( k <= n ) )
46       \textcolor{keywordflow}{if} ( abs( a(k,i) ) > eps ) \textcolor{keywordflow}{then}  \textcolor{comment}{! A valid pivot has been found}
47         found\_pivot = .true.
48       \textcolor{keywordflow}{else}                             \textcolor{comment}{! Seek a valid pivot in the next row}
49         k = k + 1
50 \textcolor{keywordflow}{      endif}
51 \textcolor{keywordflow}{    enddo} \textcolor{comment}{! end loop to find pivot}
52 
53     \textcolor{comment}{! If no pivot could be found, the system is singular.}
54     \textcolor{keywordflow}{if} ( .NOT. found\_pivot ) \textcolor{keywordflow}{then}
55       \textcolor{keyword}{write}(0,*) \textcolor{stringliteral}{' A='},a
56       \textcolor{keyword}{call }mom\_error( fatal, \textcolor{stringliteral}{'The linear system is singular !'} )
57 \textcolor{keywordflow}{    endif}
58 
59     \textcolor{comment}{! If the pivot is in a row that is different than row i, that is if}
60     \textcolor{comment}{! k is different than i, we need to swap those two rows}
61     \textcolor{keywordflow}{if} ( k /= i ) \textcolor{keywordflow}{then}
62       \textcolor{keywordflow}{do} j = 1,n
63         swap\_a = a(i,j) ; a(i,j) = a(k,j) ; a(k,j) = swap\_a
64 \textcolor{keywordflow}{      enddo}
65       swap\_b = r(i) ; r(i) = r(k) ; r(k) = swap\_b
66 \textcolor{keywordflow}{    endif}
67 
68     \textcolor{comment}{! Transform pivot to 1 by dividing the entire row (right-hand side included) by the pivot}
69     \textcolor{keywordflow}{if} (old\_answers) \textcolor{keywordflow}{then}
70       pivot = a(i,i)
71       \textcolor{keywordflow}{do} j = i,n ; a(i,j) = a(i,j) / pivot ;\textcolor{keywordflow}{ enddo}
72       r(i) = r(i) / pivot
73     \textcolor{keywordflow}{else}
74       i\_pivot = 1.0 / a(i,i)
75       a(i,i) = 1.0
76       \textcolor{keywordflow}{do} j = i+1,n ; a(i,j) = a(i,j) * i\_pivot ;\textcolor{keywordflow}{ enddo}
77       r(i) = r(i) * i\_pivot
78 \textcolor{keywordflow}{    endif}
79 
80     \textcolor{comment}{! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1}
81 
82     \textcolor{comment}{! Put zeros in column for all rows below that contain the pivot (which is row i)}
83     \textcolor{keywordflow}{do} k = i+1,n    \textcolor{comment}{! k is the row index}
84       factor = a(k,i)
85       \textcolor{comment}{! A(k,i) = 0.0  ! These elements are not used again, so this line can be skipped for speed.}
86       \textcolor{keywordflow}{do} j = i+1,n  \textcolor{comment}{! j is the column index}
87         a(k,j) = a(k,j) - factor * a(i,j)
88 \textcolor{keywordflow}{      enddo}
89       r(k) = r(k) - factor * r(i)
90 \textcolor{keywordflow}{    enddo}
91 
92 \textcolor{keywordflow}{  enddo} \textcolor{comment}{! end loop on i}
93 
94   \textcolor{comment}{! Solve system by back substituting in what is now an upper-right matrix.}
95   x(n) = r(n) / a(n,n)  \textcolor{comment}{! The last row is now trivially solved.}
96   \textcolor{keywordflow}{do} i = n-1,1,-1 \textcolor{comment}{! loop on rows, starting from second to last row}
97     x(i) = r(i)
98     \textcolor{keywordflow}{do} j = i+1,n
99       x(i) = x(i) - a(i,j) * x(j)
100 \textcolor{keywordflow}{    enddo}
101     \textcolor{keywordflow}{if} (old\_answers) x(i) = x(i) / a(i,i)
102 \textcolor{keywordflow}{  enddo}
103 
\end{DoxyCode}
\mbox{\Hypertarget{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}\label{namespaceregrid__solvers_aac4382af38975d9cfcfd6b00adafaeab}} 
\index{regrid\+\_\+solvers@{regrid\+\_\+solvers}!solve\+\_\+tridiagonal\+\_\+system@{solve\+\_\+tridiagonal\+\_\+system}}
\index{solve\+\_\+tridiagonal\+\_\+system@{solve\+\_\+tridiagonal\+\_\+system}!regrid\+\_\+solvers@{regrid\+\_\+solvers}}
\subsubsection{\texorpdfstring{solve\+\_\+tridiagonal\+\_\+system()}{solve\_tridiagonal\_system()}}
{\footnotesize\ttfamily subroutine, public regrid\+\_\+solvers\+::solve\+\_\+tridiagonal\+\_\+system (\begin{DoxyParamCaption}\item[{real, dimension(n), intent(in)}]{Al,  }\item[{real, dimension(n), intent(in)}]{Ad,  }\item[{real, dimension(n), intent(in)}]{Au,  }\item[{real, dimension(n), intent(in)}]{R,  }\item[{real, dimension(n), intent(out)}]{X,  }\item[{integer, intent(in)}]{N,  }\item[{logical, intent(in), optional}]{answers\+\_\+2018 }\end{DoxyParamCaption})}



Solve the tridiagonal system AX = R. 

This routine uses Thomas\textquotesingle{}s algorithm to solve the tridiagonal system AX = R. (A is made up of lower, middle and upper diagonals)


\begin{DoxyParams}[1]{Parameters}
\mbox{\tt in}  & {\em n} & The size of the system\\
\hline
\mbox{\tt in}  & {\em ad} & Matrix center diagonal\\
\hline
\mbox{\tt in}  & {\em al} & Matrix lower diagonal\\
\hline
\mbox{\tt in}  & {\em au} & Matrix upper diagonal\\
\hline
\mbox{\tt in}  & {\em r} & system right-\/hand side\\
\hline
\mbox{\tt out}  & {\em x} & solution vector\\
\hline
\mbox{\tt in}  & {\em answers\+\_\+2018} & If true use older, less acccurate expressions. \\
\hline
\end{DoxyParams}


Definition at line 173 of file regrid\+\_\+solvers.\+F90.


\begin{DoxyCode}
173   \textcolor{keywordtype}{integer},            \textcolor{keywordtype}{intent(in)}  :: n\textcolor{comment}{   !< The size of the system}
174   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: ad\textcolor{comment}{  !< Matrix center diagonal}
175   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: al\textcolor{comment}{  !< Matrix lower diagonal}
176   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: au\textcolor{comment}{  !< Matrix upper diagonal}
177   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(in)}  :: r\textcolor{comment}{   !< system right-hand side}
178   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)}, \textcolor{keywordtype}{intent(out)} :: x\textcolor{comment}{   !< solution vector}
179   \textcolor{keywordtype}{logical},  \textcolor{keywordtype}{optional}, \textcolor{keywordtype}{intent(in)}  :: answers\_2018\textcolor{comment}{ !< If true use older, less acccurate expressions.}
180   \textcolor{comment}{! Local variables}
181   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: pivot, al\_piv
182   \textcolor{keywordtype}{real}, \textcolor{keywordtype}{dimension(N)} :: c1       \textcolor{comment}{! Au / pivot for the backward sweep}
183   \textcolor{keywordtype}{real}    :: i\_pivot  \textcolor{comment}{! The inverse of the most recent pivot}
184   \textcolor{keywordtype}{integer} :: k        \textcolor{comment}{! Loop index}
185   \textcolor{keywordtype}{logical} :: old\_answers  \textcolor{comment}{! If true, use expressions that give the original (2008 through 2018) MOM6
       answers}
186 
187   old\_answers = .true. ; \textcolor{keywordflow}{if} (\textcolor{keyword}{present}(answers\_2018)) old\_answers = answers\_2018
188 
189   \textcolor{keywordflow}{if} (old\_answers) \textcolor{keywordflow}{then}
190     \textcolor{comment}{! This version gives the same answers as the original (2008 through 2018) MOM6 code}
191     \textcolor{comment}{! Factorization and forward sweep}
192     pivot(1) = ad(1)
193     x(1) = r(1)
194     \textcolor{keywordflow}{do} k = 2,n
195       al\_piv(k) = al(k) / pivot(k-1)
196       pivot(k) = ad(k) - al\_piv(k) * au(k-1)
197       x(k) = r(k) - al\_piv(k) * x(k-1)
198 \textcolor{keywordflow}{    enddo}
199 
200     \textcolor{comment}{! Backward sweep}
201     x(n) = r(n) / pivot(n)  \textcolor{comment}{! This should be X(N) / pivot(N), but is OK if Al(N) = 0.}
202     \textcolor{keywordflow}{do} k = n-1,1,-1
203       x(k) = ( x(k) - au(k)*x(k+1) ) / pivot(k)
204 \textcolor{keywordflow}{    enddo}
205   \textcolor{keywordflow}{else}
206     \textcolor{comment}{! This is a more typical implementation of a tridiagonal solver than the one above.}
207     \textcolor{comment}{! It is mathematically equivalent but differs at roundoff, which can cascade up to larger values.}
208 
209     \textcolor{comment}{! Factorization and forward sweep}
210     i\_pivot = 1.0 / ad(1)
211     x(1) = r(1) * i\_pivot
212     \textcolor{keywordflow}{do} k = 2,n
213       c1(k-1) = au(k-1) * i\_pivot
214       i\_pivot = 1.0 / (ad(k) - al(k) * c1(k-1))
215       x(k) = (r(k) - al(k) * x(k-1)) * i\_pivot
216 \textcolor{keywordflow}{    enddo}
217     \textcolor{comment}{! Backward sweep}
218     \textcolor{keywordflow}{do} k = n-1,1,-1
219       x(k) = x(k) - c1(k) * x(k+1)
220 \textcolor{keywordflow}{    enddo}
221 
222 \textcolor{keywordflow}{  endif}
223 
\end{DoxyCode}
