\section{Introduction and Rel. Work}
\label{sec:intro}

For new HPC installations contribution of power usage to total system cost has
been increasing steadily over the past years
\cite{koomey2011growth} and studies project this trend to continue \cite{fraunhofer}.
As a consequence energy-aware metrics have recently been gaining popularity.
Energy-to-solution, i.e. the amount of energy
consumed by a system to solve a given problem, is the most obvious of
these metrics and will be used as quality measure throughout this study.

Previous works view the application code as constant and instead focus their
energy optimization attempts to parameter tuning on either the runtime
environment, the hardware, or both \cite{journals/tpds/LiSSNC13}. The runtime
environment approach works by adjusting the number of active threads across
different parallel regions based on the regions' computational requirements
\cite{4407685}. Hardware parameter tuning involves trying to identify slack,
e.g. during MPI communication, and using \ac{DVFS} to lower energy consumption
in such sections of low computational intensity
\cite{Miyoshi:2002:CPS:514191.514200}; similar strategies can be applied to
OpenMP barriers \cite{4725146}.  The parameters employed by these algorithms
can be based on measurements made at runtime
\cite{Rountree:2009:AMD:1542275.1542340} or set statically \cite{5348808}. One
fact often ignored by \ac{DVFS} control software, however, is that hardware
delays caused by changing frequency states can often be significant
\cite{Mazouz2013} which can lead to diminishing returns in the real world.
Another form of hardware parameter optimization involves the tuning of hardware
prefetchers \cite{Wu:2011:CDM:2015039.2015516}.

In contrast to previous work, our approach focuses on increasing single-core
performance, primarily though software optimization.  Using the \ac{ECM}
performance model to guide optimizations such as SIMD vectorization, cache
blocking, non-temporal stores, and the use of \ac{COD} mode the bandwidth
consumption of a single core is maximized. This allows the chip to saturate
main memory bandwidth with the least number of active cores in the multi-core
scenario, thus reducing power as well as the energy consumed by the chip.  In a
second step, different hardware parameters such as the number of active cores
and their frequencies should be evaluated to further reduce energy consumption.
The approach is demonstrated on a 2D Jacobi solver, which acts as a proxy for
many memory-bound applications.  To produce significant results, all
experiments were performed on the four most recent server microarchitectures by
Intel, which make up about 91\% of HPC installations of the June 2016 Top 500
list.

The paper is organized as follows. Section~\ref{sec:blueprint} describes the
blueprint of and the reasoning behind our energy-efficiency optimization
approach.  Section~\ref{sec:ecm} introduces the \ac{ECM} performance model that
we use to guide optimizations.  Section~\ref{sec:testbed} contains an overview
of the systems used for benchmarking.  Core-level improvement efforts are
documented in Section~\ref{sec:scopt}, followed by a validation of single-core
results in Section~\ref{sec:res:sc}. Implications for the multi-core scenario are
discussed in Section~\ref{sec:res:mc}, followed by the conclusion in Section~\ref{sec:conc}.

\section{Optimization Approach}
\label{sec:blueprint}

Per definition, the bottleneck for bandwidth-bound applications is the
sustained bandwidth $b_\mathrm{s}$ imposed by the memory subsystem.  It is a
well-established fact that a single core of modern multi- and many-core systems
cannot saturate main memory bandwidth \cite{Treibig:2009}, necessitating the
use of multiple cores to reach peak performance. The number of cores
$n_\mathrm{s}$ required to sustain main memory bandwidth depends on single-core
performance: The faster a single-core implementation, the more bandwidth is
consumed by a single core. The more bandwidth is consumed by a single core, the
fewer cores are required to saturate main memory bandwidth.

\begin{figure}[!tb]
    \centering
    \includegraphics[width=0.47\linewidth]{blueprint-bw-crop}\hfill
    \includegraphics[width=0.47\linewidth]{blueprint-joule-crop}
    \caption{\label{fig:blueprint-example}Bandwidth and energy-to-solution for different 2D Jacobi implementations using a dataset size of 16\,GB.}
\end{figure}

This relation is illustrated in Fig.~\ref{fig:blueprint-example}, which depicts
measurement results obtained on a single socket of a standard two-socket
Broadwell-EP machine (cf. Section~\ref{sec:testbed} for hardware details).
Fig.~\ref{fig:blueprint-example}a shows memory bandwidth as function of active
cores for two different 2D Jacobi implementations.
Fig.~\ref{fig:blueprint-example}b depicts the energy-to-solution required for
both implementations.  From this real-world example two important conclusion
can be drawn:\\

\noindent 1. \emph{Memory-bound codes should be run using $n_\mathrm{s}$
cores}\\
Once the bandwidth bottleneck is hit, adding more cores no longer increases
performance.\footnote{In practise performance will decrease slightly when
using more than $n_\mathrm{s}$ cores, as prefetchers and memory subsystem
have to deal with additional memory streams which makes them less efficient;
this effect can be observed in Fig.~\ref{fig:blueprint-example}a.} Instead
using more cores increases chip power consumption resulting in a higher
energy-to-solution.  This effect is visible in
Fig.~\ref{fig:blueprint-example}b.\\

\noindent 2. \emph{Early saturation can lead to lower energy-to-solution}\\
The optimized version saturates memory bandwidth with six cores compared to the
ten cores required by the naive counterpart. Six vs. ten cores being active
typically translates into a lower power draw for version using fewer
cores.\footnote{In theory it is possible that an implementations using fewer
cores to saturate memory bandwidth is less energy-efficient than one that uses
more.  Consider, for example, an optimization that doubles single-core
performance but triples the power drawn by the core. $n_\mathrm{s}$ is halved
but energy to solution is 50\% higher nonetheless. In practise we have however
never observed such a scenario. For all optimizations described in
Section~\ref{sec:scopt} the increase in single-core power draw accompanying an
optimization is negligible.} Together with the fact that versions
($16\,\mathrm{GB}/b_\mathrm{s}$), this results in a better energy-to-solution
for the optimized version.

After establishing these facts we propose the following approach to
optimize energy consumption:
(1) Use the \ac{ECM} model to guide performance improvements for the
            single-core implementation.
(2) Attempt to lower per-core power draw by tuning hardware parameters,
            e.g. core frequency or \ac{COD} mode.
(3) Never run code with more cores than necessary ($n_\mathrm{s}$).

\begin{table*}[!tb]
\centering
\caption{\label{tab:testbed}Test machine specifications.}
\resizebox{\textwidth}{!}
\begin{tabular}{|l|c|c|c|c|}                                                                                \hline
Microarchitecture (Shorthand)   & Sandy Bridge-EP (SNB)             & Ivy Bridge-EP (IVB)                   & Haswell-EP (HSW)                      & Broadwell-EP  (BDW) \\                    \hline
Chip Model                      & Xeon E5-2680                      & Xeon E5-2690 v2                       & Xeon E5-2695 v3                       & Pre-release \\                            \hline
Release Date                    & Q1/2012                           & Q3/2013                               & Q3/2014                               & Q1/2016 \\                                \hline
non-AVX/AVX Base Freq.          & 2.7\,GHz/2.7\,GHz                 & 3.0\,GHz/3.0\,GHz                     & 2.3\,GHz/1.9\,GHz                     & 2.1\,GHz/2.0\,GHz \\                      \hline
Cores/Threads                   & 8/16                              & 10/20                                 & 14/28                                 & 18/36 \\                                  \hline
Latest SIMD Extensions          & AVX                               & AVX                                   & AVX2, FMA3                            & AVX2, FMA \\                              \hline
Core-Private L1/L2 Caches       & 8$\times$32\,kB/8$\times$256\,kB  & 10$\times$32\,kB/10$\times$256\,kB    & 14$\times$32\,kB/14$\times$256\,kB    & 18$\times$32\,kB/18$\times$256\,kB \\     \hline
Shared Last-Level Cache         & 20\,MB (8$\times$2.5\,MB)         & 25\,MB (10$\times$2.5\,MB)            & 35\,MB (14$\times$2.5\,MB)            & 45\,MB (18$\times$2.5\,MB) \\             \hline
Memory Configuration            & 4 ch. DDR3-1600                   & 4 ch. DDR3-1866                       & 4 ch. DDR4-2133                       & 4 ch. DDR4-2133 \\                        \hline
Theoretical Mem. Bandwidth      & 51.2\,GB/s                        & 59.7\,GB/s                            & 68.2\,GB/s                            & 68.2\,GB/s \\                             \hline
Sustained Copy Bandwidth        & 39.5\,GB/s (77\%)                 & 43.4\,GB/s (74\%)                     & 50.1\,GB/s (73\%)                     & 51.5 (76\%) \\                            \hline
L1$\rightarrow$Reg Bandwidth    & 2$\times$16\,B/cy                 & 2$\times$16\,B/cy                     & 2$\times$32\,B/cy                     & 2$\times$32\,B/cy \\                      \hline
Reg$\rightarrow$L1 Bandwidth    & 1$\times$16\,B/cy                 & 1$\times$16\,B/cy                     & 1$\times$32\,B/cy                     & 1$\times$32\,B/cy \\                      \hline
L1$\leftrightarrow$L2 Bandwidth & 32\,B/cy${\hat=}$2\,cy/CL         & 32\,B/cy${\hat=}$2\,cy/CL             & 64\,B/cy${\hat=}$1\,cy/CL             & 64\,B/cy${\hat=}$1\,cy/CL     \\          \hline
L2$\leftrightarrow$L3 Bandwidth & 32\,B/cy${\hat=}$2\,cy/CL         & 32\,B/cy${\hat=}$2\,cy/CL             & 32\,B/cy${\hat=}$2\,cy/CL             & 32\,B/cy${\hat=}$2\,cy/CL     \\          \hline
L3$\leftrightarrow$Mem Bandwidth (copy) & 14.5\,B/cy${\hat=}$4.4\,cy/CL & 20.0\,B/cy${\hat=}$4.4\,cy/CL     & 21.7\,B/cy${\hat=}$2.9\,cy/CL         & 24.6\,B/cy${\hat=}$2.6\,cy/CL \\          \hline
\end{tabular}}
\end{table*}

\section{The ECM Performance Model}
\label{sec:ecm}

The \ac{ECM} model \cite{Treibig:2009,hager:cpe,sthw15, Hofmann2016:1,
Hofmann2016:2} is an analytic performance model that, with the exception of
sustained memory bandwidth, works exclusively with architecture specifications
as inputs. The model estimates the numbers of CPU cycles required to execute a
number of iterations of a loop on a single core of a multi- or many-core chip.
For multi-core estimates, linear scaling of single-core performance is assumed
until a shared bottleneck, such as e.g.  main memory bandwidth, is hit.  Note
that only parts of the model relevant to this work, i.e. the single-core model,
are presented here.  Readers interested in the full ECM model can find the most
recent version for Intel Xeon, Intel Xeon Phi, and IBM POWER8 processors here
\cite{CPE:CPE3921}.

