§ INTRODUCTION

Real numbers in computer architectures are encoded in bits and represented in floating-point format. Released in 1985, the IEEE 754-1985 governs the standards for floating-point formats and arithmetic [Kahan, 1996, Zuras et al., 2008]. It underwent a revision in 2008, incorporating the addition of 16-bit half-precision, alongside 128-bit quadruple-, 64-bit double-, and 32-bit single-precision formats. The original purpose of introducing the 16-bit half-precision format was for storage. However, it has recently found utility in scientific applications such as machine learning [Björck et al., 2021, Tolliver et al., 2022] and linear algebra [Higham et al., 2019, Luszczek et al., 2019, Scott and Tuma, 2023]. This trend has increased with the emergence of new hardware that can run low-precision computations substantially faster than full-precision, e.g., Graphics Processing Units (GPUs), Intel CPUs with DL boost, ARM CPUs, Tensor Processing Units (TPUs), Field-Programmable Gate Arrays (FPGAs), and AI hardware accelerators.
While employing lower precision can potentially lead to reduced computation times, specific applications demonstrate resilience towards low-precision computations, prompting researchers to explore the creation of innovative multi- and mixed-precision algorithms. These algorithms aim to strike a balance between accuracy and computational efficiency.

In the literature, double-precision used to be the standard numerical format for a wide range of scientific computational workloads such as weather forecasting and climate modeling [Chantry et al., 2019, Lang et al., 2021], nonlinear equations [Benner et al., 2022], numerical linear algebra [Higham and Mary, 2022, Bujanović et al., 2023], ordinary differential equations [Hopkins et al., 2020], quantum chemistry and physics models [Pederson et al., 2022], quantum signal processing models [Ying, 2022], neural networks [Hrycej et al., 2022], and remote sensing [Nguyen et al., 2012, Räss et al., 2019, Li and Fotheringham, 2020]. Recently, there has been a growing interest in reducing the precision of data and computations in these applications and many others. This shift is due to the high cost of double-precision calculations, whereas lower precisions enable hardware acceleration, faster execution times, reduced memory needs, and lower energy consumption. Numerous applications have shown the significant impact of utilizing lower-precision data structures and operations, including in numerical linear algebra [Yamazaki et al., 2015, Yamazaki et al., 2015, Yang et al., 2021, Carson et al., 2022], quantum chemistry and quantum physics [Pederson et al., 2022], computer simulations of nuclear reactor models [Cherezov et al., 2023], computational fluid dynamics [Freytag et al., 2022, Lehmann et al., 2022], weather forecasting and climate modeling [Rüdisühli et al., 2013, Váňa et al., 2017, Hatfield et al., 2020, Klöwer et al., 2020], astronomy [White et al., 2023], and deep learning [Murillo et al., 2022].

Although lowering the precision offers attractive acceleration, naively reducing the precision can lead to catastrophic round-off errors as the range of the numerical representation is limited, thereby significantly changing model results. Thus, the state-of-the-art in numerical computing is moving away from traditionally specifying a notably reduced precision towards multi- and mixed-precisions to squeeze out more performance while maintaining the accuracy of results [Higham and Mary, 2022]. Multi-precision refers to using multiple different floating-point representations for different operations in a single algorithm.
Although several precisions are involved, multi-precision is still limited to prescribing a specific precision for each task of an application or algorithm. On the other hand, mixed-precision refers to using a mix of different numerical precisions within a single operation. For example, the state-of-the-art European ocean model called Nucleus for European Modelling of the Ocean (NEMO) demonstrated that from the 942 variables used for ocean simulation, 902 of them could use single-precision and the rest remain double without degrading model accuracy [Tintó Prims et al., 2019]. Other areas that started to switch to mixed-precision computing are differential equations [Ooi et al., 2020], electromagnetic field analysis [Masui and Ogino, 2019], and weather and climate [Hatfield et al., 2019]. Mixed-precision has also been used in solving high-performance computing (HPC) problems, such as those in molecular dynamics, biology, cosmology, and plasma physics. A comprehensive survey on how various HPC applications utilize mixed-precision can be found in [Netti et al., 2023].

