# Mixed-Precision HODLR Matrices: Algorithms and Software Overview

Erin Carson, Xinye Chen, Xiaobo Liu

<figure><img src="image_1.png" /></figure>

## To cite this version:

Erin Carson, Xinye Chen, Xiaobo Liu. Mixed-Precision HODLR Matrices: Algorithms and Software Overview. Doctoral. Heidelberg, France. 2025. hal-05136966

HAL Id: hal-05136966

https://hal.science/hal-05136966v1

Submitted on 8 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mixed-Precision HODLR Matrices: Algorithms and
Software Overview

Erin Carson, *Xinye Chen*, and Xiaobo Liu

Sorbonne Université, CNRS, LIP6, F-75005 Paris, France

<figure><img src="image_3.png" /></figure>

SCIENCES
SORBONNE
UNIVERSITÉ

Xinye.Chen@lip6.fr

June 30, 2025
## Outline

1. HODLR matrix

* What is HODLR matrix?

* $(p, ε)$-HODLR Matrix

2. Low precision computing routines

* Matrix-vector products

* LU factorization

* Storage reduction

3. Related materials

4. Datasets

5. Software demonstration
# What is HODLR matrix?

## Definition

An HODLR (Hierarchical Off-Diagonal Low-Rank) matrix is a structured matrix where off-diagonal blocks exhibit low numerical rank, enabling efficient storage and computation.

* **Key Property:** Hierarchical partitioning reveals rank-deficient off-diagonal sub-blocks (rank-$p$ matrix).

* **Applications:** Kernel methods, integral equations, covariance matrices.

* **Advantages:**

* Fast matrix-vector products: $O(n \log n)$ or $O(n)$.

* Efficient factorizations (e.g., LU).

* Reduced complexity vs. $O(n^3)$ for dense matrices.

(Here, $n$ denotes the matrix size.)
# What is HODLR matrix?

Convert matrix A into HODLR matrix H:

1. Partition $H^{(k)}$ into
   $$ \begin{bmatrix} H_{11}^{(k+1)} & H_{12}^{(k+1)} \\ H_{21}^{(k+1)} & H_{22}^{(k+1)} \end{bmatrix} $$

2. Low rank approximation of $H_{12}^{(k+1)}$ ($\approx U_1^{(k+1)} V(k+1)_2^T$) and $H_{21}^{(k+1)}$ ($\approx U_2^{(k+1)} V(k+1)_1^T$) until bottom level achieves or minimum block size is achieved.

<figure><img src="image_4.png" /></figure>

## Lower storage and faster arithmetic?

## Reduced precision for off-diagonal blocks.
Precision format

Table: *u*: unit roundoff. *x*<sub>min</sub>: smallest positive normalized floating-point number. *x*<sub>max</sub>: largest floating-point number. *t*: number of binary digits in the significand (including the implicit leading bit), *e*<sub>min</sub>: exponent of *x*<sub>min</sub>, and *e*<sub>max</sub>: exponent of *x*<sub>max</sub>.

<table><thead><tr><th></th><th>u</th><th>x<sub>min</sub></th><th>x<sub>max</sub></th><th>t</th><th>e<sub>min</sub></th><th>e<sub>max</sub></th></tr></thead><tbody><tr><td>quater precision (q52)</td><td>1.25 &times; 10<sup>-1</sup></td><td>6.10 &times; 10<sup>-5</sup></td><td>5.73 &times; 10<sup>4</sup></td><td>-16</td><td>-14</td><td>15</td></tr><tr><td>bfloat16 (bf16)</td><td>3.91 &times; 10<sup>-3</sup></td><td>1.18 &times; 10<sup>-38</sup></td><td>3.39 &times; 10<sup>38</sup></td><td>8</td><td>-126</td><td>127</td></tr><tr><td>half precision (fp16)</td><td>4.88 &times; 10<sup>-4</sup></td><td>6.10 &times; 10<sup>-5</sup></td><td>6.55 &times; 10<sup>4</sup></td><td>11</td><td>-14</td><td>15</td></tr><tr><td>single (fp32)</td><td>5.96 &times; 10<sup>-8</sup></td><td>1.18 &times; 10<sup>-38</sup></td><td>3.40 &times; 10<sup>38</sup></td><td>24</td><td>-126</td><td>127</td></tr><tr><td>double (fp64)</td><td>1.11 &times; 10<sup>-16</sup></td><td>2.23 &times; 10<sup>-308</sup></td><td>1.80 &times; 10<sup>308</sup></td><td>53</td><td>-1022</td><td>1023</td></tr></tbody></table>
# Fundamentals

