# Mixed precision numerical methods for solving large systems of ordinary differential equations

Mouhamad Al-Sayed Ali, Arsène Marzorati, Samuel Bernard, Jonathan Rouzaud-Cornabas

► To cite this version:

Mouhamad Al-Sayed Ali, Arsène Marzorati, Samuel Bernard, Jonathan Rouzaud-Cornabas. Mixed precision numerical methods for solving large systems of ordinary differential equations. Congrès National d'Analyse Numérique 2024, May 2024, Ile de ré, France. hal-04953553

HAL Id: hal-04953553

https://inria.hal.science/hal-04953553v1

Submitted on 18 Feb 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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Mixed precision numerical methods for solving
large systems of ordinary differential equations

M. Al Sayed Ali, S. Bernard, A. Marzorati, J. Rouzaud Cornabas

INRIA - Lyon

CANUM, Le-Bois-Plage-en-Ré, Ile de Ré, May 27-31, 2024
# Floating-point precisions

## Comparison of floating-point precisions

<table><thead><tr><th>Single precision (FP32)</th><th>Double precision (FP64)</th></tr></thead><tbody><tr><td>· 32 bit memory to represent a value.</td><td>· 64 bit memory to represent a value.</td></tr><tr><td>· 7 decimal digits of precision</td><td>· 16 decimal digits of precision</td></tr><tr><td>· Machine precision ≈ 10<sup>-8</sup></td><td>· Machine precision ≈ 10<sup>-16</sup></td></tr></tbody></table>
# Floating-point precisions

<table><thead><tr><th colspan="2">Comparison of floating-point precisions</th></tr><tr><th>Single precision (FP32)</th><th>Double precision (FP64)</th></tr></thead><tbody><tr><td>· 32 bit memory to represent a value.</td><td>· 64 bit memory to represent a value.</td></tr><tr><td>· 7 decimal digits of precision</td><td>· 16 decimal digits of precision</td></tr><tr><td>· Machine precision ≈ 10<sup>-8</sup></td><td>· Machine precision ≈ 10<sup>-16</sup></td></tr></tbody></table>

• Mixed precision = Single precision + Double precision
Floating-point precisions

Comparison of floating-point precisions

<table><thead><tr><th>Single precision (FP32)</th><th>Double precision (FP64)</th></tr></thead><tbody><tr><td>· 32 bit memory to represent a value.</td><td>· 64 bit memory to represent a value.</td></tr><tr><td>· 7 decimal digits of precision</td><td>· 16 decimal digits of precision</td></tr><tr><td>· Machine precision ≈ 10<sup>-8</sup></td><td>· Machine precision ≈ 10<sup>-16</sup></td></tr></tbody></table>

* Mixed precision = Single precision + Double precision

* Numerical methods with mixed precision run faster and use less memory but with lost of accuracy
# Floating-point precisions

<table><thead><tr><th colspan="2">Comparison of floating-point precisions</th></tr><tr><th>Single precision (FP32)</th><th>Double precision (FP64)</th></tr></thead><tbody><tr><td>· 32 bit memory to represent a value.</td><td>· 64 bit memory to represent a value.</td></tr><tr><td>· 7 decimal digits of precision</td><td>· 16 decimal digits of precision</td></tr><tr><td>· Machine precision ≈ 10<sup>-8</sup></td><td>· Machine precision ≈ 10<sup>-16</sup></td></tr></tbody></table>

* Mixed precision = Single precision + Double precision

* Numerical methods with mixed precision run faster and use less memory but with lost of accuracy

* Domain of applications : AI, Climate/weather modeling, Computational fluid dynamics, ....

* Our application : Biology and Health
• Consider the following ordinary differential equation (ODE) :

$$ \begin{cases} \dot{y}(t) = f(t, y(t)) \quad \forall t \in [t_0, T] \\ y(t_0) = y^{(0)} \end{cases} \tag{1} $$

where

• $y^{(0)}$ is a vector of $\mathbb{R}^n$.

• $f : [t_0, T] \times \mathbb{R}^n \to \mathbb{R}^n$ is a smooth function.
# Numerical methods

* Discretization in time $\Rightarrow$

$$y_{i+1} = F(y_{i-q}, y_{i-(q-1)}, \dots, y_i, f) \quad i = q, \dots, N-1.$$

where

$$\cdot y_i \approx y(t_i),$$