Computational statistics is a field that stands to benefit heavily from these innovations in numerical computing. While computational statistics has long been utilizing computational techniques to expedite the use of statistical methodologies that are computationally expensive, the use of low, multi-, and mixed-precision techniques is just gaining momentum in the literature. For instance, [Guivant et al., 2023] demonstrated that a Gaussian filter in Kalman filtering can work in low precisions. Their experiments showed that processing times can be accelerated by 50% when the covariance matrix is stored in single or half-precision. [Misra et al., 2023] showed that a compiler-inferred fixed-point precision variational inference is $8.15$X and $22.67$X faster than performing variational inference in single- and double-precision, respectively. [Maddox et al., 2022] proposed a conjugate gradients (CG) algorithm used for training Gaussian process models that utilize a mixed-precision approach wherein accuracy-critical computations are cast into full precisions while noncritical parts are performed in lower precisions. They demonstrated up to $3$X speedup compared to equivalent models at full precisions. In geostatistics, where the entries of the spatial covariance matrix quantify the dependence between measurements at two locations, [Abdulah et al., 2019] proposed a novel Cholesky factorization algorithm that works under a mixed-precision format. The algorithm prescribes double-precision for the diagonal elements while representing the off-diagonal elements in single-precision. They defined a `band' to ascertain the number of diagonal elements to be represented in a double-precision format. Their implementation has relied on the  dynamic runtime system for task scheduling on shared- and distributed-memory systems.
[Abdulah et al., 2022] extended the approach in [Abdulah et al., 2019] to include half-precision. They switched the dynamic runtime system to , which makes on-the-fly decisions on precision conversion.
[Cao et al., 2022] improved upon the approach in [Abdulah et al., 2022] by combining mixed-precision techniques with Tile Low-Rank (TLR) matrix algorithms. Furthermore, they devised a new algorithm that determines the appropriate data structure (dense or TLR) and precision (double, single, or half) for every tile at runtime.

The primary numeric format in the widely used statistical software R [R Core Team, 2023] is double-precision, thereby ensuring a substantial level of accuracy for the calculations performed. Despite the many advantages of using multi- and mixed-precision, such as lower energy usage and reduced storage and data movement costs [Higham and Mary, 2022], support for precisions lower than double-precision needs to be improved in R. Indeed in R, single-precision type can be introduced using external packages, such as the float package [Schmidt et al., 2022]. The float package can accommodate single-precision when the system has access to single-precision Basic Linear Algebra Subprograms (BLAS) or LAPACK routines. The package necessitates the underlying BLAS/LAPACK libraries, integrated with the R environment, to support single-precision computation. The package provides support for the majority of the linear algebra operations that are available in R. However, the package design does not readily facilitate easy extension to other precisions, such as half-precision. Another notable R package that supports multi-precision is the Rmpfr package which stands for Multiple Precision Floating-Point Reliable [Maechler et al., 2023]. The package serves as an interface to the MPFR library, a C library for multiple precision floating point arithmetic, which builds on the GNU Multiple Precision Arithmetic C library. The package offers extended precision floating-point representations suitable for precision-demanding tasks such as combinatorics, number theory, and cryptography, where standard floating-point arithmetic can result in unacceptable imprecision. The default precision of MPFR objects is 128-bit. Moreover, the Rmpfr package also supports arbitrary precision such that R users can specify their required precision. For example, computing $\pi$, which is equal to $4 \arctan(1)$, up to 100 decimal places, the following R command can be used: , where the number $1$ is encoded with 333-bit, the equivalent of 100 decimal places. A drawback with the Rmpfr package is that regular numeric types in R cannot be directly used with Rmpfr numbers. When interacting with functions or packages that do not support Rmpfr, an explicit conversion to numeric values is required. To interact directly with the GNU Multiple Precision Arithmetic C library [Lucas et al., 2023], the R package gmp can used.

This paper introduces a novel R package called MPCR, which stands for  Multi- and Mixed-Precision Computational Statistics in R. MPCR is an advanced package designed to provide R users with multi- and mixed-precision data structures and computations. The package harnesses the combined strength of C++ and R, empowering users with high-performance computing capabilities. Specifically tailored for researchers and data scientists working with multi- and mixed-precision arithmetic, MPCR is an invaluable tool for achieving efficient and accurate computations. The package supports three different precision types, namely, 16-bit, 32-bit, and 64-bit, as well as their combinations. Central to the package is creating an MPCR object with a user-defined precision. Leveraging the underlying BLAS/LAPACK library, MPCR facilitates low-precision calculations across various linear algebra operations. Additionally, the package supports mixed-precision calculations by incorporating a tile-based data structure within R. This novel structure enables matrices/vectors to possess different precisions, catering to scientific accuracy requirements. Furthermore, the package allows for parallel linear algebra computations on multicore systems using task-based linear algebra algorithms. The package enables users to generate matrices/vectors of four distinct types: HP-matrix, SP-matrix, DP-matrix, and MP-matrix.
HP-, SP-, and DP-matrices are matrices characterized by uniform precision entries. In an HP matrix, elements are stored in half-precision (16-bit), resulting in a smaller memory footprint and facilitating faster algorithms, albeit with a trade-off in accuracy. Each element inside a DP-matrix requires a bigger memory allocation (64-bit), resulting in slower computations but more accurate solutions. An SP-matrix, wherein each element occupies 32-bit memory, offers an intermediate balance between speed and accuracy. Beyond the trio of uniform-precision classifications, the MPCR package introduces a mixed-precision (MP-) matrix variant. An MP-matrix is a matrix subdivided into tiles (a term used in scientific computing to refer to blocks) so the tiles can be stored at different precisions. To the best of our knowledge, the MPCR package is the first package in R to introduce a tile-based data structure. The MPCR package furnishes mixed-precision functionality for linear algebra operations within the R environment.

In the first version of MPCR, we provide three tile-based operations for MP-matrices: general matrix multiplication (GEMM), triangular solve (TRSM), and tiled Cholesky factorization (POTRF) [Beaumont et al., 2020]. We highlight that MPCR has considerable potential for expansion to other tile-based matrix operations, allowing R developers to benefit from mixed-precision computation. To assess the performance of an algorithm, it is possible to create different MP-matrix structures. For instance, one can allocate higher precision to the diagonal tiles and lower precision to the off-diagonal tiles. As an illustration of the usage of MP-matrices and MP-operations in various statistical applications, this paper includes examples in Markov Chain Monte Carlo, maximum likelihood estimation in spatial statistics, principal component analysis, and Bayesian inference wherein tiled linear algebra operations are performed on appropriately structured MP-matrices to perform matrix inversion, singular value decomposition, and the like.

The subsequent sections of the paper are organized as follows: Section <ref> gives a background on floating-point arithmetic in modern technology and its realization in hardware. Section <ref> briefly reviews tile-based linear algorithms. In Section <ref>, we delve into the internal design of the package and workflow. In Section <ref>, we give an overview of the main functions in the package and how to use them with examples and code snippets. Section <ref> illustrates the efficiency, performance, and accuracy of MPCR objects and functions and compares them with native R objects and functions. In Section <ref>, we present various Frequentist and Bayesian problems (four in total) and illustrate how to plug-and-play MPCR functions to enable faster execution of computationally demanding tasks. We conclude in Section <ref>.

§ FLOATING-POINT ARITHMETIC IN MODERN TECHNOLOGY

The most widely used standard for floating-point arithmetic as of 2024 is the IEEE 754-2008 [Kahan, 1996, Zuras et al., 2008] and its correction in 2019  [IEEE, 2019]. This standard significantly improved over its predecessor, IEEE 754-1985, and includes a comprehensive set of guidelines that cover almost every aspect of the floating-point theory or simply the approximation rules for real numbers on today's digital computers. The full name of this standard is IEEE Standard 754-2019 for Floating-Point Arithmetic, which is often abbreviated as either IEEE 754-2019 or simply IEEE 754. A floating-point number system is a finite subset $\mathcal{F}$ of $\mathbb{R}$, $\mathcal{F} \subset \mathbb{R}$, which depends on its elements $(\beta, t, e_{\text{min}}, e_{\text{max}})$. A number $x$ in $\mathcal{F}$ has the form:

\begin{align*}
x = \pm m  \times \beta^{e-t+1},
\end{align*}

where $t$, $m$, and $e$ are integers known as precision, significand or mantissa, and exponent, respectively, and $\beta$ is the base which is 2 for binary (commonly used on all current computers) and 10 for decimal. Here, the significand satisfies $0 \leq m \leq \beta^t -1$,  the exponent $e$ satisfies $e_{\text{min}} \leq e \leq e_{\text{max}}$, and the IEEE standard requires that $e_{\text{min}} = 1 - e_{\text{max}}$. The number is stored in the computer with a format that consists of three fields: a sign bit, exponent bits, and significand bits.

Precision formats that take up more memory (in terms of bits) but can approximate real number values at very high accuracy are regarded as high or full precision, while those that consume smaller memory with lower accuracy are considered reduced or low precision. In the IEEE-754 standard, the double-precision 64-bit floating-point format, denoted FP64, is considered high precision, while the single-precision 32-bit floating-point format, denoted FP32, is regarded as low precision. FP32 uses 1-bit for the sign of the number, 8-bit for the exponent, and 23-bit for the significand and can represent values from $\pm (2 - 2^{-23}) \times 2^{127}$. On the other hand, FP64 uses 1-bit for the sign, 11-bit for the exponent, and 52-bit for the significand, and can represent values from $\pm (2 - 2^{-52}) \times 2^{1023}$. Another low-precision format is the half-precision 16-bit floating-point, denoted FP16. FP16 uses 1-bit for the sign, 5-bit for the exponent, and 10-bit for the significand and can represent values from $\pm (2 - 2^{-10}) \times 2^{17}$. Figure <ref> visualizes how a number $x$ takes up space in a computer's memory for each precision format.

The composition in bits of different precisions based on IEEE 754-2008 standards.

The FP32 and FP64 are the most widely used formats as a broad range of off-the-shelf general-purpose processors natively supports them. Although FP16 significantly reduces the memory requirements and increases the arithmetic throughput by 2X against FP32 and 4X against FP64, most commercial CPUs do not support FP16. However, recent chips from the leading hardware vendors, including NVIDIA, Intel, and AMD, include FP16 arithmetic in their processing units to enhance performance in applications such as machine learning, gaming, and scientific computing. NVIDIA was one of the first to offer FP16 by enabling it as a storage format in CUDA 7.5 [Sabbagh Molahosseini et al., 2012]. Since then, FP16 arithmetic has been made available for CUDA-enabled GPUs, starting with the PASCAL architecture [Abdelfattah et al., 2020]. Other notable hardware that provides support for FP16 is NVIDIA P100 (2016) and V100 (2017), A100 (2020), and H100 (2022) GPUs through tensor cores. One of the standout features of NVIDIA Tensor Cores is their ability to perform operations in mixed-precision by leveraging lower precision (like FP16) for specific calculations to increase performance and efficiency while still leveraging higher precision (like FP32 and FP64) when higher accuracy calculation is mandatory. ARM chips are another example of hardware supporting 16-bit operations through their Thumb instruction set. The instruction set reduces the memory footprint, as the smaller instruction size reduces the amount of code. New ARM chips, like the ARM Cortex series, support a mix of  32-bit (ARM instruction set) and 16-bit (Thumb instruction set) instructions, allowing running an algorithm in multi-precision mode to optimize performance and memory usage [Higham and Pranesh, 2019].

Most current hardware chips predominantly support 64-bit and 32-bit floating-point (FP) arithmetic operations. A select few, however, can support 16-bit operations, either through mixed-precision computation, as seen in NVIDIA's Tensor Cores, or exclusively, like some ARM processors. There has been a significant trend in software development towards adopting partial or full 16-bit computation, particularly in fields such as machine learning and scientific computing. This shift is driving leading hardware manufacturers and even emerging startup companies to enhance their hardware offerings with support for 16-bit FP operations. Additionally, there is a growing anticipation that support for 8-bit FP operations might become more prevalent shortly, reflecting the ongoing evolution in hardware and software domains to balance performance with computational efficiency.

§ TILE-BASED LINEAR ALGEBRA ALGORITHMS

Linear algebra algorithms underpin the linear algebra operations performed by various computing architectures. These algorithms direct the movement of data and define the tasks needed to perform the operations and arrive at the desired results. The software library LAPACK [Anderson et al., 1999], and its parallel version, ScaLAPACK [Blackford et al., 1996], offer a stable and wide range of linear algebra algorithms necessary for dense linear algebra operations. Additionally, their algorithms were made compatible to parallelism in shared-memory and distributed-memory systems by depending on MPI, OpenMP multithreading, and parallel BLAS and Basic Linear Algebra Communication Subprograms (BLACS) [Luszczek et al., 2014]. The parallelization, particularly in ScaLAPACK, is carried out using the block-cyclic decomposition technique wherein the matrix is divided into small pieces called blocks such that these blocks are then distributed to different MPI processes in a cyclical fashion [Cifani et al., 2023].  This type of parallelization scheme falls under the traditional bulk synchronous programming model wherein multiple threads will have to wait in idle until the slowest performing thread is finished [Buttari et al., 2009, Luszczek et al., 2014, Haidar et al., 2015]. Furthermore, the extent to which linear algebra tasks can be parallelized under this approach depends on the availability of parallel BLAS [Buttari et al., 2009].

Tile algorithms were developed to circumvent the limitations inherent to the block-cyclical approach. Tile algorithms propose a major rethinking of the way the matrix is decomposed in order to maximize task parallelism. The parallelization strategy of tile algorithms involves partitioning the matrix into small pieces called tiles such that their format allows for a more efficient way to access memory and is more suitable for asynchronous execution combined with dynamic scheduling of tasks across multiple platforms [Haidar et al., 2012, Dongarra et al., 2014]. Tile algorithms can be conveniently expressed using Directed
Acyclic Graphs (DAG) with nodes indicating the tasks and the edges specifying data movement and dependencies among tasks [Agullo et al., 2009, Buttari et al., 2009, Bosilca et al., 2012, Haidar et al., 2012, Haidar et al., 2015, Akbudak et al., 2017]. Several tile versions of some of the standard linear algebra algorithms, such as, Cholesky factorization [Buttari et al., 2009, Haidar et al., 2012], LU factorization [Buttari et al., 2009, Haidar et al., 2012, Dongarra et al., 2014], and QR factorization [Buttari et al., 2008, Buttari et al., 2009, Haidar et al., 2012], have been established in the literature.

§ THE MPCR PACKAGE INTERNAL DESIGN AND WORKFLOW

Historically, the R community has predominantly utilized 64-bit arithmetic for their computational tasks, aiming for a high level of accuracy that may only sometimes be necessary. In reality, arithmetic with a precision lower than 64-bit can offer the advantages of faster execution and a smaller memory footprint while maintaining the same level of accuracy as 64-bit arithmetic in numerous applications. This work introduces the MPCR package, designed to support 64-, 32-, and 16-bit arithmetic. The MPCR package also offers multi- and mixed-precision computing capabilities, aiming to optimize computational efficiency without compromising accuracy in many applications.

§.§ MPCR Internal Design

The MPCR package is designed with both performance and flexibility in mind. We rely on C++ in developing the backend of the package, and we port the code to R through the Rcpp package. We use C++ for various reasons: 1) C++ excels in scenarios demanding high performance and efficiency; 2) C++ is ideal for system-level programming that provides more flexibility in managing system resources, such as memory and processing power; 3) C++ has many tools and libraries that can help in building a robust software; 4) C++ can easily be integrated into different languages, including R and Python. A conceptual view of the software stack for the MPCR package is displayed in Figure <ref>. As can be seen in the figure, every function in MPCR is invoked from the R environment, where the Rcpp package manages such calls, effectively bridging them to the corresponding C++ functions in the C++ layer.

The MPCR package software stack.

The MPCR package offers three different types of operations: basic operations, linear algebra operations, and tile-based linear algebra operations. Detailed information about these functionalities is available in the package tutorial, accessible at <https://github.com/stsds/MPCR/blob/main/inst/doc/MPCR-manual.pdf>. The MPCR package supports numerous base functions in the R programming language, which can be readily applied to MPCR-objects. For instance, when performing the Cholesky factorization ( in base R) on a single-precision MPCR-object, the MPCR package utilizes its internal  function to execute the factorization in single-precision. Such design facilitates integrating the MPCR package into any existing R code, enabling single-precision operations with minimal modifications.

The MPCR package is designed to offer full extensibility with minimal effort. For instance, the package was built for easy support for integration of new precisions, addition of new linear operations, and inclusion of new tile-based operations. Moreover, the package uses the CMake Build system which enables easy integration for cross-platform and of external libraries. The CMake Build system also allows a full C++ testing environment, enabling a complete CI/CD pipeline for fast and stable development. Furthermore, MPCR comes with a fully documented R and C++ code, making it easy for any future development. The package was designed to be fully operational as a separate C++ module to
be extended and used in any C++ library. Additionally, MPCR uses template C++ functions offering no code redundancy. Lastly, MPCR has a fully organized code structure, offering easy and fast code navigation.

§.§ Workflow

This section outlines the operational workflow of the MPCR package, from a simple call at the R level and to its execution at the C++ level. Figure <ref> presents a high-level overview of the workflow involved in a typical MPCR operation.

The MPCR package workflow.

The process begins with specifying the appropriate precision for an MPCR object, which is then conveyed to the C++ code via the Rcpp package. The precision controller module in the C++ level assigns the correct precision to the MPCR objects and their corresponding operations. The next step involves allocating memory for the MPCR object, which depends on the structure of the object, e.g., vector or matrix, and the precision, e.g., 32-bit for single-precision and 64-bit for double-precision. Operations at different precision levels can be performed on non-tile-based MPCR objects. We refer to the non-tile-based operations collectively as native MPCR operations. On the other hand, tile-based MPCR objects utilize a different memory allocation approach, with tiles being the fundamental component of any such object. These specialized MPCR objects are equipped with tile-based linear algebra operations that can run sequentially or in parallel. Currently, three operations are supported: POTRF, GEMM, and TRSM. Once the object is allocated with particular precision, any call of a C++ function for this object is managed by the dispatcher module. This module is responsible for dispatching operations to the appropriate template function based on the precision of the inputs. To link the C++ function with the R  functions through the MPCR package, the MPCR R adapters module is used. This module acts as a mapper between the C++ and R environments. To provide a clearer description of each module, we provide examples of internal codes that illustrate the functionality of each one.

### Code snippet for an MPCR R adapters module example at the C++ level

**MPCR R adapter** is a module that allows R and C++  to interact. It acts as
a bridge between two incompatible interfaces.

std::vector <MPCR>
RSVD(MPCR *aInputA, const long    aNu, const long    aNv, const bool    aTranspose)

auto row = aInputA->GetNRow();
auto col = aInputA->GetNCol();
auto nv = aNv;
auto nu = aNu;

if (aNv < 0)
nv = std::min(row, col);

if (aNu < 0)
nu = std::min(row, col);

auto precision = aInputA->GetPrecision();
// Allocate three new objects for the SVD output.
auto d = new MPCR(precision);
auto u = new MPCR(precision);
auto v = new MPCR(precision);

SIMPLE_DISPATCH(precision, linear::SVD, *aInputA, *d, *u, *v, nu, nv,

std::vector <MPCR> output;
// Converting the output of C++ SVD funtion to a list, to match the shape of
R Return values.

return output;

### Code snippet for an MPCR dispatcher module example

**The Dispatcher** module main task is to choose the right signature for a
function according to the input and output precision. This can be done by the
help of the precision controller module.

**The Precision Controller** module main task is to decide the output precision
according to the input and the promotion strategy, and to decide the right
signature for the dispatcher.

MPCR* RRBind(MPCR *apInputA, MPCR *apInputB)

// Getting the precision for object A and B.
auto precision_a = apInputA->GetPrecision();
auto precision_b = apInputB->GetPrecision();
// Using the precision controller to decide the output precision.
auto output_precision = GetOutputPrecision(precision_a, precision_b);
// Allocating the output MPCR object with the right precision.
auto pOutput = new MPCR(output_precision);
// Using the precision controller to decide the operation signature for
the dispatcher to use.
auto operation_comb = GetOperationPrecision(precision_a, precision_b,
// Using the Operation signature to decide what template function to use.
DISPATCHER(operation_comb, basic::RowBind, *apInputA, *apInputB, *pOutput)
return pOutput;

§.§ Catch2

Catch2 is a popular C++ testing framework known for its comprehensive suite of tools that streamline writing and maintaining test cases for C++ code. As a header-only library, it enables C++ developers to easily create and manage their unit tests, offering a straightforward approach to testing.

In the MPCR package, Catch2 incorporates many unit tests, facilitating thorough testing of the functionality of the package as required. For example, integrating new features or enhancements into the package might unexpectedly impact existing code. Running these unit tests via Catch2 helps ensure software robustness throughout the development process. Furthermore, the MPCR package is designed to simplify the expansion of existing test cases, supporting future development and maintaining software reliability over time. Below is an example of using Catch2 in the MPCR package for the singular value decomposition function.

### Code snippets for C++ unit test for the SVD function using Catch2.

// Input Matrix Values.
vector <double> values = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0,
0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0,
0, 0, 0, 1, 1, 1;
// Create an MPCR Object.
MPCR a(values, FLOAT);
// Change MPCR Object to MPCR Matrix.
a.ToMatrix(9, 4);
// Validate values for testing.
vector <float> validate_values = 3.464102e+00, 1.732051e+00,
1.732051e+00, 1.922963e-16;
// Allocating three MPCR objects for the SVD output
MPCR d(FLOAT);
MPCR u(FLOAT);
MPCR v(FLOAT);
// This function will dispatch to SVD template function with the right precision.
SIMPLE_DISPATCH(FLOAT, linear::SVD, a, d, u, v, a.GetNCol(),
// Testing output d size.
REQUIRE(d.GetSize() == 4);
auto err = 0.001;
// Setting additional objects to test SVD Mathematically.
MPCR dd(9, 4, FLOAT);
for (auto i = 0; i < dd.GetSize(); i++)
dd.SetVal(i, 0);

for (auto i = 0; i < 4; i++)
dd.SetValMatrix(i, i, d.GetVal(i));

vector <double> temp_vals(81, 0);
MPCR uu(temp_vals, FLOAT);
uu.ToMatrix(9, 9);
for (auto i = 0; i < u.GetSize(); i++)
uu.SetVal(i, u.GetVal(i));

MPCR vv = v;

// Performing  A = U Sigma V^H
MPCR temp_one(FLOAT);
MPCR temp_two(FLOAT);

SIMPLE_DISPATCH(FLOAT, linear::CrossProduct, uu, dd, temp_one, false,

SIMPLE_DISPATCH(FLOAT, linear::CrossProduct, temp_one, vv, temp_two,

MPCR temp_three(FLOAT);
SIMPLE_DISPATCH(FLOAT, math::Round, temp_two, temp_three, 1);

SIMPLE_DISPATCH(FLOAT, math::PerformRoundOperation, temp_three,
temp_two, "abs");

// Checking if our SVD function values are valid mathematically
for (auto i = 0; i < a.GetSize(); i++)
REQUIRE(temp_two.GetVal(i) == a.GetVal(i));

§ OVERVIEW OF THE MPCR PACKAGE

To start using the MPCR package, it should be installed from the CRAN repository and imported into the R  environment by running these commands:

R> install.packages("MPCR")
R> library(MPCR)

While CRAN hosts the stable MPCR package, our GitHub repository contains the latest development features and fixes. To download the package from the GitHub repository:

R> library(devtools)
R> install.packages("https://github.com/stsds/MPCR")
R> library(MPCR)

§.§ Creating MPCR objects

To create an MPCR object, i.e., matrix/vector, one invokes the function  while specifying the desired  and  of the object. The default object created by the function  is a zero-vector with a number of elements equal to the  argument. To illustrate, the following codes create a single-precision MPCR zero-vector with 50 elements:

R> MPCR_object <- new(MPCR, 50, "single")
R> MPCR_object

MPCR Object: 32-Bit Precision

The function  can be used to check whether the newly created  is a vector or a matrix, e.g.:

R> MPCR_object$IsMatrix
\end{Sinput}
\begin{Soutput}
[1] FALSE
\end{Soutput}
\end{Schunk}

To display the values contained inside \verb|MPCR_object|, one can call the function \verb|PrintValues()| as follows:
\begin{Schunk}
\begin{Sinput}
R> MPCR_object$PrintValues()
\end{Sinput}
\begin{Soutput}
Vector Size : 50

[ 1 ]       0      0      0      0      0      0      0
[ 8 ]       0      0      0      0      0      0      0
[ 15 ]      0      0      0      0      0      0      0
[ 22 ]      0      0      0      0      0      0      0
[ 29 ]      0      0      0      0      0      0      0
[ 36 ]      0      0      0      0      0      0      0
[ 43 ]      0      0      0      0      0      0      0
[ 50 ]      0
\end{Soutput}
\end{Schunk}

Additionally, \verb|MPCR_object| can be reconfigured as a matrix with a certain number of rows and columns. Suppose one requires \verb|MPCR_object| to be a $5 \times 10$ matrix. The following codes perform such transformation:
\begin{Schunk}
\begin{Sinput}
R> MPCR_object$ToMatrix(5,10)
R> MPCR_object$IsMatrix
\end{Sinput}
\begin{Soutput}
[1] TRUE
\end{Soutput}
\begin{Sinput}
R> MPCR_object$PrintValues()
\end{Sinput}
\begin{Soutput}
Precision  : 32-Bit  Precision
Number of Rows : 5
Number of Columns : 10
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
\end{Soutput}
\end{Schunk}

Elements inside any \pkg{MPCR} object can be extracted and their values replaced by indicating the appropriate indices, e.g.:
\begin{Schunk}
\begin{Sinput}
R> MPCR_object[1,1] <- 1
R> MPCR_object[1,2] <- 2
R> MPCR_object[1,3] <- 3
R> MPCR_object[1,4] <- 4
R> MPCR_object[1,5] <- 5
R> MPCR_object$PrintValues()
\end{Sinput}
\begin{Soutput}
Precision  : 32-Bit  Precision
Number of Rows : 5
Number of Columns : 10
[    1    2    3    4    5    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
[    0    0    0    0    0    0    0    0    0    0    ]
\end{Soutput}
\end{Schunk}

The \pkg{MPCR} objects can also be converted to \proglang{R} objects by invoking the \pkg{MPCR} function
\verb|MPCR.ToNumericMatrix()| or \verb|MPCR.ToNumericVector()| as follows:
\begin{Schunk}
\begin{Sinput}
R> MPCR.ToNumericMatrix(MPCR_object)
\end{Sinput}
\begin{Soutput}
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]    1    2    3    4    5    0    0    0    0     0
[2,]    0    0    0    0    0    0    0    0    0     0
[3,]    0    0    0    0    0    0    0    0    0     0
[4,]    0    0    0    0    0    0    0    0    0     0
[5,]    0    0    0    0    0    0    0    0    0     0
\end{Soutput}
\begin{Sinput}
R> MPCR_vector <- new(MPCR, 5, "single")
R> MPCR.ToNumericVector(MPCR_vector)
\end{Sinput}
\begin{Soutput}
[1] 0 0 0 0 0
\end{Soutput}
\end{Schunk}

Another approach to constructing an \pkg{MPCR} object is to convert the usual \proglang{R} object to \pkg{MPCR} objects using the function \verb|as.MPCR()|. In the following, we create a $6 \times 6$ \proglang{R} matrix and convert it into a single-precision $6 \times 6$ \pkg{MPCR} matrix:
\begin{Schunk}
\begin{Sinput}
R> a <- matrix(1:36, 6, 6)
R> MPCR_matrix <- as.MPCR(a,nrow=6,ncol=6,precision="single")

R> a
\end{Sinput}
\begin{Soutput}
[,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    7   13   19   25   31
[2,]    2    8   14   20   26   32
[3,]    3    9   15   21   27   33
[4,]    4   10   16   22   28   34
[5,]    5   11   17   23   29   35
[6,]    6   12   18   24   30   36
\end{Soutput}
\begin{Sinput}
R> MPCR_matrix$PrintValues()
\end{Sinput}
\begin{Soutput}
Precision  : 32-Bit  Precision
Number of Rows : 6
Number of Columns : 6
[    1    7     13    19    25    31    ]
[    2    8     14    20    26    32    ]
[    3    9     15    21    27    33    ]
[    4    10    16    22    28    34    ]
[    5    11    17    23    29    35    ]
[    6    12    18    24    30    36    ]
\end{Soutput}
\end{Schunk}

\subsection{Creating MPCR-Tile objects}

An \pkg{MPCR}-Tile matrix can be created by calling the function \verb|new(MPCRTile, ...)| and supplying the values necessary to configure the desired tile matrix. The required arguments to the \verb|new(MPCRTile, ...)| function are listed as follows:
\begin{itemize}
\item \verb|rows|: a scalar value indicating the number of rows in the matrix;
\item \verb|cols|: a scalar value indicating the number of columns in the matrix;
\item \verb|rows_per_tile|: a scalar value indicating the number of rows in each tile;
\item \verb|cols_per_tile|: a scalar value indicating the number of columns in each tile;
\item \verb|values|: an \proglang{R} matrix or vector containing all the values that should be in the tile matrix;
\item \verb|precisions|: an \proglang{R} matrix or vector of strings, containing the precision type of each tile.
\end{itemize}

In the following, we demonstrate the usage of \verb|new(MPCRTile, ...)| by creating two tile matrices of different sizes.

Suppose we have a vector $\mathbf{z}$ with $100$ elements, and we want to construct a column tile matrix out of $\mathbf{z}$ such that it has 5 tiles, with 20 rows in each tile. Moreover, suppose we need the first two tiles to be encoded in double-precision while the rest are in single-precision. The following codes create such a column tile matrix:
\begin{Schunk}
\begin{Sinput}
R> set.seed(1234)
R> n <- 100
R> z <- runif(n, 0, 1)

R> num_tiles = 5
R> prec.z <- matrix(c(rep("double", 2), rep("single", num_tiles - 2)),
+    nrow = num_tiles, ncol = 1)

R> z.tile <- new(MPCRTile, rows = n, cols = 1, rows_per_tile = 20,
+    cols_per_tile = 1, values = z, precisions = prec.z)
R> z.tile
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 100
Number of Cols : 1
Number of Tiles : 5
Number of Tiles Per Row : 5
Number of Tiles Per Col : 1
Number of Rows Per Tile : 20
Number of Cols Per Tile : 1
\end{Soutput}
\end{Schunk}
To check that the first two tiles are in double-precision while the rest are in single-precision, we can select a specific tile and print out its attributes using the function \verb|MPCRTile.GetTile()| with three arguments, namely, \verb|matrix|, \verb|rowidx|, and \verb|colidx|. The argument \verb|matrix| is simply the tile matrix from which we want to draw out a specific tile. The arguments \verb|rowidx| and \verb|colidx| are, respectively, the row and column indices of the particular tile of interest. The following codes print out the attributes of each of the tiles in the tile matrix \verb|z.tile|:
\begin{Schunk}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,1,1)
\end{Sinput}
\begin{Soutput}
MPCR Object : 64-Bit Precision
\end{Soutput}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,2,1)
\end{Sinput}
\begin{Soutput}
MPCR Object : 64-Bit Precision
\end{Soutput}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,3,1)
\end{Sinput}
\begin{Soutput}
MPCR Object : 32-Bit Precision
\end{Soutput}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,4,1)
\end{Sinput}
\begin{Soutput}
MPCR Object : 32-Bit Precision
\end{Soutput}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,5,1)
\end{Sinput}
\begin{Soutput}
MPCR Object : 32-Bit Precision
\end{Soutput}
\end{Schunk}

To print the whole \verb|z.tile|, the command \verb|z.tile$MPCRTile.print()| or \verb|print(z.tile)| can be used. Similarly, the values of a specific tile, for instance, the third tile, can be displayed using the following command:
\begin{Schunk}
\begin{Sinput}
R> MPCRTile.GetTile(z.tile,3,1)$PrintValues()
\end{Sinput}
\end{Schunk}

Suppose we have a $100 \times 100$ matrix $\mathbf{M}$, and we want to construct a $5 \times 5$ tile matrix such that only the diagonal tiles are double-precision and the rest are single-precision. The following codes create such a tile matrix:
\begin{Schunk}
\begin{Sinput}
R> library(dplyr)
R> x <- seq(0, 1, length.out = 10)
R> locs <- expand.grid(x, x) %>% as.matrix()
R> n <- nrow(locs)
R> dist_mat <- as.matrix(dist(locs))
R> M <- exp(-dist_mat)

R> num_tiles_rowwise = 5
R> num_tiles_colwise = 5
R> num_tiles_total = num_tiles_rowwise * num_tiles_colwise
R> prec.M <- matrix(rep("single", num_tiles_total), nrow=num_tiles_rowwise,
+    ncol=num_tiles_colwise)
R> diag(prec.M) <- "double"

R> M.tile <- new(MPCRTile, rows=n, cols=n, rows_per_tile=n/num_tiles_rowwise,
+    cols_per_tile=n/num_tiles_colwise, M, prec.M)
R> M.tile
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 100
Number of Cols : 100
Number of Tiles : 25
Number of Tiles Per Row : 5
Number of Tiles Per Col : 5
Number of Rows Per Tile : 20
Number of Cols Per Tile : 20
\end{Soutput}
\end{Schunk}

\subsection{Operations on MPCR objects}
The \pkg{MPCR} package enables multi- and mixed-precision arithmetic. Arithmetic operations on \pkg{MPCR} objects can be carried out using the same symbols used to perform those operations on \proglang{R} objects. Moreover, with multi- and mixed-precision arithmetic, one can add or subtract matrices with different precisions such that the resulting \pkg{MPCR} object inherits the precision of the \pkg{MPCR} object with the highest precision. The following example adds a double-precision $2 \times 10$ matrix to a single-precision $2 \times 10$ matrix with a resulting double-precision $2 \times 10$ matrix.
\begin{Schunk}
\begin{Sinput}
R> s1 <- as.MPCR(1:20,nrow=2,ncol=10,"single")
R> s2 <- as.MPCR(21:40,nrow=2,ncol=10,"double")
R> x <- s1 + s2
R> typeof(x)
\end{Sinput}
\begin{Soutput}
MPCR Object : 64-Bit Precision
\end{Soutput}
\end{Schunk}

Other usual operations on \proglang{R} vectors and matrices can also be done on \pkg{MPCR} vectors and matrices by invoking the conventional function calls in \proglang{R}. For example, the transpose of \verb|MPCR_matrix| in the previous example can be obtained by \verb|t(MPCR_matrix)| and multiplying two \pkg{MPCR} matrices is done using the product operator \verb|%*%|, e.g., \verb|s1 %*% t(s2)|.

Linear algebra operations are also available and are called as they are in the base \proglang{R} package, such as \verb|chol()|, \verb|chol2inv()|, \verb|solve()|, and \verb|eigen()|. Furthermore, common functions used to evaluate matrices/vectors in \proglang{R} are also made available for \pkg{MPCR} objects such as \verb|cbind|, \verb|rbind|, \verb|diag|, \verb|min|, \verb|max|, etc., and can be invoked the same way as the base \proglang{R} package.
%%%%%%%%%%%% Sameh
\subsection{Operations on MPCR-Tile objects}
The \pkg{MPCR} package currently supports the three main tile-based linear algebra operations for tiled matrices, namely, tile-based matrix-matrix multiplication, Cholesky factorization, and solving of a triangular matrix equation.

\subsubsection{Tile-Based Matrix-Matrix Multiplication}
The \verb|MPCRTile.gemm()| function performs matrix-matrix multiplication of the form
\begin{equation}
\label{eqn:gemm}
\mathbf{C} = \alpha \mathbf{A} \mathbf{B} + \beta \mathbf{C}.
\end{equation}
The required arguments to the \verb|MPCRTile.gemm()| function are the following:
\begin{itemize}
\item \verb|a|: an \pkg{MPCR}-Tile matrix representing matrix $\mathbf{A}$ in Equation~(\ref{eqn:gemm});
\item \verb|b|: an \pkg{MPCR}-Tile matrix representing matrix $\mathbf{B}$ in Equation~(\ref{eqn:gemm});
\item \verb|c|: an \pkg{MPCR}-Tile matrix representing matrix $\mathbf{C}$ in Equation~(\ref{eqn:gemm});
\item \verb|transpose_a|: a flag to indicate whether the matrix in argument \verb|a| should be transposed before performing matrix multiplication; takes in the values \verb|TRUE| or \verb|FALSE| (default);
\item \verb|transpose_b|: a flag to indicate whether the matrix in argument \verb|b| should be transposed before performing matrix multiplication; takes in the values \verb|TRUE| or \verb|FALSE| (default);
\item \verb|alpha|: a scalar value representing $\alpha$ in Equation~(\ref{eqn:gemm});
\item \verb|beta|: a scalar value representing $\beta$ in Equation~(\ref{eqn:gemm});
\item \verb|num_threads|: an integer that indicates the number of threads to run using \pkg{openmp}; the default value is 1, which means serial computations with no parallelization.
\end{itemize}

The results after performing the matrix-matrix multiplication are stored in the matrix used as an input to the argument \verb|c|. This means that the original data saved into the input matrix to \verb|c| will be overwritten.

Note that the rules of basic matrix-matrix multiplication must also be followed for the tile-based version. This means that the dimensions in terms of tiles of $\mathbf{A}$, $\mathbf{B}$, and $\mathbf{C}$ must be compatible. That is, the number of tiles columnwise of $\mathbf{A}$ must be equal to the number of tiles rowwise of $\mathbf{B}$ and the resulting dimensions of their product in terms of tiles must be equal to the tile dimensions of $\mathbf{C}$.

The following codes demonstrate how to call the \verb|MPCRTile.gemm()| function:
\begin{Schunk}
\begin{Sinput}
R> A <- matrix(1:24, nrow=4, ncol=6)
R> B <- matrix(0, nrow=6, ncol=1)
R> C <- matrix(1, nrow=4, ncol=1)

R> num_tiles_rowwise_A <- 2
R> num_tiles_colwise_A <- 1
R> num_tiles_total_A <- num_tiles_rowwise_A * num_tiles_colwise_A

R> num_tiles_rowwise_B <- 1
R> num_tiles_colwise_B <- 1
R> num_tiles_total_B <- num_tiles_rowwise_B * num_tiles_colwise_B

R> num_tiles_rowwise_C <- 2
R> num_tiles_colwise_C <- 1
R> num_tiles_total_C <- num_tiles_rowwise_C * num_tiles_colwise_C

R> rows_per_tile_A <- nrow(A) / num_tiles_rowwise_A
R> cols_per_tile_A <- ncol(A) / num_tiles_colwise_A

R> rows_per_tile_B <- nrow(B) / num_tiles_rowwise_B
R> cols_per_tile_B <- ncol(B) / num_tiles_colwise_B

R> rows_per_tile_C <- nrow(C) / num_tiles_rowwise_C
R> cols_per_tile_C <- ncol(C) / num_tiles_colwise_C

R> prec.A <- matrix(rep("single", num_tiles_total_A),
+    nrow=num_tiles_rowwise_A, ncol=num_tiles_colwise_A)
R> prec.B <- matrix(rep("single", num_tiles_total_B),
+    nrow=num_tiles_rowwise_B, ncol=num_tiles_colwise_B)
R> prec.C <- matrix(rep("single", num_tiles_total_C),
+    nrow=num_tiles_rowwise_C, ncol=num_tiles_colwise_C)

R> A.tile <- new(MPCRTile, rows=nrow(A), cols=ncol(A),
+    rows_per_tile=rows_per_tile_A, cols_per_tile=cols_per_tile_A, A, prec.A)
R> B.tile <- new(MPCRTile, rows=nrow(B), cols=ncol(B),
+    rows_per_tile=rows_per_tile_B, cols_per_tile=cols_per_tile_B, B, prec.B)
R> C.tile <- new(MPCRTile, rows=nrow(C), cols=ncol(C),
+    rows_per_tile=rows_per_tile_C, cols_per_tile=cols_per_tile_C, C, prec.C)
\end{Sinput}
\end{Schunk}

Before proceeding with the matrix-matrix multiplication, we first print the values in the tile matrix \verb|C.tile| to verify that it is indeed a $4 \times 1$ column matrix of 1's.
\begin{Schunk}
\begin{Sinput}
R> print(C.tile)
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 4
Number of Cols : 1
Number of Tiles : 2
Number of Tiles Per Row : 2
Number of Tiles Per Col : 1
Number of Rows Per Tile : 2
Number of Cols Per Tile : 1

[	             1	             ]
[	             1	             ]
[	             1	             ]
[	             1	             ]
\end{Soutput}
\end{Schunk}
Suppose that $\alpha = 1$ and $\beta = 0.5$. The next line of code performs the matrix-matrix multiplication and stores the results in \verb|C.tile|:
\begin{Schunk}
\begin{Sinput}
R> MPCRTile.gemm(A.tile, B.tile, C.tile, transpose_a=F, transpose_b=F,
+    alpha=1, beta=0.5, num_threads=4)
\end{Sinput}
\end{Schunk}

Based on the values of $\mathbf{A}, \mathbf{B}, \mathbf{C}, \alpha, $ and $\beta$, we expect the answer to be a $4 \times 1$ column matrix of 0.5's. To verify that the arithmetic is correct, we check the results by printing the updated values in the tile matrix \verb|C.tile| as follows:
\begin{Schunk}
\begin{Sinput}
R> print(C.tile)
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 4
Number of Cols : 1
Number of Tiles : 2
Number of Tiles Per Row : 2
Number of Tiles Per Col : 1
Number of Rows Per Tile : 2
Number of Cols Per Tile : 1

[	           0.5	             ]
[	           0.5	             ]
[	           0.5	             ]
[	           0.5	             ]
\end{Soutput}
\end{Schunk}

\subsubsection{Tile-Based Cholesky Factorization}

The \pkg{MPCR} package has a tile-based version of the Cholesky factorization for tile matrices. That is, a positive definite tile matrix $\mathbf{A}$ can be decomposed into a product of a unique lower triangular tile matrix $\mathbf{L}$ and its transpose $\mathbf{L}^{\top}$ as follows:
\begin{equation}
\label{eqn:chol}
\mathbf{A} = \mathbf{L} \mathbf{L}^{\top}.
\end{equation}

The tile-based Cholesky factorization can be performed by calling the function \verb|chol()| and providing the following required inputs:
\begin{itemize}
\item \verb|x|: a positive definite \pkg{MPCR}-Tile matrix representing $\mathbf{A}$ in Equation~(\ref{eqn:chol});
\item \verb|overwrite_input|: a flag to indicate whether the resulting lower triangular tile matrix representing $\mathbf{L}$ in Equation~(\ref{eqn:chol}) is stored in the tile matrix used as input to the argument \verb|x| above; takes in the values \verb|TRUE| (default) or \verb|FALSE|;
\item \verb|num_threads|: an integer that indicates the number of threads to run using \pkg{openmp}; the default value is 1, which means serial computations with no parallelization.
\end{itemize}

The result after performing \verb|chol()| is a lower triangular tile matrix. Depending on the specification to the argument \verb|overwrite_input|, the resulting lower triangular tile matrix may be stored in the matrix used as the input to the argument \verb|x| or in another user-defined variable. With the $4 \times 4$ covariance matrix example below, we show how to perform a tile-based Cholesky factorization with \pkg{MPCR}:
\begin{Schunk}
\begin{Sinput}
R> library(dplyr)
R> x <- seq(0, 1, length.out = 2)
R> locs <- expand.grid(x, x) %>% as.matrix()
R> n <- nrow(locs)
R> dist_mat <- as.matrix(dist(locs))
R> M <- exp(-dist_mat)

R> num_tiles_rowwise = 2
R> num_tiles_colwise = 2
R> num_tiles_total = num_tiles_rowwise * num_tiles_colwise

R> rows_per_tile <- nrow(M) / num_tiles_rowwise
R> cols_per_tile <- ncol(M) / num_tiles_colwise

R> prec.M <- matrix(rep("single", num_tiles_total), nrow=num_tiles_rowwise,
+    ncol=num_tiles_colwise)
R> diag(prec.M) <- "double"

R> M.tile <- new(MPCRTile, rows=n, cols=n, rows_per_tile=rows_per_tile,
+    cols_per_tile=cols_per_tile, M, prec.M)
R> M.tile_chol <- chol(M.tile, overwrite=F, num_threads=4)
R> print(M.tile_chol)
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 4
Number of Cols : 4
Number of Tiles : 4
Number of Tiles Per Row : 2
Number of Tiles Per Col : 2
Number of Rows Per Tile : 2
Number of Cols Per Tile : 2

[	  1	           0	           0	           0  ]
[  0.3678794	   0.9298735	           0	           0  ]
[  0.3678795	   0.1159098	   0.9226211	           0  ]
[  0.2431167	   0.2994405	   0.2641753	   0.8839915  ]
\end{Soutput}
\end{Schunk}

\subsubsection{Tile-Based Solving of a Triangular Matrix Equation}

The \pkg{MPCR} package supports the tile-based linear algebra routine that solves for $\mathbf{X}$ in the following triangular system of equations:
\begin{equation}
\label{eqn:triangular_system_equations}
\mathbf{A} \mathbf{X} = \alpha \mathbf{B} \quad \text{or} \quad \mathbf{X} \mathbf{A} = \alpha \mathbf{B},
\end{equation}
where $\mathbf{A}$ is any triangular matrix, e.g., the matrix $\mathbf{L}$ in Equation~(\ref{eqn:chol}). The routine can be accessed by calling the function \verb|MPCRTile.trsm()|
and admitting values to the following parameters:
\begin{itemize}
\item \verb|a|: a triangular \pkg{MPCR}-Tile matrix representing $\mathbf{A}$ in Equation~(\ref{eqn:triangular_system_equations});
\item \verb|b|: an \pkg{MPCR}-Tile matrix representing $\mathbf{B}$ in Equation~(\ref{eqn:triangular_system_equations});
\item \verb|side|: a string that indicates the position of $\mathbf{A}$ relative to $\mathbf{X}$; `\verb|R|' for right side or `\verb|L|' for left side;
\item \verb|upper_triangle|: a flag that indicates which part of $\mathbf{A}$ is nonzero; \verb|TRUE| if $\mathbf{A}$ is an upper triangular matrix and \verb|FALSE| otherwise;
\item \verb|transpose|: a flag to indicate whether the matrix in argument \verb|a| should be transposed before performing the routine; takes in the values \verb|TRUE| or \verb|FALSE|;
\item \verb|alpha|: a scalar value representing $\alpha$ in Equation~(\ref{eqn:triangular_system_equations}).
\end{itemize}