The single-core prediction is made up of contributions from the in-core
execution time $T_\mathrm{core}$, i.e. the time spent executing instructions in
the core under the assumption that all data resides in the L1 cache, and the
transfer time $T_\mathrm{data}$, i.e. the time spent transferring data from its
location in the cache/memory hierarchy to the L1 cache.  As data transfers in
the cache- and memory hierarchy occur at \ac{CL} granularity we chose the
number of loop iterations $n_\mathrm{it}$ to correspond to one cache line's
"worth of work." On Intel architectures, where \ac{CL}s are 64\,B in size,
$n_\mathrm{it}=8$ when using \ac{DP} floating-point numbers, because processing
eight ``doubles'' (8\,B each) corresponds to exactly one \ac{CL} worth of work.

Superscalar core designs house multiple execution units, each dedicated to
perform certain work: loading, storing, multiplying, adding, etc.  The in-core
execution time $T_\mathrm{core}$ is determined by the unit that takes the
longest to retire the instructions allocated to it. Other constraints for the
in-core execution time may apply, e.g. the four micro-op per cycle retirement
limit of Intel Xeon cores. The model differentiates between core cycles
depending on whether data transfers in the cache hierarchy can overlap with
in-core execution time. For instance, on Intel Xeons, core cycles in which data
is moved between the L1 cache and registers, e.g. cycles in which load and/or
store instructions are retired, prohibit simultaneous transfer of data between
the L1 and L2 cache; these ``non-overlapping'' cycles contribute to
$T_\mathrm{nOL}$. Cycles in which no load or store instructions but other
instructions, such as e.g.  arithmetic instructions, retire are considered
``overlapping'' cycles and contribute to $T_\mathrm{OL}$. The in-core runtime
is the maximum of both: $T_\mathrm{core} = \max(T_\mathrm{OL},
T_\mathrm{nOL})$.

For modelling data transfers, latency effects are initially neglected, so
transfer times are exclusively a function of bandwidth. Cache bandwidths are
typically well documented and can be found in vendor data sheets. Depending on
how many \ac{CL}s have to be transferred, the contribution of each level in the
memory hierarchy ($T_\mathrm{L1L2}$, \ldots, $T_\mathrm{L3Mem}$) can be
determined.  Special care has to be taken when dealing with main memory
bandwidth, because theoretical memory bandwidth specified in the data sheet and
sustained memory bandwidth $b_\mathrm{s}$ can differ greatly. Also, in
practise $b_\mathrm{s}$ depends on the number of load and store
streams. It is therefore recommended to empirically determine $b_\mathrm{s}$
using a kernel that resembles the memory access pattern of the benchmark to be
modeled.  Once determined, the time to transfer one
\ac{CL} between the L3 cache and main memory can be derived from the CPU
frequency $f$ as $64\,\mathrm{B} \cdot f / b_\mathrm{s}$~cycles.

Starting with the Haswell-EP (HSW) microarchitecture, an empirically determined
latency penalty $T_p$ is applied to off-core transfer times.  This departure
from the bandwidth-only model has been made necessary by large core counts, the
dual-ring design, and separate clock frequencies for core(s) and Uncore all of
which increase latencies when accessing off-core data.  The penalty is added
each time the Uncore interconnect is involved in data transfers. This is the
case whenever data is transferred between the L2 and L3 caches, as data is
pseudo-randomly distributed between all last-level cache segments; and when
data is transferred between the L3 cache and memory, because the memory
controller is attached to the Uncore interconnect.  Instruction times as well
as data transfer times, e.g.  $T_\mathrm{L1L2}$ for the time required to
transfer data between L1 and L2 caches, can be summarized in shorthand
notation:
\ecmshort{T_\mathrm{OL}}{T_\mathrm{nOL}}{T_\mathrm{L1L2}}{T_\mathrm{L2L3}+T_\mathrm{p}}{T_\mathrm{L3Mem}+T_\mathrm{p}}.

To arrive at a prediction, in-core execution and data transfers times are put
together. Depending on whether there exist enough overlapping cycles to hide
all data transfers, runtime is given by either $T_\mathrm{OL}$ or the sum of
non-overlapping core cycles $T_\mathrm{nOL}$ plus contributions of data
transfers $T_\mathrm{data}$, whichever takes longer. $T_\mathrm{data}$ consists
of all necessary data transfers in the cache/memory hierarchy, plus latency
penalties if applicable, e.g. for data coming from the L3 cache:
$T_\mathrm{data}=T_\mathrm{L1L2}+T_\mathrm{L2L3}+T_\mathrm{p}$.  The prediction
is thus $T_\mathrm{ECM} = \max(T_\mathrm{OL}, T_\mathrm{nOL} +
T_\mathrm{data})$. A shorthand notation also exists for the model's prediction:
\ecmpshort{T_\mathrm{ECM}^\mathrm{core}}{T_\mathrm{ECM}^\mathrm{L2}}{T_\mathrm{ECM}^\mathrm{L3}}{T_\mathrm{ECM}^\mathrm{Mem}}{}.

Converting the prediction from time (in cycles) to performance (work per
second) is done by dividing the work per CL $W_\mathrm{CL}$ (e.g.
floating-point operations, updates, or any other relevant work metric) by the
predicted runtime in cycles and multiplying with the processor frequency $f$,
i.e.  $P_\mathrm{ECM} = W_\mathrm{CL}/T_\mathrm{ECM} \cdot f$.

\section{Experimental Testbed}
\label{sec:testbed}

All measurements were performed on standard two-socket Intel Xeon servers.  A
summary of key specifications of the four generations of processors can be
found in Table~\ref{tab:testbed}.  According to Intel's ``tick-tock'' model,
where a ``tick'' corresponds to a shrink of the manufacturing process
technology and a ``tock'' to a new microarchitecture, IVB and BDW are
``ticks''---apart from an increase in core count and a faster memory clock, no
major improvements were introduced in these microarchitectures.

HSW, which is a ``tock'', introduced AVX2, extending the already existing
256\,bit SIMD vectorization from floating-point to integer data types.
Instructions introduced by the \ac{FMA} extension are handled by two new,
AVX-capable execution units. Data paths between the L1 cache and registers as
well as the L1 and L2 caches were doubled. Due to limited scalability of a
single ring connecting the cores, HSW chips with more than eight feature a
dual-ring design.  HSW also introduces the AVX base and maximum AVX Turbo
frequencies.  The former is the minimum guaranteed frequency when running AVX
code on all cores; the latter the maximum frequency when running AVX code on
all cores (cf.  Table~3 in \cite{Intel-v3-specupdate}).  Based on workload, the
actual frequency varies between this minimum and maximum value. For a more
detailed analysis of the differences between the SNB/IVB and HSW
microarchitectures see \cite{Hofmann2016:1}.

For all measurements on SNB and IVB, the CPU frequency was fixed at the nominal
CPU frequency.  On HSW and BDW, the CPU frequency was fixed to the non-AVX base
frequency;\footnote{The guaranteed frequency when running code on all cores
that does not use AVX instructions.} in none of the single-core measurements
that feature AVX code could we observe a drop below the non-AVX base frequency.
The measured sustained main memory bandwidth displayed in
Table~\ref{tab:testbed} is that of the STREAM copy kernel \cite{McCalpin1995}
using non-temporal stores, because its memory access pattern corresponds to
that of the 2D Jacobi solver.  Energy consumption was determined by accessing
the \ac{RAPL} interface through \texttt{likwid-perfctr} \cite{Treibig:2012}.
Because node and memory power usage fluctuates depending on node configuration
and RAM manufacturer, we chose to only present results pertaining to chip power
consumption, i.e. cores, core-private caches, and Uncore (cf.~14.9 in
\cite{intel-sdm-2016}).  All code is compiled using the Intel C Compiler
version 15.0.2.

\section{Single-Core Optimizations}
\label{sec:scopt}

\begin{figure}[tb]
\lstset{

        breaklines=true,
        language=C,
        basicstyle=\small\ttfamily,
        numbers=left,
        numberstyle=\tiny,
        frame=t,
        columns=fullflexible,
        showstringspaces=false,
        escapechar=\%,
        escapebegin=\color{blue}\ttfamily\bfseries,
        keepspaces=true
}
\begin{lstlisting}
for (y=1; y<Y-1; ++y)
  for (x=1; x<X-1; ++x)
    b[y][x]=0.25 * (a[y-1][x] + a[y][x-1] +
                    a[y][x+1] + a[y+1][x]);
\end{lstlisting}
\caption{\label{fig:code:baseline} C implementation for one 2D Jacobi iteration.}
\end{figure}

Figure~\ref{fig:code:baseline} shows the source code for one 2D five-point
Jacobi sweep, i.e. the complete update of all grid points.  One grid point
update computes and stores in \texttt{b} the new state of each point from the
values of its four neighbors in \texttt{a}, which holds data from the previous
iteration.  For results to be representative, the dataset size per socket for
all measurements is 16\,GB, i.e. each of the two grids is 8\,GB in size and
made up of $32768\times32768$ double-precision numbers.  Performance is
expressed in ``lattice updates per second'' (LUP/s), i.e. scalar inner kernel
iterations per second.

\subsection{Baseline AVX Implementation}
\label{sec:scopt:baseline}

Using adequate optimization flags (\texttt{-O3 -xHost -fno-alias}) it is
trivial to generate AVX vectorized assembly for the code shown in
Figure~\ref{fig:code:baseline} with recent Intel compilers. This is why we
decided to use an AVX vectorized variant instead of scalar code as baseline.

With 256-bit AVX vectorization in place, one \ac{CL} worth of work (eight
LUP/s), consists of eight AVX loads, two AVX stores, six AVX adds, and two AVX
multiplication instructions.

To determine data transfers inside the cache hierarchy, we have to examine each
load and store in more detail. Storing the newly computed results to array
\texttt{b} involves the transfer of two \ac{CL}s to/from main memory: because
both arrays are too large to fit inside the caches, the store will miss in the
L1 cache, triggering a write-allocate of the \ac{CL} from main memory. After
the values in the \ac{CL} have been updated, the \ac{CL} will have to be
evicted from the caches eventually, triggering another main memory transfer.

The left neighbor \texttt{a[y][x-1]} can always be loaded from the L1 cache
since it was used two inner iterations before as right neighbor; access to
\texttt{a[y+1][x]} must be loaded from main memory since it was not used before
within the sweep.

Based on work by Rivera and Tseng \cite{Rivera:2000} Stengel et al. introduced
the layer condition (LC) \cite{sthw15} to help determining where data for
\texttt{a[y-1][x]} and \texttt{a[y][x+1]} is coming from.  The LC stipulates
three successive rows have to fit into a certain cache for accesses to
these data to come from this particular cache.
Assuming cache $k$ can effectively hold data up to 50\% of its nominal size
$C_k$, the LC can be formulated as
$3 \cdot N \cdot 8\,\mathrm{B} < 50\% \cdot C_k$.
For $N=32768$, three rows take up 768\,kB so on all previously introduced
machines (cf. line~8 in Table~\ref{tab:testbed}) the LC holds true for the L3
cache.\\

\noindent\textbf{ECM Model for SNB and IVB}
\label{sec:scopt:baseline:SNBIVB}

