text:
\section{Introduction}

Inverse matrices and inverse $p$-th roots, \ie $A^{-1/p}$ for a given matrix
$A$, are important for various applications in the area of scientific computing.
Examples are preconditioning, solving systems of linear equations and linear
least squares problems, non-linear optimization, Kalman filtering and solving
generalized eigenvalue problems, in particular to solve Schrödinger and Maxwell
equations.

In\label{sec:relwork:linscaling} many of these applications, the involved
matrices are very large, containing billions of entries, and also sparse. For
increasing problem sizes, it becomes more and more important to exploit this
sparsity to save both computational effort and memory resources. However, the
inverse and inverse $p$-th roots of a sparse matrix are typically not sparse
anymore, which makes exploitation of the sparsity difficult. \revision{The novel
approach proposed in this work is to only compute an approximate solution for
the inverse p-th roots of large matrices and enforcing that the result has the
same sparsity pattern as the input matrix. At the same time the method allows
massive parallelization of the required computations.} Hence, the proposed
method adopts the idea of \emph{Approximate Computing}, which denotes
the concept of sacrificing accuracy of computation results in order to increase
the efficiency of these computations~\cite{KlavikMalossiBekasEtAl2014}. \revision{We call the method proposed in this work the \emph{submatrix method}.}

There are multiple applications where approximate solutions for inverse $p$-th
roots are of use. One of these applications is preconditioning, where efficient
computations can be more important than accuracy. Another particularly important
application area are electronic structure methods to approximately solve the
electronic Schrödinger equation~\cite{Kohn1999}. Interestingly it has been shown
that, for large system sizes, the representing matrices become eventually all
sparse~\cite{Goedecker1999}. Underlying this is the concept of ``nearsightedness
of electronic matter'', which states that, at fixed chemical potential, the
electronic density depends just locally on the external
potential~\cite{prodan2005}. Consequently, with increasing system size, the
number of nonzero elements in matrices used to simulate these systems only
increases linearly, \ie the density decreases also linearly as the system size
increases. Specialized solver libraries such as CheSS~\cite{Mohr2017} exploit
this behavior \revision{but still require processing large sparse matrices}. The
method proposed in this work is especially suitable in these applications since
it exploits the increasing sparsity of the matrices \revision{and allows scaling
the parallelism with the matrix size}. Since linear-scaling electronic structure
methods entail approximation in practice~\cite{Richters2014}, it can be
sufficient to calculate an approximation for the inverse $p$-th root of the
involved matrices.

The remainder of this work is structured as follows. In
Section~\ref{sec:relwork} we give a brief overview of related work in the field
of matrix inversion, calculation of $p$-th roots and parallelization of these
operations. In Section~\ref{sec:algo} we present the submatrix method in detail.
Afterwards, we discuss the error that is introduced by using the submatrix
method in Section~\ref{sec:error} and its complexity and scalability in
Section~\ref{sec:complexity}. Finally, we present an implementation of the
proposed method using MPI~\cite{MPI} and OpenMP~\cite{OpenMP} in
Section~\ref{sec:eval} and use it to evaluate the scalability of the method in
practice using a large compute cluster. We conclude in
Section~\ref{sec:conclusion}.

\section{Foundations and Related Work}\label{sec:relwork}
Due to the importance of matrix inversions and inverse $p$-th roots
and the computational complexity of these operations, there is a large variety
of ressources on different numerical methods and on efficient implementation of
these methods. We want to give a brief overview of commonly used algorithms,
ready-to-use implementations and related work on parallelization of the required
calculations. However, an exhaustive coverage of available methods is
outside the scope of this work.

The inverse $X$ of a matrix $A$ fulfils the equation
\begin{align}
AX=I,
\end{align}
where $I$ is the identity matrix. $X$ can therefore be determined by solving
this equation using Gaussian elimination or by calculating and using an LU or
LUP decomposition of $A$. For symmetric matrices, one can compute the
singular value decomposition (SVD), and invert all singular values. A different
approach is to reduce the problem of determining the inverse of an $n\times n$
input matrix to calculating the inverse of two $n/2\times n/2$ matrices and
performing several matrix-matrix multiplications. Following this idea
recursively, inverting a matrix can be reduced to performing only matrix-matrix
multiplications~\cite[Ch. 28.2]{cormen}. Lastly, iterative methods
can be used to find the inverse of a matrix, such as the Newton-Schulz iteration
scheme~\cite{schulz1933}.

In literature, several approaches can be found to parallelize matrix inversion
or calculation of LU and SV decompositions. For example,
Van der Stappen et al.~\cite{vanderStappen1993} present an algorithm for
parallel calculation of the LU decomposition on a mesh network of transputers
where each processor holds a part of the matrix. Shen~\cite{Shen2006} evaluates
techniques for LU decomposition distributed over nodes that are connected via
slow message passing.
Dongarra et al.~\cite{Dongarra2011} demonstrate an optimized implementation of
matrix inversion on a single multicore node, focusing on the minimization of
synchronization between the different processing cores.
There are also algorithms specialized on specific applications, such as the one
described by Lin et al.~\cite{Lin2011} which can be used in 2D electronic
structure calculations to only calculate selected parts of the inverse of a
sparse matrix. For parallel calculation of the SVD, Berry et
al.~\cite{berry2006} provide an extensive overview of parallelizable methods.

For the calculation of the $p$-th root $A^{1/p}$ of a matrix $A$, a commonly
used algorithm is the one described by Higham and
Lin~\cite{higham2011,higham2013}. For the calculation of inverse $p$-th roots,
\ie $A^{-1/p}$, there are also iterative methods available, such as the ones
described by Bini et al.~\cite{bini2005algorithms} and Richters et
al.~\cite{richters17_arxiv}, which reduce the problem to repeated matrix-matrix
multiplications. For symmetric matrices, the solution can again also be computed
by building the SVD of the input matrix and applying the operation of interest
to all singular values of the matrix. \revision{For very large sparse matrices,
calculation of the SVD is typically avoided. Instead methods are used where only
the largest eigenvalues of the matrix are calculated. Iterative methods to do so
are the Lanczos algorithm~\cite{lanczos1950} and the Arnoldi
iteration~\cite{arnoldi1951}.}

Implementations for the calculation of the LU and SV decompositions and matrix
inversion are part of LAPACK~\cite{Angerson1990}, a popular software library for
numerical linear algebra. For solving large sparse systems and calculation of
singular values, ARPACK~\cite{arpack} is a well known library \revision{which is
based on the Arnoldi iteration}. There exist different implementations of these
libraries, as well as bindings for many different programming languages. With
ScaLAPACK~\cite{Choi1996} and P\_ARPACK~\cite{p_arpack}, there exist extensions
of these libraries targeting parallel execution on distributed memory systems
using MPI for message passing.

\section{Algorithm Description}\label{sec:algo}
The fundamental concept of the proposed submatrix method is to divide the large,
sparse input matrix into several smaller and dense submatrices and to apply the
desired operation, such as inversion or calculation of inverse $p$-th roots, to
all of these submatrices instead of the original matrix. Afterwards the results
of these operations are used to construct an approximate result matrix. The
overall procedure is shown in Figure~\ref{fig:overview}. In the following, we
describe the main steps of this method, in particular, building the submatrices
and assembling the final result matrix, in detail. Note that, although we only
discuss a column-based approach for building the submatrices, the method can as
well be applied in a row-based manner since we are dealing with symmetric
matrices. For ease of reading, we will often only mention the inversion of
matrices in the remainder of this section, while the same also holds for the
calculation of inverse $p$-th roots.

\begin{figure*}[t]
  \includegraphics[width=\textwidth]{construction-assembly.pdf}
  \caption{\label{fig:overview}Overview of the submatrix method}
\end{figure*}

\subsection{Building the submatrices}
The process of building the submatrices is shown in Algorithm~\ref{lst:divide}.
To simplify comprehension of the algorithm, all matrices have a dense
representation in our pseudocode. As will be discussed in
Section~\ref{sec:implnotes}, in practice a sparse representation should be used
for the sparse input and output matrices.

\begin{algorithm}[tb]
  \begin{algorithmic}

    \State $n \gets$ number of rows/columns of input matrix
    \State $A[1\dots n][1\dots n] \gets$ input matrix

    \For {$j \gets 1\dots n$}
      \State $R\gets\emptyset$
      \For {$i\gets 1\dots n$}
        \If {$A[i][j]\neq 0$}
          \State $R\gets R \cup \{i\}$
        \EndIf
      \EndFor
      \State $m\gets R$.length()
      \For {$k\gets 1\dots m$}
        \For {$l\gets 1\dots m$}
          \State submatrices$[j][k][l] \gets A[R[k]][R[l]]$
        \EndFor
      \EndFor
      \State indices$[j] \gets R$ \Comment{required later on for result assembly}
    \EndFor

  \end{algorithmic}
  \caption{\label{lst:divide}Construction of submatrices}
\end{algorithm}

To construct the $j$-th submatrix, the $j$-th column of the input matrix $A$ is
evaluated. We determine the set $R$ of row indices $i$ for which $A_{i,j}\neq
0$. The submatrix is then constructed by taking all values $A_{x,y}$ from the
input matrix where \mbox{$x,y\in R$}. For an input matrix of size $n\times n$ we
obtain a set of $n$ submatrices. The size of each submatrix is determined by the
number of nonzero elements in the corresponding column of the input matrix.

\subsection{Performing submatrix operations}
For all of the submatrices, we now perform the operation which should originally
be performed on the input matrix, \ie we either invert all submatrices or
calculate their inverse $p$-th roots. Note that the method and implementation
for these submatrix operations can be freely selected and this choice is
entirely orthogonal to the submatrix method described in this work.

\subsection{Assembling the result matrix}
After having applied the matrix operation of interest to each submatrix, we have
$n$ result submatrices. From these result submatrices we assemble an approximate
solution $X$ for the whole matrix. This procedure is shown in
Algorithm~\ref{lst:combine}.

\begin{algorithm}[tb]
  \begin{algorithmic}

    \State $n \gets$ number of rows/columns of input matrix
    \State indices$[1\dots n] \gets$ from submatrix generation stage
    \State submatrices$[1\dots n][1\dots ?][1\dots ?] \gets$ result matrices

    \State $X\gets$ zeros($n\times n$)
    \For{$j\gets 1\dots n$}
      \State $R\gets$ indices$[j]$
      \State $m\gets R$.length()
      \For{$i\gets 1\dots m$}
        \State $X[R[i]][j] \gets$ submatrices$[j][i][R$.indexof$(j)]$
      \EndFor
    \EndFor

  \end{algorithmic}
  \caption{\label{lst:combine}Assembly of result matrix}
\end{algorithm}

Similar to the construction of the submatrices, the $j$-th column of the final
result matrix is determined by the $j$-th result submatrix. We take the values
from the column of the result submatrix which was originally filled with values
from the $j$-th column of the input matrix, and copy them back to their original
position in the original matrix.