The function \verb|MPCRTile.trsm()| returns $\mathbf{X}$ and is stored in the matrix used as the input to the argument \verb|b|. To illustrate the use of the \verb|MPCRTile.trsm()| function, suppose that the resulting lower triangular tile matrix in the Cholesky factorization in the previous example is $\mathbf{A}$. In the following codes, we define $\mathbf{B}$ as an identity tile matrix and try to solve for $\mathbf{X}$:
\begin{Schunk}
\begin{Sinput}
R> B <- diag(4)
R> prec.B <- matrix(rep("single", num_tiles_total), nrow = num_tiles_rowwise,
+    ncol = num_tiles_colwise)
R> B.tile <- new(MPCRTile, rows = n, cols = n, rows_per_tile = rows_per_tile,
+    cols_per_tile = cols_per_tile, B, prec.B)
R> print(B.tile)
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 4
Number of Cols : 4
Number of Tiles : 4
Number of Tiles Per Row : 2
Number of Tiles Per Col : 2
Number of Rows Per Tile : 2
Number of Cols Per Tile : 2

[	       1	       0	       0	       0        ]
[	       0	       1	       0	       0        ]
[	       0	       0	       1	       0        ]
[	       0	       0	       0	       1        ]
\end{Soutput}
\begin{Sinput}
R> MPCRTile.trsm(a=M.tile_chol, b=B.tile, side='L', upper_triangle=F,
+    transpose=F, alpha=1)
R> print(B.tile)
\end{Sinput}
\begin{Soutput}
----------- MPCRTile Object ----------
Number of Rows : 4
Number of Cols : 4
Number of Tiles : 4
Number of Tiles Per Row : 2
Number of Tiles Per Col : 2
Number of Rows Per Tile : 2
Number of Cols Per Tile : 2