$$\cdot t_i = t_0 + ih \text{ where } h \text{ is the stepsize.}$$
# Numerical methods

* Discretization in time $\Rightarrow$

$$y_{i+1} = F(y_{i-q}, y_{i-(q-1)}, \dots, y_i, f) \quad i = q, \dots, N-1.$$

where

$$\cdot y_i \approx y(t_i),$$

$$\cdot t_i = t_0 + ih \text{ where } h \text{ is the stepsize.}$$

**Class of explicit numerical methods :**

$$y_{i+1} = y_i + h \sum_{l=1}^{q} a_l k_{l,i} \quad (2)$$

where $k_{l,i} = f(t_{l,i}, y_{l,i})$.
# Mixed precision explicit methods

Compute one/more stages $k_{\ell_0, i}$ in a lower precision.

<table><thead><tr><th>Method</th><th>Mixed method</th></tr></thead><tbody><tr><td>Explicit method</td><td>$\tilde{y}_{i+1} = y_i + h \sum_{\substack{l=1 \\ l \neq \ell_0}}^{q} a_l k_{l,i} + h a_{\ell_0} \tilde{k}_{\ell_0,i}$</td></tr></tbody></table>

where $\tilde{k}_{\ell_0, i} = f(t_{\ell_0, i}, y_{\ell_0, i}^{(L)})$ with $y_{\ell_0, i}^{(L)} = \text{Low}(y_{\ell_0, i})$.
What about the error $\max_{0 \le i \le N} \|y_i - \tilde{y}_i\|$?
Theorem

Let's assume that the function $f(t, y)$ is Lipschitz with respect to the variable $y$ with constant $c_0$. Then:

* If, for all iterations $i$, we have computed an approximation $\tilde{k}_{\ell_0,i}$ of $k_{\ell_0,i}$ in a lower precision such that $\|y_{\ell_0,i}^{(L)} - y_{\ell_0,i}\| \le \varepsilon$. Then we have

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| \le (T - t_0) |a_{\ell_0} c_0| \varepsilon + \mathcal{O}(h\varepsilon). $$

* If we compute the numerical method in single precision then we get

$$ \max_{0 \le i \le N} \|y_i - \tilde{y}_i\| = \mathcal{O}\left(\frac{\varepsilon}{h}\right). \blacksquare $$
Explicit Runge-Kutta 4

Runge-Kutta 4 (RK4) :

$$y_{i+1} = y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$$

with

$$k_1 = f(t_i, y_i),$$

$$k_2 = f(t_i + \frac{h}{2}, y_i + \frac{hk_1}{2}),$$

$$k_3 = f(t_i + \frac{h}{2}, y_i + \frac{hk_2}{2}),$$

$$k_4 = f(t_{i+1}, y_i + \frac{hk_3}{2}).$$
# Mixed precision explicit numerical methods

Compute one/more stages $k_i$ in a lower precision.

<table><caption>Table: $A_i \in \{S, D\}$ the precision of the stage $k_i$</caption><thead><tr><td>Method</td><td>Expression</td><td>Mixed method</td></tr></thead><tbody><tr><td>RK4</td><td>$\tilde{y}_{i+1} = y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$</td><td>$A_1A_2A_3A_4$</td></tr></tbody></table>
# Mixed precision explicit numerical methods

Compute one/more stages $k_i$ in a lower precision.

<table><caption>Table: $A_i \in \{S, D\}$ the precision of the stage $k_i$</caption><thead><tr><td>Method</td><td>Expression</td><td>Mixed method</td></tr></thead><tbody><tr><td>RK4</td><td>$\tilde{y}_{i+1} = y_i + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$</td><td>$A_1A_2A_3A_4$</td></tr></tbody></table>

For example, `SSDD` means that we have computed $k_1, k_2$ in single precision, $k_3, k_4$ in double precision and all other operations are performed in double precision.
## Computational environment

* **Compiler** : Fortran

* **Parallel computing** : MPI, Ethernet, 25 Gbps

* **Cluster** : Gros on Grid5000, Intel Xeon Gold 5220, 2.20GHz

* **\# processors** : 250
# Numerical Experiments

Let's consider the following benchmark model, a large nonlinear densely coupled ODEs<fnref n="1" />, ($S_i$)$_{i=1,...,N_h}$, where ($S_i$) is given by:

$$ 
\begin{aligned}
\frac{dy_1^{(i)}}{dt} &= \frac{1}{\tau} \left( \frac{\nu_{1b}(y_7^{(i)} + \psi_{1i})}{k_{1b}\left(1 + \left(\frac{y_3^{(i)}}{k_{1i}}\right)^{\rho_0}\right) + y_7^{(i)} + \psi_{1i}} - k_{1d}y_1^{(i)} \right) \\
\frac{dy_2^{(i)}}{dt} &= \frac{1}{\tau} \left( k_{2b}(y_1^{(i)})^q - k_{2d}y_2^{(i)} - k_{2t}y_2^{(i)} + k_{3t}y_3^{(i)} \right) \\
\frac{dy_3^{(i)}}{dt} &= \frac{1}{\tau} \left( k_{2t}y_2^{(i)} - k_{3t}y_3^{(i)} - k_{3d}y_3^{(i)} \right)
\end{aligned}
$$

<footnote n="1">¹R. El cheikh, S. Bernard, N. El Khatib, *A multiscale modelling approach for the regulation of the cell cycle by the circadian clock*, Journal of Theoretical Biology 426 (2017) 117-125</footnote>
$$ \frac{dy_4^{(i)}}{dt} = \frac{1}{\tau} \left( \frac{\nu_{4b}(y_3^{r_0})^{(i)}}{k_{4b}^{r_0} + (y_3^{r_0})^{(i)}} - k_{4d}y_4^{(i)} \right) $$

$$ \frac{dy_5^{(i)}}{dt} = \frac{1}{\tau} \left( k_{5b}y_5^{(i)} - k_{5d}y_5^{(i)} - k_{5t}y_5^{(i)} + k_{6t}y_6^{(i)} \right) $$

$$ \frac{dy_6^{(i)}}{dt} = \frac{1}{\tau} \left( k_{5t}y_5^{(i)} - k_{6t}y_6^{(i)} - k_{6d}y_6^{(i)} + k_{7a}y_7^{(i)} - k_{6a}y_6^{(i)} \right) $$

$$ \frac{dy_7^{(i)}}{dt} = \frac{1}{\tau} \left( k_{6a}y_6^{(i)} - k_{7a}y_7^{(i)} - k_{7d}y_7^{(i)} \right) $$

$$ \frac{dy_8^{(i)}}{dt} = \lambda \left( \frac{(k_{\text{impf}} + k_{\text{ompf}} \exp(-\eta\psi_2))k_{\text{1mpf}}^{n_0}}{(k_{\text{1mpf}}^{n_0} + (y_8^{(i)})^{n_0} + \text{sy}_{10}^{n_0})(1 - y_8^{(i)})} - d_{\text{wee}}y_9^{(i)}y_8^{(i)} \right) $$
$$ \frac{dy_9^{(i)}}{dt} = \lambda \left( \frac{k_{actw}}{k_{actw} + dw_1} (cw + C(y_7 - bbmal_0) + bbmal_0) + \dots \right) - \left( \frac{k_{actw}}{k_{actw} + dw_1} - 1 \right) \frac{kin_{actw}(y_8^{(i)})^{n_0} y_9^{(i)}}{k^{n_0}_{1\_wee_1} + (y_8^{(i)})^{n_0}} - dw_2 y_9^{(i)} $$

$$ \frac{dy_{10}^{(i)}}{dt} = \lambda (k k_{act} (y_8^{(i)} - y_{10}^{(i)})) $$

where $t \in [0, 120]$ and $\psi_{1i}, \psi_{2i}$, for $i = 1, \dots, N_h$, are given by

$$ \psi_{1i} = \frac{ks}{N_h} \sum_{j=1}^{N_h} \arctan(y_2^{(j)} - y_2^{(i)}) + \frac{\pi}{2} ks $$
$$ \psi_{2i} = N_h = 10^4. $$
* Let $y_{ref} \approx y(120)$ be the solution, in double precision, computed with the step size $h_{ref} = \frac{120}{5 \times 10^5}$.