## Lemma 1 ([Ame+22a])

Let $A = U\Sigma V^T$ where U and V have orthonormal columns, and let $\tilde{A} = \hat{U}\Sigma\hat{V}^T$ be obtained by converting U and V to precision $u_k$. Then

$$ \|\tilde{A} - A\| \le (2 + \sqrt{r_k} u_k) u_k \|\Sigma\|. $$
# HODLR matrix

To quantify the error incurred in the low-rank factorization of the off-diagonal blocks (and do error analysis), we introduced the definition.

## Definition 2 ((p, ε)-HODLR matrix)

$\tilde{H}$ is defined to be a $(p, \varepsilon)$-HODLR matrix to a $p$-HODLR matrix $H$, if every off-diagonal block at any level $k$ satisfies

$$ \|\tilde{H}_{ij}^{(k)} - H_{ij}^{(k)}\| \le \varepsilon \|H_{ij}^{(k)}\| \quad (i \ne j), \text{ where } 0 \le \varepsilon < 1. $$

* Rank-constrained → (practically) tolerance-constrained HODLR format.

## Theorem 3 (Approximation error in diagonal blocks)

Let $\tilde{H}$ be a $(p, \varepsilon)$-HODLR matrix associated with $H$. Then for the HODLR matrices $\tilde{H}_{ii}^{(k)}$, $i = 1: 2^k$, at level $k \in \{0, \dots, l\}$, it holds that

$$ \|\tilde{H}_{ii}^{(k)} - H_{ii}^{(k)}\|_F \le \varepsilon \|H_{ii}^{(k)}\|_F, \quad i = 1: 2^k. $$
# HODLR matrix

To quantify the error incurred in the low-rank factorization of the off-diagonal blocks (and do error analysis), we introduced the definition.

## Definition 2 ((p, ε)-HODLR matrix)

$\tilde{H}$ is defined to be a $(p, \varepsilon)$-HODLR matrix to a $p$-HODLR matrix $H$, if every off-diagonal block at any level $k$ satisfies

$$ \|\tilde{H}_{ij}^{(k)} - H_{ij}^{(k)}\| \le \varepsilon \|H_{ij}^{(k)}\| \quad (i \ne j), \text{ where } 0 \le \varepsilon < 1. $$

* Rank-constrained → (practically) tolerance-constrained HODLR format.

## Theorem 3 (Approximation error in diagonal blocks)

Let $\tilde{H}$ be a $(p, \varepsilon)$-HODLR matrix associated with $H$. Then for the HODLR matrices $\tilde{H}_{ii}^{(k)}$, $i = 1: 2^k$, at level $k \in \{0, \dots, \ell\}$, it holds that

$$ \|\tilde{H}_{ii}^{(k)} - H_{ii}^{(k)}\|_F \le \varepsilon \|H_{ii}^{(k)}\|_F, \quad i = 1: 2^k. $$
Questions

1. How can the approximation quality of an HODLR construction be evaluated when using mixed-precision arithmetic?

2. How is the mixed-precision scheme selected to preserve the backward stability of matrix-vector products and LU factorization in the HODLR representation?

3. How should the working precision be chosen in HODLR matrix computations, relative to the approximation error, to minimize the impact of finite-precision effects?
# HODLR matrix

1. $k \leftarrow 0, H \leftarrow H^{(0)}$

2. **while** $k < \ell$ **do**

3. **for** $i := 1$: $2^k$ **do**

4. $\hat{H}_{2i-1,2i}^{(k+1)} \leftarrow \hat{U}_{2i-1}^{(k+1)} (\hat{V}_{2i}^{(k+1)})^T$ [Compute in $u$, and store in $u_{k+1}$]

5. $\hat{H}_{2i,2i-1}^{(k+1)} \leftarrow \hat{U}_{2i}^{(k+1)} (\hat{V}_{2i-1}^{(k+1)})^T$ [Compute in $u$, and store in $u_{k+1}$]

6. **end for**

7. $k \leftarrow k + 1$

8. **end while**

9. **for** $i := 1$: $2^\ell$ **do**