[	  1	           0	           0	           0  ]
[  -0.3956231	   1.075415	           0	           0  ]
[  -0.3490305	  -0.1351055	   1.083869	           0  ]
[  -0.03670389	 -0.3239073	  -0.3239073	   1.131233  ]
\end{Soutput}
\end{Schunk}

\section{Performance Assessment} \label{sec:performance}

In this section, we aim to evaluate the performance of single-, double-, and mixed-precisions operations within the \pkg{MPCR} package at two different precision levels, i.e., 64-bit and 32-bit. This assessment includes computationally intensive tasks, i.e., non-tile and tile-based linear algebra operations. We focus on linear algebra operations that are time-intensive and those that necessitate modifications for enhanced speed.

Figure~\ref{fig:perf1} illustrates the performance of \pkg{MPCR} double- and single-precision operations compared to the \proglang{R} double-precision operations. The left part of the figure shows the speed enhancements achieved by employing \pkg{MPCR} double-precision across three different matrix dimensions, specifically $500$, $5{,}000$, and $15{,}812$, utilizing three linear algebra procedures: \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()|. The results indicate that \pkg{MPCR} double-precision operations can accelerate \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()| by $495$X, $745$X, and $612$X, respectively, compared to the \proglang{R} double-precision operations. The performance improvement is primarily attributed to using highly optimized BLAS/LAPACK libraries like Intel MKL and OpenBLAS, which are automatically used in the \pkg{MPCR} package. In contrast, by default, \proglang{R} links to RBLAS, which is not as well optimized and lacks support for single-precision in its current version. The right figure demonstrates the efficiency of \pkg{MPCR} single-precision operations relative to \proglang{R} double-precision operations. The illustration reveals that \pkg{MPCR} single-precision for \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()| significantly outperforms double-precision operations in \proglang{R}, achieving speedup of $1{,}026$X, $1{,}532$X, and $1{,}303$X, respectively.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.48\textwidth]{figures/double-prec-potrf-gemm-trsm-new.png}
\hfill
\includegraphics[width=0.48\textwidth]{figures/single-prec-potrf-gemm-trsm-new.png}
\cprotect\caption{Speedup of double-precision (left) and single-precision (right) linear algebra operations, i.e., \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()| in the \pkg{MPCR} package against double-precision operations in \proglang{R}.}
\label{fig:perf1}
\end{figure*}

The \pkg{MPCR} also provides linear scalability with different matrix sizes and consistent speedup when applying single-precision compared to double-precision. Figure~\ref{fig:perf2} shows the execution times for individual \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()| linear algebra operations with two matrix sizes, $15{,}812$ and $50{,}000$, in single-precision and double-precision formats. Consistent with expectations, a speedup of approximately $2$X is shown when precision is lowered from double to single. This increase in efficiency is scalable, as illustrated in the figure.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.48\textwidth]{figures/exec-time-potrf-gemm-trsm.png}
\cprotect\caption{Execution time of three different linear algebra operations, i.e., \verb|chol()|, \verb|crossprod()|, and \verb|backsolve()| using single-precision and double-precision operations with two different matrix sizes $15{,}812$ and $50{,}000$.}
\label{fig:perf2}
\end{figure*}

The \pkg{MPCR} package offers two distinct approaches for implementing the Cholesky factorization operation, i.e., \verb|chol()|. This operation can be executed either in the same memory space, referred to as in-place, or using new memory allocations, referred to as out-place. This is an example of performing the Cholesky factorization on a tile-based \pkg{MPCR} matrix,

\begin{Schunk}
\begin{Sinput}
R>  a <- matrix(c(1.21, 0.18, 0.13, 0.41, 0.06, 0.23,
+    0.18, 0.64, 0.10, -0.16, 0.23, 0.07,
+    0.13, 0.10, 0.36, -0.10, 0.03, 0.18,
+    0.41, -0.16, -0.10, 1.05, -0.29, -0.08,
+    0.06, 0.23, 0.03, -0.29, 1.71, -0.10,
+    0.23, 0.07, 0.18, -0.08, -0.10, 0.36), 6, 6)
R> b <- c("float", "double", "float", "float",
+    "double", "double", "float", "float",
+    "double")

R> chol_mat <- new(MPCRTile, 6, 6, 2, 2, a, b)
R> chol_values <- chol(chol_mat, overwrite_input = FALSE, num_thread = 8)
+    # out-place chol()
R> chol(chol_mat, overwrite_input = TRUE, num_thread = 8)
+    # in-place chol()
\end{Sinput}
\end{Schunk}

Figure~\ref{fig:perf3} shows the performance when using both calls of the \verb|chol()| operation. As shown, in-place \verb|chol()| outperform out-place \verb|chol()| with different matrix size. The subfigure on the left displays the performance of the single-precision \verb|chol()|, whereas the one on the right illustrates the performance when using double-precision \verb|chol()|.

\begin{figure*}[t!]
\centering
\includegraphics[width=0.48\textwidth]{figures/sp-in-out-place.png}
\hfill
\includegraphics[width=0.48\textwidth]{figures/dp-in-out-place.png}
\caption{Performance comparison between using in-place and out-place tile-based Cholesky factorization across various matrix sizes.}
\label{fig:perf3}
\end{figure*}

%Since The \pkg{MPCR} package relies on

%    \centering
%    \includegraphics[width=0.48\textwidth]{figures/multithreads.png}
%        \caption{Performance of tile-based CHOL(), GEMM(), TRSM() using different numbers of threads (up to 32 threads) for 50K matrix size.}
%    \label{fig:perf4}

\section{Applications} \label{sec:applications}

This section offers a range of practical examples demonstrating the use of the \pkg{MPCR} package for executing computationally intensive tasks in Frequentist/Bayesian statistical use cases. These examples, drawn from existing applications, illustrate the straightforward integration of different \pkg{MPCR} functions. Users dealing with demanding workflows will discover that the multi- and mixed-precision functions provided by \pkg{MPCR} significantly enhance performance while preserving the workflow of the original code. This enhancement can dramatically reduce experiment waiting times from hours to minutes or, in some cases, even seconds.

\subsection{Metropolis-Hastings Algorithm in High Dimensions}

The \pkg{MPCR} package offers efficient linear algebra operations at different precision levels, as demonstrated in the previous section. One example demonstrating the efficiency of our package is its usage in Markov Chain Monte Carlo (MCMC) simulations, which are essential for many statistical applications such as the Metropolis-Hastings (MH) algorithm.

The MH algorithm is routinely used to generate samples from a target distribution, denoted $p(\mathbf{z})$, that may otherwise be difficult to sample from [Metropolis \emph{et~al.}, 1953, Hastings, 1970]. The MH algorithm generates a Markov Chain, whose stationary distribution is the target distribution. This means that, in the long run, the samples from the Markov chain will look like the samples from the target distribution, i.e., we use a Markov chain to generate a sequence of vectors, denoted $(\mathbf{z}_0, \mathbf{z}_1, \ldots, \mathbf{z}_I)$, such that as $I \rightarrow \infty, \mathbf{z}_I \sim p(\mathbf{z})$.

The Markov Chain vectors are simulated from a distinct distribution called the proposal (conditional) distribution, represented as $q(\mathbf{z}|\mathbf{z}^{\text{current}})$. The term ``conditional'' in the context of proposal distributions in the MH algorithm is simply an indication that the mean of the proposal distribution will be set to the value of the conditioning parameter, which is the current vector in the chain, $\mathbf{z}^{\text{current}}$. The proposal distribution is chosen such that it is a distribution from which realizations are easy to simulate and that its realizations can converge to the target distribution. The two most commonly used types of proposal distributions are (a) independent and (b) random-walk. \textit{Independent} proposal distributions are those that do not depend on $\mathbf{z}^{\text{current}}$, the most current accepted sample in the chain. This means that $q(\mathbf{z} | \mathbf{z}^{\text{current}}) = q(\mathbf{z}).$ For instance, $q(\mathbf{z} | \mathbf{z}^{\text{current}}) = \mathcal{N}_{n} (\boldsymbol{\mu}, \boldsymbol{\Sigma})$. This approach is recommended when the proposal distribution is known to approximate the target distribution well. On the other hand, proposal distributions under the \textit{random walk} approach are those that depend on $\mathbf{z}^{\text{current}}$. For example, $q(\mathbf{z} | \mathbf{z}^{\text{current}}) = \mathcal{N}_{n} (\mathbf{z}^{\text{current}}, \boldsymbol{\Sigma})$ or $q(\mathbf{z} | \mathbf{z}^{\text{current}}) = \mathcal{U} (\mathbf{z}^{\text{current}} - \epsilon, \mathbf{z}^{\text{current}} + \epsilon)$, where $\mathcal{U}~(\cdot)$ denotes the uniform distribution [Hastings, 1970].

The MH algorithm often struggles in high dimensions especially when incorporating naive proposal distributions such as those listed above [Au and Beck, 2001]. A more advanced and effective alternative in such scenarios is the Metropolis-Adjusted Langevin Algorithm (MALA) [Xifara \emph{et~al.}, 2014, S{\"a}rkk{\"a} \emph{et~al.}, 2021]. The key innovation in MALA is the choice of the proposal distribution, the most critical component for the MH algorithm to be successful. The quality of the candidate samples drawn from the independent and random walk proposal distributions earlier are not great since they tend to be rejected by the algorithm. Under MALA, the suggested proposal distribution is
\begin{equation}
\label{eqn:MALA_proposal_distribution}
q(\mathbf{z} | \mathbf{z}^{\text{current}}) = \mathcal{N}_{n}\left[\mathbf{z}^{\text{current}} + \frac{h}{2} \mathbf{M} \nabla \log \{p(\mathbf{z}^{\text{current}})\}, h \mathbf{M} \right],
\end{equation}
where $h$ is a user-defined step-size, $\nabla \log \{p(\mathbf{z})\}$ denotes the gradient of $\log \{p(\mathbf{z})\}$, and $\mathbf{M}$ is a preconditioning matrix. The strength of MALA lies in incorporating the gradient of the target when drawing the next candidate sample [Durmus \emph{et~al.}, 2017].

When the number of iterations, $I$, is large, and the dimensions of the target distribution are large, the MH algorithm can be computationally burdensome to perform. The standard MH algorithm requires computing the likelihood ratio in (\ref{eqn:log_acceptance_probability}) at every iteration, incurring a computational cost of
$O(n)$ per step, which is prohibitive for large $n$ [Cornish \emph{et~al.}, 2019]. This section shows how \pkg{MPCR} has features that can ameliorate most computational costs surrounding the MH algorithm. First, \pkg{MPCR} provides multi-precision support, so vectors and matrices can be stored in double or single-precision. Users can choose precision types for particular vectors and matrices, ultimately reducing the costs of performing the MH algorithm. Second, the \pkg{MPCR} package is capable of fast inversion of a high-dimensional covariance matrix and fast vector-matrix operations involved when evaluating the Gaussian log-likelihoods in (\ref{eqn:log_acceptance_probability}).

The following shows the sequential step in MALA using \pkg{MPCR} objects and routines. Additionally, for comparison, we run the same MALA using only \proglang{R} objects.
\begin{enumerate}
\item Specify the precision. Choose among \verb|R-Double|, \verb|MPCR-Double|, and \verb|MPCR-Single|.
\begin{Schunk}
\begin{Sinput}
R> precision = 'MPCR-Single'
\end{Sinput}
\end{Schunk}

\item Specify the target distribution. \\
Suppose we want to generate a 2-dimensional (2D) spatial Gaussian random field on a $120 \times 120$ grid on the unit square with an exponential spatial covariance function model parameterized by $\boldsymbol{\theta} = (a, \sigma^2)^\top$, where $a > 0, \sigma^2>0$ are the spatial range and variance parameters, respectively. That is, our target distribution is $\mathcal{N}_{n} \{\boldsymbol{\mu}, \boldsymbol{\Sigma}(\boldsymbol{\theta})\}$, such that $n$ is the number of locations, e.g., $n = 14$,$400$, and we want to simulate an $n$-dimensional vector $\mathbf{Z} = \{ Z (\mathbf{s}_1), \ldots, Z (\mathbf{s}_{n}) \}^{\top} \in \mathbb{R}^{n}$ from the target distribution with $\boldsymbol{\mu} = \mathbf{0} \in \mathbb{R}^{n \times 1}$, and $\boldsymbol{\Sigma}(\boldsymbol{\theta}) \in \mathbb{R}^{n \times n}$ whose value at index~$(i, j)$ is
\begin{equation}
\label{eqn:univariate_exponential}
\text{cov} \{Z (\mathbf{s}_i), Z (\mathbf{s}_j)\} = \sigma^2 \exp \left( -\frac{\| \mathbf{s}_i - \mathbf{s}_j \|}{a}\right).
\end{equation}
Here, $Z (\mathbf{s})$ is the observation at location $\mathbf{s} \in \mathbb{R}^{d}$ and $d = 2$.
\begin{enumerate}[i.]
\item Code the mean vector of the target distribution.
\begin{Schunk}
\begin{Sinput}
R> M <- 120
R> n.loc <- M*M
R> locs <- cbind(rep(0:(M-1), M)/(M-1), rep(0:(M-1), each=M)/(M-1))
R> mu = rep(0, nrow(locs)) #target mean vector

R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    mu <- as.MPCR(mu, nrow=n.loc, ncol=1, precision='single')
+ }
\end{Sinput}
\end{Schunk}

\item Code the covariance matrix of the target distribution.
\begin{Schunk}
\begin{Sinput}
R> dist0 <- as.matrix(dist(locs)) # distance matrix
R> sig <- exp(-dist0/0.5) #target covariance matrix

R> if(precision == 'MPCR-Double'){
+    sig <- as.MPCR(as.matrix(sig), nrow=n.loc, ncol=n.loc,
+      precision='double')
+    }else if(precision == 'MPCR-Single'){
+    sig <- as.MPCR(as.matrix(sig), nrow=n.loc, ncol=n.loc,
+      precision='single')
+ }
\end{Sinput}
\end{Schunk}

\end{enumerate}

\item Specify the proposal distribution. \\
We code the MALA-suggested proposal distribution in Equation~(\ref{eqn:MALA_proposal_distribution}).
\begin{enumerate}[i.]
\item Set values for the parameters in MALA.
\begin{Schunk}
\begin{Sinput}
R> I = 2000
R> h = 0.01
R> pre_M = exp(-dist0/0.05)

R> if(precision == 'MPCR-Double'){
+    pre_M <- as.MPCR(as.matrix(pre_M), nrow=n.loc, ncol=n.loc,
+      precision='double')
+    }else if(precision == 'MPCR-Single'){
+      pre_M <- as.MPCR(as.matrix(pre_M), nrow=n.loc, ncol=n.loc,
+        precision='single')
+ }

R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    pre_M_scaled_inverse <- solve(pre_M$PerformMult(h))
+    L <- t(chol(pre_M$PerformMult(h)))
+    }else{
+      pre_M_scaled_inverse <- solve(h*pre_M)
+      L <- t(chol(h * pre_M))
+ }
\end{Sinput}
\end{Schunk}

