\section{Решение треугольной системы $Rx=z$ и обращение матриц $R$ и $A$\label{s6.3++}}

Как отмечено в разд.~\ref{s6.2+}, мы часто заинтересованы в решении
уравнения
\begin{equation}\label{Rx=z}
Rx=z,
\end{equation}
где $R=R(n,n)$ --- верхнетреугольная невырожденная матрица. Если нужно иметь только решение $x$, то $R^{-1}$ (для $x=R^{-1}z$) вычислять не надо. 
Следующий алгоритм обратной подстановки позволяет вычислить решение $x$
непосредственно.

\s{\large Алгоритм.}$\!\!\!\!$  Для $j=n,n-1,\ldots,1$ вычислять
\begin{equation}\label{backsubst}
x(j)=\left( z(j)-\sum_{k=j+1}^{n}R(j,k)x(k) \right)/R(j,j).
\end{equation}
По сложности этот алгоритм почти такой же, как матричное умножение. 
Он допускает записывать $x(j)$ поверх $z(j)$, что очень удобно в приложениях.

Если все же требуется иметь матрицу  $U\triangleq R^{-1}$,\label{RM1} 
то ее можно
вычислить по алгоритму окаймления, основанному на следующем легко проверяемом тождестве для треугольных матриц:
\begin{equation}\label{R^{-1}=}
\left[\begin{array}{ccc}
R_j &  & y \\
0 &  & \sigma_{j+1} 
\end{array}
\right]^{-1}=
\left[\begin{array}{ccc}
R_j^{-1} &  & -R_j^{-1}y\sigma_{j+1}^{-1} \\
 0&  & \sigma_{j+1}^{-1} 
\end{array}
\right]=R_{j+1}^{-1}\,.
\end{equation}
Это соотношение позволяет вычислять обратную матрицу $U\triangleq R^{-1}$
рекуррентно, \ie, если $R_j^{-1}=U_j$, где $R_j$ обозначает верхнюю левую часть матрицы $R$, то
\begin{equation}\label{U_{j+1}=}
U_{j+1}=\left[\begin{array}{ccc}
U_j &  & -U_j\left[R(1,j+1),\ldots,R(j,j+1)  \right]^T\sigma_{j+1} \\
0 &  & \sigma_{j+1} 
\end{array}
\right],
\end{equation}
где $\sigma_{j+1}=1/R(j+1,j+1)$.  Полагая $U=R_{-1}$, этот результат представим в алгоритмической форме.

\s{\large Алгоритм.}$\!\!\!\!$ Задать начальное значение
\begin{equation}\label{U(1,1)=}
U(1,1)=1/R(1,1).
\end{equation}
Для $j=2,\ldots, n$ вычислять по формулам (\ref{U(j,j)=})  и (\ref{U(k,j)=}):
\begin{gather}
U(j,j)=1/R(j,j),    \label{U(j,j)=} \\
U(k,j)= -\left(\sum_{i=k}^{j-1}U(k,i)R(i,j) \right)U(j,j), \quad k=1,\ldots,j-1.   \label{U(k,j)=} 
\end{gather}
\rem $R^{-1}$ вычисляется по столбцам, при этом элементы матрицы $R^{-1}$ могут записываться поверх элементов исходной матрицы $R$.\smallskip

В справочных целях приведем примеры реализации данного алгоритма на языке FORTRAN. 
Эти примеры могут помочь студентам написать свои собственные программы на других языках высокого уровня при выполнении лабораторного проекта №~6, описание которого дано ниже в подразд.~\ref{s6.8}.

\s{\large Обращение верхнетреугольной матрицы: $U:=R^{-1}$.}

Реализуются формулы (\ref{U(1,1)=}),  (\ref{U(j,j)=}) и (\ref{U(k,j)=}). Если нужно, $U$ может замещать $R$.
\[
\begin{array}{lll}
\multicolumn{3}{l}{R(N,N), \quad U(N,N), \quad R  \mbox{ и } U \mbox{ --- верхнетреугольные матрицы}.\enspace}\\\hline
&U(1,1)=1./R(1,1)&\\
&\mbox{DO }20\enspace J=2,N&\\
&U(J,J)=1./R(J,J)&\\
&JM1= J-1&\\
&\mbox{DO }20\enspace  K=1,JM1&\\
&SUM=0.&\\
&\mbox{DO }10\enspace  I=K,JM1&\\
10& SUM=SUM-U(K,I)*R(I,J)  &\\
20&U(K,J)=SUM*U(J,J)&\\\hline
\end{array}
\]
\medskip