10. $\hat{H}_{i,i}^{(\ell)} \leftarrow H_{i,i}^{(\ell)}$ [Store $\hat{H}_{i,i}^{(\ell)}$ in $u$]

11. **end for**

Define $\xi_k := \max_{|i-j|=1} \|\tilde{H}_{ij}^{(k)}\|_F / \|\tilde{H}\|_F$, $1 \le k \le \ell$, If the precisions are chosen by $u_k \le \varepsilon / (2^{k/2} \xi_k)$, we have $\|H - \tilde{H}\|_F \lesssim (2\sqrt{2\ell} + 1)\varepsilon \|H\|_F$.

* Mixed-precision representation error comparable to $\|H - \tilde{H}\|_F \le \varepsilon \|H\|_F$.
# Mixed precision HODLR construction

## Suitesparse and dense matrix!

<figure><img src="image_5.png" /><figcaption>(a) ex37</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) 1138_bus</figcaption></figure>
# Matrix-vector product

Key steps in the recursive algorithm:

1: $b \leftarrow 0 \in \mathbb{R}^n$

2: **for** $k = 1$: $\ell$ **do**

3: Partition $b$ into $2^k$ blocks, i.e., $(b_i^{(k)})_{i=1}^{2^k}$

4: Partition $x$ into $2^k$ blocks, i.e., $(x_i^{(k)})_{i=1}^{2^k}$

5: **for** $i = 1$: $2^{k-1}$ **do**

6: $b_{2i-1}^{(k)} \leftarrow b_{2i-1}^{(k)} + U_{2i-1}^{(k)}(V_{2i}^{(k)})^T x_{2i}^{(k)}$ [Compute in precision $u$]

7: $b_{2i}^{(k)} \leftarrow b_{2i}^{(k)} + U_{2i}^{(k)}(V_{2i-1}^{(k)})^T x_{2i-1}^{(k)}$ [Compute in precision $u$]

8: **end for**

9: **end for**

10: **for** $i = 1$: $2^\ell$ **do**

11: $b_i^{(\ell)} \leftarrow b_i^{(\ell)} + H_{i,i}^{(\ell)} x_i^{(\ell)}$ [Compute in precision $u$]

12: **end for**
Matrix-vector product

Lemma 4 (Working precision for matrix-vector products)

Let $A_p = X_p \Sigma_p Y_p^T = \sum_{i=1}^p \sigma_i x_i y_i^T$ be the best rank-$p$ approximation of
A. Then the error due to finite precision computation of $\hat{b} = \text{fl}(A_p x)$ will
be no larger than the error due to low-rank approximation when the
working precision has unit roundoff $u \le \sigma_{p+1}/(pn\sigma_1)$.

$\sigma_{p+1}$: the (p+1)th largest singular value.
Matrix-vector product

**Core idea:** balance the approximation error in $\hat{H}$ (mixed-precision HODLR representation) and the finite-precision computation error in $b \leftarrow \hat{H}x$.

* **Intuition:** the more inexact the low-rank representation, the lower the precision we can safely use.

**Theorem 5 (Working precision and backward error)**

If $b = \hat{H}x$ is computed (in the usual fashion of HODLR matrix-vector product) in a working precision $u \le \varepsilon/n$, then the computed $\hat{b}$ satisfies

$$ \hat{b} = \text{fl}(\hat{H}x) = (H + \Delta H)x, \quad \|\Delta H\|_F \le O(2^{\ell/2})\varepsilon\|H\|_F. $$
# Matrix-vector product – kernel matrices

Kernel matrix $K$:

$$ (i) \quad K_{ij} = \begin{cases} \frac{1}{x-y}, & \text{if } x \neq y; \\ 1, & \text{otherwise.} \end{cases} $$

$$ (ii) \quad K_{ij} = \begin{cases} \log \|x_i - x_j\|_2, & \text{if } x \neq y; \\ 0, & \text{otherwise.} \end{cases} \quad (1) $$

$$ (iii) \quad K_{ij} = \exp \left( - \frac{\|x_i - x_j\|_2^2}{2h^2} \right). $$

Evaluated at 1D and 2D point sets:

* Set $\mathbf{s}_1$: A set of uniform grid points in $[0, 1]$.