* The step size $h$ in the RK4 method is $120/10^5$.
<table><caption>Table: N = 100 000, ||y<sub>ref</sub>||<sub>∞</sub> = 2.5923</caption><thead><tr><th>Method</th><th>Time</th><th>||y-y<sub>ref</sub>||<sub>∞</sub> / ||y<sub>ref</sub>||<sub>∞</sub></th><th>Time (Double) / Time (Mixed)</th></tr></thead><tbody><tr><td>DOUBLE</td><td>7 002</td><td>1.05 × 10<sup>-7</sup></td><td>1</td></tr><tr><td>DDDS</td><td>6 237</td><td>1.26 × 10<sup>-7</sup></td><td>1.12</td></tr><tr><td>SDSD</td><td>5 547</td><td>1.02 × 10<sup>-7</sup></td><td>1.26</td></tr><tr><td>DDSS</td><td>5 400</td><td>1.02 × 10<sup>-7</sup></td><td>1.29</td></tr><tr><td>DSSD</td><td>5 389</td><td>8.99 × 10<sup>-8</sup></td><td>1.3</td></tr><tr><td>DSDS</td><td>5 432</td><td>1.17 × 10<sup>-7</sup></td><td>1.28</td></tr><tr><td>SDDS</td><td>5 464</td><td>1.41 × 10<sup>-7</sup></td><td>1.28</td></tr><tr><td>DSSS</td><td>4 563</td><td>1.01 × 10<sup>-7</sup></td><td>1.53</td></tr><tr><td>SDSS</td><td>4 607</td><td>1.21 × 10<sup>-7</sup></td><td>1.51</td></tr><tr><td>SSDS</td><td>4 574</td><td>1.87 × 10<sup>-7</sup></td><td>1.53</td></tr><tr><td>SSSD</td><td>4 551</td><td>1.01 × 10<sup>-7</sup></td><td>1.53</td></tr><tr><td>SSSS</td><td>3 683</td><td>1.13 × 10<sup>-7</sup></td><td>1.9</td></tr><tr><td>SINGLE</td><td>3 467</td><td>2.2 × 10<sup>-4</sup></td><td>2.01</td></tr></tbody></table>
<figure><img src="image_2.png" /></figure>
* We have implemented, using mixed precisions arithmetic, numerical methods (explicit and *implicit*), on CPU cluster with MPI, for solving large systems of ODE's.

* These methods *run faster* and *use less memory* with the *same quality* as methods running in "full precision".

* The numerical results show the *parallel methods running in mixed precision* are up to 1.5 – 2 times faster than those running in double precision with sufficient accuracy.
* Implemented Mixed precision's numerical methods on a Tensor Processing Unit (matrix processor).

* Tested our approach on another numerical methods.

* Use other type of precisions.
Thank you for your attention.
# Order of the Mixed method DSDS

