%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%
%
% e-mail: pjarz@ippt.gov.pl
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\documentclass[12pt]{article}
%\usepackage[polish]{babel}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
%\usepackage[dvips]{graphics}
\usepackage{graphicx}
%\usepackage[active]{srcltx}
\usepackage{epstopdf}
\usepackage{amsmath}
\usepackage{amsthm}
\usepackage[tableposition=top]{caption}
\usepackage{subcaption}
\usepackage{listings}
\usepackage{multicol}
\usepackage{fullpage}
\input{definitions.tex}
%
%% set the level of table of contents
%% 0 - only chapters,
%% 1,2,3 -chapters and sections, 4 - with paragraphs
%\setcounter{tocdepth}{4}
%
%\pagestyle{headings}
%\textheight 23.cm
%\textwidth 15.8cm
%\footskip 22.8pt
%\topmargin -1.0cm
%\oddsidemargin   0.0cm
%
%
%\setcounter{tocdepth}{4}
%\usepackage[hypertex,colorlinks=true,pagebackref]{hyperref}

% kwisn
%\usepackage[colorlinks,bookmarks,pagebackref]{hyperref}  % hyperlinks, to facilitate navigation

\begin{document}
\lstset{language=Fortran,frame=tb}
%
%\input{definitions.tex}
%\def\sh{_{\mathrm{sh}}}
%\def\ext{_{\mathrm{ext}}}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\normalsize

\begin{center}
\bigskip
\bigskip
%Version of 21.10.2016

\bigskip
\bigskip



%{\Large  \bf
%Partial Factorization Applied to Reduction of Finite Element  Matrices}

%{\Large  \bf
%Application of Partial Factorization to Reduction of FE \\ [0.5 ex]
%Matrices and to DD Equation Solvers}

{\Large  \bf
Evaluation of Partial Factorization for Reduction \\ [0.5 ex]
of Finite Element Matrices}


\vspace{5 mm} Pawel Jarzebski\textsuperscript{1}, \
Krzysztof Wisniewski\textsuperscript{1,2}

\textsuperscript{1}IFTR, Polish Academy of Sciences, Pawinskiego 5B,
02-106 Warszawa, Polska

\textsuperscript{2} e-mail: kwisn@ippt.pan.pl

% corresponding author


\end{center}


%-------------------------------------------------------------------
%\clearpage
\noindent
In this paper, we present the concept of \emph{Partial Factorization} \cite{Petra2014} and discuss
its possible applications to the Finite Element method. We consider:
(1) reduction of the element tangent matrix, which is particularly important for mixed/enhanced elements and
(2) reduction of the sub-domain matrices of
the Domain Decomposition (DD) equation solvers run either sequentially on a single machine or
in parallel on a cluster of computers.
%(2) reduction of the memory usage of the equation solver on a single machine, and
%(3) the domain decomposition equation solver run on a cluster.
We demonstrate that \emph{Partial Factorization} can be beneficial for these applications.

\bigskip
\noindent
{\bf Keywords:} \
 multi-scale models of multi-layer shells;
 mixed/enhanced finite elements; parallel computing; Domain Decomposition,  solvers;

%\tableofcontents


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\section{ Introduction}

In multi-scale computations of multi-layer shells,
the finite element models can have a size of millions of  degrees of freedom
and must be repeated many times e.g. when a  database of material properties is built
for the Design Sensitivity Analysis.
Therefore, it is necessary to improve efficiency of an FE code,
especially using concurrent capabilities of multi-core processors and multi-machine clusters.


\medskip
The minimum is to use one of the available parallel solvers,
such as \textsf{\small PARDISO}, \textsf{\small MUMPS},
\textsf{\small PaStiX},  \textsf{\small HSL} and others.
The next step is a parallelization of computations of the FE matrices and vectors, e.g. using OpenMP. \
In the third, the so-called hybrid approach, additionally
the Domain Decomposition is performed, e.g. by \textsf{\small METIS},
and computations for sub-domains are scattered over a cluster of computers,
which requires, e.g. MPI. \

