\documentclass[11pt]{article}
\usepackage{fullpage,amsmath,amsfonts,mathpazo,microtype,nicefrac,algorithm2e,graphicx,siunitx}

% Set-up for hypertext references
\usepackage{hyperref,color,textcomp}
\definecolor{webgreen}{rgb}{0,.35,0}
\definecolor{webbrown}{rgb}{.6,0,0}
\definecolor{RoyalBlue}{rgb}{0,0,0.9}
\hypersetup{
   colorlinks=true, linktocpage=true, pdfstartpage=3, pdfstartview=FitV,
   breaklinks=true, pdfpagemode=UseNone, pageanchor=true, pdfpagemode=UseOutlines,
   plainpages=false, bookmarksnumbered, bookmarksopen=true, bookmarksopenlevel=1,
   hypertexnames=true, pdfhighlight=/O,
   urlcolor=webbrown, linkcolor=RoyalBlue, citecolor=webgreen,
   pdfauthor={Chris H. Rycroft},
   pdfsubject={Harvard AM225 (Spring 2021)},
   pdfkeywords={},
   pdfcreator={pdfLaTeX},
   pdfproducer={LaTeX with hyperref}
}
\hypersetup{pdftitle={AM225: Assignment 3}}
\usepackage[font=footnotesize,skip=3pt]{caption}

\graphicspath{ {am225_hw3_files/fem_diagrams/} {am225_hw3_files/pss/} }