В случаях, когда важно или нужно экономить память компьютера, матрицы в программе объявляют как одномерные массивы (см.~подразд.~\ref{s5.2+}). 
Хотя в компьютере даже многомерно объявленные массивы всегда хранятся как одномерные, компилятор генерирует индексные выражения с операциями умножения и деления. 
Операции сложения и вычитания в компьютерах выполняются гораздо быстрее, поэтому индексы для доступа к элементам матриц следует программировать в рекуррентной инкрементной форме, экономя таким образом и время процессора  (табл.~\ref{t:effRM1}). В этой программе преобразование в треугольную форму выполняется отождествлением $J(J-1)/1+I$ с $(I,J)$. Рекуррентное инкрементное вычисление $KK$, $JJ$ и $KK$ экономит вычисления.
\medskip

Как отмечалось на с.~\pageref{RM1}, иногда требуется вычислять $R^{-1}$.
Такая ситуация появляется, если требуется найти $A^{-1}$, для которой уже выполнено преобразование $TA=R$, где $T=T^{(n-1)}$ по формуле (\ref{eq:trianHtrans}), так как в теореме
\ref{trianhaus} для этого случая $m=n$ и $A^{-1}=R^{-1}T$. 
Последнее означает, что то же самое ортогональное преобразование $T$
теперь надо применить к строкам матрицы $R^{-1}$, но уже в обратном порядке
следования элементарных преобразований, составляющих полное преобразование
$T=T^{(n-1)}$ по формуле (\ref{eq:trianHtrans}). Таким образом, возникает проблема запоминания элементарных преобразований, составляющих полное преобразование
$T=T^{(n-1)}$, чтобы позже можно было его применить в задаче отыскания $A^{-1}$ или же для решения уравнения $Ax=z$ с невырожденной матрицей $A$ после пребразования $TA=R$.


\begin{table}[h]
\caption{Эффективное обращение верхнетреугольной матрицы: $U:=R^{-1}$. Реализуются формулы (\ref{U(1,1)=}),  (\ref{U(j,j)=}) и (\ref{U(k,j)=}). 
Верхнетреугольные матрицы $R$ и $U$ хранятся как векторы размерности $N(N+1)/2$. Если нужно, $U$ может замещать $R$.}\label{t:effRM1}
\[
\begin{array}{llp{6cm}}\hline
\multicolumn{3}{c}{\rule{0pt}{14pt} R(N,N), \quad U(N,N), \quad R  \mbox{ и } U \mbox{ --- верхнетреугольные матрицы}.\enspace}\\\hline
&U(1)=1./R(1)& \copyright\quad $R(1)\equiv R(1,1)$ \\
&JJ=1  &  \\
&\mbox{DO }20\enspace J=2,N  &  \\
&\overline{JJ}=JJ  & \copyright\quad $\overline{JJ}=J(J-1)/2\hm=(J-1,J-1)$ \\
&JJ=JJ+J  & \copyright\quad $JJ=J(J+1)/2\hm=(J,J)$ \\
&U(JJ)=1./R(JJ)  &  \\
&JM1=J-1 &  \\
&KK=0  &  \\
&\mbox{DO }20\enspace K=1, JM1  &  \\
&KK=KK+K  & \copyright\quad $KK=K(K+1)/2$ \\
&KI=KK  &  \\
&SUM=0.  &  \\
&\mbox{DO }10\enspace I=K, JM1 &  \\
&SUM=SUM-U(KI)*R(\overline{JJ}+I)  & \copyright\quad $KI=(K,I),\enspace\overline{JJ}+1\hm=(I,J)$  \\
10&KI=KI+I  &\copyright\quad $KI=(K,I+1)$  \\
20&U(\overline{JJ}+K)=SUM*U(JJ)  & \copyright\quad $\overline{JJ}+K=(K,J)$ \\
\hline
\end{array}
\]
\end{table}
\medskip


Как видно из (\ref{annihilating1+}), для отражения любого вектора $y=T^{(k)}z$ от
гиперплоскости  $U_{\perp}$, заданной вектором  $u$, требуется иметь
сохраненными две величины: вектор  $u$ и скаляр $\alpha$. Поскольку
нули ниже диагонали, получающиеся в результате триангуляризации, хранить не нужно, это место можно отвести для сохранения вектора  $u$ (вместе с диагональю, поэтому диагональные элементы $s_k$ приходится хранить отдельно).  
Данное предложение схематически показано на~\ris\ref{f:TkA=}.
\begin{figure}
\input{TkA=.txt}
\caption{Сохранение преобразования $T$  и вычисление вектора $y=Tz$, $\forall y\in\Rset ^m$} \label{f:TkA=}
\end{figure}
Каким образом можно выполнить вычисление матрицы $A^{-1}$, показано
на \ris\ref{f:abcd} на конкретном примере размерностей, а именно: $m=4$, $n=4$.
\begin{figure}
\input{abcd.txt}
\caption{Вычисление матрицы $A^{-1}$ после сохранения преобразования $T$} \label{f:abcd}
\end{figure}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%