\medskip
This, however, requires several changes in the FE code and  advanced computational techniques.
In  \cite{Jarzebski2015}, we describe a parallelization of the loop over elements of FEAP while
in \cite{JarzebskiWisniewski2016}, a performance of the ompFEAP and the parallel solver HSL MA86
was established for two various composite materials.

\medskip
An objective of this paper is to describe possible applications of \emph{Partial Factorization} (PF)
to FE computations. The PF was proposed in  \cite{Petra2014} \ for the stochastic optimization problems.
We shall establish performance of this technique in two other applications:
(1) a reduction of the element tangent matrix,
which is particularly important for mixed/enhanced elements using a large number of additional parameters,  and
(2) a reduction of the sub-domain matrices to the ones expressed in terms of the interface variables
for the Domain Decomposition (DD) equation solvers, run either sequentially on a single machine or
in parallel on a cluster of computers.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\section{Partial Factorization}
\label{section:partialFactorization}

\subsection{Schur complement for  Domain Decomposition}
\label{subsubsec:SchurParallel}



It is well-known that the FE method yields the system of equations \ $\bfK \bfu= \bff$ \ of a banded structure.
If we compute these equations using the Domain Decomposition method,
i.e.  each sub-domain on a different node of a cluster of computers, then, for each sub-domain \ $i \ (i=1,...n) $, \
we have
\eqb
\begin{bmatrix}
       \bfK_i  & \bfB_{i}^{T}\\
       \bfB_{i}&  \bfC_i
     \end{bmatrix}
   \begin{bmatrix}
       \bfu_i   \\
       \bfu_i^I
     \end{bmatrix}
=        \begin{bmatrix}
       \bff_i   \\
       \bff_i^I
     \end{bmatrix},  % \qquad  i=1,...n,
\label{Eq:schurSystem}
\eqe
where we distinguish between the interface variables \ $ \bfu_i^I $ \  and the domain variables \ $ \bfu_i $. \
Besides, we assumed that the matrix is symmetric.
The domain variables must be eliminated within each domain, which yields
\eqb
\bfS_i \, \bfu_i^I = \bfs_i,
\eqe
where
\eqb
\bfS_i \doteq \bfC_i - \bfB_i \bfK_{i}^{-1} \bfB_{i}^{T}, \qquad\qquad
\bfs_i \doteq \bff_i^I-\bfB_i \bfK_{i}^{-1} \bff_{i}.
\label{Eq:ischurSystem}
\eqe
Note that the matrix \ $ \bfS_i $ \ is the Schur complement of \ $ \bfK_{i}$. \
Next step is performed for all domains taken together, for which
 we have to assemble and solve the system of equations for the interface variables,
\eqb
\bfS \, \bfu^I = \bfs,
\label{Eq:schurSystemEquationSolve}
\eqe
where
\eqb
\bfS   \doteq  \sum\limits_{i=1}^n \bfS_i  \qquad
\bfs   \doteq  \sum\limits_{i=1}^n \bfs_i  \qquad
\bfu^I \doteq  \sum\limits_{i=1}^n \bfu_i^I.
\label{Eq:schurSystemKrhs}
\eqe
Having obtained the interface variables \ $ \bfu^I $, \  we compute \ $ \bfu_i $ \ for each domain separately,
which ends the solution process.

\medskip
Note that computations of \ $ \bfS_i $, \ $ \bfs_i $ \ and \ $ \bfu_i $ \
can be performed in parallel in sub-domains, which can significantly shorten the solution time.
Most time-consuming  is computation of the Schur complement \ $ \bfS_i $ \
in sub-domains, especially of the  product \ $ \bfK_{i}^{-1} \bfB_{i}^{T}$. \
Classically this product is obtained in two steps:
(1) factorization of \ $\bfK_{i}$, \ and
(2) one back-substitution for each column of \ $\bfB^T$. \
This, however, is very inefficient and e.g.  the speedup on 12 cores is only about 3, see  \cite{JarzebskiWisniewski2016b}. \ Much faster is the method presented in the next section.

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

\subsection{Partial factorization and Schur complement}
\label{subsection:partialFactorization}

The \ LU \ decomposition algorithm for factorization of a matrix has acceptable scalability on
multicore machines and this feature is preserved when we use \emph{Partial Factorization}
to compute the Schur complement.
The PF method was proposed in  \cite{Petra2014} \ for the stochastic optimization problems;
in the current paper we propose its other applications.

