% -*- mode: Noweb; noweb-code-mode: c-mode -*-

%% It is going to be included in a section

\subsection{Gaussian Elimination}
\label{sec:gaussianeliminate}

%\begin{fullwidth}
Gaussian elimination is a method for computing the solution vector
$\vec{x}$ of a given system of linear equations.
\[A\vec{x}=\vec{b}\]
First, a reduction step transforms the augmented matrix $[A | \vec{b}]$ into
upper-triangular form.
Then, a second solving step produces the solution $\vec{x}$ by back-substitution.

%\end{fullwidth}
Let's get on with the code.
\thanks{Thanks to \hlink{Lucky's Notes}{%
  https://luckytoilet.wordpress.com/2010/02/21/solving-systems-of-linear-equations-in-haskell/}  
for the structure of the code.}
Some types and helper functions will be needed.
The reduction and solver steps must be constructed.  
And, a {\tt{}main} routine is needed as an entry point.
Here is the structure of the code.  It is written in \mynoweb style. 
\begin{fullwidth}
\nwfilename{Gauss.nw}\nwbegincode{1}\sublabel{NWGau8-Gau5-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-Gau5-1}}}\moddef{Gauss~{\nwtagstyle{}\subpageref{NWGau8-Gau5-1}}}\endmoddef
module Gauss ( gaussianReduce ) where
\LA{}Row, Column, and Matrix Types~{\nwtagstyle{}\subpageref{NWGau8-RowT-1}}\RA{}
\LA{}IO conversion functions~{\nwtagstyle{}\subpageref{NWGau8-IO*N-1}}\RA{}
\LA{}The Gaussian reduction step~{\nwtagstyle{}\subpageref{NWGau8-TheR-1}}\RA{}
\LA{}The Gaussian solver step~{\nwtagstyle{}\subpageref{NWGau8-TheO-1}}\RA{}
\LA{}Main module~{\nwtagstyle{}\subpageref{NWGau8-MaiB-1}}\RA{}
\nwnotused{Gauss}\nwendcode{}\nwbegindocs{2}\nwdocspar
Define the {\tt{}Row} and {\tt{}Column} types to be lists of double
precision floating point numbers.
Define the {\tt{}RMatrix} type to be a list indexed by rows.
Define the {\tt{}CMatrix} type to be a list indexed by columns.

\nwenddocs{}\nwbegincode{3}\sublabel{NWGau8-RowT-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-RowT-1}}}\moddef{Row, Column, and Matrix Types~{\nwtagstyle{}\subpageref{NWGau8-RowT-1}}}\endmoddef
type Row = [Double]
type Column = [Double]
type RMatrix = [Row] 
type CMatrix = [Column] 
\nwused{\\{NWGau8-Gau5-1}}\nwendcode{}\nwbegindocs{4}\nwdocspar
\subsection{The Reduction Step}
\label{subsec:reduce}

The \emph{reduction} process used in Gaussian elimination transforms a 
matrix into upper-triangular form. Pictorially, for a small example.
\[
\begin{bmatrix}
  \circled{p} & \cdot & \cdot & \cdot & |\cdot\\
  a & \cdot & \cdot & \cdot & |\cdot\\
  b & \cdot & \cdot & \cdot & |\cdot\\
  c & \cdot & \cdot & \cdot & |\cdot
\end{bmatrix}\mapsto
\begin{bmatrix}
  1 & \cdot & \cdot & \cdot & |\cdot\\
  0 & \circled{p} & \cdot & \cdot & |\cdot\\
  0 & \cdot & \cdot & \cdot & |\cdot\\
  0 & \cdot & \cdot & \cdot & |\cdot
\end{bmatrix}\mapsto
\begin{bmatrix}
  1 & \cdot & \cdot & \cdot & |\cdot\\
  0 & 1 & \cdot & \cdot & |\cdot\\
  0 & 0 & \circled{p} & \cdot & |\cdot\\
  0 & 0 & \cdot & \cdot & |\cdot
\end{bmatrix}\mapsto
\begin{bmatrix}
  1 & \cdot & \cdot & \cdot & |\cdot\\
  0 & 1 & \cdot & \cdot & |\cdot\\
  0 & 0 & 1 & \cdot & |\cdot\\
  0 & 0 & 0 & 1 & |\cdot
\end{bmatrix}\]
To reduce a matrix to upper-triangular form repeat these 
steps iteratively across all rows of the {\tt{}matrix}.
For the first row:
\begin{enumerate}
\item  Assume $p$, the {\tt{}pivot}, is not $0$ and normalize the row by
  scaling it by $1/p$.  
