libhsl.tex

\documentclass[gdweb]{geradwp}

\PassOptionsToPackage{hyphens}{url}
\usepackage[ruled,linesnumbered]{algorithm2e}
\usepackage{hyperref}

\GDcoverpagewhitespace{6.8cm}
\graphicspath{{Figures/}}
\hypersetup{colorlinks,%
citecolor={blue},
urlcolor={blue},
linkcolor={blue},
breaklinks={true}
}

\makeatletter
\ifthenelse{\isundefined{\ALG@name}}{}%
{%
\renewcommand{\ALG@name}{\sffamily\footnotesize Algorithm}
}
\makeatother
\ifthenelse{\isundefined{\SetAlCapNameFnt}}{}%
{%
\SetAlCapNameFnt{\footnotesize}
\SetAlCapFnt{\sffamily\footnotesize}
}

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%% Début - commandes de l'auteur %%%%%%%%%%%%
%% %%%%%%%% Start of author commands %%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Package verbatim
\usepackage{verbatim}

% Package natbib
\usepackage[numbers]{natbib}

% Package jlcode
\usepackage[charsperline=92,theme=default-plain]{jlcode}

% Commande LaTeX pour rajouter les liens DOI
\newcommand*{\doi}[1]{DOI: \href{http://dx.doi.org/\detokenize{#1}}{\detokenize{#1}}}
\newcommand*{\doilink}[2]{\href{http://dx.doi.org/\detokenize{#1}}{#2}}

\renewcommand{\_}{\raisebox{+0.35mm}{\textunderscore}}
\newcommand{\LibHSL}{LibHSL\xspace}
\newcommand{\HSLjll}{HSL\_jll.jl\xspace}
\newcommand{\Ipoptjll}{Ipopt\_jll.jl\xspace}
\newcommand{\LBTjll}{libblastrampoline\_jll.jl\xspace}
\newcommand{\OpenBLASjll}{OpenBLAS32\_jll.jl\xspace}

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%% Fin - commandes de l'auteur %%%%%%%%%%%%%
%% %%%%%%%%% End of author commands %%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%% Métadonnées du cahier  %%%%%%%%%%%%%%%
%% %%%%%%%%%% Working paper metadata %%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\GDtitre{\texttt{\LibHSL}: the ultimate collection for large-scale \\ scientific computation}
\GDmois{Décembre}{December}
\GDannee{2023}
\GDnumero{XX}
\GDauteursCourts{J. Fowkes, A. Lister, A. Montoison, D. Orban}
\GDauteursCopyright{Fowkes, Lister, Montoison, Orban}
\GDpostpubcitation{J. Fowkes, A. Lister, A. Montoison, D. Orban (2023). ``Un exemple de citation'', \emph{Journal of Journals}, vol. X issue Y, p. n-m}{https://www.gerad.ca/fr}
\GDsupplementname{Internet Appendix}

\begin{document}

\GDpageCouverture

\begin{GDpagetitre}

\begin{GDauthlist}
    \GDauthitem{Jaroslav Fowkes \ref{affil:1}}
    \GDauthitem{Andrew Lister \ref{affil:1}}
    \GDauthitem{Alexis Montoison \ref{affil:2}}
    \GDauthitem{Dominique Orban \ref{affil:2}}
\end{GDauthlist}

\begin{GDaffillist}
    \GDaffilitem{affil:1}{Scientific Computing Department, Rutherford Appleton Laboratory, Chilton, England}
    \GDaffilitem{affil:2}{GERAD and Department of Mathematics and Industrial Engineering, Polytechnique Montréal, Montr\'eal (Qc), Canada, H3T 2A7}
\end{GDaffillist}

\begin{GDemaillist}
    \GDemailitem{jaroslav.fowkes@stfc.ac.uk}
    \GDemailitem{andrew.lister@stfc.ac.uk}
    \GDemailitem{alexis.montoison@polymtl.ca}
    \GDemailitem{dominique.orban@gerad.ca}
\end{GDemaillist}

\end{GDpagetitre}

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%% Résumés, mots-clés, remerciements %%%%%%%%%%%
%% %%%%%%% Abstract, keywords, acknowledgements %%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\GDabstracts

\begin{GDabstract}{Abstract}
The Harwell Subroutine Library (HSL) is a renowned suite of efficient and robust numerical algorithms designed to tackle complex mathematical problems such as solving sparse linear systems or computing eigenvalues and eigenvectors of sparse matrices.
\LibHSL is a comprehensive collection of HSL packages that facilitates integration into C, Fortran or Julia projects.
One of the notable features of \LibHSL is its innovative compilation system based on Meson.
The new build system not only significantly accelerates compilation, but also ensures portability across multiple operating systems.
To further streamline the user experience, \LibHSL provides archives of precompiled libraries as well as a package called \HSLjll, specifically designed for the Julia ecosystem.
This eliminates the need for users to compile anything on their local machines, offering ready-to-use functionality right away.

\paragraph{Keywords: }
HSL, Julia, Fortran, C, sparse linear algebra, Meson, cross-compilation, portability, performance
\end{GDabstract}

\begin{GDabstract}{R\'esum\'e}
La bibliothèque de sous-programmes Harwell (HSL) est une suite renommée de méthodes numériques efficaces et robustes conçus pour résoudre des problèmes mathématiques complexes tels que la résolution de systèmes linéaires creux ou le calcul de valeurs propres et de vecteurs propres pour les matrices creuses.
\LibHSL est la collection ultime des bibliothèques HSL et facilite leur intégration dans des projets en C, Fortran ou Julia.
L'un des aspects notables de \LibHSL est son système de compilation basé sur Meson.
Ce nouveau système de compilation accélère considérablement la compilation tout en assurant la portabilité sur différents systèmes d'exploitation.
Afin de rendre l'expérience utilisateur encore plus simple, \LibHSL offre des archives contenant des bibliothèques dynamiques précompilées ainsi qu'un module \HSLjll, spécialement conçu pour l'écosystème Julia.
Les utilisateurs évitent ainsi la compilation sur leurs machines locales avec une bibliothèque prête à l'emploi.

\paragraph{Mots cl\'es\,: }
HSL, Julia, Fortran, C, alg\`ebre lin\'eaire creuse, Meson, compilation multi-plateformes, portabilit\'e, performance
\end{GDabstract}

\begin{GDacknowledgements}
Alexis Montoison is supported by an FRQNT grant and an excellence scholarship of the IVADO institute, and Dominique Orban is partially supported by an NSERC Discovery Grant.

This work began while Alexis Montoison was visiting Nick Gould, Jaroslav Fowkes and the \href{https://www.numerical.rl.ac.uk/}{Computational Mathematics Group} at the \href{https://www.ukri.org/about-us/stfc/locations/rutherford-appleton-laboratory/}{STFC Rutherford Appleton Laboratory} in December 2022 and July 2023.
He expresses its appreciation to the HSL team for their invitation and their wonderful hospitality.

The authors would like to thank:
\begin{itemize}
  \item Stefan Vigerske for the \href{https://github.com/KarypisLab/METIS}{METIS} adapter;
  \item Eli Schwartz for its advices with the \href{https://mesonbuild.com/}{Meson} build system;
  \item Geoffroy Leconte, Paul Raynaud, and Oscar Dowson for testing early versions of \HSLjll;
  \item Mosè Giordano, and Elliot Saba for their work on \href{https://github.com/JuliaPackaging/BinaryBuilder.jl}{BinaryBuilder.jl}, \href{https://github.com/JuliaLinearAlgebra/libblastrampoline}{libblastrampoline} and \href{https://github.com/JuliaPackaging/Yggdrasil}{Yggdrasil}.
\end{itemize}
\end{GDacknowledgements}

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%% Article %%%%%%%%%%%%%%%%%%%%%%%%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\GDarticlestart

\section{Introduction}

The HSL Mathematical Software Library, originally known as the Harwell Subroutine Library, is a valuable resource for researchers, scientists, and engineers working on numerical computations.
It provides a comprehensive collection of robust, well-tested, and highly efficient mathematical routines, covering areas such as sparse linear or eigenvalue problems.

\href{https://licences.stfc.ac.uk/product/libhsl}{\LibHSL} is a collection of more than 160 HSL packages.
It aims to facilitate the use of HSL in Julia \citep{bezanson-edelman-karpinski-shah-2017} as well as Fortran and C projects such as GALAHAD \citep{gould-orban-toint-2003}.
\LibHSL focuses on ease of use for users on all platforms by using the \href{https://github.com/mesonbuild/meson}{Meson} build system.
Meson enables the distribution of prebuilt binaries for the package using \href{https://github.com/JuliaPackaging/BinaryBuilder.jl}{BinaryBuilder.jl}.
Additionally, \LibHSL supports either METIS 5 or the older METIS 4 \citep{karypis-kumar-1998}.

This new package provides the source code of the included HSL packages, prebuilt binaries, and a Julia package named \HSLjll.
\HSLjll is a pre-built version of \LibHSL to be readily used in the Julia ecosystem.
Once installed, \HSLjll allows calling the shared library, containing the C and Fortran routines of the HSL packages, directly from Julia.
The HSL wrappers provided in the Julia interface HSL.jl become fully functional.
\HSLjll also offers a convenient way to employ HSL linear solvers MA27, MA57, HSL\_MA77, HSL\_MA86 and HSL\_MA97 within the Ipopt.jl interface for the IPOPT \citep{wachter-2006} nonlinear optimization solver.

Two versions of \HSLjll are included. One precompiled with OpenBLAS that requires at least Julia 1.6 and the other version precompiled with libblastrampoline (LBT) that requires at least Julia 1.9.
LBT allows one to dynamically switch the BLAS and LAPACK backends between e.g. OpenBLAS, BLIS, Intel MKL or Apple Accelerate.
\HSLjll is precompiled for various operating systems (Windows, Mac, Linux, FreeBSD) and architectures (x64, arm64, ppc64).

For information on how to use each HSL routine available through this package
please see the relevant documentation on the \href{https://www.hsl.rl.ac.uk/catalogue/}{HSL website}.
Where C interfaces are available, the C documentation should be used.

\section{Building \LibHSL with Meson}

\subsection{Meson}

To download and install \href{https://mesonbuild.com}{Meson}, we refer the user to this \href{https://mesonbuild.com/SimpleStart.html}{guide}.
By default Meson uses the \href{https://ninja-build.org/}{Ninja} build system.
After configuration, it is equivalent to build the package or install it directly with Ninja.

\subsection{Configuration}

To configure the compilation for your system, first create and initialise a build directory.
For examples in this document, we will use \jlinl{builddir} for our build directory

\begin{jllisting}
meson setup builddir [options...]
\end{jllisting}
To choose an alternative install directory you can pass the \jlinl{--prefix} argument.
Non-default compilers can be selected by setting the \jlinl{CC} and \jlinl{FC} shell variables.

\begin{jllisting}
CC=icc FC=ifort meson setup builddir --prefix=~/libhsl [options...]
\end{jllisting}
See this \href{https://mesonbuild.com/Reference-tables.html}{table} for supported compilers.
The compilation of libHSL has been tested with the Fortran compilers
\jlinl{gfortran}, \jlinl{ifort}, \jlinl{ifx} and \jlinl{nagfor}.
Note that the version 1.2.0 of Meson is required for using \jlinl{ifx} on Linux.

The list of all options supported by libHSL can be displayed with

\begin{jllisting}
meson configure
\end{jllisting}
This command can also be used from the build directory to change options,
if required, after running \jlinl{meson setup}.

To test the installation using our examples after building, configure with
the \jlinl{-Dexamples=true} option.

\begin{table}[ht]
  \label{tab:options}
  \centering
  \begin{tabular}{|l|c|l|}
    \hline
    {\bf Options}     & {\bf Default value} & {\bf Description}                                         \\ \hline
    modules           & true                & Install Fortran modules                                   \\ \hline
    examples          & false               & Generate the examples                                     \\ \hline
    libmetis\_version & 5                   & Version of the METIS library                              \\ \hline
    libblas           & blas                & BLAS library against which to link                        \\ \hline
    liblapack         & lapack              & LAPACK library against which to link                      \\ \hline
    libmetis          & metis               & METIS library against which to link                       \\ \hline
    libmpi            & mpifort             & MPI library against which to link                         \\ \hline
    libblas\_path     & $\varnothing$       & Additional directory to search for the BLAS library       \\ \hline
    liblapack\_path   & $\varnothing$       & Additional directory to search for the LAPACK library     \\ \hline
    libmetis\_path    & $\varnothing$       & Additional directory to search for the METIS library      \\ \hline
    libmpi\_path      & $\varnothing$       & Additional directory to search for the MPI library        \\ \hline
    libmetis\_include & $\varnothing$       & Additional directory to search for the METIS header files \\ \hline
    libmpi\_include   & $\varnothing$       & Additional directory to search for the MPI header files   \\ \hline
  \end{tabular}
  \caption{Project options supported by libHSL}
\end{table}

\subsection{Compilation}

After configuration, the source can be compiled in the build directory with

\begin{jllisting}
# Meson
meson compile -C builddir
\end{jllisting}

\subsection{Installation}

Finally, the library, headers, and Fortran modules can be installed to the specified path with

\begin{jllisting}
# Meson
meson install -C builddir
\end{jllisting}

This installs the library, headers, and Fortran modules in the folder specified by \jlinl{--prefix}, if the argument is used.
If the argument \jlinl{--prefix} is not used, the default installation folder is \jlinl{C:/} on Windows, and \jlinl{usr/local} otherwise.


\subsection{Running Examples}

If \LibHSL is configured to build the example programs (\jlinl{-Dexamples=true} in the setup call),
compiled versions of these will be available in your install directory as well
as a utility program \jlinl{run_example} to help run them.
By default, example programs are for the Fortran interface.
Programs postfixed with \jlinl{_c} are examples written in C.

To keep the example code simple, data is read from \jlinl{stdin} in each program.
If the example requires an input, the data can be found in the install directory
under \jlinl{examples/Fortran/data} or \jlinl{examples/C/data}.
The expected output is also available in the install directory under
\jlinl{examples/Fortran/output} or \jlinl{examples/C/output}.

To run an example on linux after installing with \jlinl{--prefix=<installdir>},
the command would look like:
\begin{jllisting}
cd <installdir>/examples/Fortran
./hsl_ma97ds < data/hsl_ma97ds.data > output.txt
\end{jllisting}
In order to simplify this and make it available to non-linux users the following is equivalent
\begin{jllisting}
cd <installdir>/examples
./run_example Fortran/hsl_ma97ds output.txt Fortran/data/hsl_ma97ds.data
\end{jllisting}

\jlinl{run_example} supports 1, 2, or 3 arguments:
\begin{itemize}
  \item \jlinl{run_example <example_program>} runs the example with no inputs and prints all outputs to stdout,
  \item \jlinl{run_example <example_program> <output_file>} runs the example with no inputs and writes all outputs to \jlinl{output_file},
  \item \jlinl{run_example <example_program> <output_file> <input_file>} runs the example with input from \jlinl{input_file} and writes all outputs to \jlinl{output_file}.
\end{itemize}

\subsection{Cross-compilation}

\LibHSL can be easily cross-compiled with the Julia package \href{https://github.com/JuliaPackaging/BinaryBuilder.jl}{BinaryBuilder.jl}.
We provide the script \jlinl{build\_tarballs.jl} to easily regenerate an \HSLjll with different options or dependencies, such as METIS 4 instead of METIS 5.
It can also be used to generate a static library \jlinl{libhsl.a} for an application outside of Julia.
Commands to cross-compile \LibHSL are available in the file \jlinl{build\_tarballs.jl}.
Note that the fortran compiler used to cross-compile \LibHSL is \jlinl{gfortran}.

\section{Use \LibHSL in C and Fortran projects}

Once \LibHSL is compiled, you can link it to your C and Fortran projects with
\begin{jllisting}
# C
-I${prefix}/include -L${prefix} -lhsl

# Fortran
-I${prefix}/modules -I${prefix}/include -L${prefix} -lhsl
\end{jllisting}

By default Meson generates a shared library but, if you prefer to do static linking,
the Meson's option \jlinl{--default-library=static} can be used to generate a
static library \jlinl{libhsl.a}.

\section{Use \HSLjll in the Julia ecosystem}

\subsection{Installation of the \HSLjll package}

To install the package \HSLjll, you only need the following commands in the Julia REPL:
\begin{jllisting}
julia> ]
pkg> dev path_to_hsl_jll
\end{jllisting}
We recommend using the version precompiled with \LBTjll if you use Julia 1.9 or future stable releases.
The version precompiled with \OpenBLASjll can be used in the long-term support (LTS) release Julia 1.6.
Due to security policy of the macOS operating system, the Mac users may need to remove the quarantine before extracting.
\begin{jllisting}
xattr -d com.apple.quarantine openblas_HSL_jll.jl-2023.11.7.zip
xattr -d com.apple.quarantine openblas_HSL_jll.jl-2023.11.7.tar.gz
xattr -d com.apple.quarantine lbt_HSL_jll.jl-2023.11.7.zip
xattr -d com.apple.quarantine lbt_HSL_jll.jl-2023.11.7.tar.gz
\end{jllisting}

\subsection{Load an LP64 BLAS and LAPACK backend}

By default, Julia ships with only an ILP64 version of OpenBLAS as BLAS and LAPACK backend.
ILP64 libraries use the 64-bit integer type (necessary for indexing large arrays, with more than $2^{31}-1$ elements), whereas the LP64 libraries index arrays with the 32-bit integer type.
\LibHSL can only be linked with LP64 versions of BLAS and LAPACK as HSL uses 32-bit integer index arrays.
Thus, one version of \HSLjll is precompiled with \OpenBLASjll and automatically loads an LP64 version of OpenBLAS with a \jlinl{using HSL_jll}.
For the version precompiled with \LBTjll, the user needs to manually load one LP64 BLAS and LAPACK backend except if \HSLjll is used with \href{https://github.com/JuliaSmoothOptimizers/HSL.jl}{HSL.jl} \citep{montoison-orban-hsl-2021} or \href{https://github.com/jump-dev/Ipopt.jl}{Ipopt.jl}.
\begin{jllisting}
# Load the LP64 version of OpenBLAS
import LinearAlgebra, OpenBLAS32_jll
LinearAlgebra.BLAS.lbt_forward(OpenBLAS32_jll.libopenblas_path)

# Load LP64 and ILP64 versions of Intel MKL
using MKL

# Load LP64 and ILP64 versions of Apple Accelerate
using AppleAccelerate
\end{jllisting}
We can verify what backends are loaded with \jlinl{LinearAlgebra.BLAS.lbt_get_config()}.
The user can also use their own LP64 version of BLAS and LAPACK.
Note that LP64 and ILP64 versions can coexist if they have different symbols, such as a suffix \jlinl{\_64} in the naming conventions.

\subsection{How to use \HSLjll in Julia packages}

A version 2.0 of \HSLjll based on this dummy \href{https://github.com/ralna/LibHSL}{\LibHSL} was precompiled with \href{https://github.com/JuliaPackaging/Yggdrasil}{Yggdrasil}.
Therefore \HSLjll is a registered Julia package and can be added as a dependency of any Julia package.
The dummy version only has one C function \jlinl{LIBHSL\_isfunctional} that returns the boolean \jlinl{false}.
The real version also exports this C function and returns the boolean \jlinl{true}.
\begin{jllisting}
using HSL_jll

function LIBHSL_isfunctional()
    @ccall libhsl.LIBHSL_isfunctional()::Bool
end

bool = LIBHSL_isfunctional()
\end{jllisting}
This function makes it possible to determine if we have a dummy version or not and consequently to call the C and Fortran routines of the genuine version.

The package \href{https://github.com/JuliaSmoothOptimizers/HSL.jl}{HSL.jl} provides wrappers for all HSL packages with a C interface or those written in Fortran 77.
Wrappers are Julia functions that call C and Fortran routines with the macro \jlinl{@ccall}.
Thus, the combination of HSL.jl and \HSLjll allows one to abstract away from the low-level C and Fortran languages and use the majority of HSL packages as if they were packages developed in Julia.
% High-level Julia functions that exploit the multiple dispatch feature are underway.
Note that HSL.jl also loads the LP64 version of OpenBLAS automatically for the user if needed.

For the Julia interface \href{https://github.com/jump-dev/Ipopt.jl}{Ipopt.jl}, \HSLjll must be used directly because it interacts with the \Ipoptjll package under the hood.

\begin{jllisting}
using JuMP
import Ipopt, HSL_jll

# An optimization problem
model = Model(Ipopt.Optimizer)
@variable(model, x)
@objective(model, Min, (x - 2)^2)

# Load the HSL solvers
set_attribute(model, "hsllib", HSL_jll.libhsl_path)

# Use the linear solver MA57
set_attribute(model, "linear_solver", "ma57")
optimize!(model)

# Use the linear solver MA97
set_attribute(model, "linear_solver", "ma97")
optimize!(model)
\end{jllisting}

The available HSL linear solvers are \jlinl{"ma27"}, \jlinl{"ma57"}, \jlinl{"ma77"}, \jlinl{"ma86"} and \jlinl{"ma97"}.
If no LP64 BLAS and LAPACK backend is loaded, Ipopt.jl loads \OpenBLASjll automatically just like HSL.jl.

\section{Future improvements}

\subsection{Support for ILP64 interfaces of BLAS and LAPACK}

Currently, libHSL only supports LP64 versions of BLAS and LAPACK.
A planned enhancement for the future is to add support for ILP64 versions of the library.
This upgrade would enable the resolution of even larger-scale problems while simplifying the integration of libHSL into optimization solvers or differential equation solvers compiled with 64-bit integers.
Extending to ILP64 would broaden the applicability of libHSL, meeting the growing demands of scientific and engineering applications.

\subsection{GPU-accelerated linear solvers}

Another anticipated improvement involves adding support for GPU versions of BLAS and LAPACK, such as CUBLAS for NVIDIA GPUs or rocBLAS for AMD GPUs.
The integration of these GPU libraries into the linear solvers could significantly accelerate computations.
However, this modification poses complex challenges, including efficient data transfers between the CPU and GPU.

\section{Bug reports}

If you think you have found a bug, feel free to open an \href{https://github.com/ralna/libHSL/issues}{issue}.
Useful suggestions and requests can also be opened as issues.
If you would like to ask a question not suitable for a bug report, feel free to start a \href{https://github.com/ralna/libHSL/discussions}{discussion}.

\bibliographystyle{abbrvnat}
\bibliography{abbrv,libhsl}

\end{document}

\endinput