\medskip
The \ LU \ decomposition of the matrix of eq.~(\ref{Eq:schurSystem}) can be written as follows:
\eqb
% \bar{K_i} =
\begin{bmatrix}
       \bfK_i   & \bfB_{i}^{T}\\
       \bfB_{i} &  \bfC_i
     \end{bmatrix}
= \begin{bmatrix}
       \bfL_{11}  & \bf0\\
       \bfL_{21}  &  \bfL_{22}
     \end{bmatrix}
     \begin{bmatrix}
       \bfU_{11}  & \bfU_{12}\\
       \bf0       &   \bfU_{22}
     \end{bmatrix}
=
     \begin{bmatrix}
       \bfL_{11}\bfU_{11}  & \bfL_{11}\bfU_{12}\\
       \bfL_{21}\bfU_{11}  & \bfL_{21}\bfU_{12}+\bfL_{22}\bfU_{22}
     \end{bmatrix}.
\label{Eq:augmentedSchurSystem}
\eqe
where \ $\bfL_{ij}$ \ and \ $ \bfU_{ij} $ \ are the block factors
of the lower and upper triangular parts, respectively.
Noting that
\ $ \bfL_{21}\bfU_{11} = \bfB_i $ \ and \  $ \bfL_{11}\bfU_{12} = \bfB^T_i $, \
we can calculate \ $\bfL_{21}$ \ and \ $\bfU_{12}$ \ first.  Then, % we obtain
\ $ \bfL_{21}\bfU_{12} = \bfB_i \bfU^{-1}_{11}\bfL^{-1}_{11}\bfB^{T}_i = \bfB_i \bfK^{-1}_i \bfB^{T}_i $ \
by \ $ \bfL_{11} \bfU_{11} = \bfK_i $. \
We note that
\eqb
\bfC_i - \bfL_{21}\bfU_{12} =  \bfC_i - \bfB_i \bfK^{-1}_i \bfB^{T}_i = \bfS_i,
\eqe
i.e. the Schur complement of eq.~(\ref{Eq:ischurSystem}).
Note that we do not need to compute  \ $\bfL_{22} $ \ and \ $  \bfU_{22}$ \ to obtain the Schur complement,
which justifies the term ``partial'' factorization.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\vspace{-5 mm}
\section{Reduction of element's matrix}
\label{section:matrixReduction}
\vspace{-2 mm}
\subsection{Set of equations for problems with additional variables}

Mixed/enhanced finite elements involve additional local parameters
which must be condensed out % on the element's level
to reduce the size of a tangent element matrix to the standard one.
The number of element parameters can be quite large,
sometimes exceeding the number of nodal variables associated with the element,
and one of the approaches is to group and condense out all parameters together,
which amounts to computation of the Schur complement.

\medskip
For the considered class of finite elements, the governing functional \ $F$ \ depends on
the nodal  displacements (and rotational parameters) designated  as  \ $\bfu_I$ \
and the elemental multipliers \ $\bfq$. \
In general, \ $ \bfq $ \ are the Lagrange multipliers (including the multipliers of stress modes)
and the multipliers of enhancing kinematical modes, used in the enhanced strain methods, such as:
Incompatible Displacements or EADG or EAS method.

\medskip
For kinematically non-linear problems,
the stationarity condition of  \ $ F(\bfu_I,\bfq) $ \
yields a system of equilibrium  equations for an element,
\ $ \bfr_u(\bfu_I,\bfq)   = {\bf0} $ \ and \ $ \bfr_q(\bfu_I,\bfq)   = {\bf0} $. \
The linearized  (Newton) form of these  equations  is  as follows:
\eqb
\left[%
\begin{array}{cc}
  \bfK & \bfL \\
  \bfL^T & \bfK_{qq} \\
\end{array}%
\right]
\left[%
\begin{array}{c}
  \Delta \bfu_I \\   \Delta \bfq \\
\end{array}%
\right]
= - \left[%
\begin{array}{c}
   \bfr_u \\      \bfr_q \\
\end{array}%
\right],
\label{set.equil.EAS.newton}
\eqe
where
\ $ \bfK \doteq {\partial \bfr_u}/{\partial \bfu_I} $, \quad $ \bfL \doteq {\partial \bfr_u}/{\partial \bfq} $ \quad
and \quad $ \bfK_{qq} \doteq {\partial \bfr_q}/{\partial \bfq} $. \
Note that \ $ \bfK $ \ and \ $ \bfK_{qq} $ \ are symmetric and, in general, indefinite.

\medskip
To eliminate the multipliers \ $ \Delta \bfq $ \ at the element level
and to reduce the size of a tangent element's matrix to the standard one,
defined by the number of nodes on the element and dofs/node,
we calculate \ $ \Delta \bfq $ \ from the second equation
and use it in the first one, which yields
\begin{equation}
\bfK^* \Delta \bfu_I = - \bfr^*, \qquad \mbox{where}\qquad
\bfK^* = \bfK - \bfL \bfK^{-1}_{qq} \bfL^T \quad \mbox{and} \quad \bfr^* = \bfr_u - \bfL\bfK^{-1}_{qq}\bfr_q.
\label{Reduced.eqs}
\end{equation}
Comparing to  eq.~(\ref{Eq:ischurSystem}), we see that
 the reduced (or condensed) matrix \ $\bfK^*$ \ is defined as the Schur complement of \ $ \bfK_{qq} $. \
Hence, the \emph{Partial Factorization} method described in Sect.~\ref{subsection:partialFactorization} \
can be applied to  speed up the above condensation.
Below, we test and compare its performance for three types of finite elements.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical results}