\item Create the functions that will compute $\nabla \log \{p(\mathbf{z}^{\text{current}})\}$ in Equation~(\ref{eqn:MALA_proposal_distribution}).
\begin{Schunk}
\begin{Sinput}
R> sig.inv <- solve(sig)
R> gradient <- function(theta){
+    sig.inv %*% theta
+ }

R> gradientStep <- function(theta, t){
+    if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+      theta - (pre_M %*% gradient(theta))$PerformMult(t)
+    }else{
+      theta - t * pre_M %*% gradient(theta)
+    }
+ }
\end{Sinput}
\end{Schunk}
\end{enumerate}

\item Create a matrix that will store the values of the accepted samples, $\mathbf{z}^{\text{current}}$.
\begin{Schunk}
\begin{Sinput}
R> z_trace <- matrix(, nrow=n.loc, ncol=I)
\end{Sinput}
\end{Schunk}

\item Initialize the chain by setting an arbitrary value to $\mathbf{z}_0 \in \mathbb{R}^n$.
\begin{Schunk}
\begin{Sinput}
R> set.seed(1234)
R> z_init <- runif(n.loc)
R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    z_init <- as.MPCR(z_init, nrow=n.loc, ncol=1, precision='single')
+ }
\end{Sinput}
\end{Schunk}

\item Set the value of $\mathbf{z}^{\text{current}}$ to $\mathbf{z}_0$.
\begin{Schunk}
\begin{Sinput}
R> z_current <- z_init
R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    z_trace[, 1] <- MPCR.ToNumericVector(z_current)
+    }else{
+      z_trace[, 1] <- z_current
+ }
\end{Sinput}
\end{Schunk}

\item For $i = 2, \ldots, I,$ iterate through the steps of MALA as follows:

\begin{enumerate}[i.]
\item Generate the next vector in the chain $\mathbf{z}_i \in \mathbb{R}^n$ by simulating a candidate vector $\mathbf{z}^{\text{candidate}}$ from (\ref{eqn:MALA_proposal_distribution}).
\begin{Schunk}
\begin{Sinput}
R> i = 2 #Loop this for i=2,...,I.
R> set.seed(i)
R> x <- rnorm(n.loc)

R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    x <- as.MPCR(x, nrow=n.loc, ncol=1, precision='single')
+ }

R> z_candidate <- gradientStep(z_current,0.5*h) + L %*% x
\end{Sinput}
\end{Schunk}

\item Compute the ``log acceptance probability'', $\log A$, using the formula:
\begin{eqnarray}
\label{eqn:log_acceptance_probability}
\log A & = & \min [0, \log \{p(\mathbf{z}^{\text{candidate}})\} - \log \{p(\mathbf{z}^{\text{current}})\} \nonumber \\
&&  + \log \{q(\mathbf{z}^{\text{current}} | \mathbf{z}^{\text{candidate}})\} - \log \{q(\mathbf{z}^{\text{candidate}} | \mathbf{z}^{\text{current}})\} ].
\end{eqnarray}
\begin{Schunk}
\begin{Sinput}
R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+    p_candidate <- MPCR.ToNumericVector((t(z_candidate-mu) %*%
+      sig.inv %*% (z_candidate-mu))$PerformMult(-0.5))
+    p_current <- MPCR.ToNumericVector((t(z_current-mu) %*%
+      sig.inv %*% (z_current-mu))$PerformMult(-0.5))

+    q_current <- MPCR.ToNumericVector((t(z_current-
+      gradientStep(z_candidate,0.5*h)) %*% pre_M_scaled_inverse %*%
+      (z_current-gradientStep(z_candidate,0.5*h)))$PerformMult(-0.5))
+    q_candidate <- MPCR.ToNumericVector((t(z_candidate-
+      gradientStep(z_current,0.5*h)) %*% pre_M_scaled_inverse %*%
+      (z_candidate-gradientStep(z_current,0.5*h)))$PerformMult(-0.5))
+    }else{
+      p_candidate <- -0.5*(t(z_candidate-mu) %*% sig.inv %*%
+      (z_candidate-mu))
+      p_current <- -0.5*(t(z_current-mu) %*% sig.inv %*%
+      (z_current-mu))

+      q_current <- -0.5*(t(z_current-gradientStep(z_candidate,0.5*h))
+      %*% pre_M_scaled_inverse %*% (z_current-gradientStep(z_candidate,
+        0.5*h)))
+      q_candidate <- -0.5*(t(z_candidate-gradientStep(z_current,0.5*h))
+      %*% pre_M_scaled_inverse %*% (z_candidate-gradientStep(z_current,
+        0.5*h)))
+ }

R> log.ratio <- min(0, p_candidate - p_current + q_current - q_candidate)
\end{Sinput}
\end{Schunk}

\item The acceptance or rejection decision of $\mathbf{z}^{\text{candidate}}$ as a new member of the chain can be performed using the following criteria:
\begin{equation*}
\mathbf{z}_i = \begin{cases}
\mathbf{z}^{\text{candidate}}, \quad \quad \; \text{ if }\log u \leq \log A, \\
\mathbf{z}^{\text{current}}, \quad \quad \quad \text{ otherwise.}
\end{cases}
\end{equation*}
Generate a uniformly distributed random number, denoted $u$, between 0 and 1.
\begin{Schunk}
\begin{Sinput}
R> if(log(runif(1)) < log.ratio) {
+      z_current <- z_candidate
+ }

R> if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+      z_trace[, i] <- MPCR.ToNumericVector(z_current)
+    }else{
+      z_trace[, i] <- z_current
+ }
\end{Sinput}
\end{Schunk}

\end{enumerate}
\end{enumerate}

\begin{figure}[t!]
\centering
\includegraphics[width=0.5\textwidth]{figures/MPCR_JSS_application_mh_result_time.pdf}
\caption{Execution time of MH algorithm using different types of precision for different number of variables.}
\label{fig:MPCR_JSS_application_mh_result_time}
\end{figure}

\begin{figure*}[t!]
\centering
\includegraphics[width=\textwidth]{figures/mh_results.png}
\caption{Initial $\mathbf{z}_0$ and the simulated values after the 2000th iteration using MALA.}
\label{fig:mh_results}
\end{figure*}

Figure~\ref{fig:MPCR_JSS_application_mh_result_time} shows the execution time when performing the naive MH algorithm with an independent proposal distribution $\mathcal{U}(\cdot)$ under different execution scenarios, i.e., \proglang{R}, \pkg{MPCR}-Double and \pkg{MPCR}-Single. The plots illustrates the overall speedup that can be obtained using \pkg{MPCR} on the simulation tasks under different numerical precisions. Figure~\ref{fig:mh_results} shows the 2D spatial fields simulated at the 2000th iteration of MALA. The simulations under different precisions are nearly identical but with a notable difference in execution time. Performing MALA with \pkg{MPCR}-Single precision proved to be~9X faster than the naive \proglang{R} implementation and~2X faster than \pkg{MPCR}-Double.