\subsection{Implementation notes}\label{sec:implnotes}
Although the inverse of a sparse matrix is typically not sparse, the approximate
solution provided by the submatrix method exhibits the exact same sparsity
pattern as the input matrix. This allows for efficient implementation of the
method based on matrices in the compressed sparse column (CSC) format which
consists of a value array (\texttt{val}), a list of row indices
(\texttt{row\_ind}) and a list of column pointers (\texttt{col\_ptr}). In
particular, the result assembly stage can be implemented by concatenating the
corresponding columns of all result submatrices to obtain the value array
\texttt{val} for the approximate result matrix. \texttt{row\_ind} and
\texttt{col\_ptr} from the input matrix can be reused for the output matrix
without any changes. If the method is applied in a row-based manner, the same
holds for matrices in the compressed sparse row (CSR) format.

\section{Applicability and Approximation Error}\label{sec:error}
Evidently, the result obtained by the submatrix method is only an approximation
of the correct result. Whether this result is still of use for an application
depends on three aspects: Is the application able to deal with results that
contain a certain error, how large can this error be to be still tolerable and
how large is the error introduced into results by using the submatrix method. In
this section, we first characterize the error based on random sparse matrices
and afterwards apply the submatrix method to \revision{two application scenarios
and show that using the submatrix method to process real-world matrices in these
application yields good results.}
Finally we discuss means to influence the approximation error.

\subsection{Error for random input matrices}\label{sec:error:generic}
To get an impression about the error introduced for arbitrary symmetric
positive-definite matrices, we generate random matrices $A$ using the
\texttt{sprandsym}\footnote{\texttt{sprandsym(size,density,1/condition,kind)} with \texttt{kind}=1} function in Matlab.
\revision{This allows us to sweep over different sizes $n$, densities $d$ and condition numbers $\kappa$ and assess the influence of these matrix properties onto the eror.}
For each set of these
parameters, we generate ten different matrices. For all of these matrices we
then use the submatrix method to obtain an approximate solution $X$ for the
inverse $p$-th root $A^{-1/p}$. To assess the error of these results, we
calculate the spectral norm of the residuals
\begin{align}
\norm{R}_2 &= \norm{X^p A - I}_2.
\end{align}
Since for a precise solution it should hold that $X^p A = I$, $R$ is a good
indicator for the introduced error. We choose the spectral norm of $R$ as a
metric because in contrast to other matrix norms like the Frobenius norm it is
relatively invariant of the matrix size. The spectral norm of a matrix is
defined as
\begin{align}
\norm{M}_2 &= \sqrt{\lambda_\text{max}(M^{H} M)},
\end{align}
where $\lambda_\text{max}(M)$ denotes the largest eigenvalue of $M$ and $M^{H}$
is the Hermitian transpose of $M$.

\begin{figure}[tb]
  \includegraphics[width=\columnwidth,trim=2cm 0.5cm 0 1cm]{surface.pdf}
  \caption{\label{fig:surface}Residual for approximately calculated inverse of random matrices using submatrix method, for different sizes and condition numbers.}
\end{figure}

Our initial evaluation has not shown a significant influence of the density of
the randomly generated matrices onto the precision of the result. We therefore
neglect this parameter in the evaluation presented here, focus on matrices with
density $d=0.05$ and discuss the influence of the size and the condition number
of the matrices. Figure~\ref{fig:surface} shows the relationship between these
matrix properties and the calculated residual for $p=1$. It shows
that the error increases for matrices with larger size and larger condition
numbers. For small matrices, the error stays relatively low even for higher
condition numbers. Similarly, for well-conditioned matrices, the error stays low
even for large matrices.

\iffalse
Both of these aspects can be explained descriptively. For a sparse matrix, its
inverse becomes less sparse with increasing condition number. Since our proposed
method retains the original sparsity pattern, we therefore neglect more
information that would be included in a precisely calculated inverse. That this
effect becomes more and more pronounced for larger matrices is caused by our
error metric which depends on the largest eigenvalue of $R^H R$, which is
expected to only grow for increasing matrix sizes.
\fi

\begin{figure}[bt]
  \includegraphics[width=\columnwidth]{error-size.pdf}
  \caption{\label{fig:error-size}Residual for approximately calculated inverse of random matrices with $\kappa=2$ using submatrix method, in relation to size of input matrix.}
\end{figure}

To demonstrate the latter, we now focus on \revision{well-conditioned} matrices with $\kappa=2$
and $d=0.05$, varying only their size. Results are shown in
Figure~\ref{fig:error-size}. It shows that for a fixed condition number, the
error introduced by using the submatrix method is limited even when further
increasing the matrix size. As shown, this not only holds for calculating the
inverse of a matrix but also for calculating inverse $p$-th roots where $p>1$.

\subsection{Applicability to other matrix operations}
Throughout this work we describe the submatrix method as a method to calculate
inverse $p$-th roots of matrices. \revision{In fact, the method can be applied
to similar matrix operations as well, such as the calculation of positive $p$-th
roots where the error behaves very similar as discussed in Section~\ref{sec:error:generic}.}
\iffalse
\begin{figure}
  \centering
  \includegraphics[width=\columnwidth]{error-other-exponents.pdf}
  \caption{Error for application of submatrix method to positive $p$-th roots.}
\label{fig:error-other-exponents}
\end{figure}

\todo[inline]{Maybe we can drop the figure and half of the following paragraph
to make up space.}
Figure~\ref{fig:error-other-exponents} exemplarily shows the error
obtained when calculating $A^{1/p}$ with positive, integer $p$, using the same
matrices as examined in Section~\ref{sec:error:generic}. However, the residual
$R$ in this case is defined as
\begin{align}
  R &= X^p A^{-1} - I.
\end{align}
\fi
\revision{In contrast to the case of inverse $p$-th roots, the
residual is calculated as $R = X^p A^{-1} - I$, which requires a matrix inversion itself.}
Hence, the residual cannot be calculated efficiently at runtime to assess the
error of an approximately computed result. Inverse $p$-th roots are therefore,
and due to their variety of target applications, presented as the main target in
this work.

\subsection{Application within electronic structure codes}\label{sec:application}
In Section~\ref{sec:error:generic} we demonstrated that using the
submatrix method \revision{for well-conditioned matrices} yields results very similar to a precisely calculated
solution. However,
whether errors are acceptable in an application depends on the the kind of
matrices used in the application and the effect that small deviations have on
the final result. In this section, we show a specific application of the
submatrix method and demonstrate its limited influence on the final result.

We choose one specific time-consuming kernel of effective
single-particle theory, the orthogonalization of a set of $n$ non-orthogonal
basis functions, and  assess the impact of using the submatrix method within
this kernel on  the band-structure energy $E_{\text{BS}}$, which corresponds to
the sum of eigenvalues of the Hamilton matrix $H$. For that purpose, we first
compute the so-called  overlap matrix $S\in \mathbb{R}^{n\times n}$, where
$S_{i,j} = \langle \varphi_i | \varphi_j \rangle$ and $\varphi_i$ are the $n$
non-orthogonal basis functions spanning the Hilbert space. Specifically, we
consider overlap matrices of bulk liquid water systems  for different numbers of
molecules and hence different dimension of the the overlap matrix, as computed
using a Daubechies Wavelet-based density functional theory (DFT)
code~\cite{mohr2014}.

\begin{figure}
  \centering
  \includegraphics[width=.7\columnwidth]{overlapMatrix.pdf}
  \caption{Structure of an examined overlap matrix.}
\label{fig:overlapMatrix}
\end{figure}

Figure~\ref{fig:overlapMatrix} shows one such overlap matrix for
a system of 128 H\textsubscript{2}O molecules. The matrix has a banded structure
and in this case a density of $d=0.25$ and a condition number of around
$\kappa=1.5$. For increasing system sizes, the density decreases lineraly with
$n$, \eg $d=0.12$ for $n=1536$, $d=0.06$ for $n=3072$ etc., while condition
numbers remain approximately constant at around $\kappa=1.5$.

In addition to the overlap matrices $S$, we also extract the density
matrix $P$ as well as the Hamilton matrix $H$ from our Wavelet-based DFT code.
This allows us to calculate the band-structure energy as
\begin{align}
  E_\text{BS} = \text{tr}(PH).
\end{align}
To orthogonalize the Hamilton matrix, we calculate
\begin{align}
  H_\text{ortho} = HS^{-1},\label{eq:ortho}
\end{align}
where $S^{-1}$ is computed using the submatrix method. From this, we again
calculate the band-structure energy as
\begin{align}
  E_\text{BS}^\text{sm} = \text{tr}(SPH_\text{ortho}),
\end{align}
and evaluate the relative error caused by the approximative inversion of $S$ as
\begin{align}
  \Delta E_\text{rel} = \abs{\frac{   E_\text{BS} - E_\text{BS}^\text{sm}   }{   E_\text{BS}   }}.
\end{align}

\begin{table}[t]
  \centering
  \caption{\label{tab:energy}Influence of using the submatrix method for orthogonalization of basis functions onto electronic band structure energy calculations.}
  \begin{tabularx}{\columnwidth}{rRRR}
    \hline\hline
    \Tstrut\Bstrut\textbf{Matrix size} $n$ & $E_\text{BS}$ & $E_\text{BS}^\text{sm}$ & $\Delta E_\text{rel}$\\
    \hline
      768 &   -372.83597 &   -372.83600 & 6.96$\times 10^{-8}$\Tstrut\\
     1536 &   -747.13928 &   -747.13933 & 7.60$\times 10^{-8}$\\
     3072 &  -1492.25282 &  -1492.25297 & 1.01$\times 10^{-7}$\\
     6144 &  -2986.25656 &  -2986.25683 & 8.76$\times 10^{-8}$\\
    12288 &  -5976.31525 &  -5976.31576 & 8.64$\times 10^{-8}$\\
    24576 & -11951.19504 & -11951.19598 & 7.85$\times 10^{-8}$\Bstrut\\
    \hline\hline
  \end{tabularx}
\end{table}

Results are shown in Table~\ref{tab:energy}. For all evaluated matrix
sizes, the relative error caused by using the submatrix method for
orthogonalization is rather small and throughout below or around $10^{-7}$.

\iffalse
\begin{table}[t]
  \centering
  \caption{\label{tab:energyI}\revision{Using the identity matrix $I$ as approximation for $S^{-1}$.}}
  \begin{tabularx}{\columnwidth}{rRRR}
    \hline\hline
    \Tstrut\Bstrut\textbf{Matrix size} $n$ & $E_\text{BS}$ & $E_\text{BS}^{I}$ & $\Delta E_\text{rel}^{I}$\\
    \hline
      768 &   -372.83597 &   -377.349354 & 1.21$\times 10^{-2}$\Tstrut\\
     1536 &   -747.13928 &   -755.763110 & 1.15$\times 10^{-2}$\\
     3072 &  -1492.25282 &  -1510.172178 & 1.20$\times 10^{-2}$\\
     6144 &  -2986.25656 &  -3022.110353 & 1.20$\times 10^{-2}$\\
    12288 &  -5976.31525 &  -6048.829821 & 1.21$\times 10^{-2}$\\
    24576 & -11951.19504 & -12093.831551 & 1.19$\times 10^{-2}$\Bstrut\\
    \hline\hline
  \end{tabularx}