To process one CL worth of data, eight AVX load, two AVX store, six AVX
addition and two AVX multiplication instructions have to executed.  Throughput
is limited by the two load units. Each load unit has a 16\,B wide data path
connecting registers and L1 cache, so retiring a 32\,B AVX load takes two
cycles. Using both load units, eight AVX loads take
$T_\mathrm{nOL}=8\,\mathrm{cy}$. Both AVX stores are retired in parallel with the
eight loads, so they do not increase $T_\mathrm{nOL}$.  Computation throughput
is limited by the single add port, which takes $T_\mathrm{OL}=6\,\mathrm{cy}$
to retire all six AVX add instructions. Both AVX multiplications can be
processed in parallel with two of the six AVX add instructions.

As established previously, three CLs have to be transferred between L3 and
memory: write-allocating and later evicting \texttt{b[y][x]} and loading
\texttt{a[y+1][x]}. On both SNB and IVB, L3-memory bandwidth is 4.4\,cy/CL (cf.
last line in Table~\ref{tab:testbed}). This results in
$T_\mathrm{L3Mem}=13.2\,\mathrm{cy}$ for both architectures.  The same three
CLs have to be transferred between L3 and L2 cache; in addition, CLs for
\texttt{a[y-1][x]} and \texttt{a[y][x+1]} have to be transferred from the L3 to
the L2 cache. Transferring five CLs at a bandwidth of 2\,cy/CL takes
$T_\mathrm{L2L3}=10\,\mathrm{cy}$ on both SNB and IVB. At a L1-L2 bandwidth of
2\,cy/CL, moving these five CLs between L2 and L1 cache takes
$T_\mathrm{L1L2}=10\,\mathrm{cy}$.  Using the ECM short notation to summarize
the inputs yields \ecm{6}{8}{10}{10}{13.2}{\mathrm{cy}} for both SNB and IVB;
the corresponding runtime prediction for SNB and IVB is
\ecmp{8}{18}{28}{41.2}{\mathrm{cy}}.

For a 2.7\,GHz SNB core the performance prediction is $P_\mathrm{ECM}^{Mem}=
\frac{8\,\mathrm{LUP/CL}}{41.2\,\mathrm{cy/CL}} \cdot 2.7\,\mathrm{GHz} =
524\,\mathrm{MLUP/s.}$ For IVB the model
predicts a performance of $582\,\mathrm{MLUP/s}$.\\

\noindent\textbf{ECM Model for HSW and BDW}
\label{sec:scopt:baseline:HSWBDW}

On HSW and BDW, the \acp{AGU} are the bottleneck for $T_\mathrm{nOL}$.  Each
load/store instruction accesses an \ac{AGU} to compute the referenced memory
address. With only two \acp{AGU} capable of performing the required addressing
operations available, retiring all ten load/store instructions takes
$T_\mathrm{nOL}=5\,\mathrm{cy}$.  HSW and BDW posses a single AVX add unit just
like SNB and IVB, so $T_\mathrm{OL}=6\,\mathrm{cy}$.

The off-core latency penalty was empirically estimated at approximately
1.6~cycles for both HSW and BDW and is applied per CL transfer that takes place
over the Uncore interconnect, i.e. all transfers between L3 and L2 caches as
well as transfers between memory and the L3 cache. The effective L3-Mem
bandwidth is thus 2.9+1.6\,cy on HSW and 2.6+1.6\,cy on BDW; the effective
L2-L3 is 2+1.6\,cy on HSW and 2+1.6\,cy on BDW.

Transferring the three required CLs then results in
$T_\mathrm{L3Mem}=8.7+4.8\,\mathrm{cy}$ on HSW resp.
$T_\mathrm{L3Mem}=7.8+4.8\,\mathrm{cy}$ on BDW. Moving five CLs between the L2
and L3 caches takes $T_\mathrm{L2L3}=10+8\,\mathrm{cy}$ on both HSW and BDW. At
a L1-L2 bandwidth of 1\,cy/CL, moving the same five CLs takes
$T_\mathrm{L1L2}=5\,\mathrm{cy}$ on both microarchitectures.  The ECM inputs
are thus \ecm{6}{5}{5}{10+8}{8.7+4.8}{\mathrm{cy}} for HSW and for BDW
\ecm{6}{5}{5}{10+8}{7.8+4.8}{\mathrm{cy}}. The corresponding runtime
predictions are \ecmp{6}{10}{28}{41.5}{\mathrm{cy}} for HSW and
\ecmp{6}{10}{28}{40.6}{\mathrm{cy}} for BDW. The performance predictions are
443\,MLUP/s for HSW and 413\,MLUP/s for BDW.

\subsection{Cache Blocking Optimization}
\label{sec:scopt:cacheblocking}

One way to increase the performance of the single-core implementation is to
reduce the amount of time spent transferring data inside the cache hierarchy.
As previously established, \texttt{a[y-1][x]} and \texttt{a[y][x+1]} are
loaded from the L3 cache, because the L1 and L2 caches are too small to
fulfill the LC for $N=32768$. Using cache blocking, it is possible to
enforce the LC in arbitrary cache levels.  This is done by partitioning the
grid into stripes along the $y$-axis; the grid is then processed stripe by
stripe.  The diameter of the stripes, also known as blocking factor $b_x$, is
chosen is such a way that the LC is met for a given cache level.
If L2 blocking is desired, the L2 cache size of 256\,kB requires that $b_x$
should be chosen smaller than 5461.

Efficient blocking for the 32\,kB L1 cache not as straightforward. Although
determining $b_x<682$ is simple using the LC, naive L1 blocking
in $x$-direction has negative side effects. With most of the data for one CL
update coming from L1, the L2 cache is less busy. This slack is detected by
hardware prefetchers, making them more aggressive, leading to data being
prefetched from main memory. With $b_x \approx 680$, the size of one stripe is
$680 \cdot 32768 \cdot 8\,\mathrm{B}=170\,\mathrm{MB}$---too large for the L3
cache, which means that data prefetched from main memory will be preempted from
the cache before it is used. This can be avoided either by disabling some of
the prefetchers or additionally blocking in $y$-direction. We used the latter,
because disabling prefetchers might degrade performance elsewhere.  To
guarantee data is used before being preempted, the size of each
chunk should be chosen smaller than 50\% of a single L3 segment,\footnote{In
the multi-core scenario, all cores might be active and store data in the shared
L3 cache; the capacity dedicated to a single core is thus that of a single L3
segment adjusted by 50\% to reflect effective cache size.} e.g. $b_y < 50\%
\cdot 2.5\,\mathrm{MB} / ( b_x \cdot 8\,\mathrm{B})$.\\

\noindent\textbf{ECM Model for SNB and IVB}
\label{sec:scopt:cacheblocking:SNBIVB}

Other than causing negligible loop overhead cache blocking does not change the
instructions that have to be retired to process one CL; thus $T_\mathrm{OL}$
and $T_\mathrm{nOL}$ remain unchanged.

L2 blocking reduces the number of CLs transferred between L3 and L2 from five
to three, lowering $T_\mathrm{L2L3}$ from ten to six cycles on both SNB and
IVB. The runtime prediction $T_\mathrm{ECM}^\mathrm{Mem}$ for both SNB
and IVB is reduced from 41.2 to 37.2\,cy, leading to a performance prediction
$P_\mathrm{ECM}^\mathrm{Mem}$ of 580\,MLUP/s for SNB and 645\,MLUP/s for IVB.

Similarly L1 blocking reduces the CL traffic between L2 and L1 caches from five
to three CLs, lowering $T_\mathrm{L1L2}$ from ten to six cycles on both SNB and
IVB. Again, the runtime prediction $T_\mathrm{ECM}^\mathrm{Mem}$ for both
microarchitectures is reduced by four cycles from 37.2 to 33.2\,cy.  The
predicted performance $P_\mathrm{ECM}^\mathrm{Mem}$ increases to 651\,MLUP/s on
SNB and 723\,MLUP/s on IVB.\\

\noindent\textbf{ECM Model for HSW and BDW}
\label{sec:scopt:cacheblocking:HSWBDW}

The effect of L2 blocking is more pronounced on HSW and BDW, because the cost
of transferring a CL is higher on these microarchitectures due to the latency
penalty. L2 blocking lowers $T_\mathrm{L2L3}$ from 10+8 to 6+4.8\,cy on HSW and
BDW. In turn, the runtime prediction $T_\mathrm{ECM}^\mathrm{Mem}$ for HSW is
reduced from 41.5 to 34.3\,cy and from 40.6 to 33.4\,cy on BDW.  The
performance prediction $P_\mathrm{ECM}^\mathrm{Mem}$ increases to 536\,MLUP/s
on HSW and 503\,MLUP/s on BDW.

Because the L1-L2 bandwidth increased from 2\,cy/CL on SNB/IVB to 1\,cy/CL on
HSW/BDW, the performance improvement offered by L1 blocking is less pronounced
than on SNB and IVB. With the number of CL transfers lowered from five to three
with L1 blocking, $T_\mathrm{L1L2}$ is reduced from five to three cycles.  The
runtime prediction $T_\mathrm{ECM}^\mathrm{Mem}$ for HSW is reduced from 34.3
to 32.3\,cy; on BDW from 33.4 to 31.4\,cy.  The predicted performance
$P_\mathrm{ECM}^\mathrm{Mem}$ is increased to 570\,MLUP/s on HSW and
535\,MLUP/s on BDW.

\subsection{Cluster-on-Die Mode}
\label{sec:scopt:cod}

As a workaround to the limited scalability of the physical ring interconnect
introduced with Westmere-EX, HSW and BDW switch to a dual-ring design. HSW uses
the so-called ``eight plus $x$'' design, in which the first eight cores of a
chip are attached to a primary ring and the remaining cores (six for the
model introduced in Section~\ref{sec:testbed}) are attached to a secondary
ring; BDW uses a symmetric design. Two queues enable data to pass between
rings. In the default (non-\ac{COD}) mode, the physical topology is hidden from the
operating system, i.e. all cores are exposed within the same \ac{NUMA} domain.

To understand the latency problems caused by the interconnect, we examine the
route data travels inside the Uncore.  Using a hashing function data is is
distributed across all L3 segments based on its memory address.  When accessing
data in the L3 cache, there is a high probability it must be fetched from
remote L3 segments.  In the worst case, this is a segment on the other ring so
there might be a high latency involved. In the case of a L3 miss the situation
gets worse. Each physical ring has attached to it a \ac{MC} and the
choice which MC to use is again based on the data's address. So a L3 miss in a
segment on one physical ring does not imply that this ring's \ac{MC} will be
used to fetch the data from memory.  That leads to cases in which a large
number of hops and multiple cross-physical ring transfers are involved when
getting data from main memory.

One way to reduce these latencies is the new \ac{COD} mode introduced together
with the dual-ring design on HSW and BDW in which cores are separated into
two physical clusters of equal size.
The latency reduction is achieved by adapting the involved hashing functions.
Data requested by a core of a cluster will only be placed in the cluster's L3
segments; in addition, all memory transfers are routed to the \ac{MC} dedicated
to the cluster.  Thus, for \ac{NUMA}-aware codes, \ac{COD} mode effectively
lowers the latency by reducing the diameter and the mean distance of the
dual-ring interconnect by restricting each cluster to its own physical ring.
For a more detailed analysis of \ac{COD} mode see \cite{Hofmann2016:1}. On the
HSW chip used for benchmarks, \ac{COD} mode lowers the interconnect latency by
0.5\,cy; on the employed BDW chip, where a single ring still has eleven cores
attached to it, the latency is only reduced by 0.3\,cy.\\