* Set $\mathbf{s}_2$: A set of uniform grid points in $[-1, 1] \times [-1, 1]$.
Matrix-vector product ($n = 2000$)

The relative backward error is evaluated as

$$ \frac{\|b - \hat{H}x\|_F}{\|K\|_F \|x\|_F} $$

* mat-1 is generated by kernel (i) evaluated at point set $s_1$;

* mat-2 is generated by kernel (ii) evaluated at point set $s_2$;

* mat-3 and mat-4 are generated by kernel (iii) (defined by the Gaussian radial basis function) choosing $h = 1$ and $h = 20$, respectively, evaluated at point set $s_2$.
Matrix-vector product ($n = 2000$)

$u$: fp64, fp32, and bf16 ($u \le \varepsilon/n$).

<figure><img src="image_7.png" /><figcaption>(a) mat-1</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) mat-2</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(c) mat-3</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(d) mat-4</figcaption></figure>
# LU factorization

## HODLR LU factorization: $H \approx LU$:

<figure><img src="image_11.png" /></figure>

### Key steps in the recursive algorithm:

1. $L_{11}, U_{11} \leftarrow \text{HODLR\_LU}(H_{11}^{(k+1)}, T_\ell, k+1)$

2. $U_{12} \leftarrow \text{Solve triangular system } L_{11}U_{12} = H_{12}^{(k+1)}$

3. $L_{21} \leftarrow \text{Solve triangular system } L_{21}U_{11} = H_{21}^{(k+1)}$

4. $H_{22}^\varepsilon \leftarrow H_{22}^{(k+1)} - L_{21}U_{12}$ (possible rank-$p$ truncation)

5. $L_{22}, U_{22} \leftarrow \text{HODLR\_LU}(H_{22}^\varepsilon, k+1)$

6. $L \leftarrow \begin{bmatrix} L_{11} & \\ & L_{21} & L_{22} \end{bmatrix}, U \leftarrow \begin{bmatrix} U_{11} & U_{12} \\ & U_{22} \end{bmatrix}$

* At the bottom level $k = \ell$, dense LU factorization of $H_{ii}^{(\ell)}$ is computed.
# LU factorization

## Theorem 6 (Working precision and backward error)

Let $\hat{H}$ be the mixed-precision $l$-level HODLR representation. If the LU decomposition of $\hat{H}$ is computed in a working precision $u \lesssim \varepsilon/n$, then the LU factorization of $\hat{H}$ satisfies

$$ \hat{L}\hat{U} = H + \Delta H, \quad \|\Delta H\|_F \lesssim \mathcal{O}(2^l)\varepsilon\|H\|_F + \mathcal{O}(2^l)\varepsilon\|\hat{L}\|_F\|\hat{U}\|_F. $$

* For a relatively large $\varepsilon$ (a coarse approximation), we can compute the LU factorization in a low precision without affecting the backward error.

* The bound remains in the same form for the LU factorization of one-precision HODLR matrices, albeit with smaller constant factors.
# LU factorization

$u$: fp64, fp32, and bf16 ($u \le \varepsilon/n$).

<figure><img src="image_12.png" /><figcaption>**Figure:** Backward error of LU factorization of ex37 for mixed precision HODLR matrix.</figcaption></figure>
# LU factorization

$u$: fp64, fp32, and bf16 ($u \le \varepsilon/n$).

<figure><img src="image_13.png" /><figcaption>**Figure:** Backward error of LU factorization of P64 for mixed precision HODLR matrix.</figcaption></figure>
Storage reduction

l = 8; purples bars correspond to ε = 10⁻⁷, green bars indicate ε = 10⁻⁴,
and blue bars indicate ε = 10⁻¹

<figure><img src="image_14.png" /><figcaption>**Figure:** Storage savings of adaptive-precision HODLR matrices relative to uniform (double) precision HODLR matrices.</figcaption></figure>
Storage reduction

ℓ = 8; purples bars correspond to ε = 10⁻⁷, green bars indicate ε = 10⁻⁴,
and blue bars indicate ε = 10⁻¹

<figure><img src="image_15.png" /><figcaption>**Figure:** Storage savings of adaptive-precision HODLR matrices relative to uniform (double) precision HODLR matrices.</figcaption></figure>
Related materials

► Software mhodlr: https://github.com/chenxinye/mhodlr

► Zenodo (experimental code): https://zenodo.org/records/13335429