<figure><img src="image_3.png" /></figure>
<table><thead><tr><th>Method</th><th>Time</th><th>||y<sub>D</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||y<sub>S</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||<sup>y<sub>D</sub> - y<sub>M</sub></sup>&frasl;<sub>y<sub>D</sub></sub>||<sub>∞</sub></th></tr></thead><tbody><tr><td>DOUBLE</td><td>7 002</td><td>0</td><td>5.8441 &times; 10<sup>-4</sup></td><td>0</td></tr><tr><td>DDDS</td><td>6 237</td><td>5.9679 &times; 10<sup>-8</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>2.5640 &times; 10<sup>-7</sup></td></tr><tr><td>DDSS</td><td>5 400</td><td>1.4211 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>4.9659 &times; 10<sup>-7</sup></td></tr><tr><td>DSDS</td><td>5 432</td><td>1.1512 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>4.8019 &times; 10<sup>-7</sup></td></tr><tr><td>DSSD</td><td>5 389</td><td>2.1592 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>9.0710 &times; 10<sup>-7</sup></td></tr><tr><td>SDDS</td><td>5 464</td><td>1.1089 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>4.9863 &times; 10<sup>-7</sup></td></tr><tr><td>SDSD</td><td>5 547</td><td>1.4250 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>4.9758 &times; 10<sup>-7</sup></td></tr><tr><td>DSSS</td><td>4 563</td><td>2.3106 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>9.3677 &times; 10<sup>-7</sup></td></tr><tr><td>SDSS</td><td>4 607</td><td>1.7184 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>5.9969 &times; 10<sup>-7</sup></td></tr><tr><td>SSDS</td><td>4 574</td><td>2.7871 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>1.0509 &times; 10<sup>-6</sup></td></tr><tr><td>SSSD</td><td>4 551</td><td>2.3202 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>9.3439 &times; 10<sup>-7</sup></td></tr><tr><td>SSSS</td><td>3 683</td><td>2.3963 &times; 10<sup>-7</sup></td><td>5.8441 &times; 10<sup>-4</sup></td><td>9.2155 &times; 10<sup>-7</sup></td></tr><tr><td>SINGLE</td><td>3 467</td><td>5.8441 &times; 10<sup>-4</sup></td><td>0</td><td>2.0789 &times; 10<sup>-3</sup></td></tr></tbody></table>
<table><caption>Table: Runge-Kutta 2, Nova, #procs = 128</caption><thead><tr><th>Method</th><th>Time</th><th><div style="text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">||y-y<sub>ref</sub>||<sub>∞</sub></div><div>||y<sub>ref</sub>||<sub>∞</sub></div></div></th><th><div style="text-align: center;"><div style="border-bottom: 1px solid black; padding-bottom: 1px;">Time(<i>Double</i>)</div><div>Time(<i>Mixed</i>)</div></div></th></tr></thead><tbody><tr><td><span style="color: blue;">DOUBLE</span></td><td>8 334</td><td>1.28 × 10<sup>-5</sup></td><td>1</td></tr><tr><td><span style="color: red;">DS</span></td><td>4 563</td><td>1.28 × 10<sup>-5</sup></td><td>1.8</td></tr><tr><td><span style="color: red;">SD</span></td><td>4 453</td><td>1.28 × 10<sup>-5</sup></td><td>1.87</td></tr><tr><td><span style="color: red;">SS</span></td><td>3 827</td><td>1.28 × 10<sup>-5</sup></td><td>2.1</td></tr><tr><td><span style="color: blue;">SINGLE</span></td><td>3 621</td><td>2.39 × 10<sup>-4</sup></td><td>2.3</td></tr></tbody></table>
<table><caption>Table: Runge-Kutta 2, Nova, #procs = 128</caption><thead><tr><th>type</th><th>Time</th><th>||y<sub>D</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||y<sub>S</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||<sup>y<sub>D</sub> - y<sub>M</sub></sup>&frasl;<sub>y<sub>D</sub></sub>||<sub>∞</sub></th></tr></thead><tbody><tr><td>DOUBLE</td><td>8 334</td><td>0</td><td>6.0503 &times; 10<sup>-4</sup></td><td>0</td></tr><tr><td>DS</td><td>4 563</td><td>2.6678 &times; 10<sup>-7</sup></td><td>6.0513 &times; 10<sup>-4</sup></td><td>9.9485 &times; 10<sup>-7</sup></td></tr><tr><td>SD</td><td>4 453</td><td>1.7214 &times; 10<sup>-7</sup></td><td>6.0503 &times; 10<sup>-4</sup></td><td>5.8557 &times; 10<sup>-7</sup></td></tr><tr><td>SS</td><td>3 827</td><td>3.3097 &times; 10<sup>-7</sup></td><td>6.0509 &times; 10<sup>-4</sup></td><td>1.1138 &times; 10<sup>-6</sup></td></tr><tr><td>SINGLE</td><td>3 621</td><td>6.0503 &times; 10<sup>-4</sup></td><td>0</td><td>2.1676 &times; 10<sup>-3</sup></td></tr></tbody></table>
<table><caption>Table: Adams Bashforth 2, Nova, #procs = 128</caption><thead><tr><th>Method</th><th>Time</th><th>||y - y<sub>ref</sub>||<sub>∞</sub><br>||y<sub>ref</sub>||<sub>∞</sub></th><th>Time(Double)<br>Time(Mixed)</th></tr></thead><tbody><tr><td>DOUBLE</td><td>3 509</td><td>3.4 × 10<sup>-5</sup></td><td>1</td></tr><tr><td>SS</td><td>1 942</td><td>3.4 × 10<sup>-5</sup></td><td>1.8</td></tr><tr><td>SINGLE</td><td>1 188</td><td>2.3 × 10<sup>-4</sup></td><td>2.3</td></tr></tbody></table>
<table><caption>Table: Adams Bashforth 2, Nova, #procs = 128</caption><thead><tr><th>Method</th><th>Time</th><th>||y<sub>D</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||y<sub>S</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||<sup>y<sub>D</sub> - y<sub>M</sub></sup>&frasl;<sub>y<sub>D</sub></sub>||<sub>∞</sub></th></tr></thead><tbody><tr><td><b>DOUBLE</b></td><td>3 509</td><td>0</td><td>5.6 × 10<sup>-4</sup></td><td>0</td></tr><tr><td><b>SS</b></td><td>1 942</td><td>3.4 × 10<sup>-7</sup></td><td>5.6 × 10<sup>-4</sup></td><td>1.331 × 10<sup>-6</sup></td></tr><tr><td><b>SINGLE</b></td><td>1 188</td><td>5.6 × 10<sup>-4</sup></td><td>0</td><td>2.026 × 10<sup>-3</sup></td></tr></tbody></table>
Thank you for your attention.
# Implicit methods for solving the ODE