\noindent\textbf{ECM Model for HSW and BDW}
\label{sec:scopt:cod:HSWBDW}

With \ac{COD} enabled, the per-CL Uncore latency penalty is reduced to 1.1\,cy on HSW resp.
1.3\,cy on BDW.  This leads to $T_\mathrm{L3Mem}=8.7+3.3$\,cy on HSW resp.
7.8+3.9\,cy on BDW.  Because the Uncore is also involved in L2-L3 transfers,
the latency improvement caused by enabling \ac{COD} also positively influences
$T_\mathrm{L2L3}$. To transfer three CLs, 6+3.3\,cy are required on HSW resp.
6+3.9\,cy on BDW.  The resulting ECM inputs are \ecmshort{6}{5}{3}{6+3.3}{8.7+3.3}
for HSW and \ecmshort{6}{5}{3}{6+3.9}{7.8+3.9} for BDW.  The corresponding
runtime prediction is \ecmpshort{6}{8}{17.3}{29.3} on HSW and
\ecmpshort{6}{8}{17.9}{29.6} on BDW. The performance predicted by the ECM model
is 628\,MLUP/s for HSW and 567\,MLUP/s on BDW.

\subsection{Non-Temporal Stores}
\label{sec:scopt:ntstores}

Streaming or \ac{NT} stores are special instructions that
avoid write-allocates on modern Intel microarchitectures.
Without NT stores, storing the newly computed result \texttt{b[y][x]} triggers
a write-allocate. The old data is brought in from memory and travels through
the whole cache hierarchy.  Thus the first benefit of using NT stores is that
the unnecessary transfer of \texttt{b[y][x]} from memory to the L1 cache no
longer takes place. In addition, NT stores will also strip some cycles off the
time involved getting the new result to main memory. Using regular stores, the
newly computed result is written to the L1 cache, from where the data has to be
evicted at some point through the whole cache hierarchy into memory.  Using NT
stores, \acp{CL} are sent via the L1 cache to the \acp{LFB}; from there, they
are transfered directly to memory and do not pass through the L2 and L3 caches.
Although the benefits should apply equally to all microarchitectures, there are
shortcomings in SNB and IVB that make single-core implementations using NT
stores slower than their regular stores counterpart. The positive impact of NT
stores can only be leveraged in multi-core scenarios on these
microarchitectures, which is why we chose to omit ECM models and measurements
for the NT store implementation for SNB and IVB.\\

\begin{table*}[!tb]
\centering
\caption{\label{tab:single-core-results}Summary of ECM inputs, ECM predictions, as well as measured performance, power consumption, and energy-to-solution for one 2D Jacobi iteration with a dataset size of 16\,GB.}
\resizebox{\textwidth}{!}
\begin{tabular}{|l|l|c|c|c|c|c|c|}                                                                                \hline
                &               & ECM                                                           & ECM pre-                                                  & $P_\mathrm{ECM}^\mathrm{Mem}$ & Measured  & Chip      & Chip Energy-\\
\textmu{}arch   & Version       & input [cy]                                                    & diction [cy]                                              & [MLUP/s]      & [MLUP/s]  & Power [W] & to-Solution [J] \\                \hline
SNB             & Baseline      & \ecmshort{6}{8}{10}{10}{13.2}                                 & \ecmpshort{8}{18}{28}{41.2}                               & 524           & 514       & 35.9      & 75.0 \\
                & L2 blocked    & \ecmshort{6}{8}{10}{\textbf{6}}{13.2}                         & \ecmpshort{8}{18}{\textbf{24}}{\textbf{37.2}}             & 580           & 623       & 36.4      & 62.7 \\
                & L1 blocked    & \ecmshort{6}{8}{\textbf{6}}{6}{13.2}                          & \ecmpshort{8}{\textbf{14}}{\textbf{20}}{\textbf{33.2}}    & 651           & 672       & 36.6      & 58.4 \\                           \hline
IVB             & Baseline      & \ecmshort{6}{8}{10}{10}{13.2}                                 & \ecmpshort{8}{18}{28}{41.2}                               & 582           & 539       & 32.7      & 65.3 \\
                & L2 blocked    & \ecmshort{6}{8}{10}{\textbf{6}}{13.2}                         & \ecmpshort{8}{18}{\textbf{24}}{\textbf{37.2}}             & 645           & 651       & 34.1      & 56.1 \\
                & L1 blocked    & \ecmshort{6}{8}{\textbf{6}}{6}{13.2}                          & \ecmpshort{8}{\textbf{14}}{\textbf{20}}{\textbf{33.2}}    & 722           & 714       & 33.0      & 49.6 \\                           \hline
HSW             & Baseline      & \ecmshort{6}{5}{5}{10+8}{8.7+4.8}                             & \ecmpshort{6}{10}{28}{41.5}                               & 443           & 435       & 50.9      & 125.6\\
                & L2 blocked    & \ecmshort{6}{5}{5}{\textbf{6+4.8}}{8.7+4.8}                   & \ecmpshort{6}{10}{\textbf{20.8}}{\textbf{34.3}}           & 536           & 529       & 51.0      & 103.5\\
                & L1 blocked    & \ecmshort{6}{5}{\textbf{3}}{6+4.8}{8.7+4.8}                   & \ecmpshort{6}{\textbf{8}}{\textbf{18.8}}{\textbf{32.3}}   & 570           & 579       & 52.1      & 96.6 \\
                & L1 b.+CoD     & \ecmshort{6}{5}{3}{6+\textbf{3.3}}{8.7+\textbf{3.3}}          & \ecmpshort{6}{8}{\textbf{17.3}}{\textbf{29.3}}            & 628           & 625       & 51.3      & 88.0 \\
                & L1 b.+CoD+nt  & \ecmshort{6}{5}{\textbf{2}}{\textbf{2+1.1}}{\textbf{5.8+2.2}} & \ecmpshort{6}{\textbf{7}}{\textbf{10.1}}{\textbf{18.1}}   & 1016          & 951       & 51.0      & 57.6 \\                           \hline
BDW             & Baseline      & \ecmshort{6}{5}{5}{10+8}{7.8+4.8}                             & \ecmpshort{6}{10}{28}{40.6}                               & 413           & 407       & 44.2      & 116.5\\
                & L2 blocked    & \ecmshort{6}{5}{5}{\textbf{6+4.8}}{7.8+4.8}                   & \ecmpshort{6}{10}{\textbf{20.8}}{\textbf{33.4}}           & 503           & 489       & 44.3      & 97.4 \\
                & L1 blocked    & \ecmshort{6}{5}{\textbf{3}}{6+4.8}{7.8+4.8}                   & \ecmpshort{6}{\textbf{8}}{\textbf{18.8}}{\textbf{31.4}}   & 535           & 509       & 44.5      & 93.7 \\
                & L1 b.+CoD     & \ecmshort{6}{5}{3}{6+\textbf{3.9}}{7.8+\textbf{3.9}}          & \ecmpshort{6}{8}{\textbf{17.9}}{\textbf{29.6}}            & 567           & 561       & 44.6      & 86.0 \\
                & L1 b.+CoD+nt  & \ecmshort{6}{5}{\textbf{2}}{\textbf{2+1.3}}{\textbf{5.2+2.6}} & \ecmpshort{6}{\textbf{7}}{\textbf{10.3}}{\textbf{18.1}}   & 928           & 862       & 45.4      & 56.6 \\                           \hline
\end{tabular}}
\end{table*}

\noindent\textbf{ECM Model for HSW and BDW}
\label{sec:scopt:ntstores:HSWBDW}

Transferring \texttt{a[y+1][x]} between the L1 and L2 caches takes 1\,cy.
Although the transfer is not strictly between the L1 and L2 cache, the cycle
spent transferring the CL for \texttt{b[y][x]} from the L1 cache to the
\acf{LFB} is booked in $T_\mathrm{L1L2}$ as well, making for a total L1-L2
transfer time of 2\,cy.  Loading \texttt{a[y+1][x]} from the L3 to the L2 cache
takes $T_\mathrm{L2L3}=2+1.1$\,cy on HSW and 2+1.3\,cy on BDW. Getting the CL
containing \texttt{a[y+1][x]} from memory takes 2.9+1.1\,cy on HSW and
2.6+1.3\,cy on BDW; again, although the transfer of \texttt{b[y][x]} is
strictly not between the L3 cache and memory, the transfer time to send the CL
from the \ac{LFB} to memory is booked in $T_\mathrm{L3Mem}$, making for a total
L3-Mem transfer time of 5.8+2.2\,cy on HSW and 5.2+2.6\,cy on BDW.  In summary,
the full ECM inputs are \ecmshort{6}{5}{2}{2+1.1}{5.8+2.2}\,cy on HSW and
\ecmshort{6}{5}{2}{2+1.3}{5.2+2.6}\,cy on BDW; the corresponding ECM runtime
prediction is \ecmpshort{6}{7}{10.1}{18.1} on HSW and
\ecmpshort{6}{7}{10.3}{18.1}.  The in-memory performance prediction by the ECM
model is $P_\mathrm{ECM}^\mathrm{Mem}=1016\,\mathrm{MLUP/s}$ on HSW and
928\,MLUP/s on BDW.

\section{Single-Core Results}
\label{sec:res:sc}

Table~\ref{tab:single-core-results} contains a summary of the ECM inputs and
predictions discussed in Section~\ref{sec:scopt}, as well as measurements of
performance, power, and energy consumption for one 2D Jacobi iteration using a
16\,GB dataset.  The model correctly predicts performance with a mean error of
3\% and a maximum error of 7\%, which indicates that all single-core
performance engineering measures work as intended. On SNB and IVB performance
increases of around 1.3$\times$ are achieved; on HSW with 2.2$\times$ resp. BDW
with 2.1$\times$ improvement the increases are even more pronounced.
Measurements obtained via the \ac{RAPL} interface indicate increases in power
consumption due to optimizations are negligible (in the range of 2\%). With
power draw almost constant, this means that the performance gains directly
translate into energy improvements.

An interesting observation regarding single-core power consumption\footnote{The
\ac{RAPL} counters can not report power consumption of individual cores but
only that of the whole package. Thus reported values also include power drawn
by Uncore facilities, e.g. all L3 segments and the interconnect.} can be made
when comparing the different microarchitectures that can be explained using
Intel's ``tick-tock'' model.  Power decreases with ``ticks,'' i.e. a shrink in
manufacturing size and the accompanying decreases in dynamic power; power
increases with ``tocks,'' i.e.  major improvements in microarchitecture: The
``tick'' from SNB which uses 32\,nm to IVB which uses 22\,nm technology
corresponds to a 10\% decrease in power consumption. HSW, the only ``tock'' in
Table~\ref{tab:single-core-results} uses the same 22\,nm process as IVB but
introduced major improvements in the microarchitecture (cf.
Section~\ref{sec:testbed}) which lead to a more than 50\% higher power draw.
BDW, a ``tick'' which uses a
14\,nm process uses 11\% less power than HSW.  While it is tempting to
generalize from these results, note that the reported numbers are specific to the
2D Jacobi application and the used chip models. We also believe that while a
surge in power consumption occurs with the HSW ``tock'' the cause to be a
combination of microarchitectural improvements and the higher core
count.\footnote{Despite only one core being active in the measurements, all L3
cache segments are active and draw power; a chip with more cores will thus draw
more power in single-core use cases.}