Tests were performed using only 1 core of a multi-core machine
(2 processors Xeon X5650 2.66GHz with 6 cores each, running under Linux). %  and 24GB DDR3 1333MHz RAM memory,
%Test code was compiled with Intel Fortran Compiler 11.0 with compiler flag O3.
This is in accord with parallelization of a loop over elements in FEAP,
see \cite{Jarzebski2015}  for details,
by which each core processes a different finite element.

\medskip
The stiffness matrices used in computations were obtained for the central (irregular) elements of the
patch tests described in [MacNeal, Harder, 1985].
The  computations were repeated 1 million times for each matrix. % to obtain a good estimate of time.
The speedups are presented in Table~\ref{tab:results1}, \ where the best results for each matrix are boldfaced.

%--------------------------------------------------------------------------------------------
\begin{table}[htbp]
	 \caption{Speedup for particular solvers and methods.
	 	% related to the time for scheme of eq.~(\ref{Reduced.eqs})  	and the LUDCMP  routine.
	 }
 %\scriptsize
 \small
  \centerline{
    %\renewcommand{\arraystretch}{1.3}
    \begin{tabular}{l|ccc|ccc|ccc} \hline
            &  \multicolumn{3}{c|}{Shell elements} &  \multicolumn{3}{c|}{Solid-shell elements}  &  \multicolumn{3}{c}{3D elements} \\
   Solver \qquad Method   &     EAS10  &   HW35  &  HW43  & EAS10 &  HW29   & HW47 & EADG12 &   EAS30  &   HW60  \\
   \hline
%      Crout \cite{Press1999}        & 1.00      &  1.00    & 1.00   &  1.00     &   1.00  & 1.00  &  1.00     &   1.00  & 1.00   & 1.00  &  1.00   & 1.00  \\
      DSYTRF \hfill 1RHS           & 0.26  &  0.48   & 0.68 &   0.30     &   0.56  & 0.60    & 0.29      &  0.45    & 0.75 \\
      MA64 \hfill 1RHS            & 0.57  &  1.02   & 1.36  &   0.72     &   1.37  & 1.51   & 0.84      &  1.23    & 1.74 \\ [0.7 ex]
%      &  \multicolumn{3}{c|}{} &  \multicolumn{3}{c|}{}  &  \multicolumn{3}{c}{} \\
      DSYTRF \hfill mRHS      & 0.97  &  1.72   & 2.06 &    1.04     &   1.93  & 1.82   & 1.12      & 1.49    & 2.20 \\
      MA64 \hfill mRHS        & 0.62  &  1.04   & 1.37  &    0.77     &   1.39  & 1.82 & 0.90      &  1.26    & 1.75 \\ [0.7 ex]
%      &  \multicolumn{3}{c|}{} &  \multicolumn{3}{c|}{}  &  \multicolumn{3}{c}{} \\
      DSYTRF \hfill ownPF    & \textbf{1.27}  &  \textbf{2.47}   & \textbf{3.37}   &  \textbf{1.37}  &  \textbf{2.64}  & \textbf{2.52} & \textbf{1.57} & \textbf{1.86} & \textbf{2.41}  \\
      MA64 \hfill PF     & 0.58  &  1.41   & 1.91  & 0.72 &    1.90  & 2.10  & 0.85      &  1.49    & 2.31   \\
       \hline
     Matrix density [\%] & 91.20     &   100.00 &  58.70 &   98.40    &   48.20 &  29.20 & 40.60 & 22.20   & 13.70 \\
     Reference time [secs]        & 6.77      &  49.96     & 129.08  &  10.46 & 75.19 & 108.19  & 13.34 &  53.88  & 192.54  \\ \hline
    \end{tabular}
    }

  \label{tab:results1}