\item Repeatedly, multiply the normalized row by $a,\,b,\,c$ and
  subtract the result row from second, third, and fourth rows, respectively.
\end{enumerate}
\nwenddocs{}\nwbegindocs{5}\nwdocspar
\nwenddocs{}\nwbegincode{6}\sublabel{NWGau8-TheR-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-TheR-1}}}\moddef{The Gaussian reduction step~{\nwtagstyle{}\subpageref{NWGau8-TheR-1}}}\endmoddef
gaussianReduce :: RMatrix -> RMatrix
gaussianReduce matrix = foldl reducerow matrix [0..length matrix-1] where
  reducerow :: RMatrix -> Int -> RMatrix
  reducerow m r = let 
    \LA{}Pick the \code{}row\edoc{} to reduce by, its \code{}pivot\edoc{}, and normalize this pivot \code{}row\edoc{}~{\nwtagstyle{}\subpageref{NWGau8-Pic1G-1}}\RA{}
    \LA{}Construct a function that reduces other rows~{\nwtagstyle{}\subpageref{NWGau8-Coni-1}}\RA{}
    \LA{}Apply the reduction function to rows below the \code{}pivot\edoc{}~{\nwtagstyle{}\subpageref{NWGau8-Appu-1}}\RA{}
    \LA{}Piece the matrix back together~{\nwtagstyle{}\subpageref{NWGau8-PieU-1}}\RA{}
\nwused{\\{NWGau8-Gau5-1}}\nwendcode{}\nwbegindocs{7}\nwdocspar
\begin{itemize}
\item Pick out the {\tt{}r}-th {\tt{}row} of matrix {\tt{}m}, using {\tt{}!!},
  Haskell's indexing operator.
\item Pick out the {\tt{}r}-th value in the {\tt{}r}-th row and call it $p$,
  the \emph{pivot}.
  To keep things simple, assume these pivot elements along the main
  diagonal never equal $0$. 