\section{Chip-Level Observations}
\label{sec:res:mc}

Although the \ac{ECM} model can be used to predict multi-core scaling behavior
\cite{CPE:CPE3921,Hofmann2016:1,sthw15}, due to space constraints chip-level discussion is restricted
to empirical results.
The graphs in Figure~\ref{fig:chip-perf} show the improvements
discussed in Section~\ref{sec:scopt} and relate measured performance and
energy-to-solution for different core counts.\footnote{The
leftmost measuring point of each graph corresponds to one core. Following the
line attached to a point to the next corresponds to one more core being active.
For demonstration purposes the purple graph in Fig.~\ref{fig:chip-perf}a has
some core counts annotated.} IVB (Fig.~{\ref{fig:chip-perf}a) and
HSW (Fig.~\ref{fig:chip-perf}b) are chosen as representatives to demonstrate
that the different optimizations have different impacts on performance and
energy consumption depending on the microarchitecture.

On IVB, running on $n_\mathrm{s} = 5$ instead of ten cores reduces energy
consumption by 13\% 44.4\,J to 34.8\,J; on HSW, running on seven instead of 14
cores amounts to a reduction of 16\% from 57.0\,J to 47.6\,J. The positive
impacts of L1 blocking can be observed for both microarchitectures,
corresponding to a further decrease of energy consumption by 16\% to 29.1\,J on
IVB resp. 10\% to 43.0\,J on HSW.  \ac{COD} mode
brings energy consumption on HSW further down by 6\% to 40.5\,J.

When it comes to NT stores the differences in the microarchitectures
become visible. On IVB, the per-core performance of the implementation using
L1 blocking and NT stores (dark blue line in
Fig.~\ref{fig:chip-perf}a) is \textit{lower} than that of the version using
regular stores.  A version using NT stores without L1 blocking
(bright blue line) does not even manage to sustain memory bandwidth. Due to this bad
per-core performance, NT stores on SNB/IVB have almost no
positive impact on energy consumption!  On HSW, per-core performance as
expected is exactly $1.5\times$ faster with NT stores, enabling a
lowering in energy consumption by 33\% to 27.0\,J.

Another difference in microarchitecture surfaces when examining
frequency-tuning for potential energy savings.  Before HSW, i.e. on SNB and
IVB, the chip's Uncore frequency was set to match the frequency of the fastest
active core. Because the Uncore contains L3 segments, ring interconnect, and
memory controllers, L3 and memory bandwidth is a function of core frequency on
SNB/IVB. This effect can be observed when setting the frequency to 1.2\,GHz on
IVB (magenta line in Fig~\ref{fig:chip-perf}a): Despite providing the best
energy-to-solution, performance is hurt badly. On HSW, the Uncore is clocked
independently. This leads to a situation in which (due to frequency-induced lower
per-core performance) more cores are required to saturate main memory bandwidth
(cf. magenta line in Fig.\ref{fig:chip-perf}b); however, the lower per-core
power consumption translates into overall energy savings of 11\%, lowering
energy-to-solution to 24.0\,J. These results indicate that energy consumption
could be further reduced if the chips offered frequencies below their current 1.2\,GHz
floor.

Table~\ref{tab:multi-core-results} contains a summary of the final result
pertaining to energy saving for each microarchitecture. The ``reference'' value
corresponds to energy consumption of the baseline implementation running on all
cores clocked at the chip's nominal frequency. Energy-to-solution of the most
energy-efficient version is listed in the ``optimized'' column along with the
number of active cores and their frequency that was used to obtain the results
in the ``configuration'' column.

\begin{figure}[!tb]
    \centering
    \includegraphics[width=0.49\linewidth]{ivb-annotated-crop}\hfill
    \includegraphics[width=0.49\linewidth]{hsw-crop}
    \caption{\label{fig:chip-perf}Performance vs. energy-to-solution for different core counts on IVB (a) and HSW (b).}
\end{figure}

\section{Conclusion}
\label{sec:conc}

We have applied a new energy-optimization approach to a 2D Jacobi solver and
analyzed its effects on a range of recent Intel multi-core chips. Using the
execution-cache-memory (ECM) model single-core software improvements were
described and their accuracy validated by measurements. For the first time, the
ECM model has been (a) extended to incorporate non-temporal (NT) stores and (b)
applied to a Broadwell-EP chip. We found energy consumption can reduces by a
factor of 2.1$\times$ on Sandy Bridge-EP, 2.0$\times$ on Ivy Bridge-EP,
2.4$\times$ on Haswell-EP, and 2.3$\times$ on Broadwell-EP. Further, we found
that while NT stores can increase performance on Sandy- and Ivy Bridge-based E5
processors, a direct positive impact on energy consumption could not be
observed; only in combination with frequency tuning do NT stores offer a better
energy-to-solutions on these architectures.  Measurements indicate that this
problem has been solved on Haswell-EP and Broadwell-EP. Moreover, our results
indicate that future microarchitectures that keep core and Uncore frequencies
decoupled could offer improved energy-efficiency if core frequencies below the
current 1.2\,GHz floor were available.  Beyond these immediate results we have
demonstrated the viability of our energy-optimization approach.

\bibliographystyle{abbrv}
\footnotesize{


\section{Introduction}
Isometric embeddings of metric spaces into infinite-dimensional Banach spaces have a long tradition.
Classical  results are due to
Fr\'echet, Urysohn,  Kuratowski,  Banach,
 \cite{
  Frechet1910,
 Urysohn1927,
Kuratowski1935,
Banach1955}.
For an introduction to the subject we refer to Heinonen's survey \cite{MR2014506}.
The case of  embeddings into
finite-dimensional Banach spaces is much harder, even when one considers bi-Lipschitz embeddings in place of isometric embeddings.
It is a wide open problem to give an intrinsic characterizations of those metric spaces which admit bi-Lipschitz embeddings into some Euclidean space. See for example
\cite{MR1744622,
MR1387955,
Lang-Plaut,
MR2915474,
MR3253787}.

The situation is quite different for quasisymmetric maps (see \cite[Chapter~10-12]{Heinonenbook} for an introduction to the theory of quasisymmetric embeddings).
A metric space quasisymmetrically embeds into some Euclidean spaces if and only if it is doubling (see \cite[Theorem~12.1]{Heinonenbook}).
More specifically, Assouad proved the following result (see \cite{Assouad83}, and also \cite{MR2990137, MR3108866}):
if $(X,d)$ is a doubling metric space and $\alpha\in (0,1)$ then the metric space $(X,d^\alpha)$ admits a bi-Lipschitz embedding into some Euclidean space.
If $d_E$ is the Euclidean distance and $\alpha=\log 2 / \log 3$ then the metric space
$([0,1], d_E^\alpha )$ is bi-Lipschitz equivalent to the von Koch snowflake curve, so
 $(X,d^\alpha)$  is said to be the {\em $\alpha$-snowflake} of  $(X,d)$.

The Assouad Embedding Theorem  is sharp in that there
  are examples (none of which are trivial) of doubling spaces which do not admit bi-Lipschiz embeddings into any Euclidean space,
  even though each of their $\alpha$-snowflakes do.
See
\cite{
Semmes,
Semmes2,
Laakso,
Cheeger-Kleiner2}.
We also  stress that it has been known that snowflakes of doubling   spaces
in general do not isometrically  embed  in any Euclidean space.
Indeed,  the space
  $([0,1], d_E^{1/2})$ does not, see
  \cite[Remark~3.16(b)]{MR2014506}.

  The main aim of this paper
is to show that if some $\alpha$-snowflake of a metric space
isometrically embedds into a finite dimensional Banach space then the
metric space in question is finite.

  Our main result is the following.

\begin{theorem}\label{thm:main}
For any $n \in \mathbb N$ and $\alpha \in (0, 1)$ there is an $N \in \mathbb N$ such that if a metric space $(X,d)$ has cardinality at least $N$ then    $(X,d^\alpha)$ does not admit an isometric embedding into any $n$-dimensional normed linear space.
\end{theorem}

The techniques that we use in Section~\ref{sec:thm:main1} for the proof of Theorem~\ref{thm:main} can also be used to study more general notions of snowflakes.
For this purpose, we introduce general   snowflaking functions.
We say that a function $h: \mathbb{R}_{\geq} \to \mathbb{R}_{\geq}$ is a  {\em snowflaking function} if the following hold:
\begin{enumerate}
\item[(S1)] $h(0) = 0$.
\item[(S2)] $h$ is concave.
\item[(S3)] $ \frac{h(t)}{t} \to \infty$, as ${t \to 0}$.
\item[(S4)] $ \frac{h(t)}{t} \to 0$, as ${t \to \infty} $.
\end{enumerate}
Let $h$ be a snowflaking function.
Then function $h$ is  weakly increasing and, if $d$ is a metric on a set $X$ then $h \circ d$ is also a metric on $X$.
Given a snowflaking function $h$ and a metric space $(X,d)$
we say that the metric space $(X, h \circ d)$ is the  {\em $h$-snowflake} of $(X,d)$.
If $h(t) = t^\alpha$ for some $\alpha \in (0,1)$ then $(X, h \circ d)$ is the $\alpha$-snowflake of $(X,d)$.
In Section~\ref{sec:thm:main} we prove the following.
\begin{theorem}\label{thm:main1}
Let $h$ be a snowflaking function and $(X,d)$ a metric space.
If the $h$-snowflake of $(X,d)$  admits an isometric embedding into some finite-dimensional Banach space then $X$ is finite.
\end{theorem}
\begin{remark}
Note that for general snowflaking functions there may not be any bound on the number of points one can embed, see Remark \ref{rmk:nobound}.
If one removes either of the requirements (S3) or (S4), then we say that $(X, h \circ d)$ is a degenerate snowflake (at zero or at infinity, respectively). Indeed, in such cases the conclusion of Theorem~\ref{thm:main1} does not hold, in general, see Proposition \ref{prop:necessity}.
\end{remark}

We conclude the introduction with few other simple observations about embeddings into Euclidean spaces.
Every $\alpha$-snowflake of $\mathbb{R}^n$, $\alpha\in (0,1)$, isometrically embeds
 into the Hilbert space $\ell^2$ of square summable sequences, see \cite[Remark~3.16(d)]{MR2014506}.
 For any $n\in \mathbb N$, there is a metric space of cardinality $n$ such that for any $\alpha \in (0,1)$ its $\alpha$-snowflake
can be isometrically embedded into $\mathbb{R}^{n-1}$ (just take the vertices of the standard simplex).
There is a 4-point metric space which has an $\alpha$-snowflake which cannot be isometrically embedded into $\mathbb{R}^4$,
and so cannot be isometrically embedded into any Euclidean space, (just take the vertices of the (3,1) complete bipartite graph).
Every finite metric space has an $\alpha$-snowflake which admits an isometric embedding into some Euclidean space, see Proposition~\ref{prop:existence} below.

\section{Corollaries}

In this section we record a few small results and corollaries.
First of all, note that as a corollary of Theorem \ref{thm:main1},  by Ascoli-Arzel\`a Theorem we immediately obtain:
\begin{corollary}\label{lipschitz}
Suppose that $(X,d)$ is infinite.
For any snowflaking function $h$ and $n \in \mathbb{N}$ there is a $\delta > 0$ such that $(X, h \circ d)$ does not admit a $(1 + \delta)$-bilipschitz embedding into any normed linear space of dimension $n$.
\end{corollary}

Theorem \ref{thm:main} shows that there is a bound (depending on $\alpha$ and $n$) on the cardinality of a metric space whose $\alpha$-snowflake admits an isometric embedding into $\mathbb R^n$.
The next easy result says that any finite metric space can be isometrically embedded in some Euclidean space after some $\alpha$-snowflaking.
\begin{proposition}\label{prop:existence}
If $(X,d)$ is a finite metric space of cardinality $n$, then some $\alpha$-snowflake of $(X,d)$ admits an isometric embedding into the Euclidean space
 $\mathbb R^n$.
\end{proposition}

Combining Theorem~\ref{thm:main} and Proposition~\ref{prop:existence} we have:\begin{corollary}
Given a metric space $(X,d)$,
there is an $\alpha$-snowflake of $(X,d)$ that isometrically embedds into some finite-dimensional normed space if and only if $X$ is finite.
\end{corollary}

\begin{proof}[Proof of Proposition \ref{prop:existence}]
We show that $(X,d^\alpha)$ admits an isometric embedding into $\mathbb R^n$ when $\alpha$ is sufficiently close to $1$.
In fact we show that there is an $\epsilon > 0$ such that if $d'$ is any metric on $X$ satisfying $1 - \epsilon < d'(x,y) < 1 +\epsilon$ for all distinct $x,y \in X$ then $(X,d')$ admits an isometric embedding into $\mathbb R^n$.
The proposition follows as $d(x,y)^\alpha \to 1$ as $\alpha \to 0$ for any distinct $x,y \in X$.

Consider the points $p_j:= \tfrac{1}{\sqrt{2}} e_j$ where $e_1, \ldots, e_n$ is the canonical basis of $\mathbb R^n$.
Thus $d(q_i,q_j) = 1$ for all distinct $1 \leq i,j \leq n$.
The proposition is thus a direct consequence of the following claim:

{\bf Claim A.} There is an $\epsilon>0$ such that if $\rho_{ij}\in \mathbb R$, with $1\leq i<j\leq n$ are such that
$|\rho_{ij} - 1 | < \epsilon$, then there exist $q_1, \ldots, q_n\in \mathbb R^n$ such that
$\norm{q_i-q_j} = \rho_{ij}$ for all $1 \leq i < j \leq n$.

To show the claim, consider
the vector subspace $U$ of $\mathbb R^{n\times n}
$
defined by the upper-triangular matrices, i.e.,
the elements of $ U$ are elements of the form $  (a_{ij})_{ij}$ with
$1\leq i<j\leq n$ and $a_{ij}\in \mathbb R$.
Let
$F:\mathbb R^{n\times n} \to U$ be given by
$$F(q_1, \ldots, q_n) := (\norm{q_i-q_j}^2)_{ij} \quad \text{ for }\,q_1, \ldots, q_n\in \mathbb R^n.$$
Denote by
${\bf p}$ the vector of vectors
${\bf p}:= (\tfrac{1}{\sqrt{2}} e_1, \ldots, \tfrac{1}{\sqrt{2}} e_n)$.
Notice that
$F({\bf p})=(1)_{ij}$.

Claim A holds if and only if $F(\bf p)$ lies in the interior of the image of $F$.
Thus it suffices to show that the differential $(dF)_{\bf p}$ has maximal rank and apply the inverse function theorem.

Denote by $q_l^k$ the $k$-th component of an $n$-vector $q_l$.
Notice that $e_l^k= \delta_l^k$, where $\delta_l^k$ is the
Kronecker symbol.
Set $E_{lk} = (\delta_i^l \delta_j^k)_{ij}$ Notice that
$\{E_{lk}, 1\leq i<j\leq n\}$ span the set $U$.
Let us show that any $E_{lk}$ is in the image of $(dF)_{\bf p}$.
Fix $l$ and $k$ so that
$1\leq l<k\leq n$.
By direct calculation,
\begin{eqnarray*}
\dfrac{\partial F( {\bf q})}{\partial q_l^k}  =
\left( 2 \langle q_i -q_j ,
\delta_i^l e_k - \delta_j^l e_k\rangle \right)_{ij}
 =
\left( 2 ( q_i^k -q_j^k )
(\delta_i^l - \delta_j^l ) \right)_{ij}.
\end{eqnarray*}
Evaluating at ${\bf p}$, we get
$$
\left.\dfrac{\partial F( {\bf q})}{\partial q_l^k} \right|_{{\bf q}={\bf p}} =
\left( 2 \frac{1}{\sqrt{2}}( \delta_i^k -\delta_j^k )
(\delta_i^l - \delta_j^l ) \right)_{ij}.
$$
Since $l\neq k$, we have
$ \delta_i^k\delta_i^l\equiv 0$ and
$\delta_j^k\delta_j^l \equiv 0$.
Since $l< k$, and $i<j$, we have
$ \delta_i^k\delta_j^l \equiv 0$.
Hence
$$
\left.\dfrac{\partial F( {\bf q})}{\partial q_l^k} \right|_{{\bf q}={\bf p}} =
\left(  \sqrt{2} ( -\delta_j^k )
\delta_i^l \right)_{ij} = - \sqrt{2}  E_{lk}.
$$
This concludes the proof of Claim A.
\end{proof}

\section{Proofs of the main results}

In the next subsection we first give two simple geometric lemmas which allow us to prove our results for non-euclidean normed vector spaces. The proofs of Theorem~\ref{thm:main} and Theorem~\ref{thm:main1} rely on the Ramsey-theoretic fact that given sufficiently many points in $\mathbb R^n$
there must be a triple that forms as large angle as we want. We recall this result at the end of the subsection in Lemma~\ref{lma:angles}.
In the following two subsections we then prove Theorem~\ref{thm:main} and Theorem~\ref{thm:main1}.
In the final subsection we prove that if $h$ is a degenerate snowflaking function then there exists an infinite metric space whose $h$-snowflake isometrically embeds into $2$-dimensional Euclidean space.

\subsection{Geometric lemmas}

Throughout this section $V$ is an $n$-dimensional normed linear space with norm $\|\cdot\|$.
For the remainder of this section we fix   an inner product $\langle\cdot,\cdot\rangle$ on $V$ for which John ellipsoid property holds. Namely, we have
\begin{equation}\label{John:eq}
B \subset B_{V} \subset \sqrt{n}B,
\end{equation}
where $B_{V}$ is the $\|\cdot\|$-unit ball and $B$ is the the unit ball in the $l_2$ metric associated to the inner product $\langle\cdot,\cdot\rangle$. For $u,v \in V$ we denote the length given by the inner product by $uv := \sqrt{\langle u-v,u-v\rangle}$.

Given two vectors $u,v \in V$ let $\angle(u,v)$ be their angle with respect to the above inner product $\langle\cdot,\cdot\rangle$:
\[\angle(u,v) := \arccos \dfrac{ \langle u,v\rangle}{ \langle u,u\rangle  \langle v,v\rangle}.\]
Given three points $x,y,z
\in V$ we set
$\angle_y(x,z) := \angle (x-y, z-y)$.

Let $x,y,z \in V$ and let  $z' := x + \langle z-x,y-x \rangle(y-x)$ be the projection of $z$ onto the line through $x$ and $y$.

\begin{lemma}\label{lma:closetoEuclidean}
There exist constants $    \varepsilon , C \in (0, \infty)$
 depending only on $n$ such that

 if $\angle_x(y,z) < \varepsilon$ then
 \[
  \big|\|x-z\| - \|x-z'\|\big| \le C\|x-z'\|\angle_x(y,z).
 \]
\end{lemma}
\begin{proof}
By the John Ellipsoid Theorem, \eqref{John:eq}, there is an angle $\theta$ and a height $h$ only depending on $n$ such that
for any point $p$ on the sphere with respect to $\norm{\cdot}$.
The cone $C^- (p, \theta, \ell, 0)$ at $p$ with opening $\theta$ and  height $\ell$ in the direction of $0$ is inside the ball, whereas
the  cone $C^+ (p, \theta, \ell, p)$ at $p$ with opening $\theta$ and  height $\ell$ in the direction of $p$ is outside the ball.

The inequality that we need to prove is translation and dilation invariant. So we assume that $x=0$ and that $z'$ lies on the unit sphere.
Let $\beta := \angle_x(y,z)$. We shall require $\beta  < \varepsilon < \frac\pi2$ for a small enough $\epsilon$.
Let $t, s, t_1, t_2, t_3$ be as in the picture.
\begin{figure}[!h]
  \centering
  \setlength{\unitlength}{0.1\textwidth}
  \begin{picture}(5,1.9)

   \put(0,2){\line(3,-1){6}}
       \put(0,0){\line(1,0){6}}
        \put(0,2){\line(0,-2){2}}
                        \put(0,2){\line(1,0){2}}
                            \put(2,2){\line(0,-2){2}}
                            \put(1,2){\line(0,-2){2}}
                                \put(6,-0.2){$x$}
                                                                \put(4.9,0.1){$\beta$}
           \put(1,0){\line(1,2){1}}
           \put(1,0){\line(-1,2){1}}
    \put(2,-0.22){$t_2$}
        \put(1,-0.22){$z'$}
                \put(0.7,0.05){$\theta$}
                \put(0,-0.22){$t_3$}
                    \put(-0.2,2){$t$}
                                        \put(2.1,2){$t_1$}
    \put(0.8,1.5){$z$}
        \put(1.6,1.5){$s$}
          \end{picture}
  \caption{The point $w$ is between $t$ and $s$.}
\end{figure}

If $\epsilon$ is small enough, then
$\sin \theta \cos \beta - 2 \cos \theta \sin \beta > {1 \over 2}\sin\theta.$
Setting $r:= tz'$ we have
$$t_2y = t_3y- t_3t_2 = r \dfrac{\sin\theta}{\sin \beta} \cos \beta - 2 r \cos \theta> \dfrac{r\sin\theta}{2\sin\beta}.$$
Also, if $\epsilon$ is small enough, then the point $t$ is in $C^+ (p, \theta, \ell, p)$ and so $s\in C^- (p, \theta, \ell, 0)$.
In particular, between $s$ and $t$ there is a point $w$ such that
$\|x-w\| = \|x-z'\|$.

We conclude using, in order, the definition of $w$, the fact that on a line any two norms are a multiple of each other, the triangle inequality, the properties of $s$ and $t$, and the previous bound:
\begin{eqnarray*}
1& \leq &\dfrac{\|x-z\| }{ \|x-z'\|} =
\dfrac{\|x-z\| }{ \|x-w\|}
=
\dfrac{xz }{ xw}
\leq
\dfrac{xw +wz }{ xw}
=
1+ \dfrac{wz }{ xw}
\leq  1+ \dfrac{ts }{ sy} \\
&\leq & 1+ \dfrac{tt_1 + tz' }{ t_2y}
\leq  1+ \dfrac{ 2r \cos \theta + r }{  \frac{r}{2}\sin\theta / \sin \beta}
\leq  1+   C_\theta \sin \beta, \\
\end{eqnarray*}
where $C_\theta  =  ( 4  \cos \theta +2) /  \sin\theta    $.
\end{proof}

The following bound is another easy consequence of the John Ellipsoid theorem.
\begin{lemma}\label{lemma234124}
There exists a constant $K>0$ depending only on $n$ such that
when $  \angle_x(y,z),\angle_y(x,z) \in(0, \frac{\pi}{4})$ we have
 \begin{equation}\label{eq:anglecomparison}
  K^{-1}\frac{\|z'-y\|}{\|z'-x\|} \le \frac{\angle_x(y,z)}{\angle_y(x,z)} \le K \frac{\|z'-y\|}{\|z'-x\|}.
 \end{equation}

\end{lemma}
\begin{proof}
 Set $\alpha =   {\angle_x(y,z)}$, and $ \beta:={\angle_y(x,z)}$. By \eqref{John:eq}
 \begin{eqnarray*}
 \frac{\alpha}{\beta} &\le& \frac{\tan\alpha}{{1 \over 2}\tan\beta}
  =2 \frac{ zz'/xz' }{  zz'/yz'}
  =2 \frac{  yz'} { xz' }
 \le 2\frac{\sqrt{n} \|z'-y\|}{\|z'-x\|}.
 \end{eqnarray*}
By symmetry, we get the other inequality.
 \end{proof}

Our arguments rely on the following Ramsey-theoretic result. Explicit bounds on the number of points that one can have in $\mathbb R^n$ without forming an angle larger than a given bound can be found in \cite{MR841305}. A proof of Lemma~\ref{lma:angles} can also be found in \cite{MR2737262}.

 \begin{lemma}\label{lma:angles}
 For any $n \in \mathbb N$ and $0 < \beta < \pi$ there is an $N \in \mathbb N$ such that if $S \subseteq \mathbb{R}^n$ has cardinality at least $N$ then there are distinct $x,y,z \in S$ such that $\beta \leq \angle(xyz) \leq \pi$.
\end{lemma}

\subsection{Proof of Theorem~\ref{thm:main}}\label{sec:thm:main}

The proof of Theorem~\ref{thm:main} as well as the proof of Theorem~\ref{thm:main1} combines two observations: snowflaking forbids the formation of large angles whereas the fact that we have many points forces such angles to exist.
In the proof of Theorem~\ref{thm:main} the special form of the snowflaking function allows us to directly prove a bound on the cardinality of the snowflaked space that can be embedded in the normed vector space.

\begin{proof}[Proof of Theorem~\ref{thm:main}]
Fix $n$ and $\alpha$ and let $C,\epsilon, K$ be given by Lemma~\ref{lma:closetoEuclidean} and Lemma~\ref{lemma234124}.
 Now take $\theta < \min\{\epsilon,\frac{\pi}{4},\frac{2-2^\alpha}{3C(2K+2^\alpha)}\}$.

We let $N \in \mathbb N$ be the constant in Lemma~\ref{lma:angles} with $\beta = \pi - \theta$.
 We suppose towards a contradiction that $(X,d)$ has cardinality $N$ and that there is an isometric embedding $\iota$ of $(X,d^\alpha)$ into an $n$-dimensional normed vector space $(V,\|\cdot\|)$.
 By Lemma~\ref{lma:angles} there exist three isometrically embedded points $x,y,z \in \iota(X)$ such that
 $\pi - \theta < \angle_z(x,y) \leq \pi$.
 We declare $\delta_z := \angle_z(x,y).$
 We may assume that $\delta_z \ne \pi$.
We also declare $\delta_x := \angle_x(y,z)$ and $\delta_y := \angle_y(x,z)$.
 Note that $0 < \delta_x,\delta_y < \theta$.

 Let $z' \in \mathbb R^n$ be the orthogonal projection of $z$ on the line passing through $x$ and $y$, i.e. $z' := x + \langle z-x,y-x \rangle(y-x)$.
Lemma~\ref{lma:closetoEuclidean} yields
 \begin{equation}\label{eq:dy}
  \left|\|x-z\| - \|x-z'\|\right| \le C\|x-z'\|\delta_x
 \end{equation}
 and
 \begin{equation}\label{eq:dx}
  \left|\|y-z\| - \|y-z'\|\right| \le C\|y-z'\|\delta_y.
  \end{equation}

 We now estimate
 \[
  (\|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha})^{2\alpha}
 \]
 from above.
 Subadditivity of $t \mapsto t^\alpha$ yields:
 \begin{align*}
   (\|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha})^{2\alpha} &= ((\|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha})^2)^{\alpha}\\
   &= (\|y-z\|^\frac2\alpha + \|z-x\|^\frac2\alpha + 2\|y-z\|^\frac1\alpha\|z-x\|^\frac1\alpha)^\alpha\\
   & \le \|y-z\|^2 + \|z-x\|^2 + 2^\alpha\|y-z\|\|z-x\|.
 \end{align*}
 Estimating the obtained terms from above using \eqref{eq:dx} and \eqref{eq:dy} we have
 \begin{align*}
  (\|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha})^{2\alpha} & \le
 \|y-z'\|^2(1+C\delta_y)^2 + \|z'-x\|^2(1+C\delta_x)^2 \\
 & \quad\, + 2^\alpha \|y-z'\|\|z'-x\|(1+C\delta_x)(1+C\delta_y).
 \end{align*}
 Now we use the fact that $C\delta_x < 1$ and Lemma~\ref{lemma234124} to obtain
 \begin{align*}
  \|y-z'\|^2(1+C\delta_y)^2 & \le \|y-z'\|^2(1+3C\delta_y)
  \le \|y-z'\|^2\left(1+3CK\delta_x\frac{\|z'-x\|}{\|z'-y\|}\right)\\
  & \le \|y-z'\|^2 + 3CK\theta\|y-z'\|\|x-z'\|
 \end{align*}
 and similarly
 \[
  \|x-z'\|^2(1+C\delta_x)^2 \le \|x-z'\|^2 + 3CK\theta\|y-z'\|\|x-z'\|.
 \]
 Again by $C\delta_x < 1$ we get $(1+C\delta_x)(1+C\delta_y) < 1+3C\theta$.
 Collecting the estimates together and using the fact that $6CK\theta + 2^\alpha(1+3C\theta) < 2$ we have:
 \begin{align*}
  (\|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha})^{2\alpha} & \le \|y-z'\|^2 + 3CK\theta\|y-z'\|\|x-z'\|\\
  & \quad\, + \|x-z'\|^2 + 3CK\theta\|y-z'\|\|x-z'\|\\
  & \quad\, + 2^\alpha(1+3C\theta) \|y-z'\|\|z'-x\| \\
  & = \|y-z'\|^2 + \|x-z'\|^2\\
  & \quad\, + (6CK\theta + 2^\alpha(1+3C\theta))\|y-z'\|\|z'-x\|\\
  & < \|x-z'\|^2 + \|z'-y\|^2 + 2\|x-z'\|\|z'-y\|\\
   & = (\|x-z'\|+\|z'-y\|)^2 = \|x-y\|^2.
 \end{align*}

 Therefore
 \begin{align*}
  d(\iota^{-1}(x),\iota^{-1}(y)) & =  \|x-y\|^\frac{1}{\alpha} > \|y-z\|^\frac{1}{\alpha} + \|z-x\|^\frac{1}{\alpha}\\
  & =  d(\iota^{-1}(z),\iota^{-1}(y)) +  d(\iota^{-1}(x),\iota^{-1}(z)).
 \end{align*}
This contradicts the triangle inequality in $(X,d)$.
\end{proof}

\subsection{Proof of Theorem~\ref{thm:main1}}\label{sec:thm:main1}

In the proof of Theorem~\ref{thm:main1} we use the same geometric lemmas (Lemma~\ref{lma:closetoEuclidean} and Lemma~\ref{lemma234124}) as in the proof of Theorem~\ref{thm:main}.
However, in the proof of Theorem~\ref{thm:main1} the choice of a sequence of points giving the contradiction depends not only on the snowflaking function $h$, but also on the first element of the sequence. Therefore no upper bound (depending on $h$ and $n$) on the number of points that can be snowflake embedded can in general be obtained in Theorem \ref{thm:main1}.

\begin{proof}[Proof of Theorem~\ref{thm:main1}]

Suppose to the contrary that $X$ is infinite and that there exists an isometric embedding $\iota \colon (X,h\circ d) \to V$ where $V$ is an $n$-dimensional normed vector space. We divide our proof into two cases.
An infinite bounded subset of $\mathbb{R}^n$ is not discrete, so one of the following holds:
\begin{enumerate}
 \item[(i)] $\iota(X)$ is unbounded;
   \item[(ii)] $\iota(X)$ is not discrete.
\end{enumerate}

If (i) holds we will arrive at a contradiction with the condition (S4) of a snowflaking function. If (ii) holds, a contradiction follows with (S3).

\textbf{Case (i): Suppose $\iota(X)$ is unbounded}\\
Observe that (S4) implies that the existence of a function $T \colon \mathbb R_{\ge} \to \mathbb R_{\ge}$ such
that for any $t > 0$ and $S \ge T(t)$ we have
\begin{equation}\label{tapio1}
 \frac{t}{h(t)}\frac{h(S)}{S} \le \frac12.
\end{equation}
Combining \eqref{tapio1} with (S1) and (S2) we get
\begin{equation}\label{tapio2}
 h(S+t) \le h(S) + t\frac{h(S)}{S} = h(S) + \frac{t}{h(t)}\frac{h(S)}{S} h(t)
 \le h(S) + \frac12 h(t).
\end{equation}
Now fix $x_0,x_1 \in X$, $x_0\ne x_1$. Since $(X,d)$ is unbounded, there exists a point $x_2 \in X$ with
 $\angle_{\iota(x_2)}(\iota(x_0), \iota(x_1) )\leq \pi/4$
 and
$d(x_2,x_i) > T(d(x_0,x_1))$ for $i = 0,1$. We continue inductively. Suppose $(x_i)_{i=0}^{N-1} \subset X$ have been chosen.
Now we select $x_N \in X$ satisfying
\begin{equation}\label{iota_x}
\angle_{\iota(x_N)}(\iota(x_i), \iota(x_j) )\leq \pi/4
\text{  and }
 d(x_i,x_N) > T(d(x_i,x_j)),\quad\forall i,j <N .
 \end{equation}

Let $\epsilon, C , K$ be the constants from Lemma~\ref{lma:closetoEuclidean} and Lemma~\ref{lemma234124}.
Set
\begin{equation}\label{eq:delta}
 \delta=\min\{ \epsilon, \pi/4, \frac{1}{2C(1+K)} \}.
\end{equation}
By Lemma~\ref{lma:angles} there exist
$x,y,z$ in $\{\iota(x_\ell)\}_{\ell\in \mathbb N}$ such that $\angle_{z}(x,y) > \pi - \delta$.
By the condition \eqref{iota_x}, there exist
$ i, j, k \in \mathbb N$ with $k > \max\{i,j\}$ such that
$x=\iota(x_j)$,
$y=\iota(x_k)$,
$z= \iota(x_i)$.
Let $z'$ be the orthogonal projection of $z$ to the line passing through
$x$ and $y$. On the one hand, we have that
\begin{align*}
 h(d(x_j,x_k)) & = \|x-y\| \\
 & = \|x-z'\| + \|z'-y\|\\
(\text{Using  Lemma~\ref{lma:closetoEuclidean}})\quad & \ge \|x - z\| - C\|x-z'\|\angle_{x}(z,y)\\
 & \hspace{2cm} +  \|z - y\| - C\|y-z'\|\angle_{y}(z,x)\\
 (\text{Using  Lemma~\ref{lemma234124}})\quad & \ge \|x - z\| - C\|x-z'\|\angle_{x}(z,y)\\
 & \hspace{2cm} +  \|z - y\| - CK\|x-z'\|\angle_{x}(z,y)\\
 & = \|x - z\| + \|z - y\| - C(1+K)\|x - z'\|\angle_{x}(z,y)\\
  & \ge \|x - z\| + \|z - y\| - C(1+K)\|x - z\|\angle_{x}(z,y)\\
 & = h(d(x_i,x_k)) + (1-C(1+K) \angle_{x}(z,y)) h(d(x_i,x_j))\\
 (\text{Using  that } \angle_{x}(z,y)< \delta)   \quad & > h(d(x_i,x_k)) + \frac12 h(d(x_i,x_j)).
\end{align*}
On the other hand,
first notice that since
$h$ is concave and positive, then $h$ has to be weakly increasing.
Secondly, by the definition of the function $T$ we then have that for any three $ i, j, k \in \mathbb N$ with $k > \max\{i,j\}$
\begin{equation}\label{eq:a_bit_shorter}
 h(d(x_j,x_k)) \le h(d(x_i,x_k)+d(x_i,x_j)) \le h(d(x_i,x_k)) + \frac12 h(d(x_i,x_j)).
\end{equation}
Therefore we have a contradiction.

\textbf{Case (ii): Suppose $\iota(X)$ is not discrete}\\
This time we observe that (S3) implies the existence of a function $\tilde T \colon \mathbb R_{\ge} \to \mathbb R_{\ge}$ such
that $\tilde T(r) \le r$ for all $r$ and for any $S > 0$ and $0 < t \le \tilde T(S)$ we have
\eqref{tapio1}, and hence \eqref{tapio2}, using (S1) and (S2).

Let $y$ be an accumulation point of $X$. First we select $x_0 \in X \setminus \{y\}$.
Next we take a radius $r_0 > 0$ so that for all $y_1,y_2 \in B(y,r_0)$ we have both
$\angle_{\iota(x_0)}(\iota(y_1), \iota(y_2) )\leq \pi/4$
 and
$d(y_1,y_2) < \tilde T(d(x_0,y_i))$ for $i = 0,1$. Now we select a point $x_1 \in B(y,r_1) \setminus \{y\}$.
We continue inductively. Suppose $(x_i)_{i=0}^{N-1} \subset X$ have been chosen.
Now we take a radius $r_{N-1} < r_{N-2}$ such that for all $y_1,y_2 \in B(y,r_0)$ we have
\begin{equation}\label{iota_x_new}
\angle_{\iota(x_i)}(\iota(y_1), \iota(y_2) )\leq \pi/4
\text{  and }
 d(y_1,y_2) < \tilde T(d(x_i,y_j)),\quad\forall i <N, j = 1,2.
 \end{equation}
Then we select a point $x_N \in B(y,r_{N-1}) \setminus \{y\}$.

With the points $\{x_i\}$ we now arrive at a contradiction with the same argument as in the case (i).
Let $\delta$ be as in \eqref{eq:delta}. Again by Lemma~\ref{lma:angles} there exist
$x,y,z$ in $\{\iota(x_\ell)\}_{\ell\in \mathbb N}$ such that $\angle_{z}(x,y) > \pi - \delta$,
but this time by the condition \eqref{iota_x_new}, there exist
$ i, j, k \in \mathbb N$ with $k < \min\{i,j\}$ such that
$x=\iota(x_j)$,
$y=\iota(x_k)$,
$z= \iota(x_i)$. Now we continue verbatim the proof in the case (i).
\end{proof}

\subsection{Necessity of (S3) and (S4)}

We end this paper by showing that the conditions (S3) and (S4) of generalized snowflakes are indeed needed for Theorem \ref{thm:main1} to hold.

\begin{proposition}\label{prop:necessity}
 Suppose $h \colon \mathbb R_\ge \to \mathbb R_\ge$ satisfies (S1) and (S2) but fails to satisfy (S3) or (S4). Then there is an
 infinite metric space $(X,d)$ such that $(X,h \circ d)$ admits an isometric embeddeding into $2$-dimensional Euclidean space.
\end{proposition}
\begin{proof}
 We only treat the case where (S4) fails. The case where (S3) fails follows in a similar way.

We construct a sequence of points $\{x_i\}_{i = 0}^{\infty}$ in $\mathbb R^2$ such that $(X,d) := (\{x_i\}_{i = 1}^{\infty},h^{-1}\circ d_E)$
 is a metric space. For this purpose we fix a sequence $\{\alpha_i\}_{i=1}^\infty$ of positive angles such that $\sum_{i=1}^\infty \alpha_i < \frac{\pi}{2}$.
 Depending on the sequence $\{\alpha_i\}_{i = 1}^{\infty}$ and the function $h$ we construct an increasing sequence $\{t_i\}_{i=1}^\infty$ of positive real numbers which determine the Euclidean distance between $x_{i-1}$ and $x_i$.
 For notational convenience we let $c(t) = \frac{h(t)}t$. Notice that by assumption $c(t) \searrow c > 0$ as $t \to \infty$.
 This allows us to select for every $i \in \mathbb N$ a real number $t_i > 0$ such that for all $s,t \ge t_i$ we have
 \begin{equation}\label{eq:tichoice}
  c(s)(c(s+t)-c) + c(t)(c(s+t)-c) \le 2(c(s)c(t) \cos(\pi-\alpha_i) + c(s+t)^2).
  \end{equation}

 Now, using the sequences $\{\alpha_i\}_{i = 1}^{\infty}$ and $\{t_i\}_{i = 1}^{\infty}$ we define the sequence $\{x_i\}_{i \in \mathbb{N}}$ as follows.
 We set
 $x_0 := (0,0)$, $x_1 := (t_1,0)$, and  inductively for $n \geq 2$ declare
 \[
  x_n := x_{n-1} + \left(t_n\sin(\sum_{j=1}^{n-1}a_j),t_n\cos(\sum_{j=1}^{n-1}a_j)\right).
 \]

 In order to see that $(X,d)$ is a metric space we need to check that the triangle inequality holds. For this purpose let $0 \le i < j < k$ be three integers. Let $d_E$ be the Euclidean metric.
 The only nontrivial inequality that we have to verify is
 \[
  h^{-1}(d_E(x_i,x_k)) \le   h^{-1}(d_E(x_j,x_k)) +   h^{-1}(d_E(x_j,x_k)).
 \]
 Denoting $s :=   h^{-1}(d_E(x_i,x_j))$ and $t :=   h^{-1}(d_E(x_j,x_k))$ the above inequality is equivalent to
 \begin{equation}\label{eq:wewant}
  d_E(x_i,x_k) \le h(s+t).
 \end{equation}
 By (S1) and (S2) we can estimate
 \begin{equation}\label{eq:cestim}
  c(t) = \frac{h(t)}t \le \frac{h(t+s) - cs}{t} = c(t+s) + (c(t+s) - c)\frac{s}{t}.
 \end{equation}
 Since $s,t \ge t_j$, by applying the law of cosines,
 \eqref{eq:cestim} and \eqref{eq:tichoice} we obtain
 \begin{align*}
  d_E(x_i,x_k)^2 - h(s+t)^2 & = d_E(x_i,x_j)^2 + d_E(x_j,x_k)^2 \\
  & \qquad - 2d_E(x_i,x_j)d_E(x_j,x_k)\cos(\angle_{x_j}(x_i,x_k)) - h(s+t)^2\\
  & \le h(s)^2 + h(t)^2 - 2h(s)h(t)\cos(\pi-\alpha_j) - h(s+t)^2\\
  & = s^2c(s)^2 + t^2c(t)^2 - 2stc(s)c(t) \cos(\pi-\alpha_j) - (s+t)^2c(s+t)^2\\
  & = s^2(c(s)^2-c(s+t)^2) + t^2(c(t)^2 - c(s+t)^2)\\
 & \qquad - 2st(c(s)c(t) \cos(\pi-\alpha_j) + c(s+t)^2)\\
 \text{(Using \eqref{eq:cestim})} \qquad  & \le st\big(c(s)(c(s+t)-c) + c(t)(c(s+t)-c)\\
  & \qquad- 2(c(s)c(t) \cos(\pi-\alpha_j) + c(s+t)^2)\big)\\
  \text{(Using \eqref{eq:tichoice})} \qquad & \le 0
 \end{align*}
 and thus \eqref{eq:wewant} holds.
\end{proof}

\begin{remark}\label{rmk:nobound}
 The proof of Proposition \ref{prop:necessity} can be modified to show that there is a snowflake function $h$ (satisfying all the conditions (S1)--(S4)) such that for every $n \in \mathbb N$ there exists a metric space $(X_n,d_n)$ with cardinality $n$ so that $(X_n,h\circ d_n)$ embeds
 isometrically into $2$-dimensional Euclidean space.

 Indeed, suppose we are given $n \in \mathbb N$ and we have already defined $h$ on $[0,T_{n-1}]$ and that the slope of $h$ at $T_{n - 1}$ is $c_{n-1} > 0$. Then we can define $h$ on an arbitrary long interval $[T_{n-1},T_n]$ as $h(t) = H(T_{n-1}) + c_nt$, where $0 < c_n < c_{n-1}$. By the proof of Proposition \ref{prop:necessity},
 by taking $S$ large enough, there exist $n$ points in $\mathbb R^2$ with the distance between any two of them between $h(T)$ and $h(S)$ so that
 they are an $h$-snowflake of some metric space. We can define $h$ near $0$ so that it satisfies (S3) and if we take $c_n \searrow 0$ then (S4) is also satisfied. Therefore $h$ will have all the required properties.
\end{remark}