\end{table}
\fi

\revision{
The band structure of the considered overlap matrices may suggest that the low deviation between $E_\text{BS}$ and
$E_\text{BS}^\text{sm}$ comes from the similarity between $S$ and the identity
matrix $I$, and therefore $I$ could be used in Equation~(\ref{eq:ortho}) as an
approximation for $S^{-1}$. However, doing so leads to relative errors between
$1.15\times 10^{-2}$ and $1.21\times 10^{-2}$, \ie five orders of magnite higher
than when using the approximate inverse provided by the submatrix method.}

\iffalse
However, as shown in Table~\ref{tab:energyI}, the relative errors of the
calculated energies are around five orders of magnitude higher in this case,
which demonstrates the quality of the approximation provided by the submatrix
method.
\fi

\subsection{\revision{Application as preconditioner}}\label{sec:preconditioning}
\revision{In the last section we showed an application where the input matrices are
well conditioned and the results provided by the submatrix method are good enough
for the application. Now we demonstrate using the submatrix method to process ill-conditioned matrices in order to obtain a
preconditioner that can be used to iteratively solve systems of linear equations.}

\begin{table}[t]
  \centering
  \caption{\label{tab:preconditioning}\revision{Number of iterations required to solve Equation~(\ref{eq:lineq}) for different matrices $A$ using CG with different preconditioners.}}
  \begin{tabularx}{\columnwidth}{lRRRRR}
    \hline\hline
    \Tstrut\Bstrut\textbf{Matrix} & $n$ & $\kappa$ & \textbf{None} & \textbf{SM} & \textbf{ILU(0)} \\
    \hline

    1138\_bus & 1138 & 8.5$\times10^{06}$ & 2120 & 151 & 139\Tstrut\\
    bcsstk08 & 1074 & 2.6$\times10^{07}$ & --- & 41 & 27\\
    bcsstk09 & 1083 & 9.5$\times10^{03}$ & 194 & 56 & ---\\
    bcsstk10 & 1086 & 5.2$\times10^{05}$ & --- & 85 & 182\\
    bcsstk11 & 1473 & 2.2$\times10^{08}$ & --- & 273 & 477\\
    bcsstk12 & 1473 & 2.2$\times10^{08}$ & --- & 273 & 477\\
    bcsstk13 & 2003 & 1.1$\times10^{10}$ & --- & 409 & ---\\
    bcsstk14 & 1806 & 1.2$\times10^{10}$ & --- & 54 & 262\\
    bcsstk15 & 3948 & 6.5$\times10^{09}$ & --- & 177 & 591\\
    bcsstk16 & 4884 & 4.9$\times10^{09}$ & 464 & 32 & 35\\
    bcsstk21 & 3600 & 1.8$\times10^{07}$ & --- & 224 & ---\\
    bcsstk23 & 3134 & 2.7$\times10^{12}$ & --- & 1269 & ---\\
    bcsstk24 & 3562 & 2.0$\times10^{11}$ & --- & 300 & 244\\
    bcsstk26 & 1922 & 1.7$\times10^{08}$ & --- & 325 & 337\\
    bcsstk27 & 1224 & 2.4$\times10^{04}$ & 907 & 66 & 19\\
    bcsstk28 & 4410 & 9.5$\times10^{08}$ & --- & 668 & 755\\
    bcsstm12 & 1473 & 6.3$\times10^{05}$ & 2790 & 7 & 12\\
    Chem97ZtZ & 2541 & 2.5$\times10^{02}$ & 86 & 10 & 1\\
    crystm01 & 4875 & 2.3$\times10^{02}$ & 70 & 8 & 2\\
    ex10hs & 2548 & 5.5$\times10^{11}$ & --- & --- & ---\\
    ex10 & 2410 & 9.1$\times10^{11}$ & --- & --- & ---\\
    ex13 & 2568 & 1.1$\times10^{15}$ & --- & --- & ---\\
    ex33 & 1733 & 7.0$\times10^{12}$ & --- & 1052 & ---\\
    ex3 & 1821 & 1.7$\times10^{10}$ & --- & --- & ---\\
    ex9 & 3363 & 1.2$\times10^{13}$ & --- & --- & ---\\
    mhd3200b & 3200 & 1.6$\times10^{13}$ & --- & 6 & 3\\
    mhd4800b & 4800 & 8.2$\times10^{13}$ & --- & 6 & 2\\
    msc01050 & 1050 & 4.6$\times10^{15}$ & --- & --- & ---\\
    msc01440 & 1440 & 3.3$\times10^{06}$ & --- & 89 & 155\\
    msc04515 & 4515 & 2.3$\times10^{06}$ & 4411 & 357 & ---\\
    nasa1824 & 1824 & 1.9$\times10^{06}$ & --- & 275 & 264\\
    nasa2146 & 2146 & 1.7$\times10^{03}$ & 282 & 67 & 12\\
    nasa2910 & 2910 & 6.0$\times10^{06}$ & --- & 282 & 760\\
    nasa4704 & 4704 & 4.2$\times10^{07}$ & --- & 1100 & 570\\
    plat1919 & 1919 & 1.2$\times10^{17}$ & --- & --- & ---\\
    plbuckle & 1282 & 1.3$\times10^{06}$ & 1965 & 76 & 69\\
    sts4098 & 4098 & 2.2$\times10^{08}$ & --- & 67 & 119\\
    Trefethen\_2000 & 2000 & 1.6$\times10^{04}$ & 435 & 6 & 5\Bstrut\\

    \hline\hline
  \end{tabularx}
\end{table}

\revision{To demonstrate this application scenario, we obtain sparse, symmetric,
positive definite matrices from the SuiteSparse matrix library~\cite{florida}.
We select all matrices $A$ with size $1000\leq n \leq 5000$ that fulfill these
requirements. For these matrices we solve the system
\begin{align}
  Ax=b,\quad b=[1,1,\dots,1]^T\label{eq:lineq}
\end{align}
using the Conjugate Gradient (CG) method. We set the threshold for the residual
to $10^{-6}$ and limit the number of iterations by $2n$.
Table~\ref{tab:preconditioning} shows the number of iterations required for CG
to converge towards a solution.}
\revision{For preconditioning, we use the submatrix method to obtain
\begin{align}
  K \approx A^{-1/2}.
\end{align}
Instead of solving Equation~(\ref{eq:lineq}), we now solve the system given by
\begin{align}
  K^T A K y = K^T b
\end{align}
using the CG method. The solution $x$ for Equation~(\ref{eq:lineq}) can then be
computed as
\begin{align}
  x = K y.
\end{align}
Again, results are shown in Table~\ref{tab:preconditioning}. For comparison, we
also include the number of iterations required when using an
ILU(0)\footnote{incomplete LU decomposition with zero fill-in} preconditioner.
The results show that using the submatrix method for preconditioning is not only
competitive to the use of ILU(0) but enables CG to converge in more of the
cases.}

\subsection{Controlling the approximation error}
We have demonstrated the use of the submatrix method in \revision{two
applications that can highly benefit from the speedup and the additional
parallelism and still yields good results}. However, there may be applications
that are less tolerant to errors but still can benefit from using the submatrix
method.

If an application requires a lower error than what is provided by the
solution calculated using the submatrix method, an iterative method like those
described by Bini et al.~\cite{bini2005algorithms} and Richters et
al.~\cite{richters17_arxiv} can be used to refine the solution obtained from the
submatrix method. The result obtained by using the submatrix method then acts as
an initial guess for these iterative methods. While we validated that such a
refinement of a solution generated by the submatrix method works in principle
and converges within very few iterations, a detailed evaluation of combining the
submatrix method with iterative methods remains for future work.

In the contrary case, if the application has a particularly high
resiliency against errors in the inverse matrix, the submatrix method can also
be combined with other approximation techniques to achieve further performance
gains. Since using the submatrix method is orthogonal to the implementation of
the operations performed on the single submatrices, these submatrix calculations
can be performed in an approximate manner as well. Using an iterative method,
precision can be scaled by the number of iterations. Additionally, calculations
can be performed using low precision arithmetic or other Approximate Computing
techniques~\cite{lass17_esl}.

\section{Complexity and Scalibility}\label{sec:complexity}
We now want to discuss the time complexity and scalability of the submatrix
method and show that, although for an $n\times n$ matrix $n$ submatrices need to
be processed, it can still provide a significant reduction in time required for
determining a matrix inverse or its inverse $p$-th root.

Note that the considerations in this section only hold if the density
of the input matrix shows in each of its columns (or rows). An obvious
counterexample are arrowhead matrices, where using the submatrix method cannot
provide any speedup since the first submatrix has the same size as the original
input matrix.

\subsection{Single-threaded scenario}
We first want to discuss the general complexity of matrix inversion, both using
conventional methods and using our proposed submatrix method. While from a
theoretical standpoint, inversion of matrices is not harder than multiplication,
and therefore $\mathcal{O}(n^{2.81})$ using Strassen's algorithm~\cite[pp. 79,
829]{cormen}, or even $\mathcal{O}(n^{2.373})$ using Coppersmith and Winograd's
algorithm~\cite{LeGall2014}, in practice methods such as Gaussian elimination or
building and using the LU decomposition for inversion which have time complexity
$\mathcal{O}(n^3)$  are commonly used. In the following, we define $I(n)$ as the
time required for a precise matrix inversion, abstracting from a concrete
implementation.

For a sparse $n\times n$ input matrix, using the submatrix method requires
performing $n$ matrix inversions for smaller but dense matrices. To be more
efficient in a single-threaded application scenario, these submatrices have to
be significantly smaller than the original input matrix. In the
following, we assume a uniformly filled, sparse input matrix. Let $d$ be the
density of this matrix, then the average size $m\times m$ of the submatrices is
determined by:
\begin{align}
m &= d\cdot n.
\end{align}
If the density $d$ is small enough, such that
\begin{align}
n\cdot I(d\cdot n) &< I(n),\label{eq:breakeven}
\end{align}
then the submatrix method has lower run time than a precise inversion, even in a
single-threaded scenario.