► Software mhodlr documentations:
https://mhodlr.readthedocs.io/en/latest/? badge=stable

► E. Carson, X. Chen, X. Liu. Mixed precision HODLR matrices.
https://arxiv.org/abs/2407.21637v3. SIAM J. Sci. Comput., To appear.
Datasets

<table><caption>Table: Summary of test matrices from [Ame+22b] and the SuiteSparse collection [DH11].</caption><thead><tr><td>Dataset</td><td>Size</td><td>Nonzeros</td><td>Description</td></tr></thead><tbody><tr><td>1138_bus</td><td>1,138</td><td>4,054</td><td>Power network problem</td></tr><tr><td>bcsstk08</td><td>1,074</td><td>12,960</td><td>Structural problem</td></tr><tr><td>cavity18</td><td>4,562</td><td>138,040</td><td>Computational fluid dynamics problem</td></tr><tr><td>ex37</td><td>3,565</td><td>67,591</td><td>Computational fluid dynamics problem</td></tr><tr><td>LeGresley_2508</td><td>2,508</td><td>16,727</td><td>Power network problem</td></tr><tr><td>P64</td><td>4,096</td><td>16,777,216</td><td>The discretization of a Poisson equation</td></tr><tr><td>psmigr_1</td><td>3,140</td><td>543,160</td><td>Economic problem</td></tr><tr><td>saylr3</td><td>1,000</td><td>3,750</td><td>Computational fluid dynamics problem</td></tr></tbody></table>
Print mphA

The construction of the vanilla HODLR, mixed-precision, and
adaptive-precision HODLR matrices is performed via the following:

1 u1=precision('d'); u2=precision('s'); u3=precision('h');

2 u4=precision('b'); u5=precision('q52');

3 % define the precisions, build the collection of precisions

4 u_chain = prec_chain(u1, u2, u3, u4, u5); % the order
matters!

6 hA=hodlr(A,depth,min_block_size,'svd',epsilon);

7 rA=hA.dense; % recover from the HODLR format

9 % build the precisions according to u_chain (each layer uses
the corresponding precison)

10 mphA=mphodlr(u_chain,A,depth,min_block_size,'svd',epsilon);

11 mprA=mphA.dense; % recover from the HODLR format

13 % build the precisions automatically

14 aphA=amphodlr(u_chain,A,depth,min_block_size,'svd',epsilon);

15 aprA=aphA.dense; % recover from the HODLR format
Construct HODLR matrix

hodlr with properties:

U1: [250×249 double]

V2: [249×250 double]

U2: [250×249 double]

V1: [249×250 double]

D: []

A11: [1×1 hodlr]

A22: [1×1 hodlr]

level: 1

shape: [500 500]

max_level: 5

bottom_level: 5

vareps: 1.0000e-04

min_block_size: 2
# Print mphA

mphodlr with properties:

U1: [250×249 double]

V2: [249×250 double]

U2: [250×249 double]

V1: [249×250 double]

D: []

A11: [1×1 mphodlr]

A22: [1×1 mphodlr]

level: 1

precSettings: {1×5 cell}

shape: [500 500]

max_level: 5

bottom_level: 5

vareps: 1.0000e-04

min_block_size: 2
Print aphA

amphodlr with properties:
U1: [250x249 double]
V2: [249x250 double]
U2: [250x249 double]
V1: [249x250 double]
D: []
A11: [1x1 amphodlr]
A22: [1x1 amphodlr]
level: 1
shape: [500 500]
max_level: 5
bottom_level: 5
normOrder: [8.3423e+04 2.0891e+04 ... ] (1x6 double)
precIndex: [2 2 2 2 2]
unitRoundOff: [1.1102e-16 1.1102e-16 ... ] (1x6 double)
min_block_size: 2
vareps: 1.0000e-04
prec_sets: {1x6 cell}
Visualize structure

After converting the dense matrix into the HODLR representation, one can easily use `mhodlr` to visualize the structure of the HODLR matrix via:

1. `visualize_hodlr(hA)`

The visualization to a HODLR matrix `hA`.
Visualize structure

<figure><img src="image_16.png" /><figcaption>Figure: HODLR matrix structure to a 20 × 20 matrix with depth of 3.</figcaption></figure>
# Matrix computations