\subsection{Maximum Likelihood Estimation of Mat\'ern Covariance Function}

Another example demonstrating the advantages of the \pkg{MPCR} package is its application in evaluating high-dimensional likelihood functions, which are commonly needed in spatial statistics. Suppose we have a 2D spatial Gaussian random field with $14{,}400$ locations, each associated with a single measurement as shown in Figure~\ref{fig:MPCR_JSS_application_mle_simulated_values}. Our goal is to model this Gaussian field using the Mat\'{e}rn spatial covariance function:
\begin{equation}
\label{eqn:univariate_matern}
\text{cov} \{Z (\mathbf{s}_i), Z (\mathbf{s}_j)\} = \frac{\sigma^2}{2^{\nu - 1} \Gamma \left( \nu \right)} \mathcal{M}_{\nu} \left( \frac{\| \mathbf{s}_i - \mathbf{s}_j \|}{a}\right),
\end{equation}
parameterized by $\boldsymbol{\theta} = (\nu, a, \sigma^2)^\top$, where $\nu > 0, a > 0, \sigma^2>0$ are the smoothness, spatial range, and variance parameters, respectively [Guttorp and Gneiting, 2006, Mat{\'e}rn, 2013]. Here, $\mathcal{M}_{\nu}(x) = x^{\nu} \mathcal{K}_{\nu}(x)$, $\mathcal{K}_{\nu}(\cdot)$ is the modified Bessel function of the second kind of order $\nu$, and $\Gamma(\cdot)$ is the gamma function.

\begin{figure}[b!]
\centering
\includegraphics[width=7cm]{figures/MPCR_JSS_application_mle_simulated_values.pdf}
\caption{Sample spatial field on a 120 $\times$ 120 grid in a unit square generated from a Mat\'{e}rn spatial covariance function model with $\boldsymbol{\theta} = (1, 0.05, 1)^\top$.}
\label{fig:MPCR_JSS_application_mle_simulated_values}
\end{figure}

We can estimate the values of the parameters in the model in (\ref{eqn:univariate_matern}) via maximum likelihood estimation (MLE), i.e., the best estimates for the model parameters are those that maximize the Gaussian log-likelihood function:
\begin{equation} \label{eqn:loglikelihood_function}
l(\boldsymbol{\theta}) = -\frac{n}{2} \log (2 \pi) - \frac{1}{2} \log |\boldsymbol{\Sigma}(\boldsymbol{\theta})| - \frac{1}{2} \mathbf{Z}^{\top} \boldsymbol{\Sigma}(\boldsymbol{\theta})^{-1} \mathbf{Z},
\end{equation}
where $|\cdot|$ takes the determinant of $\boldsymbol{\Sigma}(\boldsymbol{\theta})$. This can be achieved by using non-linear optimizers such as those built-in in \proglang{R}, namely, \verb|optim| and \verb|nlm|. The \proglang{R} package \pkg{nloptr}, which serves as the \proglang{R} interface to the \pkg{NLopt} library, also provides a suite of different optimization routines that can also be used here. Conventionally, these non-linear optimizers operate in the opposite direction. That is, they work on finding parameters that minimize a given function. Hence, MLE is usually done by minimizing the negative Gaussian log-likelihood.

When $n$ is large,  the MLE becomes impractical due to the extensive execution time and the significant memory required to store the covariance matrix. This is because evaluating the Gaussian log-likelihood function in (\ref{eqn:loglikelihood_function}) costs $O(n^2)$ for storage and $O(n^3)$ for computation [Abdulah \emph{et~al.}, 2018] and such costs are incurred at each iteration of the algorithm. This section demonstrates how the MLE operation can be sped up using \pkg{MPCR} objects rather than native \proglang{R} objects. Similar to the example in the preceding section, we cast the covariance matrix $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ as an \proglang{R}-Double matrix and as an \pkg{MPCR}-Double matrix and compare their corresponding MLE runtimes. MLE can also be performed in lower precisions such as single or half. However, the Cholesky factorization, a key linear algebra operation involved in MLE, is known to suffer badly from round-off errors [Higham and Pranesh, 2021, Maddox \emph{et~al.}, 2022]. Hence, casting $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ as \pkg{MPCR}-Single may result in $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ being computationally non-positive definite. Nevertheless, despite the sensitivity of some linear algebra operations to lower precisions, we can still reduce the computational costs of MLE by pursuing a mixed-precision arithmetic on $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ such as those proposed in [Abdulah \emph{et~al.}, 2022] and [Cao \emph{et~al.}, 2022]. Essentially, we can store portions of $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ in varying precisions such that $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ is kept positive definite but with reduced computational overhead when performing the entire MLE. For our mixed-precision strategy, we cast $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ as an \pkg{MPCR}-Tile object and employ tile-based linear algebra operations, i.e., tile-based Cholesky factorization. Moreover, we employ the (i) band strategy [Abdulah \emph{et~al.}, 2022] and the (ii) adaptive precision-aware runtime decision strategy [Cao \emph{et~al.}, 2022] in determining the appropriate precision for each tile. In the band strategy, we specify the number of tiles closest to the diagonal that will be cast as double-precision, and the rest will be in single-precision. Meanwhile, the adaptive precision-aware runtime decision strategy will rely on the values of the global covariance matrix's Frobenius norm and each tile's Frobenius norm to determine the most appropriate precision for every single tile. Algorithm~\ref{alg:adaptive_precision_aware} outlines the steps to performing the adaptive precision-aware runtime decision strategy.

\begin{algorithm}
\caption{Adaptive Precision-Aware Runtime Decision by [Cao \emph{et~al.}, 2022].}
\label{alg:adaptive_precision_aware}
\begin{algorithmic}[1]
\State Compute the Frobenius norm of $\boldsymbol{\Sigma}(\boldsymbol{\theta})$, i.e., $\| \boldsymbol{\Sigma}(\boldsymbol{\theta})\|_F = \sqrt{\sum_{i = 1}^n \sum_{j = 1}^n | \Sigma_{ij} (\boldsymbol{\theta})|^2}$.

\State Set the value of $NT$, the number of tiles in one dimension.

\State Divide $\boldsymbol{\Sigma}(\boldsymbol{\theta})$ into tiles, i.e., $\boldsymbol{\Sigma}_{ij}(\boldsymbol{\theta}), i = 1, \ldots, M_1$ and $j = 1, \ldots, M_2$, such that $NT = M_1 M_2$.

\State Choose two precisions for the mixed-precision strategy and set the values of the variables $u_{\text{low}}$ and $u_{\text{high}}$ to their equivalent machine epsilons, e.g., $u_{\text{low}} = 2^{-24}$ and $u_{\text{high}} = 1e^{-8}$.

\For{$i = 1:{M_1}$}
\For{$j = 1:{M_2}$}
\State Compute the Frobenius norm of tile $\boldsymbol{\Sigma}_{ij}(\boldsymbol{\theta})$.
\If {$\| \boldsymbol{\Sigma}_{ij}(\boldsymbol{\theta})\|_F \leq u_{\text{high}} \| \boldsymbol{\Sigma}(\boldsymbol{\theta})\|_F / (u_{\text{low}} NT)$}
\State Store tile $\boldsymbol{\Sigma}_{ij}(\boldsymbol{\theta})$ in the lower (single) precision.
\Else
\State Store tile $\boldsymbol{\Sigma}_{ij}(\boldsymbol{\theta})$ in the higher (double) precision.
\EndIf

\EndFor
\EndFor

\end{algorithmic}
\end{algorithm}

In the following, we outline the steps to performing (A) double-precision MLE with \proglang{R}-Double and \pkg{MPCR}-Double covariance matrices and (B) mixed-precision MLE involving covariance matrices as \pkg{MPCR}-Tile matrices with double- and single-precision.

\begin{enumerate}[A.]
\item \textit{double-precision}
\begin{enumerate}[1.]
\item Simulate a sample 2D spatial field from a Mat\'{e}rn spatial covariance function model.
\begin{Schunk}
\begin{Sinput}
R> cov.matern <- function(x, nu, a, sigma_sq){
+    if(nu == 0.5) return(sigma_sq*exp( -x / a))
+    ismatrix <- is.matrix(x)
+    if(ismatrix){nr=nrow(x); nc=ncol(x)}
+    x <- c(x / a)
+    output <- rep(1, length(x))
+    n <- sum(x > 0)
+    if(n > 0) {
+      x1 <- x[x > 0]
+      output[x > 0] <-
+        (1/((2^(nu - 1)) * gamma(nu))) * (x1^nu) * besselK(x1, nu)
+    }
+    if(ismatrix){
+      output <- matrix(output, nr, nc)
+    }
+    return(sigma_sq * output)
+ }

R> M <- 120
R> n.loc <- M*M
R> locs <- cbind(rep(0:(M-1), M)/(M-1), rep(0:(M-1), each=M)/(M-1))
R> x <- as.matrix(dist(locs)) # distance matrix
R> theta <- c(1, 0.05, 1) # true parameters
R> cov.R <- cov.matern(x, theta[1], theta[2], theta[3])

R> # Simulate the spatial field.
R> library(mgcv)
R> set.seed(4)
R> z.R <- rmvn(1, rep(0,M*M), cov.R)
\end{Sinput}
\end{Schunk}

\item Create a function that computes the value of the negative Gaussian log-likelihood function. The function \verb|nll| below takes in the argument \verb|pars| which is an \proglang{R} vector object that represents $\boldsymbol{\theta}$. Moreover, the parameter values contained in the vector \verb|pars| are transformed to ensure that their values remain within their valid ranges.
\begin{Schunk}
\begin{Sinput}
R> nll <- function(pars, type, precision){
+    nu_param = 2 * 1 / (1 + exp(-pars[1]))
+    a_param = exp(pars[2])
+    sigma_param = exp(pars[3])

+    cat(c(nu_param, a_param, sigma_param), '\n')

+    #Creating the covariance matrix as R-double-precision
+    V <- cov.matern(x, nu_param, a_param, sigma_param)
+    if(type == 'FULL'){
+      if(precision == 'MPCR-Double'){
+        #Casting the covariance matrix as MPCR double-precision
+        cov_mpcr <- as.MPCR(V, nrow=M*M, ncol=M*M, precision='double')
+        L <- chol(cov_mpcr)
+        d <- log(diag(L))
+        log.det.cov <- 2*d$Sum()
+        z.mpcr <- as.MPCR(z.R, nrow=M*M, ncol=1, precision='single')
+        inner_product <- MPCR.trsm(a=L, b=z.mpcr, side='L',
+          upper_triangle=T, transpose=T, alpha=1)
+        inner.prod <- inner_product$SquareSum()
+      }else if(precision == 'R-Double'){
+        L <- t(chol(V))
+        log.det.cov <- 2*sum(log(diag(L)))
+        z.new <- forwardsolve(L, z.R)
+        inner.prod <- sum(z.new^2)
+      }
+    }else if(type == 'BANDED'){
+      cov.tile <- new(MPCRTile, nr, nc, tr, tc, V, prec.cov.banded)
+      # NOTE: the result of chol is a lower triangular matrix
+      cov.tile_chol <- chol(cov.tile, overwrite = F, num_threads = 4)
+      d <- log(cov.tile_chol$Diag())
+      log.det.cov <- 2*d$Sum()
+      z.tile <- new(MPCRTile, nr, 1, tc, 1, z.R, prec.z)
+      MPCRTile.trsm(a=cov.tile_chol, b=z.tile, side='L',
+        upper_triangle=F, transpose=F, alpha=1)
+      inner.prod <- z.tile$SquareSum()
+    }else if(type == 'ADAPTIVE'){
+      cov.tile <- new(MPCRTile, nr, nc, tr, tc, V, prec.cov)
+      V_norm <- sqrt(sum(V^2))
+      upper_bound <- u_high * V_norm / (nt * u_low)
+      for(ii in 1:tr_total){
+        for(jj in 1:tc_total){
+          test <- MPCRTile.GetTile(cov.tile, ii, jj)
+          tile_norm <- sqrt(test$SquareSum())
+          if(tile_norm < upper_bound){
+            cov.tile$ChangeTilePrecision(ii, jj, "single")
+          }
+        }
+      }
+      # NOTE: the result of chol is a lower triangular matrix
+      cov.tile_chol <- chol(cov.tile, overwrite = F, num_threads = 4)
+      d <- log(cov.tile_chol$Diag())
+      log.det.cov <- 2*d$Sum()
+      z.tile <- new(MPCRTile, nr, 1, tc, 1, z.R, prec.z)
+      MPCRTile.trsm(a=cov.tile_chol, b=z.tile, side='L',
+        upper_triangle=F, transpose=F, alpha=1)
+      inner.prod <- z.tile$SquareSum()
+    }
+    # Computing the negative Gaussian log-likelihood
+    nll <- 0.5*inner.prod + 0.5*log.det.cov + 0.5*n.loc*log(2*pi)
+    return(nll)
+ }
\end{Sinput}
\end{Schunk}

\item Perform non-linear optimization. In this example, we use the \verb|nloptr| function with the \verb|BOBYQA| subroutine [Powell \emph{et~al.}, 2009] to minimize the log-likelihood function.
\begin{Schunk}
\begin{Sinput}
R> library(nloptr)
R> opts <- list("algorithm" = "NLOPT_LN_BOBYQA", "xtol_rel" = 1e-8,
+    "maxeval" = 1000)
R> init <- c(-0.3, -1.5, -0.3)
R> fit_R <- nloptr(x0 = init, eval_f = nll, type = 'FULL',
+    precision = 'R-Double', opts = opts)
R> fit_mpcr_double <- nloptr(x0 = init, eval_f = nll,
+    type = 'FULL', precision = 'MPCR-Double', opts = opts)
\end{Sinput}
\end{Schunk}
\end{enumerate}

\item \textit{Mixed-Precision}
\begin{enumerate}[1.]
\item Define the tile sizes and the low and high precisions involved in the mixed-precision MLE.
\begin{Schunk}
\begin{Sinput}
R> u_low = 2^(-24)
R> u_high = 1e-8

R> nr <- nc <- n.loc
R> tr <- 2400
R> tc <- 2400
R> tr_total <- n.loc / tr
R> tc_total <- n.loc / tc
R> nt <- tr_total * tc_total

R> # Precisions for z
R> prec.z <- matrix(rep("single", tr_total), tr_total, 1)

R> # Precisions for cov
R> create_prec_banded_mat <- function (prec_matrix, bandwidth) {
+    d <- dim(prec_matrix)
+    outside_band <- .row(d) - bandwidth >= .col(d) |
+      .col(d) - bandwidth >= .row(d)
+    prec_matrix[outside_band] <- 'single'
+    return(prec_matrix)
+ }
R> prec.cov <- matrix(rep("double", nt), n.loc/tr, n.loc/tc)
R> prec.cov.banded <- create_prec_banded_mat(prec.cov, bandwidth=2)
\end{Sinput}
\end{Schunk}

\item Perform non-linear optimization.
\begin{Schunk}
\begin{Sinput}
R> fit_tile_adaptive <- nloptr(x0 = init, eval_f = nll,
+    type = 'ADAPTIVE', precision = NULL, opts = opts)
R> fit_tile_banded <- nloptr(x0 = init, eval_f = nll,
+    type = 'BANDED', precision = NULL, opts = opts)
\end{Sinput}
\end{Schunk}

\end{enumerate}
\end{enumerate}

\begin{table}[h!]
\caption{Summary of MLE results under different types of precision for $n = 14$,$400$. Here, the \texttt{nll} column reports the value of the negative Gaussian log-likelihood function when evaluated with the corresponding parameter estimates.}
\centering
\scalebox{0.85}{
\begin{tabular}{c c c c c c c c}
\toprule
Precision & \texttt{nll} & Parameter Estimates & Execution Time \\[0.5ex] \hline
\proglang{R}-Double & $-7,077$ & $\hat{\boldsymbol{\theta}} = (0.9862619, 0.0513225, 0.9893821)^\top$ & 48.96 hours \\
\pkg{MPCR}-Double & $-7,077$ & $\hat{\boldsymbol{\theta}} = (0.9862613, 0.05132348, 0.9894147)^\top$ & 3.62 hours \\
\pkg{MPCR}-Mixed (adaptive) & $-7,077$ & $\hat{\boldsymbol{\theta}} = (0.9862476, 0.05126803, 0.9873119)^\top$ & 3.44 hours \\
\pkg{MPCR}-Mixed (banded) & $-7,077$ & $\hat{\boldsymbol{\theta}} = (0.9870746, 0.05101309, 0.9812639 )^\top$ & 2.30 hours \\
\hline
\end{tabular}
\label{tab:MPCR_JSS_application_mle_results}
\end{table}

Table~\ref{tab:MPCR_JSS_application_mle_results} summarizes the runtimes of different approaches to performing MLE. The results show significant improvement of the mixed-precision strategy over those full precision approaches, with a reduction in the execution time alongside a minimal loss of accuracy.

\subsection{Principal Component Analysis (PCA)}

Another application example is the Principal Component Analysis (PCA) method. PCA is one of the most popular dimensionality-reduction methods in statistics. The goal of PCA is to construct a new dataset with a relatively small number of independent variables or predictors from a high-dimensional dataset containing a large number of variables or predictors. PCA proceeds by finding a linear combination of the original variables to capture most of the variation observed in the original dataset.

Performing PCA on a space-time dataset is more widely referred to as Empirical Orthogonal Function (EOF) analysis. In EOF analysis, a space-time process $Z(\mathbf{s}, t)$ can be represented as a linear combination of empirical orthogonal functions (EOFs), i.e.,
\begin{equation*}
Z(\mathbf{s}, t) = \sum_{k = 1}^K \text{PC}_k (t) \text{EOF}_k(\mathbf{s}),
\end{equation*}
where $\text{EOF}_k(\mathbf{s})$, the $k$-th EOF evaluated at spatial location $\mathbf{s}$, describes the spatial patterns of variability, while $\text{PC}_k$, the $k$-th principal component evaluated at time $t$, describes how the amplitude of the $k$-th EOF changes with time. The index $k$ orders the EOFs such that $\text{EOF}_1 (\mathbf{s})$ corresponds to the spatial pattern that explains the most variance relative to the total variance observed in the data. Likewise, $\text{EOF}_2 (\mathbf{s})$ characterizes the spatial pattern that accounts for the most variance relative to the remaining total variance, and so on.

The EOFs are obtained from the data by computing the spectral decomposition of an empirical covariance matrix or a singular value decomposition (SVD) of the centered data matrix. A matrix $\mathbf{Z} \in \mathbb{R}^{T \times N}$, where $N$ is the number of spatial locations and $T$ is the number of time points where the observations were made, can be factored using its SVD, i.e., $\mathbf{Z} = \mathbf{U} \mathbf{D} \mathbf{V}^{\top},$ where $\mathbf{U}$ is a $T \times T$ orthonormal matrix containing the left singular vectors and $\mathbf{V}$ is an $N \times N$ orthonormal matrix containing the right singular vectors. The matrix $\mathbf{D}$ is a $T \times N$ diagonal matrix containing the singular values of decreasing magnitude along the main diagonal. The cost of such SVD on a dense space-time data matrix of size $T \times N$ is $O(T^2 N + N^3)$, thereby imposing limitations when performing EOF analyses on datasets with large $T$ and $N$ [Li \emph{et~al.}, 2019].

\begin{figure}[t!]
\centering
\includegraphics[width=\textwidth]{figures/MPCR_JSS_application_pca_result_eof.pdf}
\caption{Plots of the first three EOFs of the zonal wind measurements with the corresponding execution time of the SVD using different execution scenarios.}
\label{fig:MPCR_JSS_application_pca_result_eof}
\end{figure}

In the following, we show how to use the \pkg{MPCR} package to perform EOF analysis on a high-dimensional dataset of 3 hourly averaged zonal wind measurements (in $m/s$) on a regular grid with a nominal $1^{\circ}$ horizontal resolution from November 23, 2009 to November 22, 2014, produced by the US National Center for Atmospheric Research (NCAR). We focus the analysis within the region $[7^\circ N, 85^\circ N] \times [180^\circ W, 20^\circ W]$. The dataset has $10{,}707$ spatial locations and $14{,}601$ time points.
\begin{enumerate}
\item Download the NetCDF file\\
\texttt{b.e21.BHISTcmip6.f09\_g17.LE2-1001.001.cam.h3.UBOT.2010010100-2014123100.nc} \\
from the link: \\
\textcolor{blue}{\url{https://www.earthsystemgrid.org/dataset/ucar.cgd.cesm2le.atm.proc.3hourly_ave.UBOT.html}}. \\
Alternatively, the dataset can be downloaded from this Google Drive link: \\
\textcolor{blue}{\url{https://drive.google.com/file/d/1rwBUeh8OY_wGyhkNKTltdADfsvxudGvP}}.

\item Import and process the NetCDF file.
\begin{Schunk}
\begin{Sinput}
R> library(ncdf4)
R> file_name <- "b.e21.BHISTcmip6.f09_g17.LE2-1001.001.cam.h3.UBOT.
+    2010010100-2014123100.nc"
R> nc <- nc_open(file_name)

R> lat <- ncvar_get(nc,"lat")
R> lon0 <- ncvar_get(nc, "lon")
R> lon <- (lon0+180) %% 360 - 180

R> lat_rng <- c(7,85)
R> lon_rng <- c(-180,-20)

R> lat_indx <- which(lat >= lat_rng[1] & lat <= lat_rng[2])
R> lon_indx <-which(lon >= lon_rng[1] & lon <= lon_rng[2])

R> #Extract the variable within the region of interest
R> anom_bob <- ncvar_get(nc,"UBOT",start = c(lon_indx[1],lat_indx[1],1),
+    count = c(length(lon_indx),length(lat_indx),-1))
\end{Sinput}
\end{Schunk}

\item Arrange the data in a matrix

\begin{Schunk}
\begin{Sinput}
R> s2 <- which(!is.na(anom_bob[,,1]))

R> #Vectorize the array.
R> ubot <- matrix(0,nrow = dim(anom_bob)[3], ncol = length(s2))
R> for(i in 1:dim(anom_bob)[3])
+    ubot[i,] <- anom_bob[,,i]
\end{Sinput}
\end{Schunk}

\item Compute the SVD of the data matrix for 3 precision types.
\begin{Schunk}
\begin{Sinput}
R> #Perform the native R svd() function on an R-Double data matrix.
R> svd_R <- svd(ubot)

R> #Perform the MPCR svd() function on an MPCR-Double data matrix.
R> ubot_mpcr_double <- as.MPCR(ubot, nrow = nrow(ubot), ncol = ncol(ubot),
+    precision = 'double')
R> svd_mpcr_double <- svd(ubot_mpcr_double)

R> #Perform the MPCR svd() function on an MPCR-Single data matrix.
R> ubot_mpcr_single <- as.MPCR(ubot, nrow = nrow(ubot), ncol = ncol(ubot),
+    precision = 'single')
R> svd_mpcr_single <- svd(ubot_mpcr_single)
\end{Sinput}
\end{Schunk}

\item Extract the $\mathbf{U}$, $\mathbf{V}$, and $\mathbf{D}$ vectors.
\begin{Schunk}
\begin{Sinput}
R> eof_R_u <- svd_R$u
R> eof_R_v <- svd_R$v
R> eof_R_d <- svd_R$d

R> eof_mpcr_double_u <- MPCR.ToNumericMatrix(svd_mpcr_double$u)
R> eof_mpcr_double_v <- MPCR.ToNumericMatrix(svd_mpcr_double$v)
R> eof_mpcr_double_d <- MPCR.ToNumericVector(svd_mpcr_double$d)

R> eof_mpcr_single_u <- MPCR.ToNumericMatrix(svd_mpcr_single$u)
R> eof_mpcr_single_v <- MPCR.ToNumericMatrix(svd_mpcr_single$v)
R> eof_mpcr_single_d <- MPCR.ToNumericVector(svd_mpcr_single$d)
\end{Sinput}
\end{Schunk}

\item Visualize the EOFs.
\begin{Schunk}
\begin{Sinput}
R> library(RColorBrewer)
R> library(fields)
R> library(maps)

R> loc <- as.matrix(expand.grid(x=lon[lon_indx],y=lat[lat_indx]))[s2,]
R> coltab <- colorRampPalette(brewer.pal(9,"BrBG"))(2048)

R> #Specify k or the order of the EOF you wish to plot.
R> EOF_INDEX = 1

R> #Map the k-th EOF for R-Double.
R> quilt.plot(loc,eof_R_v[, EOF_INDEX],nx=length(lon[lon_indx]),
+    ny=length(lat[lat_indx]),xlab="",ylab="",col = coltab,
+    zlim = c(-0.04, 0.04))
R> map('world',fill=F, wrap = c(-180, 180), add = TRUE)

R> #Map the k-th EOF for MPCR-Double.
R> quilt.plot(loc,eof_mpcr_double_v[, EOF_INDEX],nx=length(lon[lon_indx]),
+    ny=length(lat[lat_indx]),xlab="",ylab="",col = coltab,
+    zlim = c(-0.04, 0.04))
R> map('world',fill=F, wrap = c(-180, 180), add = TRUE)

R> #Map the k-th EOF for MPCR-Single.
R> quilt.plot(loc,eof_mpcr_single_v[, EOF_INDEX],nx=length(lon[lon_indx]),
+    ny=length(lat[lat_indx]),xlab="",ylab="",col = coltab,
+    zlim = c(-0.04, 0.04))
R> map('world',fill=F, wrap = c(-180, 180), add = TRUE)
\end{Sinput}
\end{Schunk}

\item Compute the percent of variance explained by each EOF.
\begin{Schunk}
\begin{Sinput}
R> plot(100*(eof_R_d)^2/sum(eof_R_d^2), pch = 16,
+    ylab="Percentage of variance [%]",
+    xlab="", xlim = c(1,20), col = "red")
R> lines(100*(eof_mpcr_double_d)^2/sum(eof_mpcr_double_d^2),
+    col = "green", lwd = 2, lty = 1)
R> lines(100*(eof_mpcr_single_d)^2/sum(eof_mpcr_single_d^2),
+    col = "blue", lwd = 4, lty = 3)
\end{Sinput}
\end{Schunk}

\end{enumerate}

The EOFs are the column vectors of $\mathbf{V}$, and the PCs are column vectors of $\mathbf{U}$. Figure~\ref{fig:MPCR_JSS_application_pca_result_eof} shows a map of the first three EOFs. It can be seen that the EOFs obtained when performing SVD under the three different scenarios are visually identical.

The squares of the singular values in the main diagonal of $\mathbf{D}$ specify the variance that each EOF explains.  Figure~\ref{fig:MPCR_JSS_application_pca_results_eigenvalue} plots the percentage of the variance attributable to each of the first twenty EOFs. Like the EOF results, the percent of explained variances are visually identical across different precision types.

\begin{figure}[t!]
\centering
\includegraphics[width=0.8\textwidth]
\caption{Plots of the percentage of variance explained by the first 20 ordered eigenvalues using different execution scenarios. Under all three scenarios, the first three EOFs explain $\sim$36\%, $\sim$6\%, and $\sim$4\% of the total variance, respectively.}
\label{fig:MPCR_JSS_application_pca_results_eigenvalue}
\end{figure}

\subsection{Integrated Nested Single-Half-Arithmetic Laplace Approximation (INSHALA)}

Another example to show the capabilities of the \pkg{MPCR} package is large-scale Bayesian inference problems, which operate on large precision matrices.

The integrated nested Laplace approximations (INLA) method developed by [Rue \emph{et~al.}, 2009] is a widely known and applied framework for approximate Bayesian inference. INLA is formulated for latent Gaussian models, including but not limited to linear mixed models, spatial, and spatio-temporal models [Rue \emph{et~al.}, 2017, Bakka \emph{et~al.}, 2018]. INLA has a well-maintained dedicated \proglang{R} package called \pkg{R-INLA} (\url{http://www.r-inla.org}).

A detailed discussion of each step in the INLA method, complete with \proglang{R} code examples, is accessible at the website \url{https://stefansiegert.net/inla-project/inla-from-scratch}. In the following demonstration, we closely align with their presentation, highlighting how the \pkg{MPCR} functions can be seamlessly incorporated into their coding framework.

Consider the following hierarchical model:
\begin{itemize}
\item The observations $y_1, \ldots, y_n$ are from a Bernoulli process with parameters $p_1, \ldots p_n$, i.e.,
\begin{equation*}
p(y_t | p_t) = \begin{cases}
p_t, & y_t = 1,\\
1 - p_t, & y_t = 0,
\end{cases}
\end{equation*}
where $t = 1, 2, \ldots, n$.
\item The parameters $p_1, \ldots, p_n$ are assumed to change over time smoothly, and they depend on a latent process $x_1, \ldots, x_n$ through a logit link:
\begin{equation*}
p_t = \frac{e^{\beta x_t}}{1 + e^{\beta x_t}}.
\end{equation*}

\item The process $x_1, \ldots, x_n$ is a Gaussian process, i.e., $(x_1, \ldots, x_n)^{\top} = \{X(t_1), \ldots, X(t_n)\}^{\top} \sim \mathcal{N}_{n} \{\mathbf{0}, \boldsymbol{\Sigma}(\alpha)\}$, where $\Sigma_{ij}(\alpha) = \text{cov} \{X (t_i), X (t_j)\} = \sigma^2 \exp ( - | t_i - t_j | / \alpha)$.

\item The hyperparameter $\alpha$ has a uniform prior $U(0, 1)$.
\end{itemize}

The joint log-density of the observations, latents, and hyperparameters is given by
\begin{eqnarray}
\label{eqn:joint_log_density}
\log\{ p(\mathbf{y}, \mathbf{x}, \alpha)\} & = & \log\{ p(\mathbf{y} | \mathbf{x}, \alpha)\} + \log\{ p(\mathbf{x} | \alpha)\} + \log\{p(\alpha)\} \nonumber \\
& \propto & \sum_t^n \{\beta x_t y_t - \log (1 + e^{\beta x_t})\} + \frac{1}{2} \log |\mathbf{Q}| - \frac{1}{2} \mathbf{x}^{\top} \mathbf{Q} \mathbf{x} + \log\{p(\alpha)\},
\end{eqnarray}
where $\mathbf{Q} = \boldsymbol{\Sigma}^{-1}(\alpha)$ is the inverse covariance matrix or more commonly termed precision matrix. The term ``precision'' in the context of the precision matrix $\mathbf{Q}$ should not be confused with the term ``precision'' used to describe numerical formats, which is the subject of this paper.

Fitting a model using the INLA method has many operations based on linear algebra where the INLA method and the \pkg{R-INLA} package proved to be computationally efficient. However, this efficiency results from working only with sparse matrices [Abdul-Fattah \emph{et~al.}, 2023]. In this regard, \pkg{MPCR} can be used to overcome the obstacles INLA encounters with dense matrices. Herein, we introduce the Integrated Nested Single-Half-Arithmetic Laplace Approximation (INSHALA) powered by the \pkg{MPCR} package. INSHALA is a more computationally efficient version of INLA that operates in low precisions such as 32-bit (single) precision formats. In the model in (\ref{eqn:joint_log_density}) above, for example, the INSHALA method would store the precision matrix $\mathbf{Q}$ as an \pkg{MPCR}-Single-precision matrix to improve computational performance.

In the following, we show how the INLA-inspired INSHALA methodology works. In particular, we illustrate how to approximate the posterior distribution of $\alpha$ with \pkg{MPCR}-Double and \pkg{MPCR}-Single-precision computations. Furthermore, the classical INLA modeling implemented using native \proglang{R} objects is also shown for comparison.

Note that the log posterior distribution of the hyperparameter $\alpha$ has the following expression:
\begin{eqnarray}
\label{eqn:posterior_distribution}
\log \{ p(\alpha | \mathbf{y}) \} \propto \log\{ p(\mathbf{y}, \mathbf{x}, \alpha)\} - \log \{ p(\mathbf{x} | \alpha, \mathbf{y}) \},
\end{eqnarray}
which is derived from another factorization of the joint log-density in (\ref{eqn:joint_log_density}). Here the value of $\log\{ p(\mathbf{y}, \mathbf{x}, \alpha)\}$ can be computed by a straightforward evaluation of the expression in (\ref{eqn:joint_log_density}). The value of the conditional log distribution of the latent variables $\log \{ p(\mathbf{x} | \alpha, \mathbf{y}) \}$, on the other hand, is obtained using Laplace approximation. A key step in the INLA method is to define a function $f(\mathbf{x})$ such that by applying the definition of conditional probability, we can express $p(\mathbf{x} | \alpha, \mathbf{y})$ as:
\begin{eqnarray}
\label{eqn:latent_posterior_distribution}
p(\mathbf{x} | \alpha, \mathbf{y}) = \frac{e^{f(\mathbf{x})}}{\int e^{f(\mathbf{x})} d \mathbf{x}},
\end{eqnarray}
and then apply a Laplace approximation to the denominator $\int e^{f(\mathbf{x})} d \mathbf{x}$, i.e.,
\begin{eqnarray}
\label{eqn:laplace_approximation}
\int e^{f(\mathbf{x})} d \mathbf{x} \approx e^{f(\mathbf{x}_0)} (2 \pi)^{n / 2} | - \mathbf{H}f(\mathbf{x}_0) |^{-1/2}.
\end{eqnarray}
Here, $\mathbf{x}_0$ is the mode of $f(\mathbf{x})$ and $\mathbf{H}f(\mathbf{x}_0)$ is the Hessian matrix of $f(\mathbf{x})$ evaluated at the mode. Note that the form in (\ref{eqn:laplace_approximation}) appears as a result of a second-order Taylor expansion of $f(\mathbf{x})$ around the mode (see \cite {opitz2017latent, flagg2023integrated} for some derivations).

Applying (\ref{eqn:laplace_approximation}) to (\ref{eqn:latent_posterior_distribution}), we get the Laplace approximation of the conditional log distribution $\log \{p(\mathbf{x} | \alpha, \mathbf{y})\}$ evaluated at mode $\mathbf{x}_0$ as follows:
\begin{eqnarray}
\label{eqn:laplace_approximation_conditional_distribution}
\log \{p(\mathbf{x}_0 | \alpha, \mathbf{y})\} \approx -\frac{n}{2} \log(2 \pi) + \frac{1}{2} \log | - \mathbf{H}f(\mathbf{x}_0) |.
\end{eqnarray}
By subtituting (\ref{eqn:laplace_approximation_conditional_distribution}) to (\ref{eqn:posterior_distribution}), the approximation of the log posterior distribution in (\ref{eqn:posterior_distribution}) is:
\begin{eqnarray}
\label{eqn:posterior_distribution_approximation}
\log \{ p(\alpha | \mathbf{y}) \} \propto \log\{ p(\mathbf{y}, \mathbf{x}_0, \alpha)\} - \frac{1}{2} \log | - \mathbf{H}f(\mathbf{x}_0) |.
\end{eqnarray}

In order to compute the posterior distribution of $\alpha$ using the formula in (\ref{eqn:posterior_distribution_approximation}), we need to know the form of $f(\mathbf{x})$ and consequently its Hessian matrix $\mathbf{H}f(\mathbf{x})$. In the example above, a great choice for $f(\mathbf{x})$ in (\ref{eqn:latent_posterior_distribution}) is such that $f(\mathbf{x}) \propto \log p(\mathbf{y}, \mathbf{x}, \alpha)$. From the formulation of $\log p(\mathbf{y}, \mathbf{x}, \alpha)$ in (\ref{eqn:joint_log_density}), we write $f(\mathbf{x})$ as a function of the components that depend on $\mathbf{x}$ as follows:
\begin{eqnarray}
\label{eqn:f}
f(\mathbf{x}) = \sum_t^n \{\beta x_t y_t - \log (1 + e^{\beta x_t})\} - \frac{1}{2} \mathbf{x}^{\top} \mathbf{Q} \mathbf{x}.
\end{eqnarray}
By the above form of $f(\mathbf{x})$, the gradient vector and Hessian matrix of $f(\mathbf{x})$ are the following:
\begin{eqnarray}
\label{eqn:gradient_vector}
\nabla f(\mathbf{x}) & = & \mbox{vec} \left\{\beta y_t - \frac{\beta e^{\beta x_t}}{1 + e^{\beta x_t}} \right\} - \mathbf{Q} \mathbf{x}, \text{ and} \\
\label{eqn:hessian_matrix}
\mathbf{H} f(\mathbf{x}) & = & -\mathbf{Q} - \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_t}}{(1 + e^{\beta x_t})^2} \right\},
\end{eqnarray}
where $\nabla f(\mathbf{x})$ is a column vector of size $t \times 1$ and $\mathbf{H} f(\mathbf{x})$ is a $t \times t$ matrix. Here, $\mbox{vec()}$ constructs a column vector while $\mbox{diag()}$ builds a diagonal matrix.

Finally, the mode $\mathbf{x}_0$ of $f(\mathbf{x})$ can be found iteratively using the steps outlined in Algorithm~\ref{alg:solving_for_mode}.
\begin{algorithm} [h!]
\caption{Solving for the mode $\mathbf{x}_0$ iteratively. \\
Source: \url{https://stefansiegert.net/inla-project/inla-from-scratch}}
\label{alg:solving_for_mode}
\begin{algorithmic}[1]
\State Set an initial value for $\mathbf{x}_0$.

\State Taylor expand $f(\mathbf{x})$ around $\mathbf{x}_0$ and keep only the terms that depend on $\mathbf{x}$, i.e.,
\begin{eqnarray*}
f(\mathbf{x}) & \approx & f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)^{\top} (\mathbf{x} - \mathbf{x}_0) + \frac{1}{2} (\mathbf{x} - \mathbf{x}_0)^{\top} \mathbf{H}f(\mathbf{x}_0) (\mathbf{x} - \mathbf{x}_0) \\
& \propto & \{ \nabla f(\mathbf{x}_0)^{\top} - \mathbf{x}_0^{\top} \mathbf{H}f(\mathbf{x}_0) \} \mathbf{x} + \frac{1}{2} \mathbf{x}^{\top} \mathbf{H}f(\mathbf{x}_0) \mathbf{x}.
\end{eqnarray*}

\State Find the mode of the second-order Taylor expansion above. This mode serves as a new and improved estimate of the true mode of $f(\mathbf{x})$. This can be done by setting its first derivative to zero and solving for $\mathbf{x}$, i.e.,
\begin{eqnarray*}
\nabla f(\mathbf{x}) & \propto & \nabla f(\mathbf{x}_0) - \mathbf{H}f(\mathbf{x}_0) \mathbf{x}_0 + \mathbf{H}f(\mathbf{x}_0) \mathbf{x} \\
& = & \mbox{vec} \left\{\beta y_t - \frac{\beta e^{\beta x_{0,t}}}{1 + e^{\beta x_{0,t}}} \right\} - \mathbf{Q} \mathbf{x}_0 - \Bigg[-\mathbf{Q} - \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \Bigg] \mathbf{x}_0 \\
&& + \Bigg[-\mathbf{Q} - \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \Bigg] \mathbf{x}  \\
& = & \mbox{vec} \left\{\beta y_t - \frac{\beta e^{\beta x_{0,t}}}{1 + e^{\beta x_{0,t}}} \right\} + \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \mathbf{x}_0 \\
&& - \left[\mathbf{Q} + \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \right] \mathbf{x} = \mathbf{0},
\end{eqnarray*}
where the mode of the second-order Taylor expansion of $f(\mathbf{x})$ is the solution to the equation
\begin{eqnarray*}
\mathbf{x} = \left[\mathbf{Q} + \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \right]^{-1} \left[\mbox{vec} \left\{\beta y_t - \frac{\beta e^{\beta x_{0,t}}}{1 + e^{\beta x_{0,t}}} \right\} + \mbox{diag} \left\{\frac{\beta^2 e^{\beta x_{0,t}}}{(1 + e^{\beta x_{0,t}})^2} \right\} \mathbf{x}_0\right].
\end{eqnarray*}

\State Update the value of $\mathbf{x}_0$ to the value of $\mathbf{x}$ above, i.e., $\mathbf{x}_0 \leftarrow  \mathbf{x}$.

\State Repeat until the distance between $\mathbf{x}_0$ and $\mathbf{x}$ is less than a user-defined tolerance, \verb|tol|, i.e., $\| \mathbf{x}_0 - \mathbf{x} \|^2 < \verb|tol|$.
\end{algorithmic}
\end{algorithm}

In the following, we outline how to perform the approximation of $\log \{ p(\alpha | \mathbf{y}) \}$ in (\ref{eqn:posterior_distribution_approximation}) via the INSHALA method using the \pkg{MPCR} package.

\begin{enumerate}
\item Simulate the data.
\begin{enumerate}
\item Simulate $x_1, \ldots, x_n$ from an exponential covariance function model.
\begin{Schunk}
\begin{Sinput}
R> library(mgcv)
R> library(Matrix)
R> library(dplyr)

R> n <- 50^2
R> alpha_true = 0.6
R> sigma_true = 0.1

R> x <- seq(0, n, length.out = n)
R> locs <- cbind(x, 0)
R> dist0 <- as.matrix(dist(locs)) # distance matrix in time
R> cov.R <- sigma_true * exp(-dist0 / alpha_true)

R> nr <- nrow(cov.R)
R> nc <- ncol(cov.R)

R> set.seed(4)
R> # Simulate x_1, x_2, ..., x_n
R> x_true <- rmvn(1, rep(0, n), cov.R)
\end{Sinput}
\end{Schunk}
\item Simulate $p_1, \ldots, p_n$ using the logit link.
\begin{Schunk}
\begin{Sinput}
R> beta_true = 10
R> p_true <- 1 / (1 + exp(-beta_true * x_true))
\end{Sinput}
\end{Schunk}

\item Simulate $y_t$ as Bernoulli$(p_t)$.
\begin{Schunk}
\begin{Sinput}
R> y <- rbinom(n, 1, p_true)
\end{Sinput}
\end{Schunk}
\end{enumerate}

\item Create a function that computes $\mathbf{Q}$.
\begin{Schunk}
\begin{Sinput}
R> calc_Q = function(alpha, precision) {
+    cov_R <- cov <- sigma_true * exp( - dist0 / alpha)
+    if(precision == 'MPCR-Double'){
+      cov <- as.MPCR(as.matrix(cov_R), nrow=nrow(cov_R), ncol=ncol(cov_R),
+        precision='double')
+    }else if(precision == 'MPCR-Single'){
+      cov <- as.MPCR(as.matrix(cov_R), nrow=nrow(cov_R), ncol=ncol(cov_R),
+        precision='single')
+    }
+    Q <- solve(cov)
+    return(Q)
+ }
\end{Sinput}
\end{Schunk}

\item Create a function that computes the log prior for $\alpha$.
\begin{Schunk}
\begin{Sinput}
R> calc_lprior = function(alpha, a=1, b=1) {
+    (a-1) * log(alpha) + (b-1) * log(1 - alpha)
+ }
\end{Sinput}
\end{Schunk}

\item Create a function that computes the log joint distribution in (\ref{eqn:joint_log_density}).
\begin{Schunk}
\begin{Sinput}
R> calc_ljoint = function(y, x, alpha, a=1, b=1, precision) {
+    Q <- calc_Q(alpha, precision)
+    chol_Q <- chol(Q)
+    if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+      d <- log(diag(chol_Q))
+      logdet_Q_half <- d$Sum()
+      x.mpcr <- as.MPCR(x, nrow=length(x), ncol=1, precision='single')
+      inner_product <- chol_Q %*% x.mpcr
+      quad_form <- inner_product$SquareSum()
+    }else{
+      logdet_Q_half <- chol_Q %>% diag %>% log %>% sum
+      quad_form <- crossprod(chol_Q %*% x) %>% drop
+    }
+    res <- sum(beta_true * x * y - log1p(exp(beta_true * x))) +
+      logdet_Q_half - 0.5 * quad_form + calc_lprior(alpha, a, b)
+    return(res)
+ }
\end{Sinput}
\end{Schunk}

\item Create a function that computes $f(\mathbf{x})$ in (\ref{eqn:f}), its gradient (\ref{eqn:gradient_vector}) and its negative Hessian~(\ref{eqn:hessian_matrix}).
\begin{Schunk}
\begin{Sinput}
R> calc_ff = function(x, alpha, precision) {
+    sum(beta_true * x * y - log1p(exp(beta_true * x))) -
+      0.5 * drop(as.matrix(x %*% calc_Q(alpha, precision) %*% x))
+ }

R> calc_grad_ff = function(x, alpha, precision) {
+    beta_true * y -
+      beta_true * exp(beta_true * x) / (1 + exp(beta_true * x)) -
+      drop(as.matrix(calc_Q(alpha, precision) %*% x))
+ }

R> calc_neg_hess_ff = function(x, alpha, precision) {
+    Q <- calc_Q(alpha, precision)
+    other_term <- diag(beta_true^2 * exp(beta_true * x) /
+      (1 + exp(beta_true * x))^2)
+    if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+      other_term.mpcr <- as.MPCR(other_term, nrow=nrow(other_term),
+        ncol=ncol(other_term), precision='single')
+      return(Q + other_term.mpcr)
+    }else{
+      return(Q + other_term)
+    }
+ }
\end{Sinput}
\end{Schunk}

\item Create a function that computes the mode of $f(\mathbf{x})$.
\begin{Schunk}
\begin{Sinput}
R> # the function g(x) = log(p(y | x,theta)), its gradient and hessian
R> calc_g0 = function(x) {
+    sum(beta_true * x * y - log1p(exp(beta_true*x)))

R> calc_g1 = function(x) {
+    beta_true * y - beta_true * exp(beta_true*x) / (1 + exp(beta_true*x))
+ }

R> calc_g2 = function(x) {
+    (-1) * beta_true^2 * exp(beta_true*x) / (1 + exp(beta_true*x))^2
+ }

R> calc_x0 = function(alpha, precision, tol=1e-12) {
+    Q <- calc_Q(alpha, precision)
+    x <- x0 <- rep(0, n)
+    while(1) {
+      g1 <- calc_g1(x)
+      g2 <- calc_g2(x)
+      diag_g <- bandSparse(n=n, k=0, diagonals=list(g2))
+      mode <- (g1 - x0 * g2)
+      if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+        diag_g.mpcr <- as.MPCR(as.matrix(diag_g), nrow=nrow(diag_g),
+          ncol=ncol(diag_g), precision='double')
+        mode.mpcr <- as.MPCR(mode, nrow=length(x0), ncol=1,
+          precision='single')
+        x <- MPCR.ToNumericVector(solve(Q - diag_g.mpcr) %*% mode.mpcr)
+      }else{
+        x <- drop(solve(Q - diag_g) %*% mode)
+      }
+      if(mean((x-x0)^2 < tol)) {
+        break
+      } else {
+        x0 <- x
+      }
+    }
+    return(x)
+ }
\end{Sinput}
\end{Schunk}

\item Create a function that approximates the log posterior in Equation~(\ref{eqn:posterior_distribution_approximation}) up to an additive constant.
\begin{Schunk}
\begin{Sinput}
R> calc_lpost = function(alpha, precision) {
+    x0 <- calc_x0(alpha, precision)
+    chol_h <- chol(calc_neg_hess_ff(x0, alpha, precision))
+    if(precision %in% c('MPCR-Double', 'MPCR-Single')){
+      d <- log(diag(chol_h))
+      logdet_h_half <- d$Sum()
+    }else{
+      logdet_h_half <- sum(log(diag(chol_h)))
+    }
+    calc_ljoint(y, x0, alpha, precision = precision) - logdet_h_half
+ }
\end{Sinput}
\end{Schunk}

\item For each precision type, compute the approximation of the log posterior via the iterative method.
\begin{Schunk}
\begin{Sinput}
R> alpha_vec = seq(0.05, 0.95, len=21)
R> lpost_R = sapply(alpha_vec, calc_lpost, precision = 'R-Double')
R> lpost_mpcr_d = sapply(alpha_vec, calc_lpost, precision='MPCR-Double')
R> lpost_mpcr_s = sapply(alpha_vec, calc_lpost, precision='MPCR-Single')
\end{Sinput}
\end{Schunk}

\item Create a function that normalizes the posterior via numerical integration.
\begin{Schunk}
\begin{Sinput}
R> calc_Z = function(alpha_vec, lpost_vec) {
+    nn <- length(alpha_vec)
+    hh <- alpha_vec[2] - alpha_vec[1]
+    ww <- c(1, rep(c(4,2), (nn-3)/2), c(4,1))
+    return(sum(ww * exp(lpost_vec)) * hh / 3)
+ }
\end{Sinput}
\end{Schunk}

\item Normalize the posterior.
\begin{Schunk}
\begin{Sinput}
R> lpost_R = lpost_R - mean(lpost_R)
R> Z_R = calc_Z(alpha_vec, lpost_R)

R> lpost_mpcr_d = lpost_mpcr_d - mean(lpost_mpcr_d)
R> Z_mpcr_double = calc_Z(alpha_vec, lpost_mpcr_d)

R> lpost_mpcr_s = lpost_mpcr_s - mean(lpost_mpcr_s)
R> Z_mpcr_single = calc_Z(alpha_vec, lpost_mpcr_s)
\end{Sinput}
\end{Schunk}
\end{enumerate}

\begin{figure*}[t!]
\centering
\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth]{figures/MPCR_JSS_application_inla_result_posterior.pdf}
\caption{Posterior distribution of $\alpha$}
\label{fig:MPCR_JSS_application_inla_result_posterior}
\end{subfigure}%
\begin{subfigure}[t]{0.5\textwidth}
\centering
\includegraphics[width=\textwidth, height=7.5cm]{figures/MPCR_JSS_application_inla_result_time.pdf}
\caption{Execution Time}
\label{fig:MPCR_JSS_application_inla_result_time}
\end{subfigure}
\caption{Results of the INSHALA experiments: (a) Approximation of the posterior distribution of $\alpha$ for $n = 2500$. The bold vertical line marks the true value of $\alpha = 0.6$. (b)~Execution time of computing the posterior distribution for various $n$ under different types of precision.}
\label{fig:INSHALA_results}
\end{figure*}

We run the aforementioned for different sizes of the precision matrix $\mathbf{Q}$. Figure~\ref{fig:MPCR_JSS_application_inla_result_posterior} plots the resulting posterior distribution of $\alpha$ when $n = 2,500$ for each precision type. From the figure, it can be seen that the INSHALA method yields identical results to that of the INLA method which was obtained using \proglang{R}-double-precision. Figure~\ref{fig:MPCR_JSS_application_inla_result_time} plots the execution time in log scale when running the INLA and INSHALA procedures. From the figure, it can be seen that for reasonably large values of $n$, our proposed INSHALA method yields a non-trivial decrease in execution time while maintaining results identical to those of the INLA method.

%\section{Comparison with existing package} \label{sec:comparison}

\section{Discussion and Future Work} \label{sec:discuss}

Employing multi- and mixed-precision computation in scientific applications can significantly speedup the processing of extensive, high-dimensional problems. In this work, we introduce \pkg{MPCR}, a package designed to facilitate multi- and mixed-precision computation within the \proglang{R} environment. This package is integrated with well-tuned BLAS/LAPACK libraries, e.g., Intel MKL and OpenBLAS, ensuring optimal performance on available hardware. Herein, we demonstrate the usage of the package in \proglang{R} through various examples and show how it can enhance the efficiency of existing applications from the execution time perspective by reducing the computation precision while preserving the same accuracy of double precision. The package also provides mixed-precision support of matrix operations that enable faster linear algebra algorithms, e.g., tile-based linear algebra using parallel execution and a mix of precisions, i.e., 32-bit and 64-bit, through leverages OpenMP for parallelization on multicore systems. \pkg{MPCR} also offers a robust framework that can support additional precisions, like 128-bit and 16-bit, provided the underlying BLAS/LAPACK library and hardware architecture are compatible. The seamless integration with the optimized BLAS/LAPACK libraries positions \pkg{MPCR} as a highly effective alternative to traditional \proglang{R} operations, mainly when \proglang{R} is linked to a less optimized library such as RBLAS.

We plan to extend the \pkg{MPCR} package in our future work by integrating GPU support as a high-performance backend. This upgrade will empower \proglang{R} users to execute their linear algebra operations more efficiently. Additionally, we aim to incorporate support for 16-bit precision through tensor core technology available in modern GPU hardware, such as V100, A100, and H100. These advanced technologies can offer a significant speed boost compared to higher precision operations.

%\section{Conclusion and future work} \label{sec:conclusion}

%As computational workloads in scientific applications increase, there is a demand for methods in computational statistics to provide high compute power in low precision formats. The current state-of-the-art in numerical computing is gearing towards lower numerical precision.
%Multi- and mixed-precision computing deliver the promise of higher computational performance with no significant impact on accuracy.

%We introduced the \proglang{R} package \pkg{MPCR} which enables users to carry out computations on multi-precision vectors, matrices, and tile matrices layouts

%We have shown the effectiveness
%of the package in various tasks in statistics. \textcolor{blue}{providing more than an order of magnitude higher performance than what is available in IEEE double-precision}. The package leverages both the user-friendly interface of \proglang{R} and state-of-the-art linear algebra libraries in high-performance computing. aiming to maintain the operation’s accuracy while simultaneously enhancing its speed

%\textcolor{purple}{Aside from new developments in hardware architectures, carefully designed numerical algorithms that remain efficient and stable given the varying precisions of data and computations are important forces propelling multi- and mixed-precision computing. }

%% -- Optional special unnumbered sections -------------------------------------

\section*{Computational details}

The results in this paper were obtained using
\proglang{R}~4.3.1 with the
\pkg{MPCR}~1.0.0 package. \proglang{R} itself
and all packages used are available from the Comprehensive
\proglang{R} Archive Network (CRAN) at
\url{https://CRAN.R-project.org/}.

The experiments were conducted on an aarch64-apple-darwin20 (64-bit) platform running under macOS Sonoma 14.0.

\section*{Acknowledgment}
%This research received support from the King Abdullah University of Science and Technology (KAUST) in Saudi Arabia. Our gratitude extends to the team at the KAUST Supercomputing Laboratory (KSL) and the Extreme Computing Research Center (ECRC) for providing the computational resources that were essential for the experiments conducted in this study.

\begin{leftbar}
This study was funded by the King Abdullah University of Science and Technology (KAUST) in Saudi Arabia. We would like to thank the team at BrightSkies for their technical assistance in the development of the \pkg{MPCR} package and to Omar Marzouk and Merna Moawad for their support.
\end{leftbar}

%% -- Bibliography -------------------------------------------------------------
%% - References need to be provided in a .bib BibTeX database.
%% - All references should be made with \cite, \citet, \citep, \citealp etc.
%%   (and never hard-coded). See the FAQ for details.
%% - JSS-specific markup (\proglang, \pkg, \code) should be used in the .bib.
%% - Titles in the .bib should be in title case.
%% - DOIs should be included where available.

\begin{thebibliography}{84}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Abdelfattah \emph{et~al.}, 2020]
Abdelfattah A, Anzt H, Boman EG, Carson E, Cojean T, Dongarra J, Gates M,
Gr{\"u}tzmacher T, Higham NJ, Li S, \emph{et~al.} (2020).
\newblock \enquote{A survey of numerical methods utilizing mixed precision
\newblock \emph{arXiv preprint arXiv:2007.06674}.

[Abdul-Fattah \emph{et~al.}, 2023]
Abdul-Fattah E, Van~Niekerk J, Rue H (2023).
\newblock \enquote{INLA+--Approximate Bayesian inference for non-sparse models
using HPC.}
\newblock \emph{arXiv preprint arXiv:2311.08050}.

[Abdulah \emph{et~al.}, 2022]
Abdulah S, Cao Q, Pei Y, Bosilca G, Dongarra J, Genton MG, Keyes DE, Ltaief H,
Sun Y (2022).
\newblock \enquote{Accelerating geostatistical modeling and prediction with
mixed-precision computations: A high-productivity approach with parsec.}
\newblock \emph{IEEE Transactions on Parallel and Distributed Systems},
\textbf{33}(4), 964--976.

[Abdulah \emph{et~al.}, 2018]
Abdulah S, Ltaief H, Sun Y, Genton MG, Keyes DE (2018).
\newblock \enquote{Exageostat: A high performance unified software for
geostatistics on manycore systems.}
\newblock \emph{IEEE Transactions on Parallel and Distributed Systems},
\textbf{29}(12), 2771--2784.

[Abdulah \emph{et~al.}, 2019]
Abdulah S, Ltaief H, Sun Y, Genton MG, Keyes DE (2019).
\newblock \enquote{Geostatistical modeling and prediction using mixed precision
tile Cholesky factorization.}
\newblock In \emph{2019 IEEE 26th international conference on high performance
computing, data, and analytics (HiPC)}, pp. 152--162. IEEE.

[Agullo \emph{et~al.}, 2009]
Agullo E, Demmel J, Dongarra J, Hadri B, Kurzak J, Langou J, Ltaief H, Luszczek
P, Tomov S (2009).
\newblock \enquote{Numerical linear algebra on emerging architectures: The
PLASMA and MAGMA projects.}
\newblock In \emph{Journal of Physics: Conference Series}, volume 180, p.
012037. IOP Publishing.

[Akbudak \emph{et~al.}, 2017]
Akbudak K, Ltaief H, Mikhalev A, Keyes D (2017).
\newblock \enquote{Tile low rank Cholesky factorization for climate/weather
modeling applications on manycore architectures.}
\newblock In \emph{International Conference on High Performance Computing}, pp.
22--40. Springer.

[Anderson \emph{et~al.}, 1999]
Anderson E, Bai Z, Bischof C, Blackford S, Dongarra JDJ, Croz JD, Greenbaum A,
Hammarling S, McKenney A, Sorensen D (1999).
\newblock \emph{LAPACK Users' Guide}.
\newblock Third edition. SIAM, Philadelphia, Pennsylvania, USA.

[Au and Beck, 2001]
Au SK, Beck JL (2001).
\newblock \enquote{Estimation of small failure probabilities in high dimensions
by subset simulation.}
\newblock \emph{Probabilistic Engineering Mechanics}, \textbf{16}(4), 263--277.

[Bakka \emph{et~al.}, 2018]
Bakka H, Rue H, Fuglstad GA, Riebler A, Bolin D, Illian J, Krainski E, Simpson
D, Lindgren F (2018).
\newblock \enquote{Spatial modeling with R-INLA: A review.}
\newblock \emph{Wiley Interdisciplinary Reviews: Computational Statistics},
\textbf{10}(6), e1443.

[Beaumont \emph{et~al.}, 2020]
Beaumont O, Langou J, Quach W, Shilova A (2020).
\newblock \enquote{A makespan lower bound for the tiled cholesky factorization
based on alap schedule.}
\newblock In \emph{Euro-Par 2020: Parallel Processing: 26th International
Conference on Parallel and Distributed Computing, Warsaw, Poland, August
24--28, 2020, Proceedings}, pp. 134--150. Springer.

[Benner \emph{et~al.}, 2022]
Benner P, Iannazzo B, Meini B, Palitta D (2022).
\newblock \enquote{Palindromic linearization and numerical solution of
nonsymmetric algebraic T-Riccati equations.}
\newblock \emph{BIT Numerical Mathematics}, \textbf{62}(4), 1649--1672.

[Bj{\"o}rck \emph{et~al.}, 2021]
Bj{\"o}rck J, Chen X, De~Sa C, Gomes CP, Weinberger K (2021).
\newblock \enquote{Low-precision reinforcement learning: running soft
actor-critic in half precision.}
\newblock In \emph{International Conference on Machine Learning}, pp. 980--991.

[Blackford \emph{et~al.}, 1996]
Blackford LS, Choi J, Cleary A, Petitet A, Whaley RC, Demmel J, Dhillon I,
Stanley K, Dongarra J, Hammarling S, Henry G, Walker D (1996).
\newblock \enquote{{ScaLAPACK}: {A} Portable Linear Algebra Library for
Distributed Memory Computers --- Design Issues and Performance.}
\newblock pp. 95--106.
\newblock