We have

$$y_{i+1} = a_i + h\beta_0 f(t_{i+1}, y_{i+1}) \quad i = q, \dots, N-1,$$

where $a_i = F(y_i, y_{i-1}, \dots, y_{i-q+1})$.

Let $x_i = \frac{y_{i+1} - a_i}{\beta_0}$ then the numerical method becomes

$$G_i(x_i) = 0, \quad (3)$$

where

$$G_i(x) = x - hf(t_{i+1}, a_i + \beta_0 x).$$
<figure><img src="image_4.png" /></figure>

The nonlinear system $G_i(x_i) = 0$ is solved by the Newton method.
<figure><img src="image_5.png" /></figure>

The nonlinear system $G_i(x_i) = 0$ is solved by the Newton method.

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

The Newton method computes a sequence $(x_i^{(m)})$ such that
$$ \lim_{m \to +\infty} x_i^{(m)} = x_i \text{ when } x_i^{(0)} \approx x_i. $$
<figure><img src="image_7.png" /></figure>

The nonlinear system $G_i(x_i) = 0$ is solved by the Newton method.

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

The Newton method computes a sequence $(x_i^{(m)})$ such that
$$ \lim_{m \to +\infty} x_i^{(m)} = x_i \quad \text{when} \quad x_i^{(0)} \approx x_i. $$

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

This sequence is given by $x_i^{(m+1)} = x_i^{(m)} + z_i^{(m)}$ where $z_i^{(m)}$ is the solution of

$$ G_i'(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)}). $$
The nonlinear system $G_i(x_i) = 0$ is solved by the Newton method. ●

The Newton method computes a sequence $(x_i^{(m)})$ such that $\lim_{m \to +\infty} x_i^{(m)} = x_i$ when $x_i^{(0)} \approx x_i$. ●

This sequence is given by $x_i^{(m+1)} = x_i^{(m)} + z_i^{(m)}$ where $z_i^{(m)}$ is the solution of

$$G_i'(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)}).$$

To compute good initial guess $x_i^{(0)}$ we use the Line-search algorithm. ●
Line-search algorithm for computing $x_i^{(0)}$

* The goal of this algorithm is to compute a local minimum of the function $H_i$ defined by $H_i(x) = \frac{1}{2} \|G_i(x)\|^2$.
Line-search algorithm for computing $x_i^{(0)}$

* The goal of this algorithm is to compute a local minimum of the function $H_i$ defined by $H_i(x) = \frac{1}{2} \|G_i(x)\|^2$.

* It computes a sequence $(u_i^{(k)})_{k \ge 0}$ s.t. $\nabla H_i(u_i^{(k)}) \underset{k \to +\infty}{\to} 0$.
Line-search algorithm for computing $x_i^{(0)}$

• The goal of this algorithm is to compute a local minimum of the function $H_i$ defined by $H_i(x) = \frac{1}{2} \|G_i(x)\|^2$.

• It computes a sequence $(u_i^{(k)})_{k \ge 0}$ s.t. $\nabla H_i(u_i^{(k)}) \xrightarrow{k \to +\infty} 0$.

• $u_i^{(k+1)} = u_i^{(k)} + \lambda_k p_i^{(k)}$ where $p_i^{(k)}$ solves

$$ (G_i'(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)}). $$
Line-search algorithm for computing $x_i^{(0)}$

• The goal of this algorithm is to compute a local minimum of the function $H_i$ defined by $H_i(x) = \frac{1}{2} \|G_i(x)\|^2$.

• It computes a sequence $(u_i^{(k)})_{k \ge 0}$ s.t. $\nabla H_i(u_i^{(k)}) \xrightarrow{k \to +\infty} 0$.

• $u_i^{(k+1)} = u_i^{(k)} + \lambda_k p_i^{(k)}$ where $p_i^{(k)}$ solves

$$ (G_i'(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)}). $$

• $\lambda_k$ is a **positive** scalar chosen to satisfy the Goldstein-Armijo condition:

$$ h_i(u_i^{(k+1)}) \le h_i(u_i^{(k)}) + \alpha \lambda_k \nabla H_i(u_i^{(k)})^T p_i^{(k)} \quad (4) $$

where $\alpha \in (0, 1/2)$ is a parameter typically set to $10^{-4}$.
Solution of the linear systems $G'_i(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)})$ and $(G'_i(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)})$.
Solution of the linear systems $G_i'(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)})$ and $(G_i(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)})$.