We now want to determine, by what rate the density $d$ has to decrease so that
for increasing matrix sizes the asymptotic run time does not grow faster than
using conventional methods for matrix inversion. We therefore assume that matrix
inversion has at least time complexity $n^2$, \ie $I(n) = \Omega(n^2)$. To
fulfill Equation~(\ref{eq:breakeven}), it then needs to hold that for
sufficiently large $n$
\begin{align}
n \cdot {(d \cdot n)}^2 &< n^2\nonumber\\
d^2 &< n^{-1}\nonumber\\
d &< n^{-0.5}.
\end{align}
From this we can deduce the following asymptotic relation:
\begin{align}
d = \mathcal{O}(n^{-0.5}) &\Rightarrow S(n,d) = \mathcal{O}(I(n)),\label{eq:o-notation}
\end{align}
where $S(n,d)$ is the time required to calculate an approximate inverse using
the submatrix method. If $d$ decreases faster than with rate $n^{-0.5}$, then
$S(n,d)$ increases slower than $I(n)$ for larger $n$. Using methods where
$I(n)=\Omega(n^3)$ further relaxes the requirements on $d$, in that
$d=\mathcal{O}(n^{-1/3})$ suffices to fulfill Equation~(\ref{eq:breakeven}) for
large $n$.

Note that we neglect the time required for building the submatrices and
assembling the final result matrix, as their influence on execution time is
negligible compared to the involved matrix inversions in asymptotic
considerations.

\subsection{Parallel execution of submatrix operations}
Although using the submatrix method can reduce execution time even in a
single-threaded environment for very sparse matrices, its strength is to allow
massively parallel execution. All submatrix operations are entirely independent
from each other such that the inversion of an $n\times n$ matrix can be
distributed over $n$ compute nodes. Each compute node can construct its own
submatrix from the input matrix. The final result matrix has to be assembled on
a single node but as described in Section~\ref{sec:implnotes}, this step
consists of a simple concatenation of $n$ arrays. Communication between nodes is
only required for initial data distribution and for the final collection of all
results. Provided that $n$ compute nodes can be used for execution of the
algorithm, a speedup is already achievable if all submatrices are significantly
smaller than the original input matrix, \ie each column of the original matrix
contains a significant fraction of zero-elements.

\subsection{Application to Electronic Structure Methods}\label{sec:complexity:linscaling}
In Section~\ref{sec:relwork:linscaling}, we motivated our method with the
principal of linear scaling techniques in density functional theory, based on
the nearsightedness of electronic matter. With respect to the matrices for which
an inverse or an inverse $p$-th root needs to be calculated, this means that
while for growing systems the total number of matrix elements increases
quadratically, the density of the matrix decreases linearly with $n^{-1}$.
Consequently, the number of nonzero elements in the matrix increases only
linearly with $n$.

Based on this fact, the submatrix method is particularly suitable for solving
these problems. In particular, since
\begin{align}
  n^{-1} < n^{-0.5},
\end{align}
the density of matrices decreases faster than required in
Equation~(\ref{eq:o-notation}). From that it follows that the asymptotic run
time of the submatrix method in a single-threaded environment is limited by that
of a precise inversion for the appliations discussed here.

Again, the strong advantage of the submatrix method is the possibility of
parallel execution on many compute nodes. In the case of linear scaling methods
in density functional theory, this means that for growing systems the execution
can be parallelized onto more and more nodes while the size of the single
submatrices stays constant. As long as the number of compute nodes can be scaled
with $n$ as well, the overall execution time can even be held constant.

\section{Performance Evaluation}\label{sec:eval}
To evaluate the performance and scalability of the proposed method, we built a
distributed implementation using MPI and OpenMP\@. We run this implementation on
a compute cluster comprised of 65 nodes. Each node features two Intel Xeon
E5-2670 CPUs with a total of 16 CPU cores and 64~GByte of memory. All nodes are
connected via 40~Gbit/s QDR InfiniBand. We use one node as a control node,
leaving the remaining 64 nodes with a total of 1024 CPU cores for handling the
workload.
In the following, we first describe details of our implementation, and
then present results obtained from running our implementation on our compute
cluster.

\subsection{Details of our Implementation}
Our implementation makes use of Intel MPI~\cite{IntelMPI} to distribute work
over a large number of compute nodes and to collect all results in order to
build up the final result matrix. The MPI rank 0, in the following called main
process, reads the input matrix stored in CSC format from persistent storage
into memory. Metadata such as an identifyer for the matrix, as well as its total
size and number of nonzero elements, are then sent via \texttt{MPI\_Bcast} to
all nodes.

\subsubsection{Data distribution and work assignment}
There are different possible ways to make the input matrix available to all
other MPI ranks, which we call worker processes in the following. In principal,
it would be sufficient to send single submatrices to the workers which then
perform the inversion. In this case, all submatrices would have to be
constructed within the main process, which would clearly present a bottleneck.
Instead, we make the whole input matrix available to all worker processes which
then autonomously construct their submatrices. In our environment all systems
have access to a shared file system which allows all processes to read in the
input matrix from persistent storage. Since this scenario cannot generally be
assumed, we additionally implemented distribution of data via
\texttt{MPI\_Bcast} to all worker processes. We found that, for the data we use
in our evaluation, \revision{both variants provide comparable performance}.

Assuming we have $w$ workers, each worker needs to process \mbox{$x=n/w$}
submatrices. In our implementation, each worker processes a contiguous set of
submatrices, \ie the worker with rank $k$ is responsible for submatrices
$(k-1)x$ to $kx-1$. Therefore, depending on its rank and the total number of
ranks, each worker process can determine autonomously, which of the submatrices
it has to process.\footnote{Note that if the number of worker processes does not
divide the size of the matrix, some workers need to process one additional
submatrix.} It builds the submatrix according to Algorithm~\ref{lst:divide} and
calls the LAPACK functions \texttt{dgetrf} to obtain an LU decomposition and
\texttt{dgetri} to calculate the inverse of the submatrix. In our evaluation we
use Intel MKL~\cite{IntelMKL} as a highly optimized implementation for these
LAPACK routines. After inversion, the worker selects the section of the result
matrix which is relevant for the final result matrix and stores it in a buffer.
Since the main process just needs to concatenate these buffers to create the
final result matrix, a single call to \texttt{MPI\_Gatherv} is sufficient to
perform data collection and assembly after all submatrices have been processed.

\subsubsection{Multi-threading using OpenMP}
The described implementation already allows to distribute the load over many
nodes and CPU cores. To utilize multiple CPU cores on a single node,
multiple MPI ranks could be placed on a node, or multiple cores could be used for
processing a single submatrix by using a multi-threaded LAPACK implementation.
Having multiple MPI ranks on the same node comes at the cost of data duplication
in memory and overall increased MPI communication load. Using multiple cores for
a single submatrix operation has shown to provide lower speedup than additional
parallelization of submatrix operations.

In our implementation, we therefore use OpenMP to process $c$ submatrices in
parallel on a node featuring $c$ CPU cores. We do so by calling all submatrix
operations within an OpenMP parallel for loop:

\begin{center}
\texttt{\small \#pragma omp parallel for schedule(dynamic)}
\end{center}

To allow OpenMP to fully utilize the available CPU cores, we explicitly disable
the multi-threading functionality provided by Intel MKL by calling
\texttt{mkl\_set\_num\_threads(1)}.

\subsubsection{Limitations of our Implementation}\label{sec:limitations}
As described, each worker process is responsible for a contigous set of
submatrices and all workers are responsible for the same number ($\pm 1$) of
submatrices. This can lead to workload imbalance between the different workers,
if the input matrix exhibits a pattern such that certain sets of columns contain
significantly more or significantly fewer nonzero values than other sets of
columns. It is important to note that this is a limitation of our implementation
and not a conceptual issue of the proposed submatrix method. In practice, there
are different ways to deal with this issue in order to create an optimized
implementation that does not exhibit this load imbalance:

\emph{Shuffling input matrices:} To balance the load between all worker
processes, the mapping between submatrices and workers can be shuffled randomly.
Clusters of full colums which result in larger submatrices would then not be
assigned to a single worker but distributed over all workers. This can, for
example, be implemented using a pseudo-random but deterministic permutation, so
that each worker can still autonomously determine the submatrices it is
responsible for.

In our implementation, each worker concatenates the results of its submatrix
operations in a buffer which is then sent as a whole to the main process. The
main process therefore only needs to collect and concatenate $w$ arrays for $w$
worker processes. If submatrices are shuffled, collection and cocatenation of
$n$ arrays would be required in the main process instead. \revision{Apart from
this, there is no additional computational effort required for this load
balancing technique.}

\emph{Dynamic work scheduling:} Instead of assigning a fixed set of submatrices
to a worker process, work can be scheduled dynamically. Each worker could
request work packages from the main process using MPI\@. The size of these work
packages can be chosen in the range from one single submatrix up to $n/w$
submatrices in order to trade off load balancing and additional communication
effort. \revision{Concepts similar to OpenMP's \emph{guided} scheduling could also be implemented to
minimize scheduling overhead.}
Note that in our implementation, we already use dynamic work scheduling for the
parallel processing of multiple submatrices on a single node by using OpenMP's
\emph{dynamic} scheduler.

\subsubsection{Availability}
Our prototype implementation as well as scripts used in our evaluation are
published under MIT license and can be found on
\url{https://github.com/pc2/SubmatrixMethod}.

\subsection{Results}

\subsubsection{Scalability for increasing number of CPU cores}\label{sec:res:scalability}
We use our implementation of the submatrix method to calculate an approximate
inverse of multiple random matrices with size $n=32768$, condition number
$\kappa =2$ and density $d=0.01$. We vary the number of utilized CPU cores in
the range from 1 to 1024 and measure the total wall clock time required to
obtain a result. \revision{We consider a set of balanced matrices whose columns have
roughly the same number of nonzero
elements\footnote{generated using \texttt{sprandsym(size,density,1/condition,kind)} with
\texttt{kind}=2} and a set of unbalanced matrices which exhibit visible patterns
in the distribution of
values\footnote{generated using \texttt{sprandsym(size,density,1/condition,kind)} with
\texttt{kind}=1}.}
Results
are shown in Table~\ref{tab:scaling-d1} and Figure~\ref{fig:scaling-d1}. \revision{As a reference}, we also show the time required for a precise matrix inversion using
Intel MKL's implementation of the \texttt{dgetrf} and \texttt{dgetri} routines,
utilizing up to 16 CPU cores on a single node.

\begin{figure}
  \includegraphics[width=\columnwidth]{scaling-d1.pdf}
  \caption{\label{fig:scaling-d1}\revision{Scalability of the submatrix method for a random matrix of size $n=32768$ and density $d=0.01$.}}
\end{figure}

\begin{table}[b]
  \centering
  \caption{\label{tab:scaling-d1}\revision{Time in ms required for inversion of a matrix with size $n=32768$ and density $d=0.01$ using the submatrix method.}}
  \begin{tabularx}{\columnwidth}{rRRRR}
    \hline\hline
    \Tstrut & \multicolumn{2}{c}{\textbf{Balanced matrix}} & \multicolumn{2}{c}{\textbf{Unbalanced matrix}}\\
    \Bstrut\textbf{Cores} & \textbf{Wall time} & \textbf{Speedup} & \textbf{Wall time} & \textbf{Speedup}\\
    \hline
         1 & 578,140 & 1.0 & 1,150,366 & 1.0\Tstrut\\
         2 & 255,081 & 2.3 & 586,778 & 2.0\\
         4 & 118,534 & 4.9 & 304,941 & 3.8\\
         8 & 63,644 & 9.1 & 162,792 & 7.1\\
        16 & 32,405 & 17.8 & 82,571 & 13.9\\
        32 & 16,216 & 35.7 & 42,760 & 26.9\\
        64 & 8,485 & 68.1 & 22,692 & 50.7\\
       128 & 4,242 & 136.3 & 12,447 & 92.4\\
       256 & 2,339 & 247.2 & 7,402 & 155.4\\
       512 & 1,293 & 447.1 & 5,447 & 211.2\\
      1024 & 870 & 664.5 & 4,765 & 241.4\Bstrut\\
    \hline\hline
  \end{tabularx}
\end{table}

\iffalse
It shows that for the generated input matrices, our implementation of the
submatrix method scales linearly up to 256 CPU cores and after that the curve
slowly starts flattening out. The reason for this is the workload imbalance
described in Section~\ref{sec:limitations}, \ie some workers have already
finished their calculations while others are still running. This implementation
weakness shows up for the matrices we generated with Matlab's \texttt{sprandsym}
function because these exhibit structural patterns, causing contiguous sets of
columns having significantly more nonzero elements than others. The less
submatrices each worker has to solve, the more pronounced is the effect of
workload imbalance.

We model the run time of an enhanced implementation of the submatrix method that
balances the load by taking the average run time of all workers and adding the
time required for communication. This idealized model is also shown in
Table~\ref{tab:scaling-d1} and Figure~\ref{fig:scaling-d1}. It shows that
without load imbalance, the submatrix method keeps scaling linearly also above
256 CPU cores.

To verify that the difference between our measured run times and the idealized
model is in fact caused by patterns in the generated input matrices and can be
defeated with additional implementation effort, we repeated the measurements
using the largest overlap matrix from our application scenario described in
Section~\ref{sec:application}. \revision{In this case, the measured run time on
1024 cores exceeds that of the idealized model by only 10\%.}
\fi

\revision{The results show that the submatrix method overall scales well
over a large number of processors. Comparing the data for
balanced and unbalanced matrices, two distinct effects can be observed:
\begin{enumerate}
\item Even for a low number of CPU cores, the submatrix method performs better for balanced matrices. The reason for this is that if some columns contain significantly more nonzero elements than others, the resulting submatrices are larger in size and the time to process them increases cubically with their size.
\item For the imbalanced matrices in our scenario, the curve starts to flatten at around 256 cores and scaling beyond 512 cores provides diminishing returns. The reason for this is that the number of submatrices per worker becomes small enough such that load imbalance between workers has an increasing effect. This effect could be countered by implementing some form of load balancing, as discussed in Section~\ref{sec:limitations}.
\end{enumerate}
For over 512 cores, even for the balanced matrices in our scenario the additional speedup is limited. This is caused by the overall short runtime of the algorithm and therefore increased influence of communication time (around 32\%).
}

\iffalse
As shown in Figure~\ref{fig:balanced-matrix}, in this case our implementation
achieves nearly the same performance as the idealized model. The slight drop in
scalability above 256 cores even for the idealized model is caused by the
overall low run time for this smaller input matrix (below 200ms) and therefore
increasing influence of communication delays.

\begin{figure}
  \includegraphics[width=\columnwidth]{balanced-matrix.pdf}
  \caption{\label{fig:balanced-matrix}Comparison between our implementation and the idealized model for an overlap matrix from Section~\ref{sec:application} (size $n=24576$, density $d=0.007$).}
\end{figure}
\fi

On a single node, Intel MKL as well nearly scales linearly with the number of
CPU cores. Only for 16 cores there is a slight efficiency drop\revision{, likely caused by the NUMA architecture of our compute nodes}. Using ScaLAPACK
to distribute execution of the utilized library functions over multiple nodes
may allow to further increase the number of CPU cores. However, due to
increasing communication overhead, the potential for scaling is limited in this
case. Related work that uses ScaLAPACK for matrix inversion, describes
decreasing performance for execution on more than 64 CPU
cores~\cite{Kolberg2008}.

\subsubsection{Run time for growing matrices}
We now evaluate how the total execution time develops for increasing matrix
sizes, given a fixed number of CPU cores.
We consider two different scenarios: a fixed density of $d=0.01$ \revision{and matrix sizes ranging from $2^{11}$ to $2^{18}$} and a density that decreases
linearly with $n$ as encountered in applications like electronic structure
methods. For the latter, we set $d=0.16 \cdot 1024 / n$ \revision{and consider sizes from $2^{10}$ to $2^{20}$}.

The results of this evaluation are shown in Table~\ref{tab:matrix-sizes} and
Figure~\ref{fig:matrix-sizes}.
\revision{In the table we also show the fraction of the total wall clock time spent on communication and the fraction of compute time spent on building the submatrices. Note that the assembly of the result matrix is performed implicitly by \texttt{MPI\_Gatherv} and therefore accounted as communication time.}
It clearly shows that for matrices
with linearly decreasing density, the required run time \revision{only increases linearly with the matrix size}, as expected based on the discussion in
Section~\ref{sec:complexity:linscaling}. Combining this result with the
possibility for linear performance scaling with the number of CPU cores, the run
time can be held constant by increasing the number of cores with $n$ for growing
matrices. \revision{The data also shows that for increasing size of the submatrices, as shown in the upper half of Table~\ref{tab:matrix-sizes}, the overhead required for communication and for building the submatrices decreases. For fixed-size submatrices, the overhead stays relatively constant. Note that times for communication fluctuate in our measurements due to shared usage of the underlying InfiniBand network.}

\begin{figure}
  \includegraphics[width=\columnwidth]{matrix-sizes.pdf}
  \caption{\label{fig:matrix-sizes}\revision{Required time for inversion of a matrix using the submatrix method (1024 cores) and Intel MKL (16 cores).}}
\end{figure}

\begin{table}[t]
  \centering
  \caption{\label{tab:matrix-sizes}\revision{Time in ms required for inversion of a matrix using the submatrix method on 1024 cores.}}
  \begin{tabularx}{\columnwidth}{rRRRR}
    \hline\hline
    \Tstrut\Bstrut\textbf{Size} & \textbf{Density} & \textbf{Wall time} & \textbf{MPI Comm.} & \textbf{Submat. Constr.}\\
    \hline
        2,048 & $1.0\times 10^{-2}$ & 1 & --- & ---\Tstrut\\
        4,096 & $1.0\times 10^{-2}$ & 2 & --- & ---\\
        8,192 & $1.0\times 10^{-2}$ & 14 & 57.1\% & 58.8\%\\
       16,384 & $1.0\times 10^{-2}$ & 87 & 39.1\% & 60.5\%\\
       32,768 & $1.0\times 10^{-2}$ & 868 & 31.7\% & 57.8\%\\
       65,536 & $1.0\times 10^{-2}$ & 8,380 & 13.5\% & 56.8\%\\
      131,072 & $1.0\times 10^{-2}$ & 87,977 & 5.0\% & 47.0\%\\
      262,144 & $1.0\times 10^{-2}$ & 1,085,176 & 1.5\% & 36.8\%\Bstrut\\
    \hline
      1,024 & $1.6\times 10^{-1}$ & 9 & 22.2\% & 61.3\%\Tstrut\\
      2,048 & $8.0\times 10^{-2}$ & 15 & 26.7\% & 61.0\%\\
      4,096 & $4.0\times 10^{-2}$ & 32 & 37.5\% & 60.5\%\\
      8,192 & $2.0\times 10^{-2}$ & 48 & 33.3\% & 60.3\%\\
      16,384 & $1.0\times 10^{-2}$ & 96 & 38.5\% & 60.4\%\\
      32,768 & $5.0\times 10^{-3}$ & 220 & 51.7\% & 60.9\%\\
      65,536 & $2.5\times 10^{-3}$ & 482 & 54.4\% & 61.4\%\\
      131,072 & $1.3\times 10^{-3}$ & 990 & 54.7\% & 62.2\%\\
      262,144 & $6.3\times 10^{-4}$ & 1977 & 55.2\% & 63.0\%\\
      524,288 & $3.1\times 10^{-4}$ & 4020 & 56.8\% & 63.7\%\\
      1,048,576 & $1.6\times 10^{-4}$ & 7609 & 49.9\% & 64.0\%\Bstrut\\
    \hline\hline
  \end{tabularx}
\end{table}

\section{Conclusion}\label{sec:conclusion}
In this work we presented the submatrix method, which can be used to calculate
an approximate inverse of matrices, as well as inverse $p$-th roots. Following
the idea of Approximate Computing, it allows the result to deviate from an
exactly calculated solution in order to utilize the sparsity of the input matrix
and to allow massively parallel execution of the involved calculations. For an
$n\times n$ matrix, the workload can be distributed over $n$ nodes.

A particularly well suited application for the submatrix method are electronic
structure methods in density functional theory. In these applications, for
growing matrices
their density decreases linearly at the same time. In this
case, the submatrix method exhibits a \revision{linear} increase in execution time for
growing systems. As long as the number of available CPU cores can be scaled with
the same rate, execution time can even be held constant.

We showed that the error introduced by using the submatrix method is limited for
well-conditioned input matrices \revision{and demonstrated its use for
preconditioning of ill-conditioned matrices}. We discussed the scalability of
the algorithm both theoretically and in a practical evaluation on a large
compute cluster.

\revision{An emerging question for future work is how this method can be
combined with iterative methods to refine the solution or to increase efficiency
of the submatrix operations. Additionally, the dense nature of the submatrices
make the submatrix operations well suited for the use of accelerator hardware
which may be explored in future to develop high-performance implementations of
the submatrix operations.}

\begin{acks}
This project has received funding from the \grantsponsor{erc}{European Research
Council (ERC)}{https://erc.europa.eu/} under the European Union's Horizon 2020
research and innovation programme (grant agreement No~\grantnum{erc}{716142})
and from the \grantsponsor{dfg}{German Research Foundation
(DFG)}{http://www.dfg.de/en/index.jsp} under the project PerficienCC
(grant agreement No~\grantnum{dfg}{PL 595/2-1}). The evaluation was
performed on resources provided by the Paderborn Center for Parallel Computing
(PC\textsuperscript{2}).
\end{acks}

\bibliographystyle{ACM-Reference-Format}


BSF Recruitment 2019: The online process for 1072 Head Constable Post has been started at the official website of BSF i.e. bsf.nic.in. Eligible candidates can apply to the posts by 12 June 2019.
The online applications for BSF Recruitment 2019 for 1072 Head Constable Posts were started from 14 May and will continue till 12 June 2019. This is the biggest opportunity for the candidates who wish to make their carrier in BSF.
Candidates can check a brief detail about BSF Recruitment 2019 for 1072 Head Constable Posts by scrolling down.
Interested candidates can apply to the post through the online mode only. The online application for BSF Head constable Jobs 2019 can be submitted upto 12 June 2019. Candidates can check the details in the provided PDF.


text:
Линетчина () — село на Украине, находится в Оратовском районе Винницкой области.

Код КОАТУУ — 0523180403. Население по переписи 2001 года составляет 8 человек. Почтовый индекс — 22653. Телефонный код — 4330.
Занимает площадь 0,183 км².

Адрес местного совета
22653, Винницкая область, Оратовский р-н, с. Балабановка, ул. Ленина, 25

Ссылки
 Линетчина на сайте Верховной рады Украины

Населённые пункты Оратовского района

redpajama_set_name: RedPajamaWikipedia


MHEC is located in the Nancy Grasmick Building, Maryland State Department of Education (MSDE), next to the east side of the parking facility. Please enter from the Liberty Street entrance to access elevators directly to the tenth floor. Planning, Academic Affairs and Research and Finance Policy departments are on the tenth floor. The Office of Student Financial Aid and Information Systems departments are on the ground floor.
The parking facility mentioned above is operated by Arrow Garage and is next to the Nancy Grasmick Building (MSDE). Rates are $8 for one hour, $11 for two hours and $16 for all day. There is additional parking located within 2 blocks of the MSDE building on Fayette, Howard, Liberty and Lombard Streets.
The Light Rail travels from Baltimore County's Hunt Valley through the heart of Downtown Baltimore to Cromwell Station/Glen Burnie in Anne Arundel County. Seven light rail stops are in Downtown, Baltimore. There is a stop on Howard Street, half a block from the entrance to MHEC. Hours of operation: Mon - Fri 6am-11pm; Sat 8am-11pm; Sun 11am-7pm.
The Baltimore Metro runs between Owings Mills and Johns Hopkins Hospital. It operates Mon-Fri - 5am-midnight and Sat, Sun and Holidays - 6am-midnight. The nearest stop to the MHEC entrance is the Metro-Charles Center Station, 2 blocks east.
Location of MHEC entrance on Liberty St. in Baltimore is shown on the right side of picture as you look south on Liberty St.
For additional direction information please call MHEC at (410) 767-3301.


For over four decades, MAC Jets has provided a variety of aviation services including aircraft charter, aircraft sales and management, aircraft maintenance and parts, and flight training. With only one call we are able to provide an extensive range of services in a safe, reliable manner with truly professional and courteous assistance.
Maine Aviation Sales provides Worldwide Aircraft Sales, Buying & Brokerage along with Management. Our maintenance department, Aircraft Maintenance of Maine provides Quality Professional Services to include Annual Inspections, Engine Changes, General Maintenance and Parts. Our Charter Department, Portland Air Charter provides on Demand Service 24/7 in our Pressurized Executive Cabin Class Cessna 421B. Visit our Website to learn more about our Business, www.maineaviationsales.com.
People ship millions of pounds of cargo through the Jetport annually. From birthday presents for the kids to lobsters for friends, we have cargo carriers that can transport your packages to whatever destination you choose.


Concerned about global warming, Dr. Peter Kalmus cut his emissions to below a tenth of the US average, and found that this made his life more interesting, satisfying, and joyful -- not less. He started riding a bike, driving on vegetable oil, and engaging in "slow travel" without airplanes. He switched to "freegan vegetarianism" and grows much of his own food. He has simplified and decarbonized his life, and he has never looked back.
He now does everything he can to spread this message and to help humanity transition away from fossil fuels. He believes that the problem is not only technological or scientific, but also spiritual: it calls humans to set aside the notion that we are exceptional within the web of life, and instead to live, gratefully, in alignment with the biosphere.
Dr. Peter Kalmus is a scientist studying clouds and climate change at NASA's Jet Propulsion Laboratory (speaking on his own behalf). He holds a PhD in physics from Columbia University. Peter is a contributing editor and columnist at YES! Magazine.
Peter has been active in Transition Pasadena for several years, and his book was developed through discussions, community actions, and presentations with them. Peter is now interested in exploring fossil fuel reduction in community settings, such as local church and neighborhood groups. In this TeleSeminar, in addition to introducing his book, Peter will explore how his work could complement Transition Streets. Please join us!
Consider picking up a copy of Peter's book. Peter is donating all profits to groups connecting individual agency to collective climate action.


A luxurious Organic Moisturising Lotion base blended using 99% natural ingredients. This product is formulated to accept the addition of essential oils.
Base Formula's Organic Moisturising Lotion is formulated with 99% natural ingredients and 90.85% organically certified ingredients. It is suitable for use as a facial, hand and bodycare lotion and is ideal for dry and sensitive skin. It has been formulated with high quality ingredients to gently replenish moisture and oil lost by the skin and can be used therapeutically to soothe and relieve the discomfort caused by rough and dry skin. The lotion has been formulated to be used alone or blended with essential oils, carrier oils and other cosmetic ingredients to produce a range of tailor made skincare products.
Professional Tip: The experts at Base Formula recommend that no more than 1% or 10ml of any essential oil and 2% (20ml) of liquid ingredients such as aloe vera or hydrolats be added to a litre of our cosmetic bases. Exceeding this amount may destabilise the blend. Over time this may cause the bases to lose some of their luxurious properties including viscosity and silky texture.


text:
Gambino mobster wins big at Atlantic City and Las Vegas poker tournaments
October 26, 2012 Dapper_Don
A Gambino wiseguy charged with conspiring to distribute marijuana went after a different kind of pot in Atlantic City and Las Vegas, law enforcement sources say.
While Drug Enforcement Administration agents probed his alleged buying and selling of several pounds of marijuana, Lee Fama, 44, won thousands in casino poker tournaments.
Fama hauled in $22,000 in three poker tourneys at the Borgata in Atlantic City this year, and $666 at a tourney during the summer at Las Vegas's Rio hotel, poker Web sites report.
Fama, of Bensonhurst, Brooklyn, is not on the list of 266 people excluded from New Jersey casinos. Nor is he on the much shorter list of people excluded from Nevada casinos.
Casino exclusion rules differ in the two states. New Jersey can bar "a career or professional offender."
Nevada's rules are more lenient — a spokesman for the Nevada Gaming Control Board said the state probably wouldn't bar a convicted criminal from gambling, but might bar that person from operating or working in a casino.
Fama was arrested on the conspiracy charge Tuesday, and is expected to be released on bail later today.
In August 2011, Fama was caught on wiretaps buying two pounds of "high quality" marijuana from the cooperating witness, says a DEA affidavit.
In November, Fama sold another pound of marijuana to the same witness, the affidavit said.
His lawyer, Fred Sosinsky, said he's going to fight the drug charges, and denied his client is a made member of the Gambino family. A law-enforcement source say he is.
http://www.nypost.com/p/news/local/tourneys_vegas_gambino_mobster_under_tIu68AqkG5fYQWmCjo1FEI
Gambino Crime Family / Lee Fama

redpajama_set_name: RedPajamaCommonCrawl


text:
2018 MIAA Softball State Championship
Taunton Softball 2018 MIAA State Champions with a 13-2 Win over Wachusett.
Tigers End Long Wait, Bring Home State Championship
WORCESTER, Mass. – In the bottom of the fourth inning, Wachusett started to hit the ball hard, scoring a run on back-to-back doubles and cutting into the Taunton lead. Although the Tigers still led by three runs, momentum appeared to be shifting to the Central champs and Taunton went to the plate looking for an instant response to take back control.
With one run home, the Tigers had the bases loaded with two outs and up stepped junior Alexa White. Wachusett starter Courtney Lanpher got two strikes on Taunton's second baseman, but then White turned on a pitch, lifting it to deep left. The Taunton bench was almost holding its breath watching the ball sail towards the fence and almost exploded when the ball cleared it.
Click here for a photo gallery from this game.
White's grand slam broke Saturday night's Div. 1 state championship game wide open, lifting Taunton (26-2) to a convincing 13-2 win over the Mountaineers and ending an 11-year wait for the program's fifth state title...................................
Taunton Seniors Cap Career With State Title
WORCESTER, Mass. – Quiet, inspirational, and crazy — the perfect combination.
That's how Taunton head coach Dave Lewry described each of his three senior leaders — Liana Duarte, Kya Enos, and Mackenzie Handrahan.
All three of those Tigers played a vital role this season, and unsurprisingly, were a big part of Taunton winning the D1 State Championship with a 13-2 decision over previously unbeaten Wachusett.
"They are the nucleus of the team," Lewry said. "What a great combination of seniors. I never had to use the word 'senioritis' this season. Usually, we have to talk about it with them being out of school at the end of the season. But we never had to mention it all year and that's a tremendous credit to those three...............................


text:
<!doctype html>
<script src="/resources/testharness.js"></script>
<script src="/resources/testharnessreport.js"></script>
<script>
promise_test(async t => {
  let abort_signal;
  let onabort_called = false;
  appHistory.onnavigate = t.step_func(e => {
    abort_signal = e.signal;
    abort_signal.onabort = () => onabort_called = true;
    e.transitionWhile(Promise.resolve());
  });

  await appHistory.navigate("?1").finished;
  assert_false(abort_signal.aborted);
  assert_false(onabort_called);
}, "AppHistoryNavigateEvent's transitionWhile() does not signal AppHistoryNavigateEvent's signal");
</script>

redpajama_set_name: RedPajamaGithub


Brother UK Ltd (BUK) are the primary governing authority for the Partner Portal and any user by activating this portal gives their implied authority to be bound by these terms and conditions.
The Partner Portal is owned and sponsored by BUK who have full rights in all copyright and intellectual property associated with the portal. Any user has viewing rights only and materials (unless explicitly instructed on the portal to do so) may not be downloaded or copied without the express written permission of BUK. Unauthorised use of any materials will be subject to the provisions of clause 3b.
To become a Registered Member of the Partner Portal you must be an employee of a Brother Authorised Reseller (either retail or e-selling), before participation. Details of participating resellers can be obtained direct from BUK.
b) By participating the Registered Member warrants full responsibility for their own legal and ethical behaviour. Any actions that Brother consider inappropriate will result in that Registered Member and it's employer being forbidden from further participation in any BUK initiatives; the Registered Member by accepting these terms and conditions of participation accepts Brother's right to take any appropriate legal action for theft or fraud.
BUK reserves the right at any time to change the Partner Portal without any liability to any user under any circumstances.
If the Registered Member ceases employment with the Brother Authorised Reseller, access to Partner Portal will not be permitted. Any person who continues access to BrotherZone after their employment has been terminated will be subject to action as stated in clause 3b.
In the event of any disputes the final decision will be solely at the discretion of BUK.
Personal data which the Registered Member provides will be used by BUK in accordance with the Data Protection Act 1998; this does not affect your statutory rights.


Artificial intelligence has been around for half a century, but it has only recently become a real player in the healthcare space. Harjinder Sandhu, the CEO and founder of SayKara, says that's all down to data.
Sandhu was joined on stage by Michael Calhoun, CEO and co-founder of medical imaging analysis company Mindshare Medical, and Anisha Sood, a partner at Echo Health Ventures. Sandhu said the data barrier for AI has tumbled down in the past few years.
"Now every device is capturing data. Everybody's using electronic medical records, all the data is getting captured in electronic form," he said.
And the applications for AI could be just as varied and far-reaching as the vast swaths of data involved in the healthcare system, everything from predicting which patients will get sick to assisting doctors in the operating room.
But Sandhu and Calhoun agreed that, for now, AI isn't anywhere near replacing doctors or other humans in the healthcare world.
"It may happen eventually, for a lot of use cases. It's not happening right now. So the way to go about it, if you're a startup in this space, is find ways to add value to the doctors today," Sandhu said.
Watch the full panel discussion of AI's role in healthcare above and read all our coverage from the 2018 GeekWire Summit.


Mercedes-Benz Actros Gigaspace refrigerated box Trailer "Herpa Weltgeschichte Nr. 5.1, Die Dinosaurier"
After more than 13 years and more than 1.8 million kilometers, the legendary Dino truck from the Herpa world history series received a new tractor unit. The Herpa model comes as a unique and limited edition.


text:
ArticleBio
Younis Mahmoud
Younis Mahmoud Biography
Biography by Clarence Published on 18 Aug, 2016 Updated on 18 Aug, 2016
Facts of Younis Mahmoud
1979 , February-18
Birth Nation:
6 Feet 1 Inch
Younis Mahmoud Khalaf
Birth Place/City
Working For
Al-Talaba S
Married to
View more / View Less Facts of Younis Mahmoud
Younis Mahmoud Khalaf, currently a free agent striker, was born on February 18, 1979, in Dibis, Kirkuk, Iraq. As a result, he is Iraqi by nationality and ethnicity.
Younis Mahmoud started his career as a shooting guard for Kahraba Al-Dibis basketball team. He decided to play football when the team's football coach Muwafag Nouraddin explain the profit of playing football in Irag as there is more money involved. So, he switches to Kahraba Al-Dibis football team and regularly played for the side.
Younis Mahmoud, in 1999, left Kahraba Al-Dibis and moved to Kirkuk FC. There, he scored his first goal against Al-Kadhimiya in Round 8 of the 1999-2000 Iraqi Premier League. Later, he was downgraded to lower division for a while but managed to get on the top division before the end of the season. During his stay at Kirkuk Mahmoud was discovered by Iraqi's largest clubs both Al-Shorta and Al-Talaba and was invited to go on trial by Al-Shorta but was shortly traded to Al-Talaba.
During Younis Mahmoud's first match for Al-Talaba at 2001-2002 Iraqi Elite League, he scored hat-trick within 10 minutes against his former club Kirkuk and won the match by 8-0. During his career at Al-Talaba, he won his first ever trophy after they defeated Al-Zawraa with a goal score 2-1 in Baghdad Cup 2001. During this season Mahmoud scored 13 goals for Al-Talaba. The team also won the Iraq FA Cup, in which Mahmoud scored 6 goals and assisted the winning goal in the final match. The same season he scored 4 goals in the Umm Al Ma'arik Cup, Al-Talab managed to get to final but ironically lose against Al-Shorta.
Younis Mahmoud signed a contract with Al-Khor of Qatar Stars League in 2004. During the 2004-05 Qatar Stars League, Mahmoud scored 19 goals and the team finished in 3rd place. In 2006, he was transferred to Al-Gharafa and played here for straight five years. In 2006-07 Qatar Stars League, Al-Gharafa was the runner-up and Mahmoud was honored for the first time with the Golden Boot after he scored 19 goals in the league. He scored 2 goals in the final, both in extra time and won the 2007 Sheikh Jassim Cup. Mahmoud scored 16 goals and won the Qatar Stars League title for the first time in his career in the 2007-08 season.
Younis Mahmoud, during the 2008-09 season,played for Al-Arabi club on loan and won the Sheikh Jassim Cup . With this club, he played the 2008-09 Qatar Stars League and scored 2 goals in the league. In the middle of the 2008-09 season, he returned to Al-Gharafa and the team won the Qatar Stars League, in which Mahmoud scored 1 goal for the team. In 2009, for the second time, he earned the Golden Boot for scoring 21 goals for Al-Gharafa and won the Emir of Qatar Cup for the third time.
In 2011, Younis Mahmoud moved to Al-Wakrah and played in the 2011-12 Qatar Stars League, in which the team finishes at 7th place and Mahmoud scored 8 goals. During his career in Al-Wakrah, he managed to win the 2011 Qatar Stars Cup in which he scored 1 goal in the semi-final. Soon after Al-Wakrah team were knocked out from 2012-13 Qatari Stars Cup, Mahmoud signed with Al-Saad on 8 February 2013 and stayed here till the end of 2012-13 season. During his stayed on this team, he managed to win the league title for the fourth time and was equalled Nasser Kamile's record goal in the Qatar Crown Prince Cup on 26 April of this season.
On September 16, 2013, Mahmoud signed a four-month contract with Al-Ahli club of Saudi Professional League. He ended the contract before two weeks of the official end. He then again returned to Irag after the 2015 AFC Asian Cup and played for Erbil and played only in AFC Cup games , because according to him he did not want to compete against his former team Al-Talaba in the Iraqi Premier League. However, on 12 September 2015, Mahmoud returned to his previous club Al-Talaba and was honored with a captain's armband and jersey number 9. At the end of 2015-16 season, Mahmoud said that he will be taking a break from football for the coming season but denied that he had retired.
Mahmoud has married three times but his wives named is not disclosed by any media. However, his net worth and salary are not yet disclosed, he is living in millions. Although he is rich, Mahmoud lives a simple life.
Al-Talaba SC
Footaballer
Scarlett Byrne 1990 -
11 Mar, 2020 Biography by Arsad
Mikey Foster 1987 -
05 Mar, 2020 Biography by Shally
Ewan McGregor 1971 -
04 Mar, 2020 Biography by Basanta
Fisher Stevens
Emma Stone Married To Fiance Dave McCary? Pregnancy Rumors?
23 Sep, 2020 News by Luna
Carly Steel Gets Engaged To Beau Jacob Andreou In Front Of Eiffel Tower
02 Dec, 2019 News by Joey Jordan
Lyssa Chapman's Relationship With Girlfriend Leiana Evensen
16 Oct, 2019 News by Joey Jordan
Why Lyssa Chapman "Baby Lyssa" Divorced Her First Husband Brahman Galanti
Trending Biography
Natalie Hanby 1990 -
Jen Wilson 1982 -
Lauren Hammersley 1971 -
Daniel Bellomy
Margot Rooker
Jermaine Jakes 1980 -
Al Pacino: Personal And Professional Life-Net Worth, Age, Height, Movies & TV Shows, Wife & Children
20 Oct, 2020 News by Luna
How Much Is Peter Berg's Net Worth? His Sources Of Income, Salary, And Career
05 Aug, 2020 News by Luna
Legendary Comedian And Actor Jerry Stiller Dies At 92: His Career, Personal Life, And Net Worth
17 May, 2020 News by Dabria
Jennifer Carpenter 1979 -
18 Jan, 2021 Biography by Shrijan
Susan Waren
18 Jan, 2021 Biography by Sammy
Michael C. Hall 1971 -
Melanie Bright
18 Jan, 2021 Biography by Basanta
Erin Payne 1971 -
Sam Caird
Using Articlebio.com you are agree in our Terms of Use. Contents published by users are under Creative Commons License.


A big congratulations to the Welcome Baby program at Miller Children's & Women's Hospital Long Beach for receiving a special proclamation from Long Beach Mayor Robert Garcia!
Parenting education, breastfeeding support, information on infant development and bonding and attachment, plus connections to community services and resources — those are some of the key benefits Welcome Baby provides. And, as Silvia and Abraham explain in this video, those are all the things these parents happily received from the program.


text:
Villa Santo Sospir, photograph by Loïc Thébaud, loicthebaud.com
The villa where Jean Cocteau communed with the Gods
Our new book Interiors: The Greatest Rooms of the Century features a Godly work created by the poet and filmmaker
Our new book, Interiors: The Greatest Rooms of the Century, is a beautiful, global celebration of residential interior design and decorating. With 400 rooms organised by designer, the book goes beyond decorators, interior pros and architects, and highlights also the exquisite homes created by artists, style icons, fashion designers, and film actors, all of whom have made a unique contribution to the way we arrange and enjoy our homes.
Some of these creations were made relatively quickly, as formal commissions, yet others were the fruits of long-standing friendships, and say as much about the room's creator as the home's owner.
Jean Cocteau's work at his friend Francine Weisweiller's house, Villa Santo Sospir, on France's Mediterranean coast, certainly falls into the second category.
"In 1949 Jean Cocteau was invited by art patron Francine Weisweiller to stay at her home on the French Riviera," our new book reveals. "Wesweiller had met Cocteau on the set of Les Enfants Terribles and persuaded her husband, Alec, to invest in the film.
"Cocteau came for a week, but asked his hostess if he could paint a mural in the living room, and, over twelve years, Cocteau decorated most of the house in murals and mosaics, sometimes inviting other artists, such as Pablo Picasso, to join in. In characteristic faux-naïf style, Cocteau painted Apollo, the Greek god of sun, philosophy, poetry, and music, with piercing eyes and blazing hair."
The midnight blue edition of Interiors
The poet, artist and film director's images might have been applied in a knowingly naïve style, but both the images, and the act of painting meant a lot to him.
"Myth for Cocteau was a vehicle for exploring the role of the contemporary poet, and he evoked Greek gods, goddesses, and goats around the villa, transforming it from bourgeois beach house into a salon for elevated thinking, creativity, and sensuous lounging enjoyed by other guests, including Marlene Dietrich and composer Francis Poulenc," writes Giles. "For Cocteau, it was all about self-expression: "When I was working at Santo Sospir, I became myself a wall and these walls spoke for me." The walls, he said, were like a skin that he tattooed with brushstrokes, creating a mark of enchantment."
For more on interior design from the beginning of the 20th century up until the present day, order a copy of Interiors: The Greatest Rooms of the Century here; you can even choose your own version. There are four to choose from: saffron yellow, platinum gray, merlot red and midnight blue. Pick the right one for your own interior here.


text:
Barr & Stroud Limited was a pioneering Glasgow optical engineering firm. They played a leading role in the development of modern optics, including rangefinders, for the Royal Navy and for other branches of British Armed Forces during the 20th century. There was a non-military arm of the company which made medical equipment, like photocoagulators and electronic filters, some of which were used by the BBC. The company and its intellectual property passed through Pilkington group to Thales Optronics. The Barr and Stroud name was sold on to an importer of optical equipment, who use the trademarked name for a line of binoculars and similar instruments.

History

Archibald Barr and William Stroud had been associated from as early as 1888 when the two men were professors of, respectively, engineering and physics at the Yorkshire College (now the University of Leeds).  In 1891, they were approached by the Admiralty to submit a design for a short-base rangefinder for trial.  By this time, Barr had returned to Scotland and taken the post of Regius Professor of Civil Engineering and Mechanics at the University of Glasgow. Although apart, Barr and Stroud kept in close touch and in 1892 they were awarded with a contract for six of their rangefinders.

In 1895, Barr & Stroud's Patents Ltd was renting workshop space near the university, at 250 Byres Road, Glasgow, but demand for the product soon necessitated a move to larger premises in Ashton Lane, Glasgow.  By 1904, 100 men were working for the company in a new purpose-built factory in Anniesland, Glasgow designed by Campbell Douglas.  Shortly thereafter, in 1909, Stroud resigned his chair at the University of Leeds and moved to Glasgow to work for the company full-time.  Barr, in spite of a distinguished teaching career at Glasgow University, followed his example in 1913.  Together they formed Barr & Stroud Ltd. that year.

In 1914, they began extensions to the Anniesland works in order to meet the sharp increase in demand for their rangefinders that followed on the outbreak of the First World War. The war years saw the development of other products, including a torpedo depth recorder, a periscope rangefinder, fire-control systems and a dome sight for aircraft. During World War I there was a problem with the supply of binoculars to the armed forces and apart from the British makers, binoculars were bought in from various sources. So it was only natural that the military would be looking to rectify this problem. At the time Zeiss was probably the leading makers and the military thought this may cause a problem in the future.

It was c.1919 when the company started producing their first binoculars which were supplied to the British Navy and from then on the company continued to operate independently until c.1977 when they were taken over by the  Pilkington Group.  In 1992 operations moved from the original factory in Anniesland to a new plant in Linthouse on the site of the former Alexander Stephen and Sons shipyard.  In 2000 the company became a subsidiary of the French company, Thales Group, and in 2001 Barr & Stroud Ltd became Thales Optronics Ltd.

The Barr and Stroud brand name was then bought by Eastleigh-based Optical Distribution Services Ltd, who re-registered as Barr and Stroud Ltd in 2008. The new company has developed a new range of binoculars and telescopes.
The new range of Barr & Stroud binoculars are currently made in China (Nov. 2011) and distributed by Optical Vision Ltd.

Engines
In the 1920s Barr & Stroud Ltd started offering sleeve valve motorcycle engines based on a design by Peter Burt and L.J. McCollum. In a half-page advert in Motor Cycle magazine in 1922 they encouraged readers to contact them at Anniesland, Glasgow, for more information, or contact the manufacturers currently offering their 350cc single cylinder engine, naming Beardmore-Precision, Royal Scot, Edmund and Diamond motorcycles. In 1923 a new motorcycle manufacturer, Grindlay-Peerless offered a 999cc V-twin Barr & Stroud engined motorcycle, later adding a 500cc and 350cc single. They also offered J.A.P. engines, and it appears the racing success of the J.A.P. engines made them the preferred choice, and so the sleeve valve engines were dropped in 1927.

Computers
Barr and Stroud constructed the first computer to be built in Scotland, the pioneering SOLIDAC minicomputer for the University of Glasgow, assembled between 1958 and 1963 as an attempt to expand into electronics.

Notes

External links

 Short History of Barr and Stroud
 Post WW1 trials of Barr & Stroud rangefinders
 Current website
 Details of their WW1 era rangefinders
 Details of their WW1 era fire control instruments

Naval artillery
Anti-aircraft weapons
Military computers
Royal Navy
Firearm sights
Armoured fighting vehicle vision and sighting equipment
Manufacturing companies based in Glasgow
Engineering companies of the United Kingdom
Optics manufacturing companies

fi:Stereoetäisyysmittari


Through the York County Beekeepers Association, students will also have the opportunity to visit local hives for observation and hands-on experience. Beginner Bee School instructor is Larry Peiffer, Master Beekeeper and former Vice President of the Maine State Beekeepers Association.
Course Fee: $95, includes beginning beekeeping textbooks, other materials, and a 1-year membership in the York County Beekeepers Association. A reduced fee of $45 applies to one of a pair of students who choose to share the text and materials ($140/pair).
Sponsored by the Maine State Beekeepers Association and University of Maine Cooperative Extension York County.
For more information, contact UMaine Extension in York County at 1.800.287.1535 (in Maine), 207.324.2814 or (TDD) 1.800.287.8957 or e-mail elizabeth.clock@maine.edu.
If you are a person with a disability and need an accommodation to participate in this program, please call Frank Wertheim at 800.287.1535 to discuss your needs. Receiving requests for accommodations at least 7 days before the program provides a reasonable amount of time to meet the request; however, all requests will be considered.


1. How do I purchase a print from your site?
To purchase a print you must browse to the category of your choice from the navigation menu in the header at the top of the pages. From there you will look over the images until you find one or two or…that you want to purchase. Next you will choose a print size (mandatory) from the dropdown menu under the print. If all you want is the print you will make sure the mat and frame options are set to none and click add to cart. To purchase a print with a mat, choose the print size then the print size mat in either single or double from the drop down menu options and click add to cart. If after choosing your print size and type of mat, you also want to add a frame to your print choose the frame size and colour then click add to cart.
Keep in mind that if you choose a 5×7 print you will need a 8×10 mat and frame, if you choose a 8.5×11 print you will need a 11×14 mat and frame, if you choose a 11×14 print you will need a 16×20 mat and frame, if you choose a 13×19 print you will need a 18×24 mat and print, if you choose a 16×20 print you will need a 20×24 mat and print. If you want a 24×36 print you will have to order a custom mat and frame as we don't stock precut mats or frames for that larger size. To view our mat colors click on this link. Here are images of our picture frames.
Print sizes larger than 11×14 will be shipped without glass as damage is very common with sizes over 11×14. Contact your local glass shop for a price for your print/mat/frame.
After making your choice and adding to cart, your price (in blue) will change from my base price (15.00) to the price for your selection. Your shopping cart will show up at the bottom of the page with your order details or in the side bar if you choose to keep browsing my website. Once you decide to check out please follow the prompts after making payment choice. Once your funds have cleared, your order will be shipped.
2. Do you ship glass with your framed prints?
Only up to print sizes of 11×14 everything above this size will not be shipped because of consistent damage occurring during transit. Check with your local glass shop for prices on frame size you are looking at.


text:
The iSelect Approach for Accredited Investors
There is a disconnect in the startup funding ecosystem. Most ordinary Americans only learn about exciting new venture opportunities after it is too late to invest.
Perhaps one of the reasons for this is that investing in venture capital is hard. It takes time and expertise to adequately source, identify, and monitor a portfolio of venture investments. Diversifying across several deals typically requires $1.0 million or more of investable capital. Monitoring venture investments requires constant vigilance. The task is unmanageable for most people despite the potential returns.
As the first build-your-own venture fund, iSelect solves this problem by managing the diligence burden for investors, providing diversification by allowing minimum investments in a single company for as low as $5,000 and continually monitoring investments after.
A Portfolio of Venture Investments
Strict Selection Criteria
Informed Decisions
iSelect is the first investment vehicle that affords accredited investors access to an evergreen portfolio of venture investment opportunities that have been thoroughly vetted by a committee of successful entrepreneurs, domain experts, and professional venture investors. We allow accredited investors to realize the potential to earn venture capital returns via a product that is diversified and risk-managed.
Each company in the iSelect portfolio is subjected to an extensive three-tiered due diligence regime, including review by a FINRA registered broker-dealer. This takes the hard work of vetting potential investments off of fund investors. Once selected for review by our Venture Team, potential companies are studied by an industry-specific Selection Committee that's composed of independent entrepreneurs, business executives, world-renowned scientists and researchers, and accomplished venture investors in market niches including Healthcare, Agriculture, Energy, Technology and more, before being submitted to iSelect's Investment Committee for final review.
As part of the diligence process, iSelect's Selection Committees dig into the fundamentals underlying each company: (1) Is the venture commercially ready for business? (2) Does the venture have the ability to create the barriers of entry necessary for a sustainable competitive advantage? (3) Is the management team in place capable of growing the business responsibly? (4) Has the company taken the appropriate steps to identify and address the risks of the venture to the extent possible for the relevant stage of the company. Companies that fall short on any of these criteria are not accepted to the iSelect platform.
Our investors have access to all of this research as well, in the form of a summarized and packaged diligence file that is made available to investors and their financial advisors for each venture investment in the iSelect portfolio. Post-investment, companies provide quarterly reports so that investors can monitor the progress of their portfolio, all in an effort to help accredited investors to more easily evaluate, understand, and manage their personalized series of early-stage venture securities.
iSelect works with investors and their financial advisors to help create a diversified portfolio of venture investments, and is an effective investment strategy for growth, retirement and wealth preservation. By identifying and investing in only the most promising early-stage companies, we are able to not only support the work of smart, innovative entrepreneurs, but also reap the reward of potential investment returns over time.
Must be an accredited investor who is comfortable with a hold period that can be as long as 5-7 years or more.
The minimum investment in iSelect is $50k, no less than $5k per company.
Each investor, in conjunction with their financial advisor, can select which companies to invest in, and how much to invest (subject to the $5,000 per company minimum).
iSelect investments can be held in a retirement account, a normal brokerage account, or a trust, and the Fund is structured for investment by individual investors, family offices, retirement plans, and institutional investors.


Is Android the next Symbian?
Home Android Is Android the next Symbian?
The dust has settled and the time has come for the reality to kick in. Whatever threats we have seen with Android have now been quelled, but new threats are emerging. Running a business is like playing chess. You move a piece to support another piece, only to find that you have weakened the position of yet another piece.
I just learnt that Symbian was started as a open handset alliance, much like Android. In case of Symbian, the shares were not equally distributed and Nokia had majority share. This made other members uncomfortable who sold their shares and left. After a while, Nokia was the sole Licensee and the licensor. There were 11 participants in the Symbian alliance, which also included Motorola.
Android's open handset alliance has 84 members in it. The alliance has mobile operators, software companies, commercialization companies, semiconductor companies and handset manufacturers. With so much of diversity and so many vested interests, it is hard to think that Android might go the Symbian way. But what's stopping Android from going the Symbian way?
How much of openness will Google let Android have?
Sriram Vadlamani is the Editor and co-founder of TheGadgetFan and a columnist at Asian Correspondent. You can follow him on twitter @indianomics (Email: [email protected] | Bio). Add him to your Google circle.
Samsung Has a New Problem : S7 Camera Lens Shattering Without Any Impact!
Best 8 Custom ROMs for Your Android Device!
How to Unlock the True Potential of Your Android Smartphone, Without Rooting!