**Table:** Summary of typical HODLR matrix operations and their complexity in `mhodlr`. Note: Fused multiply-add complexity assumes optimized implementation; matrix-matrix multiplication may be $O(p^2n \log n)$ in some optimized cases.

<table><thead><tr><th>Operation</th><th>Complexity</th><th>Function call</th><th>Mixed-precision call</th></tr></thead><tbody><tr><td>Matrix-vector mult.</td><td>O(pn log n)</td><td>dot(A, v, ...)</td><td>mdot(A, v)</td></tr><tr><td>Matrix addition</td><td>O(p<sup>2</sup>n log n)</td><td>add(A, B, ...)</td><td>madd(A, B, ...)</td></tr><tr><td>Matrix-matrix mult.</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>matmul(A, B, ...)</td><td>mmatmul(A, B, ...)</td></tr><tr><td>Fused multiply-add</td><td>O(p<sup>2</sup>n log n)</td><td>fma(A, B, C, ...)</td><td>mfma(A, B, C, ...)</td></tr><tr><td>Matrix inversion</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>inverse(A)</td><td>minverse(A)</td></tr><tr><td>LU decomposition</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>hlu(A)</td><td>mhlu(A)</td></tr><tr><td>Cholesky decomp.</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>h chol(A)</td><td>mh chol(A)</td></tr><tr><td>QR decomposition</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>hqr(method, A)</td><td>mhqr(method, A)</td></tr><tr><td>Linear solve</td><td>O(p<sup>2</sup>n log<sup>2</sup> n)</td><td>lrsolve(A, B)</td><td>mlrsolve(A, B)</td></tr></tbody></table>
# Matrix computations

The precision setting is performed by

1 $u = \text{precision('h'); % specify the precision}$

2 $\text{set_precision(u)}$;

Then all the HODLR arithmetic function starting with “m” (see Table 3 for references) will be performed in the specified precision:

1 $C1 = \text{add}(hA, B, hA);$

2 $C2 = A + B + A;$

3 $\text{disp}(\text{norm}(C1 \cdot \text{dense} - C2, 'fro'));$

5 $C1 = \text{sub}(hA, B, hA);$

6 $C2 = A - B - A;$

7 $\text{disp}(\text{norm}(C1 \cdot \text{dense} - C2, 'fro'));$

9 $\text{set_precision(u); % turns on a global precision}$

10 $C1 = \text{madd}(hA, B, hA);$

11 $C2 = A + B + A;$

12 $\text{disp}(\text{norm}(C1 \cdot \text{dense} - C2, 'fro'));$
# Matrix computations

The LU factorization is performed by the routine

1. % Output the factors L and U are hodlr format as default

2. [L, U] = hlu(hA);

3. err = norm(hdot(L, U, 'dense') - A, 'fro');

4. disp(err);

5.

6. % Set the factors L and U to the dense matrix format.

7. [L, U] = hlu(hA, 'dense');

8. err = norm(L * U - A, 'fro');

9. disp(err); % print error

The working precision for LU factorization can be specified by the method `mhlu`:

1. u = precision('h');

2. [L, U] = mhlu(hA, 'hodlr');

3. err = norm(hdot(L, U, 'dense') - A, 'fro');

4. disp(err); % print error
# Matrix computations

We provide three implementations for QR factorizations.

## Working precision

1. `[Q, R] = hqr(hA, 'lt');`

2. `[Q, R] = hqr(hA, 'lt2');`

3. `[Q, R] = hqr(hA, 'dk');`

## Multiple precision

1. `[Q, R] = mhqr(hA, 'lt');`

2. `[Q, R] = mhqr(hA, 'lt2');`

3. `[Q, R] = mhqr(hA, 'dk');`
References I

[DH11] Timothy A. Davis and Yifan Hu. "The University of Florida Sparse Matrix Collection". In: 38.1 (2011). DOI: 10.1145/2049662.2049663.

[Ame+22a] Patrick Amestoy et al. "Mixed precision low-rank approximations and their application to block low-rank LU factorization". In: *IMA Journal of Numerical Analysis* 43.4 (Aug. 2022), pp. 2198-2227. DOI: 10.1093/imanum/drac037.

[Ame+22b] Patrick Amestoy et al. "Mixed precision low-rank approximations and their application to block low-rank LU factorization". In: 43.4 (2022), pp. 2198-2227. DOI: 10.1093/imanum/drac037.