\item Normalize the row by dividing each element in it by the
  pivot $p$.
  The Haskell \emph{idiom} that does this is to map the 
  \emph{divide by p)} function {\tt{}(/\ p)} across the row.
  Call the normalized row {\tt{}row'}.
\end{itemize}
\nwenddocs{}\nwbegincode{8}\sublabel{NWGau8-Pic1G-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-Pic1G-1}}}\moddef{Pick the \code{}row\edoc{} to reduce by, its \code{}pivot\edoc{}, and normalize this pivot \code{}row\edoc{}~{\nwtagstyle{}\subpageref{NWGau8-Pic1G-1}}}\endmoddef
row = m !! r
p = row !! r  
row' = map (/ p) row  
\nwused{\\{NWGau8-TheR-1}}\nwendcode{}\nwbegindocs{9}\nwdocspar
Now, to reduce another row, say {\tt{}nrow}, by {\tt{}row'} apply the function
\[nr * a - b\tag{where {\tt{}nr} is the {\tt{}r}-th element in {\tt{}nrow}}\]
to each entry $a$ in {\tt{}row} and $b$ in {\tt{}nrow}.
The Haskell \emph{idiom} for this is to \emph{zip} the function with the
values in {\tt{}row'} and {\tt{}nrow}.
\nwenddocs{}\nwbegincode{10}\sublabel{NWGau8-Coni-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-Coni-1}}}\moddef{Construct a function that reduces other rows~{\nwtagstyle{}\subpageref{NWGau8-Coni-1}}}\endmoddef
reduceonerow nrow = let nr = nrow !! r in zipWith (\\a b -> nr*a - b) row' nrow
\nwused{\\{NWGau8-TheR-1}}\nwendcode{}\nwbegindocs{11}\nwdocspar
Next, map {\tt{}reduceRow} across all rows in
{\tt{}matrix} below the {\tt{}pivot} row.
\nwenddocs{}\nwbegincode{12}\sublabel{NWGau8-Appu-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-Appu-1}}}\moddef{Apply the reduction function to rows below the \code{}pivot\edoc{}~{\nwtagstyle{}\subpageref{NWGau8-Appu-1}}}\endmoddef
nextrows = map reduceonerow (drop (r+1) m)
\nwused{\\{NWGau8-TheR-1}}\nwendcode{}\nwbegindocs{13}\nwdocspar
And finally, piece the results back together:
Concatenate the first {\tt{}r} rows from {\tt{}m},
row {\tt{}row'} and the reduces {\tt{}nextrows}.
\nwenddocs{}\nwbegincode{14}\sublabel{NWGau8-PieU-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-PieU-1}}}\moddef{Piece the matrix back together~{\nwtagstyle{}\subpageref{NWGau8-PieU-1}}}\endmoddef
in take r m ++ [row'] ++ nextrows
\nwused{\\{NWGau8-TheR-1}}\nwendcode{}\nwbegindocs{15}\nwdocspar

To turn {\tt{}gaussianReduce} into a standalone program, a {\tt{}main}
module, an \emph{entry} point, must be established.
IO can be tricky.  The structure of the input must be known, and that
must be translated into the structure of the Haskell function the
implements the program.

Let's agree, for the purpose of these notes, that the program is
executed with input redirected from \emph{standard input}, like so
\begin{quote}
  \texttt{./Gauss < gauss.dat}
\end{quote}
Assume the contents of \texttt{gauss.dat} is a string of lines separated by
newline characters like so:
\[\begin{array}{@{}ccccccc@{}}
    a_{00} & a_{01} & a_{02} & \ldots & a_{0(n-1)} & b_{0} & \textbackslash n \\
    a_{10} & a_{11} & a_{12} & \ldots & a_{1(n-1)} & b_{1} & \textbackslash n \\
   \vdots  & \vdots & \vdots & \vdots & \vdots & \vdots  \\
    a_{n-1,0} & a_{n-1,1} & a_{n-1,2} & \ldots & a_{n-1,n-1} & b_{n-1} & \textbackslash n 
  \end{array}\]
Let's agree to use the Haskell \emph{idiom} {\tt{}raw\ <-\ getContents}
to input all of \texttt{gauss.dat} into one {\tt{}String} called {\tt{}raw}.
\begin{itemize}
\item The {\tt{}lines} function breaks {\tt{}raw} up into a list of strings
  {\tt{}[String]}, separated at the newline. Call the result {\tt{}rows}. 
\item Mapping the {\tt{}words}function over {\tt{}rows} breaks
  each string in {\tt{}rows} list of {\tt{}String}s. These values need to be
  converted from {\tt{}String} to {\tt{}Double}.
\end{itemize}

\nwenddocs{}\nwbegincode{16}\sublabel{NWGau8-MaiB-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-MaiB-1}}}\moddef{Main module~{\nwtagstyle{}\subpageref{NWGau8-MaiB-1}}}\endmoddef
main :: IO () 
main = do
  raw <- getContents
  let rows = lines raw
  let rows' = map words rows
  let matrix = map stringsToDoubles rows'
  print $ gaussianReduce matrix
\nwused{\\{NWGau8-Gau5-1}}\nwendcode{}\nwbegindocs{17}\nwdocspar

\nwenddocs{}\nwbegincode{18}\sublabel{NWGau8-IO*N-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-IO*N-1}}}\moddef{IO conversion functions~{\nwtagstyle{}\subpageref{NWGau8-IO*N-1}}}\endmoddef
stringToDouble = read :: String -> Double

stringsToDoubles :: [String] -> [Double]
stringsToDoubles = map stringToDouble 
\nwused{\\{NWGau8-Gau5-1}}\nwendcode{}\nwbegindocs{19}\nwdocspar

\end{fullwidth}

\subsection{Test the Reduction Step}

\label{subsec:testreduce}

\todo[noline]{Assignment, Part 1: Complete the steps in this section. (25 points)}