* These linear systems are solved by an iterative method.
Solution of the linear systems $G_i'(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)})$ and $(G_i'(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)})$.

* These linear systems are solved by an iterative method.

* This method computes an approximation $\tilde{x}$ of a linear system, of the form $Ax = b$, such that

$$ \|A\tilde{x} - b\| < \varepsilon_1 \|b\|, \quad (5) $$

where $\varepsilon_1 > 0$ is some tolerance threshold.
Solution of the linear systems $G_i'(x_i^{(m)})z_i^{(m)} = -G_i(x_i^{(m)})$ and $(G_i(u_i^{(k)}))p_i^{(k)} = -G_i(u_i^{(k)})$.

• These linear systems are solved by an iterative method.

• This method computes an approximation $\tilde{x}$ of a linear system, of the form $Ax = b$, such that

$$ \|A\tilde{x} - b\| < \varepsilon_1 \|b\|, \quad (5) $$

where $\varepsilon_1 > 0$ is some tolerance threshold.

• We will use the following mixed precision iterative method.
**Algorithm 1** Mixed Precision GMRES (Generalized Minimal RESidual)
for $Az = b$

Input: $z_0$ initial guess

1: for $i = 0, \dots, imax$ or until converged

2: Compute $r_i = b - Az_i$ in **double** precision.

3: Solve $Au = r_i$ in **single** precision.

4: Update $z_{i+1} = z_i + u$ in **double** precision.

5: end
## Theorem

Let's suppose that for any vector $v$, we have $||v - v^{(L)}|| = O(\epsilon)$. Let's $x^{(M)}$ and $x^{(D)}$ denote, respectively, the solution computed by Mixed precision GMRES and double precision GMRES. Then we have:

$$||b - Ax^{(M)}|| = ||b - Ax^{(D)}|| + O(\epsilon).$$

## Remark

If $\epsilon$ is small enough then the mixed solution $x^{(M)}$ satisfies

$$||b - Ax^{(M)}|| \le c||b||$$

where $c < 1$, and therefore $\tilde{p}_i^{(k)}$ is a descent direction and the Newton method will be locally convergent.
# Mixed precision implicit numerical method (MPIN)

(MPIN1) Set $i = 0$. For $i \le NT$, compute the numerical method $y_i$ as follows:

A) Line-search (LS) algorithm for computing $x_i^{(0)}$:

(LS1) Let $u_i^{(0)} = 0$. While $||G_i(u_i^{(k)})|| > \varepsilon$, do:

(LS2) Compute an approximation $\tilde{p}_i^{(k)} \approx p_i^{(k)}$ by **MP GMRES**.

(LS3) Compute $u_i^{(k+1)} = u_i^{(k)} + \lambda_k \tilde{p}_i^{(k)}$.

(LS4) Set $x_i^{(0)} = u_i^{(k)}$.
## B) Inexact Newton (IN):

(IN1) Let $m = 0$ and repeat until $||G_i(x_i^{(m)})|| \le \varepsilon$:

(IN2) Compute an approximation $\tilde{z}_i^{(m)} \approx z_i^{(m)}$ by **MP GMRES**.

(IN3) $x_i^{(m+1)} = x_i^{(m)} + \tilde{z}_i^{(m)}$.

(MPIN2) Set $x_i = x_i^{(m)}$.

(MPIN3) Compute $y_{i+1} = a_i + \beta_0 x_i$.

(MPIN4) $i = i + 1$.

(MPIN5) $y_M = y_{NT}$.
* Let $T_M$ be the total time in seconds for computing $(y_i)$ by the above algorithm.

* $(T_D, Y_D)$ for Double precision algorithm.

* $(T_S, Y_S)$ for Single precision algorithm.
In the table below, we present the numerical results when the circadian clock's problem is solved by the BDF method of order 2 that is given by