% Macro definitions
\newcommand{\N}{\mathbb{N}}
\newcommand{\Z}{\mathbb{Z}}
\newcommand{\Q}{\mathbb{Q}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\p}{\partial}
\renewcommand{\vec}[1]{\mathbf{#1}}
\newcommand{\vx}{\vec{x}}
\newcommand{\vv}{\vec{v}}
\newcommand{\vp}{\vec{p}}
\newcommand{\rs}{\sqrt{6}}
\newcommand{\Atol}{\textit{Atol}}
\newcommand{\Rtol}{\textit{Rtol}}
\newcommand{\Trans}{\mathsf{T}}

\DeclareMathOperator{\atantwo}{atan2}

\renewcommand{\topfraction}{0.9}
\renewcommand{\bottomfraction}{0.8}
\renewcommand{\textfraction}{0.2}
\renewcommand{\floatpagefraction}{0.7}
\renewcommand{\dblfloatpagefraction}{0.7}

\begin{document}
\section*{AM225: Assignment 3 (due 5~PM, March 25)}
\noindent \textit{Complete \textbf{at least two} questions. If you submit
answers for more, your grade will be calculated using the best scores. The
questions are roughly ordered by difficulty, from easiest to hardest.}

\begin{enumerate}
  \item \textbf{Strassen's algorithm.} The standard approach for multiplying two
    $n\times n$ matrices together, $AB=C$, requires $O(n^3)$ operations. In 1969,
    Strassen discovered an algorithm that requires asymptotically fewer operations.
    Suppose $n$ is even, and write the matrix product as
    \begin{equation}
      \left(
      \begin{array}{cc}
	A_{00} & A_{01} \\
	A_{10} & A_{11}
      \end{array}
      \right)
      \left(
      \begin{array}{cc}
	B_{00} & B_{01} \\
	B_{10} & B_{11}
      \end{array}
      \right)
      =
      \left(
      \begin{array}{cc}
	C_{00} & C_{01} \\
	C_{10} & C_{11}
      \end{array}
      \right)
    \end{equation}
    where each term is a submatrix of size \smash{$\tfrac{n}{2} \times \tfrac{n}{2}$}.
    If the seven quantities
    \begin{align}
      Q_0 &= (A_{00}+A_{11})(B_{00}+B_{11}), \\
      Q_1 &= (A_{10}+A_{11})B_{00}, \\
      Q_2 &= A_{00} (B_{01}-B_{11}), \\
      Q_3 &= A_{11} (B_{10}-B_{00}), \\
      Q_4 &= (A_{00}+A_{01}) B_{11}, \\
      Q_5 &= (A_{10}-A_{00})(B_{00}+B_{01}), \\
      Q_6 &= (A_{01}-A_{11})(B_{10}+B_{11})
    \end{align}
    are computed, then the entries in $C$ are given by
    \begin{align}
      C_{00} &= Q_0 + Q_3 - Q_4 + Q_6, \\
      C_{10} &= Q_1 + Q_3, \\
      C_{01} &= Q_2 + Q_4, \\
      C_{11} &= Q_0 + Q_2 - Q_1 + Q_5.
    \end{align}
    This procedure can be applied recursively. Since the original problem of
    dimension $n$ is replaced by seven problems of dimension
    \smash{$\tfrac{n}{2}$}, the total number of operations required is
    \smash{$O(n^{\log_2 7}) = O(n^{2.807})$}. While Strassen's method requires
    more operations than the standard implementation, the improvement in
    exponent will eventually make it worthwhile for large matrices.
    \begin{enumerate}
      \item Write a program to implement the standard matrix multiplication
	algorithm. Time your code on random matrices of dimensions that are
	powers of two from $n=16$ upward.\footnote{Try to make $n$ as large
	as your system will allow. For example, for $n=8192$, you will need at
	least $3 \times 8192^2 \times (8 \text{~bytes}) = 403\text{~MB}$ of
	memory to represent the three matrices.}
      \item Implement Strassen's algorithm, and time your code on random
	matrices of dimensions that are powers of two from $n=16$ upward.
	Compare your results to part (a). Fit your timing results to a power
	law $\alpha n^\beta$ and report the values of $\alpha$ and
	$\beta$.\footnote{You should find $\beta<3$ as expected.}
      \item Time the \href{http://www.netlib.org/blas/}{BLAS} matrix--matrix
	multiply routine, \texttt{dgemm}, and compare your results to the
	previous sections. You should find that the BLAS routine is faster, but
	has $O(n^3)$ scaling.
      \item \textbf{Optional.} Extend your codes from parts (a) and (b) to use
	multithreading and repeat the timing comparison.
      \item \textbf{Optional.} Extend your code from part (b) to work on
	matrices that are not powers of two.
      \item \textbf{Optional.} It is difficult to beat BLAS for matrices of
	practical size, due to its fine tuning of caching performance. Could a
	combination of Strassen and BLAS be faster than BLAS itself?
    \end{enumerate}
  \item \textbf{Fractal-based preconditioning} In lecture 9, the preconditioned
    conjugate gradient (PCG) algorithm was introduced, and tested on a linear
    system to find radial basis function weights. The program files contain a
    program called \texttt{rbf\_time.cc}, which compares the speed of solving
    the system with LAPACK and PCG for $k$ points varying from 10 to 10000.
    This code uses the compact Wendland radial function with cutoff radius set
    to \smash{$r(k)=5/\sqrt{k}$}, which results in approximately 15 entries in
    each line of the linear system. The PCG uses a block Jacobi preconditioner
    with blocks of size \smash{$b(k)=\sqrt{k}$}.\footnote{In many cases $b(k)$
    will not exactly divide $k$. In this case a final block is used with size
    less than $b(k)$.} Each block is inverted using dense linear algebra via
    LAPACK.

    Ideally, we want the blocks to be as dense as possible, and thus it is
    beneficial to order the points in the linear system so that they are
    physically nearby. The current code orders the points by the $y$ coordinate
    (see \texttt{rbf::init\_random}), which helps, but points can still be
    delocalized in the $x$ direction, resulting in sparse Jacobi blocks.
    \begin{enumerate}
      \item Write a program to calculate the total non-zero matrix entries
	$T(k)$ and the non-zero matrix entries in the Jacobi blocks
	$P(k)$.\footnote{You could do this by adding your own function to the
	\texttt{rbf} class.} Make a plot of the proportion of non-zero matrix
	entries in the Jacobi blocks $P(k)/T(k)$ as a function of $k$, for a
	variety of values from $k=10$ to $k=10000$.
      \item The Hilbert curve is a space-filling fractal that bijectively maps
	a one-dimensional coordinate $q$ continuously into a square with
	coordinates $\vec{x}=(x,y)$. Detailed information about the construction of
	Hilbert curves is available online
	(\textit{e.g.}~\href{https://en.wikipedia.org/wiki/Hilbert_curve}{Wikipedia}).
	For the given example, the sample points $\vec{x}_i$ lie in the square
	$[-1,1]^2$. Write a program that computes the 1D coordinate $q_i$
	for each point, and then re-sorts the points based on these $q$ values,
	thereby ensuring that points are physically close together in the
	indexing. Make a corresponding plot of $P(k)/T(k)$ like part (a)
	using this point initialization method.
      \item Use the program \texttt{rbf\_time.cc} to compare the time of PCG
	for the original initialization method and the fractal-based method,
	over the range from $k=10$ to $k=10000$.
      \item \textbf{Optional.} Try varying the functions $r(k)$ and $b(k)$ to
	see how they affect your results in the previous parts.
      \item \textbf{Optional.} Are there better ways to organize the points
	than the Hilbert curve?
    \end{enumerate}
  \item \textbf{Alternative cubic elements.} In the lectures, we considered
    solving the elliptic problem
    \begin{equation}
      -\frac{d}{dx} \left( x \frac{du}{dx} \right) = f(x) \label{eq:exprob1}
    \end{equation}
    for a function $u(x)$ on $\Omega=(1,2)$ with the boundary conditions
    \begin{equation}
      u(1)=0, \qquad \left. \frac{du}{dx} \right|_{x=2} = g \label{eq:exprob2}
    \end{equation}
    where $g$ is a real constant. The interval was divided into $N$
    subintervals of equal length, and the finite elements were constructed from
    cubics with a nodal basis. Following the notation introduced the lectures,
    this corresponds to a standard element shown below. \vspace{-1em}
    \begin{center}
      \includegraphics[width=0.3\textwidth]{fe_1d_cubic}
    \end{center}\vspace{-12em}
    Modify this example code to solve the problem using the alternative cubic
    basis whereby the function values and first derivatives are specified,
    which corresponds to the standard element shown below.
    \vspace{-1.5em}
    \begin{center}
      \includegraphics[width=0.3\textwidth]{fe_1d_cubic2}
    \end{center}\vspace{-12em}
    With this basis, your finite element functions will be in $C^1(\Omega)$.
    \begin{enumerate}
      \item For $N=3$, make a plot of the finite element basis functions
	that is similar to Fig.~1 in the
	\href{https://courses.seas.harvard.edu/courses/am225/notes/am225_fe_example.pdf}{example
	notes}. In the notes, the finite element problem on $N$ intervals had a
	basis with $3N$ functions.\footnote{In the notes, there was one
	additional function $\psi_0$ which was included for completeness, but
	was not considered the finite-element linear system since it was
	constrained to zero by the Dirichlet condition at $x=1$.} For this
	alternative basis, what is the dimension of the finite element function
	space when $N$ intervals are used?
      \item Write a program to solve the elliptic problem using the alternative
	cubic basis.\footnote{You can either write your own code, or modify the
	one from the example.} Using the manufactured solution considered in
	the notes, make a log--log plot of the $L_2$ error of the numerical
	solution as a function of $N$ for a range of values from $N=10$ to
	$N=1000$. On your plot, include also the $L_2$ errors for the original
	code with the nodal basis.
      \item Make another log--log plot comparing the $L_2$ error for the two
	methods, but for this case, put the wall clock time for the computation
	on the horizontal axis. Which basis is more efficient for achieving a
	desired level of accuracy?
    \end{enumerate}
  \item \textbf{Simulating an MRI.\footnote{This question was written by Dan
    Fortunato (TF, Spring 2018).}} Magnetic resonance imaging (MRI) is a
    noninvasive medical imaging technique in which a body is subject to a
    magnetic field operating at a frequency $f$. A simple model of an MRI is
    given by the Helmholtz equation
    \begin{equation}\label{eq:helmholtz}
      \nabla^2 u(\vx) + k^2 u(\vx) = v(\vx),
    \end{equation}
    where $u(\vx)$ is the electric field strength, $k$ is the wave number, and
    $v(\vx)$ is the electric excitation. In general, both $u(\vx)$ and $k$ can
    be complex-valued. Note that $k^2 = \omega^2 \mu \epsilon$, where $\omega =
    2\pi f$, $\mu$ is the magnetic permeability, and $\epsilon$ is the electric
    permittivity. In free space, $k^2 = k_0^2 = \omega^2 \mu_0 \epsilon_0$,
    where $\mu_0 = 4\pi\times\SI{e-7}{\henry\per\meter}$ and $\epsilon_0 =
    \SI{8.8542e-12}{\farad\per\meter}$.
    \begin{enumerate}
      \item Implement a second-order accurate finite difference scheme for
	Eq.~\ref{eq:helmholtz} in free space on the domain $[0,1]^2$ meters with
	homogeneous Dirichlet boundary conditions. Set $v(\vx)$ to be the
	impulse $v(\vx) = \delta(\vx - (0.5,0.5))$ and use a frequency of $f =
	\SI{21.3}{\mega\hertz}$. Assume $u(\vx)$ is complex-valued.
	Solve the resulting system using (i) a direct banded solver in
	LAPACK\footnote{If you index the grid in the
	conventional way, \textit{i.e.} $(i,j)$ maps to $i+(N+1)j$, then the matrix
	you need to solve is banded with $N+1$ lower and $N+1$ upper diagonals.
	In this case, there is a LAPACK routine \texttt{zgbsv} that provides a
	huge saving in memory and execution time over a full dense solver. See
	\texttt{banded\_test.cc} for an example.} and (ii) a fast
	solver based on the discrete sine transform\footnote{Use FFTW
	(\url{www.fftw.org}) or another library to perform the FFTs.} (DST).
	Make a single log--log plot of the timings of both methods using $N+1$
	gridpoints in each dimension for $N =
	16,\,32,\,\ldots,\,1024$.\footnote{You should have no problem in using
	$N=1024$ with the DST method. However, the direct solver may become too
	expensive, so you may skip some of the larger grid sizes for this
	case.}
      \item Human tissue can be modeled as a material with $\mu=\mu_0$ and
	$\epsilon = \epsilon(\vx)$, causing $k = k(\vx)$ to vary in space:
        \begin{equation}\label{eq:helmholtz_variable}
          \nabla^2 u(\vx) + k(\vx)^2 u(\vx) = v(\vx).
        \end{equation}
	Use method (i) to simulate an MRI of a human head with $N=256$ and
	$v(\vx) = \delta(\vx - (0.6,0.7))$, and show the resulting MRIs by plotting $|u|$ for
	$f=\SI{21.3}{\mega\hertz}$ and $f=\SI{298.3}{\mega\hertz}$. To model
	the head, use the provided binary data \texttt{mri\_data\_256.dat},
	which gives values for $\epsilon(\vx) / \epsilon_0$ as an $(N+1) \times
	(N+1)$ matrix of complex numbers stored in row-major order.
      \item Does the DST-based method work for Eq.~\ref{eq:helmholtz_variable}?
        Explain why or why not.
      \item As $N$ gets large, method (i) becomes computationally infeasible,
	and iterative methods are used instead. A simple iterative method for
	solving Eq.~\ref{eq:helmholtz_variable} that can utilize the DST-based
	solver is
	\begin{equation}
	  \nabla^2 u^{(i+1)} + k_0^2 u^{(i+1)} = v - \left(k(\vx)^2 - k_0^2\right) u^{(i)},
        \end{equation}
	where $u^{(i)}$ denotes the solution at iteration $i$ and $u^{(0)} =
	0$. Implement this method and produce a log--log plot of the
	$L^\infty$-norm relative error at each iteration, taking the
	solutions from (b) as the true solutions for $f=\SI{21.3}{\mega\hertz}$
	and $f=\SI{298.3}{\mega\hertz}$. Does convergence depend on $f$?
      \item \textbf{Optional.} Define the error at iteration $i$ to be $\delta
	u^{(i)} = u^{(i)} - u$. Show that
        \begin{equation}
          \delta u^{(i+1)} = \left(\nabla^2 + k_0^2\right)^{-1} \left(k_0^2 - k(\vx)^2\right) \delta u^{(i)}
        \end{equation}
        and derive a condition on $k(\vx)$ that guarantees convergence.
      \item \textbf{Optional.} Solve Eq.~\ref{eq:helmholtz_variable} using the
	iterative method with $N=1024$ and $f=\SI{21.3}{\mega\hertz}$ and
	record the time to convergence. Estimate the time it would take for
	method (i) to solve this problem. Use the provided data
	\texttt{mri\_data\_1024.dat}.
    \end{enumerate}
  \item \textbf{Testing the Schur complement for the Poisson equation.}
    Figure~\ref{fig:pss} shows an example of a \textit{perfect squared square},
    whereby a square of side length 112 is tessellated by twenty one smaller
    squares with unique integer side lengths. A text file \texttt{pss.txt} is
    provided that lists the square side lengths and $(x,y)$ offsets from the
    bottom left corner.

    In the lectures, we considered the model problem of solving the Poisson
    equation $-\nabla^2 u=f$ on the domain $[0,1]^2$ with zero Dirichlet
    boundary conditions. A C++ class to solve this problem using the FFT
    is provided in the course examples (\texttt{poisson\_fft.cc}).

    Suppose now that the domain is discretized into $112 \times 112$ grid.
    Specifically, the domain is divided into squares with side length
    $h=\nicefrac1{112}$, with gridpoints located on the corners of the squares.
    Consider the following two equivalent methods to solve this problem:
    \begin{itemize}
      \item Applying the FFT solver to the entire grid. Following the the
	notation of lecture 9, this grid corresponds to using $N=111$, since
	there is a matrix of $111\times 111$ unknown entries in the grid
	interior.
      \item Using the Schur complement method to solve the twenty-one
	non-overlapping square domains numbered $k=1,\ldots 21$ with grid sizes
	shown. In this case, the final domain $k=22$ corresponds to all
	gridpoints that lie on the internal square edges in Fig.~\ref{fig:pss}.
	In the AM225 code examples Git repository, the directory
	\texttt{2c\_fft+ddecomp/pss} contains a code \texttt{pss\_test.cc} that
	reads in the perfect squared square information, and prints ASCII art
	showing how the grid is subdivided. It also prints out a table of
	adjacencies for all gridpoints in the final domain with $k=22$.
    \end{itemize}
    Write codes to solve this problem using the two methods above, for the
    source term of $f(x,y)=e^{x-y}$. Show that your two codes give identical
    results up to numerical precision.
\end{enumerate}

\begin{figure}[b]
  \begin{center}
    \includegraphics[height=4in]{pss}
    \caption{A perfect squared square.\label{fig:pss}}
  \end{center}
\end{figure}

\end{document}