\end{table}

\emph{Solver and method.}
For reference, we used the scheme of eq.~(\ref{Reduced.eqs})  and the LUDCMP  routine of
Numerical Recipes [Press et al., 1999], which performs the $LU$ factorization using the Crout method;
the times obtained are designated ``Reference time'' in Table~\ref{tab:results1}. \ 
The speedup is computed as a quotient $ {\rm time}_{method}/{\rm time}_{ref} $. \
Besides, we tested two routines based on the Gauss elimination: \
DSYTRF of  [LAPACK ver. 3.2] and MA64 of [HSL, 2013]. \
To obtain \ $\bfK^{-1}_{qq} \bfL^T$, \ the back-substitution routine is called in two ways:
either for each column of \ $ \bfL^T$ separately (``1RHS''),
or for all columns of \ $ \bfL^T$ together  (``mRHS''). \
Besides, ``ownPF'' indicates our modification of the code to perform \emph{Partial Factorization}.
(Note that the PF is not suitable for the Crout method.)

\medskip
\emph{Finite elements.}
Three types of elements were tested: 4-node shells,  8-node solid-shells and 8-node 3D elements.
Their formulations are designated by:
HW - based on the Hu-Washizu functional and enhanced for shell and solid-shell elements \cite{Wisniewski2010, Wisniewski2012},
EAS  - based on the potential energy with the Enhanced Assumed Strain, and
EADG - based on the potential energy   with Enhanced Assumed Displacement Gradient.
The number of additional parameters follows these letters.
%\medskip
The results presented in Table~\ref{tab:results1} \ indicate that
\emph{Partial Factorization} is beneficial for all elements
and that the solver ``DSYTRF'' and the method ``ownPF'' provide the best speedup. % for all tested elements.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\vspace{-5 mm}
\section{Domain Decomposition solvers using PF}
\vspace{-2 mm}
\subsection{Sequential solver for single machine}
\label{section:memoryReduction}

The parallel algorithm for solving a system of equations of Sect.~\ref{section:partialFactorization}
can also be performed sequentially on a single machine.
This reduces the memory usage, thus allowing to compute bigger examples.

\medskip
We divide the domain into two sub-domains.
The memory usage is reduced because the fill-in of the sub-domain matrix of eq.~(\ref{Eq:schurSystem})
is less than half of that for the whole matrix.
However, after factorization of the first sub-domain matrix
we have to deallocate the memory to compute the second sub-domain matrix.
Hence, we need to factorize this matrix once more
after solving  eq.~(\ref{Eq:schurSystemEquationSolve}),
which  increases the time of computation.
We used the above method with two solvers: either Pardiso or modMA86,
which is the HSL MA86 solver with our implementation of \emph{Partial Factorization}.