\hypertarget{subsec:testreduce}{The \mynoweb source file \texttt{Gauss.nw} is
\hlink{here}{http://cs.fit.edu/~wds/classes/aa/Code/Gauss.nw}.}
\begin{itemize}
\item Running \texttt{notangle Gauss.nw > Gauss.hs} generates the Haskell
  code. 
\item Running \texttt{noweave -index -delay Gauss.nw > Gauss.tex} generates
  a \LaTeX~file that is included in a \emph{wrapper} \LaTeX~file,
  that produces this document.
\end{itemize}
You may not have the \mynoweb tools. You can retrieve the Haskell code
from a link in the first step of this part of the assignment.
\begin{fullwidth}
\begin{enumerate}
\item  Download \hlink{Gauss.hs}{http://cs.fit.edu/~wds/classes/aa/Code/Gauss.hs}.
\item Install the Glasgow Haskell Compiler and use its interpreter
  \texttt{ghci} to load \texttt{Gauss.hs}. Check the reduction
  code is correct on some simple cases, for instance,
  \begin{itemize}
  \item {\tt{}*Gauss>\ gaussianReduce\ [[]]}  -- the empty matrix
  \item {\tt{}*Gauss>\ gaussianReduce\ [[1]]} -- a $1\times1$ matrix that needs no normalization
  \item {\tt{}*Gauss>\ gaussianReduce\ [[2]]} -- a $1\times1$ matrix that needs normalization
  \item {\tt{}*Gauss>\ gaussianReduce\ [[1,2]]} -- the equation $2x=1$
  \item {\tt{}*Gauss>\ gaussianReduce\ [[1,2],[2,3]]} -- an inconsistent system
  \item {\tt{}*Gauss>\ gaussianReduce\ [[1,-2,1,4],[2,3,-1,5],[3,1,4,7]]} -- a picked out of the air example
  \end{itemize}
\item Write code to generate $n$ lists of random {\tt{}Double}s of length
  $n+1$ that represents  a linear system $Ax=b$, where $A$ is an
  $n\times n$ matrix and $b$ is an $n\times1$ vector.  
\item Run your data generation code to generate several data files of
  varying sizes $n\times n+1$.
\item Compile the Haskell source
  \hlink{Gauss.hs}{http://cs.fit.edu/~wds/classes/aa/Code/Gauss.hs} 
  with \texttt{ghc} profiling options, see the \texttt{ghc}
  \hlink{User's Guide on Profiling}{%
    http://downloads.haskell.org/~ghc/latest/docs/html/users_guide/profiling.html}
  for instructions on this. 
\item Execute {\tt{}Gauss} on your data files, and collect running times.
\item Plot the running times.  Find a curve that approximates
  the data.  The $R^{2}$ value for the approximation should be close
  to $1$.
\end{enumerate}

\end{fullwidth}
\subsection{Analyze the Reduction step}
\label{subsec:analyzereduce}

\todo[noline]{Assignment, Part 2: Perform a mathematical analysis of the
  reduction algorithm. (25 points)}

A big-$O$ time complexity can be computed for each function in the
code.  For instance, the \emph{anonymous} function {\tt{}({\char92}a\ b\ ->\ nr*a\ -\ b)}
has time complexity $O(1)$.  The time cost to map it over a list of
length $n+1$ is $O(n)$.

\begin{fullwidth}
What are the big-$O$ time complexities for the functions below.
Explain your reasoning for each function.
In particular, identify the size and type for the input and output of
each function.
\begin{itemize}
\item {\tt{}stringsToDoubles}
\item {\tt{}map\ stringsToDoubles}
\item {\tt{}reduceonerow}
\item {\tt{}reduceRow}
\item {\tt{}take\ r\ m\ ++\ [row']\ ++\ nextrows}
\item {\tt{}gaussianReduce}
\end{itemize}
Does the empirical run time data from experiments in
you performed when 
\hyperlink{subsec:testreduce}{testing {\tt{}gaussianReduce}} 
agree with the analytic analysis?
  
\end{fullwidth}
\subsection{The Gaussian Solver}
\label{subsec:solver}

\todo[noline]{Assignment, Part 3: Complete the Gaussian elimination algorithm
  by implementing the solver step. (25 points)}

The solving step in Gaussian elimination uses \emph{back substitution}
to solve for the values in $\vec{x}$ in order
$x_{n-1},\,x_{n-2},\ldots,\,x_{0}$. Write Haskell code that implements
the solver step.

\begin{fullwidth}
\nwenddocs{}\nwbegincode{20}\sublabel{NWGau8-TheO-1}\nwmargintag{{\nwtagstyle{}\subpageref{NWGau8-TheO-1}}}\moddef{The Gaussian solver step~{\nwtagstyle{}\subpageref{NWGau8-TheO-1}}}\endmoddef
-- gaussianSolve :: RMatrix -> [Double]
-- Your code goes here

\nwused{\\{NWGau8-Gau5-1}}\nwendcode{}

\nwixlogsorted{c}{{Apply the reduction function to rows below the \code{}pivot\edoc{}}{NWGau8-Appu-1}{\nwixu{NWGau8-TheR-1}\nwixd{NWGau8-Appu-1}}}%
\nwixlogsorted{c}{{Construct a function that reduces other rows}{NWGau8-Coni-1}{\nwixu{NWGau8-TheR-1}\nwixd{NWGau8-Coni-1}}}%
\nwixlogsorted{c}{{Gauss}{NWGau8-Gau5-1}{\nwixd{NWGau8-Gau5-1}}}%
\nwixlogsorted{c}{{IO conversion functions}{NWGau8-IO*N-1}{\nwixu{NWGau8-Gau5-1}\nwixd{NWGau8-IO*N-1}}}%
\nwixlogsorted{c}{{Main module}{NWGau8-MaiB-1}{\nwixu{NWGau8-Gau5-1}\nwixd{NWGau8-MaiB-1}}}%
\nwixlogsorted{c}{{Pick the \code{}row\edoc{} to reduce by, its \code{}pivot\edoc{}, and normalize this pivot \code{}row\edoc{}}{NWGau8-Pic1G-1}{\nwixu{NWGau8-TheR-1}\nwixd{NWGau8-Pic1G-1}}}%
\nwixlogsorted{c}{{Piece the matrix back together}{NWGau8-PieU-1}{\nwixu{NWGau8-TheR-1}\nwixd{NWGau8-PieU-1}}}%
\nwixlogsorted{c}{{Row, Column, and Matrix Types}{NWGau8-RowT-1}{\nwixu{NWGau8-Gau5-1}\nwixd{NWGau8-RowT-1}}}%
\nwixlogsorted{c}{{The Gaussian reduction step}{NWGau8-TheR-1}{\nwixu{NWGau8-Gau5-1}\nwixd{NWGau8-TheR-1}}}%
\nwixlogsorted{c}{{The Gaussian solver step}{NWGau8-TheO-1}{\nwixu{NWGau8-Gau5-1}\nwixd{NWGau8-TheO-1}}}%
\nwbegindocs{21}\nwdocspar
\subsection{Test the Solver Step}
\label{subsec:testsolve}

Test your {\tt{}gaussianSolve} code following steps similar to those
outlined in \hyperlink{subsec:testreduce}{testing {\tt{}gaussianReduce}}.

\end{fullwidth}
\subsection{Analyze the Solver step}
\label{subsec:analyzesolve}

\todo[noline]{Assignment, Part 4: Perform empirical and mathematical analyses of the
  solver algorithm. (25 points)}

\begin{fullwidth}
\begin{enumerate}
\item What is the big-$O$ time complexity of your {\tt{}gaussianSolve} code?
\item Write a function that tests the accuracy of the computed
  solution $\vec{x}'$.
  \begin{enumerate}
  \item Use the $L_{\infty}$ norm to measure and report solution accuracy.
    \[\|b'-b\|_{\infty} = \max\Set{\abs{b_{0}'-b_{0}},\,\abs{b_{1}'-b_{1}},\ldots,\,\abs{b_{n-1}'-b_{n-1}}}\]
    where $b'=A\vec{x}'$.
  \item Hilbert matrices $H_{n}$ are famous examples of \emph{ill-conditioned}
    matrices.  Informally, this means solving $H_{n}\vec{x}=b$
    accurately is difficult. 
    The Hilbert $H_{5}$ is
    \[H_{5}={\begin{bmatrix}
    1&{\frac {1}{2}}&{\frac{1}{3}}&{\frac {1}{4}}&{\frac {1}{5}}\\
    {\frac{1}{2}}&{\frac {1}{3}}&{\frac {1}{4}}&{\frac {1}{5}}&{\frac{1}{6}}\\
    {\frac {1}{3}}&{\frac {1}{4}}&{\frac{1}{5}}&{\frac {1}{6}}&{\frac {1}{7}}\\
    {\frac{1}{4}}&{\frac {1}{5}}&{\frac {1}{6}}&{\frac {1}{7}}&{\frac{1}{8}}\\
    {\frac {1}{5}}&{\frac {1}{6}}&{\frac{1}{7}}&{\frac {1}{8}}&{\frac {1}{9}}\end{bmatrix}}\]
   In general, the value in row $i$, column $j$ is 
   \[(H_{n})_{ij}=\frac{1}{i+j+1},\quad i,\,j=0,\ldots,n-1\]
   Test your code on Hilbert matrices.
  \end{enumerate}
\end{enumerate}

\bibliography{cs}
\end{fullwidth}

%In particular, identify the size and type for the input and output of
%{\tt{}gaussianSolve} and internal functions it calls.

\nocite{Lipo:11}
\nocite{Haskell98}
\nocite{ghcu:14}
\nwenddocs{}