$$y_{i+2} = \frac{4}{3}y_{i+1} - \frac{1}{3}y_i + \frac{2}{3}hf(t_{i+2}, y_{i+2}).$$
<table><caption>Table: BDF2 scheme, with $n = 100 000$, cluster = nova, nbprocs = 128.</caption><thead><tr><th>NT</th><th>T<sub>S</sub></th><th>T<sub>M</sub></th><th>T<sub>D</sub></th><th><sup>||y<sub>ref</sub> - y<sub>S</sub>||<sub>∞</sub></sup> / ||y<sub>ref</sub>||<sub>∞</sub></th><th><sup>||y<sub>ref</sub> - y<sub>M</sub>||<sub>∞</sub></sup> / ||y<sub>ref</sub>||<sub>∞</sub></th><th><sup>||y<sub>ref</sub> - y<sub>D</sub>||<sub>∞</sub></sup> / ||y<sub>ref</sub>||<sub>∞</sub></th></tr></thead><tbody><tr><td>1000</td><td>2 394</td><td>2 900</td><td>4 125</td><td>8.3 × 10<sup>-3</sup></td><td>8.2 × 10<sup>-3</sup></td><td>8.3 × 10<sup>-3</sup></td></tr><tr><td>10000</td><td>5 565</td><td>8 387</td><td>10 051</td><td>1.3 × 10<sup>-3</sup></td><td>3.9 × 10<sup>-4</sup></td><td>3.9 × 10<sup>-4</sup></td></tr><tr><td>20000</td><td>15 190</td><td>22 943</td><td>26 773</td><td>6.4 × 10<sup>-3</sup></td><td>5.5 × 10<sup>-4</sup></td><td>5.5 × 10<sup>-4</sup></td></tr></tbody></table>
<table><caption>Table: BDF2 scheme, with $n = 100 000$, cluster = nova, nbprocs = 128.</caption><thead><tr><th>NT</th><th>||y<sub>D</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||y<sub>S</sub> - y<sub>M</sub>||<sub>∞</sub></th><th>||y<sub>D</sub> - y<sub>S</sub>||<sub>∞</sub></th><th>||<sup>y<sub>D</sub> - y<sub>M</sub></sup>/<sub>y<sub>D</sub></sub>||<sub>∞</sub></th><th>||<sup>y<sub>D</sub> - y<sub>S</sub></sup>/<sub>y<sub>D</sub></sub>||<sub>∞</sub></th></tr></thead><tbody><tr><td>1000</td><td>9.8 × 10<sup>-14</sup></td><td>5 × 10<sup>-4</sup></td><td>5 × 10<sup>-4</sup></td><td>5 × 10<sup>-13</sup></td><td>1.9 × 10<sup>-3</sup></td></tr><tr><td>10000</td><td>1.1 × 10<sup>-4</sup></td><td>3.2 × 10<sup>-3</sup></td><td>3.2 × 10<sup>-3</sup></td><td>3.4 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-2</sup></td></tr><tr><td>20000</td><td>4.4 × 10<sup>-16</sup></td><td>1.7 × 10<sup>-2</sup></td><td>1.7 × 10<sup>-2</sup></td><td>6.2 × 10<sup>-16</sup></td><td>6.1 × 10<sup>-2</sup></td></tr></tbody></table>
<table><caption>Table: BDF2 scheme, with $n = 100 000$, cluster = nova, nbprocs = 128.</caption><thead><tr><th>NT</th><th>1/N<sub>h</sub> &sum;<sub>l=1</sub><sup>N<sub>h</sub></sup> &radic; (y<sub>D</sub>(l) - y<sub>M</sub>(l)) / y<sub>D</sub>(l)</th><th>1/N<sub>h</sub> &sum;<sub>l=1</sub><sup>N<sub>h</sub></sup> &radic; (y<sub>D</sub>(l) - y<sub>S</sub>(l)) / y<sub>D</sub>(l)</th></tr></thead><tbody><tr><td>1000</td><td>1.0 &times; 10<sup>-13</sup></td><td>2.4 &times; 10<sup>-4</sup></td></tr><tr><td>10000</td><td>5.8 &times; 10<sup>-8</sup></td><td>1.5 &times; 10<sup>-3</sup></td></tr><tr><td>20000</td><td>4.4 &times; 10<sup>-17</sup></td><td>8.8 &times; 10<sup>-3</sup></td></tr></tbody></table>
Thank you for your attention.