\medskip
To test this method we computed a linear example of an elastic cube,
loaded by uniform unit forces on the top face and fixed at the bottom face.
A mesh of \ $ N \times  N \times N $ \ 3D 8-node standard elements was used with \ $ N=64$. \
%Computations were performed on the machine with
%2 processors Xeon X5650 2.66GHz (6 cores each) and 24GB DDR3 1333MHz RAM memory.
The results presented in Table~\ref{Tab:MemoryReduction} \ indicate that
the new method saves around 32-38\%  memory but is about 37-40\% slower than the standard one.
We see that each of the solvers tested has its merits.
%-----------------------------------------------------------------------
\begin{table}[h]
 \footnotesize
  \centering
  \caption{Memory and time usage by a sequential DD solver.} %Unit cube test $N=64$.}
  \label{Tab:MemoryReduction}
\begin{tabular}{c|ccc|ccc}
   & \multicolumn{3}{c|}{Memory [GB]}   & \multicolumn{3}{c}{Time [sec]}  \\
  Solver        & Standard    & New       &  Change [\%]    & Standard    & New       & Change [\%]  \\
 \hline
 Pardiso &  14.16       & 8.71      & -38   & 210.98      & 288.90     & +37   \\
 modMA86 &  14.47       & 9.91      & -32   & 156.07      & 218.83     & +40   \\
 \hline
\end{tabular}
\end{table}
%-----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\vspace{-10 mm}
\subsection{Parallel solver for cluster of computers}
\label{section:distributedSolver}


The parallel algorithm for solving a system of equations of Sect.~\ref{section:partialFactorization}
is naturally suited for parallel solution on a cluster of computers.
In each sub-domain we use the  modMA86 solver, which is the HSL MA86 solver with
our implementation of \emph{Partial Factorization}.\
Each sub-domain is computed by different computational node of a cluster.
On each node we use multithreaded version of modMA86 solver with use of 12 threads.
The equation for interfaces (\ref{Eq:schurSystemEquationSolve}) is solved using
the  solver for dense matrices MA64 of HSL.
Additional special treatment is necessary for the cases of 4 and 8 nodes;
more details on our implementation can be found in \cite{JarzebskiWisniewski2016b}. \

\medskip
Our solver is compared to the distributed version of WSMP solver of IBM, see \cite{Gupta2015},
by computations on 2,4, and 8 computational nodes of the type characterized in Sect.~\ref{section:memoryReduction}. \
The speedup is shown according to run on 1 node with 12 threads.
Numerical results for the cube example (the same as in  Sect.~\ref{section:memoryReduction}) \
are presented in Figs.~\ref{Fig:TimeWSMPmy} \  and \ref{Fig:ScalWSMPmy}. \
We note that our solver compares favorably with  WSMP (it uses also  14\% less memory).
%and the speedup is around 2.7 on 8 nodes, much below the ideal one.
%Nonetheless,
An additional speedup is provided by the OMP parallelization of the  MA86 solver
which, on the 12-core machine, is about 9.1.
Hence, the total speedup provided by our cluster solver is about 24.4.
%---------------------------------------------------------------------------------------------
\hspace{-0.5cm}
\begin{figure}[h]
\graphicspath{ {chap5/} }
    \begin{subfigure}[b]{0.5\textwidth}
     \scalebox{0.6}{\includegraphics{time_wsmp_my.eps}}
     \caption{}
     \label{Fig:TimeWSMPmy}
    \end{subfigure}
    \begin{subfigure}[b]{0.5\textwidth}
     \centering
     \scalebox{0.6}{\includegraphics{scal_wsmp_my.eps}}
     \caption{}
     \label{Fig:ScalWSMPmy}
    \end{subfigure}
    \caption{Comparison of time and speedup of cluster solvers.}
 \end{figure}
%-----------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\vspace{-10 mm}
\section{Conclusions}
\label{section:conclusions}

We applied the PF technique to two problems involving the matrix reduction 
%(1) the elemental tangent matrix, and 
%(2) the sub-domain matrices for the DD equation solvers, 
%run either on a single machine or on a cluster.
and for both  the use of this technique has proven beneficial.
For the first problem (Sect.~3), it saves time and the best speedup is provided 
when it is combined with the solver ``DSYTRF''.
For the second problem (Sect.~4), it either 
saves  32-38\% of memory, when run on a single machine,
or provides the total speedup  24.4 on the cluster.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\section{Discussion and Conclusions}
%\label{section:conclusions}
%
%In the current paper, we describe the possible applications of \emph{Partial Factorization}.
%We implemented concept of partial factorization in HSL MA86 code.
%With implementation of this concept we could derive two algorithms:
%(1) algorithm for reduction of memory usage in the solver.
%We achieve 32\% less memory usage.
%(2) domain decomposition solver.
%We achieve around 2.68 speedup.
%Additionally we proved that
%the implementation of the \emph{Partial Factorization} is beneficial for
%condensation of elements matrices.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{-5 mm}
\small
\begin{thebibliography}{99}

\bibitem{Gupta2015}  %[Gupta 2015]
Gupta A.,
\it WSMP: Watson Sparse Matrix Package, %Version 15.01,
\rm IBM Research Report, Watson, 2015


%\bibitem{Hogg2010}  %[Hogg 2010]
%Hogg J.D., Scott J.A.:
%\it An indefinite sparse direct solver for multicore machines.
%\rm Technical Report TR-RAL-2010-011 (2010)

%\bibitem{Hogg2010b}  %[Hogg 2010b]
%Hogg J.D., Scott J.A.,
%\it A note on the solve phase of a multicore solver,
%\rm Technical Report TR-RAL-2010-007, 2010.


\bibitem{Jarzebski2015}
Jarzebski P., Wisniewski K., Taylor R.L.,
\it   On paralelization of the loop over elements in FEAP,
\rm Computational Mechanics, \textbf{56}(1): 77-86, 2015.

\bibitem{JarzebskiWisniewski2016}  %[Jarzebski 2016]
Jarzebski P., Wisniewski K.,
\it Performance of the parallel FEAP in calculations of effective material properties using RVE,
\rm  Advances in Mechanics, Kleiber et al. (Eds), Taylor \& Francis, London, 241-244, 2016.

\bibitem{JarzebskiWisniewski2016b}  %[Jarzebski 2016b]
Jarzebski P., Wisniewski K.,
\it Application of partial factorization for domain decomposition solver,
\rm  In preparation, 2016.

%\bibitem{KarypisKumar95}  %[Karypis, Kumar, 1995]
%Karypis G., Kumar V.:
%\it A fast and high quality multilevel scheme for partitioning irregular graphs.
%\rm International Conference on Parallel Processing, 113-122 (1995)
%
%\bibitem{METIS} %[METIS]
%METIS, http://glaros.dtc.umn.edu/gkhome/metis/metis/overview
%
%\bibitem{MKL} %[IntelMKL]
%Intel MKL, http://software.intel.com/en-us/intel-mkl
%
%\bibitem{MUMPS} %[MUMPS]
%MUMPS: a MUltifrontal Massively Parallel sparse direct Solver, http://http://mumps.enseeiht.fr/
%
%\bibitem{MPI} %[MPI]
%MPI, http://www.mpi-forum.org
%
%\bibitem{OpenMP} %[OpenMP]
%OpenMP, http://openmp.org
%
%\bibitem{PARDISO} %[PARDISO]
%Intel PARDISO, http://software.intel.com/sites/products/documentation
%/hpc/mkl/mklman/index.htm\#GUID-7E829836-0FEF-46B2-8943-86A022193462.htm

%\bibitem{Press1999}
% Press W. H. et al.,
%\emph{Numerical Recipes in Fortran 77}, Cambridge Univeristy Press, 1999.

\bibitem{Petra2014}   % [Petra, 2014]
Petra C.G. et al.,
\it An augmented incomplete factorization approach for computing the Schur complement in stochastic optimization,
\rm SIAM J. Sci. Comput., \textbf{36}(2): C139-C162, 2014.

\bibitem{Wisniewski2010}  % [Wisniewski, 2010]
Wisniewski K.,
\it Finite Rotation Shells. Basic Equations and Finite Elements for Reissner Kinematics.,
\rm Springer, 2010.

\bibitem{Wisniewski2012}  %[Wisniewski, Turska 2012]
Wisniewski K., Turska E.,
\it Four-node mixed Hu-Washizu shell element with drilling rotation,
\rm \IJNME, \textbf{90}: 506-536, 2012.

\end{thebibliography}
\normalsize

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\end{document}
