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

HAL
open science

CGM-based parallel solutions for a class of non-serial
polyadic dynamic-programming problems

Jerry LACMOU ZEUTOUO

► To cite this version:

Jerry LACMOU ZEUTOUO. CGM-based parallel solutions for a class of non-serial polyadic dynamic-programming problems. Distributed, Parallel, and Cluster Computing [cs.DC]. Université de Dschang (Cameroun), 2022. English. NNT: . tel-03936771

HAL Id: tel-03936771

https://hal.science/tel-03936771

Submitted on 12 Jan 2023

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.

Public Domain
REPUBLIC OF CAMEROON
PEACE-WORK-FATHERLAND

UNIVERSITY OF DSCHANG

POST GRADUATE SCHOOL

RÉPUBLIQUE DU CAMEROUN
PAIX-TRAVAIL-PATRIE

UNIVERSITÉ DE DSCHANG

ÉCOLE DOCTORALE

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

DSCHANG SCHOOL OF SCIENCE AND TECHNOLOGY

FUNDAMENTAL COMPUTER SCIENCE, ENGINEERING, AND
APPLICATIONS RESEARCH UNIT (URIFIA)

TOPIC :

CGM-BASED PARALLEL SOLUTIONS FOR A CLASS OF NON-SERIAL
POLYADIC DYNAMIC-PROGRAMMING PROBLEMS

Thesis publicly defended for the award of a Doctorat/PhD degree in
Computer Science

Option : Networks and Distributed Services

Specialty : Parallel Computing

By

LACMOU ZEUTOUO Jerry

Registration Number : CM-UDS-11SCI0267

Master of Science in Computer Science

Under the direction of

KENGNE TCHENDJI Vianney

Associate Professor

LÉLÉ Célestin

Professor

On December 21, 2022 in front of the examination panel consisting of

President

TAYOU DJAMEGNI Clémentin Professor University of Dschang

Reporters

LÉLÉ Célestin Professor University of Dschang

KENGNE TCHENDJI Vianney Associate Professor University of Dschang

Examinators

TIEUDJO Daniel Professor University of Ngaoundere

NDOUNDAM René Associate Professor University of Yaounde I

DJOTIO NDIE Thomas Associate Professor University of Yaounde I
# Dedication

To my parents :
LACMOU Jean and KUELA SEGNOU Antoinette
Acknowledgements

First of all, I would like to thank the Almighty Lord for the graces and mercies He has granted upon me throughout this work.

I have lived a fascinating adventure during these years of thesis work. It made me grow in both my professional and personal lives. I faced difficult moments that almost made me give up. But thanks to my supervisors, Pr. LÉLÉ Célestin and Pr. KENGNE TCHENDJI Vianney, who encouraged me, took care of me and gave me all the necessary materials to carry out this work, I did not give up. You have made me the man I am today. I also thank Pr. MYOUPO Jean Frédéric who took a keen regard over my work. It is through him that we were able to carry out our experimentations on the Dolphin cluster of the MatriCS platform of the University of Picardie Jules Verne.

I thank all the jury members for the honor you have bestowed upon me in judging this work. I also express my gratitude to the teaching staff of the Department of Mathematics and Computer Science. Their invaluable advice and expertise have allowed me to enrich my knowledge in many fields since I entered the University of Dschang in 2011.

What about my parents, LACMOU Jean and KUELA SEGNOU Antoinette, to whom I dedicate this work? Since I was born, you have made a lot of effort to ensure my education. You have labored and patiently waited until the completion of this thesis. I can only say thank you. I cannot forget the full support of my brothers and sisters, Olivier, Marlène, Dimitry, and Rhina. May the Almighty Lord infinitely bless the DIFO family, the FOGAING family, the SOPFOSSI family, the LOBE family, the SOUMANA family, the NIMPA family, the ESSOUNGUE family, the NGOUEGNI family, the MEBANN family, and the DJOUMESSI family.

My beloved friend, MEBANN Elodie, has always taken care of me throughout this work. I will always remember all the sacrifices and risks you took for me. My fan club did not abandon me and did not stop encouraging me : NGUEMKAP Romuald, GUIFO Yvan, DJOUMESSI Kerol, NOUKELA Christian, TESSA Xavier, DJOUMESSI Gina, NANKEM Ida, and LAMGOAH Marie-Noël.
I thank my research team. Although we are in different fields, we were able to make ourselves understood to each other. My elders, Dr. YANKAM Yannick and Dr. NKONJOH Armel, have contributed a lot to the completion of this work. The same goes for my cadets MVAH Fabrice, BOGNING Hermann, KAMGA Ingrid, TESSA Colette, KOMBOU Carole, and DJEUFACK Vadèle. The glances I have had at your respective works have significantly enhanced the quality of mine.

Thanks to the board of the Koossery Technology company, especially to Mr. BAKENEGHE Jean-Claude. Without him, this work would be hopeless because he allowed me to survive financially. He taught me rigor in professional work, which is obviously reflected in this thesis.

I will be eternally grateful to all the people who contributed in any way to the writing of this thesis, especially to all those who have proofread this document.
# Table of Contents

<table><tr><td>Dedication</td><td>i</td></tr><tr><td>Acknowledgements</td><td>ii</td></tr><tr><td>Abstract</td><td>vii</td></tr><tr><td>Résumé</td><td>ix</td></tr><tr><td>List of Acronyms</td><td>x</td></tr><tr><td>List of Tables</td><td>xi</td></tr><tr><td>List of Figures</td><td>xiii</td></tr><tr><td>List of Algorithms</td><td>xvii</td></tr><tr><td>General Introduction</td><td>1</td></tr><tr><td>Background</td><td>1</td></tr><tr><td>Research problem</td><td>5</td></tr><tr><td>Research aim</td><td>6</td></tr><tr><td>Our contributions</td><td>7</td></tr><tr><td>Thesis outline</td><td>8</td></tr><tr><td>Chapter 1 • Parallel Computing and Dynamic Programming</td><td>9</td></tr><tr><td>1.1 Introduction</td><td>9</td></tr><tr><td>1.2 High-performance computing</td><td>9</td></tr><tr><td>1.3 Taxonomies of parallel computer architectures</td><td>11</td></tr><tr><td>1.3.1 Classification according to the number of instruction streams and data streams</td><td>11</td></tr><tr><td>1.3.2 Classification according to the memory</td><td>12</td></tr><tr><td>1.3.3 Classification according to the network topology</td><td>14</td></tr><tr><td>1.3.4 Classification according to the granularity</td><td>14</td></tr></table>
<table><tr><td>1.4</td><td>Designing parallel algorithms</td><td>16</td></tr><tr><td>1.4.1</td><td>Manual parallelization versus automatic parallelization</td><td>16</td></tr><tr><td>1.4.2</td><td>Control parallelism versus data parallelism</td><td>18</td></tr><tr><td>1.4.3</td><td>Parallel programming models</td><td>18</td></tr><tr><td>1.4.4</td><td>Performance of a parallel algorithm</td><td>22</td></tr></table>

<table><tr><td>1.5</td><td>Parallel computing models</td><td>24</td></tr><tr><td>1.5.1</td><td>PRAM model</td><td>25</td></tr><tr><td>1.5.2</td><td>Systolic model</td><td>26</td></tr><tr><td>1.5.3</td><td>Hypercube model</td><td>27</td></tr><tr><td>1.5.4</td><td>BSP model</td><td>29</td></tr><tr><td>1.5.5</td><td>CGM model and motivation behind the choice of this model</td><td>31</td></tr></table>

<table><tr><td>1.6</td><td>Dynamic programming</td><td>32</td></tr><tr><td>1.6.1</td><td>Recalling the divide-and-conquer technique</td><td>33</td></tr><tr><td>1.6.2</td><td>Building a dynamic-programming solution</td><td>35</td></tr><tr><td>1.6.3</td><td>Principles of dynamic programming</td><td>38</td></tr></table>

<table><tr><td>1.7</td><td>Taxonomy of dynamic-programming formulations</td><td>39</td></tr></table>

<table><tr><td>1.8</td><td>General dynamic-programming formulation of the studied problems</td><td>41</td></tr></table>

<table><tr><td>1.9</td><td>Summary</td><td>43</td></tr></table>

## Chapter 2 • Parallelization of the Studied Problems : State of the Art 45

<table><tr><td>2.1</td><td>Introduction</td><td>45</td></tr></table>

<table><tr><td>2.2</td><td>Minimum cost parenthesizing problem</td><td>45</td></tr></table>

<table><tr><td>2.2.1</td><td>Overview</td><td>45</td></tr></table>

<table><tr><td>2.2.2</td><td>Sequential algorithm of Godbole (1973)</td><td>47</td></tr></table>

<table><tr><td>2.2.3</td><td>Dynamic graph model of Bradford (1994)</td><td>51</td></tr></table>

<table><tr><td>2.2.4</td><td>CGM-based parallel solution of Kechid and Myoupo (2009)</td><td>53</td></tr></table>

<table><tr><td>2.2.5</td><td>CGM-based parallel solution of Kengne and Myoupo (2012)</td><td>60</td></tr></table>

<table><tr><td>2.3</td><td>Optimal binary search tree problem</td><td>62</td></tr></table>

<table><tr><td>2.3.1</td><td>Overview</td><td>62</td></tr></table>

<table><tr><td>2.3.2</td><td>Sequential algorithm of Knuth (1971)</td><td>66</td></tr></table>

<table><tr><td>2.3.3</td><td>CGM-based parallel solution of Kengne et al. (2016)</td><td>69</td></tr></table>

<table><tr><td>2.4</td><td>Triangulation of a convex polygon problem</td><td>76</td></tr></table>

<table><tr><td>2.4.1</td><td>Overview</td><td>76</td></tr></table>

<table><tr><td>2.4.2</td><td>Sequential algorithm of Yao (1982)</td><td>77</td></tr></table>

<table><tr><td>2.4.3</td><td>CGM-based parallel solution of Kechid and Myoupo (2008b)</td><td>81</td></tr></table>

<table><tr><td>2.4.4</td><td>CGM-based parallel solution of Myoupo and Kengne (2014a)</td><td>85</td></tr></table>

<table><tr><td>2.4.5</td><td>Drawbacks of sequential and CGM-based parallel solutions</td><td>89</td></tr></table>
<table><tr><td>2.4.6 Our fast sequential algorithm</td><td>90</td></tr><tr><td>2.4.7 Experimental results</td><td>93</td></tr><tr><td>2.5 Summary</td><td>96</td></tr></table>

## Chapter 3 • Reconciliation of the Minimization of the Number of Communication Rounds and the Load-Balancing of Processors 98

3.1 Introduction 98

3.2 Dynamic graph model of the OBST problem 99

3.3 First dynamic graph partitioning: irregular partitioning technique 103

3.3.1 Blocks' dependency analysis 108

3.3.2 Mapping blocks onto processors 110

3.3.3 CGM-based parallel algorithm for solving the MPP 111

3.3.4 CGM-based parallel algorithm for solving the OBST problem 113

3.3.5 Experimental results 114

3.3.6 Drawback of the irregular partitioning technique 120

3.4 Second dynamic graph partitioning: k-block splitting technique 121

3.4.1 Blocks' dependency analysis of the MPP 123

3.4.2 CGM-based parallel algorithms to solve the MPP 124

3.4.3 Experimental results 127

3.4.4 Drawback of the k-block splitting technique 131

3.5 Third dynamic graph partitioning: four-splitting technique 133

3.5.1 CGM-based parallel algorithm to solve the MPP 136

3.5.2 CGM-based parallel algorithm to solve the OBST problem 137

3.5.3 Experimental results 139

3.6 Summary 149

## General Conclusion 151

Re-stating the research problem 151

Results obtained and critical analysis 152

Further work 155

## Bibliography 157

## Appendix A • List of Publications a
Abstract

We are interested in the parallelization of a class of non-serial polyadic dynamic-
programming problems in this thesis. These problems are characterized by a strong
dependency between subproblems. To design efficient and portable parallel solu-
tions to solve these problems, the CGM (coarse-grained multicomputer) model is
the suitable choice because of its simplicity and its compatibility with most super-
computers.

A CGM-based parallel algorithm is a succession of computation and communica-
tion rounds. The solutions proposed in the literature give the end-user the possi-
bility of minimizing the number of communication rounds or balancing the load
between processors because both objectives are conflicting. Moreover, their main
drawback is to foster the latency time of processors, which accounts for most of
the global communication time.

In this work, we propose an irregular partitioning technique of the dependency
graph to tackle these conflicting objectives. It consists in subdividing the depen-
dency graph into subgraphs (or blocks) of variable size. It ensures that the blocks of
the first steps (or diagonals) are of large sizes to minimize the number of commu-
nication rounds. Thereafter, it decreases these sizes along the diagonals to increase
the number of blocks in these diagonals and allow processors to stay active as long
as possible. These blocks are fairly distributed among processors to minimize their
idle time and balance the load between them. Nevertheless, this strategy induces a
high latency time of processors. Indeed, varying the blocks' sizes does not enable
them to start evaluating some blocks as soon as the data they need are available. To
get over this shortcoming, we propose strategies to evaluate a block as a sequence
of computation and communication steps of a set of small-size blocks. The experi-
mental results obtained show a significant performance gain compared to the most
efficient solutions proposed in the literature.

**Keywords:** Dynamic Programming, Parallel Algorithm, CGM Model, Dependency Graph, Irregular Partitioning
# Solutions parallèles sur le modèle CGM pour une classe de problèmes de programmation dynamique de type polyadique non-serial
# Résumé

Nous nous intéressons à la parallélisation d'une classe de problèmes de programmation dynamique de type polyadique non-sérial dans cette thèse. Ces problèmes se caractérisent par une forte dépendance entre les sous-problèmes. Pour concevoir des solutions parallèles efficaces et portables à la classe de problèmes sus-citée, le modèle de calcul CGM (coarse-grained multicomputer) est le choix idéal en raison de sa simplicité et de sa compatibilité avec la plupart des superordinateurs.

Un algorithme CGM est une succession de rondes de calcul et de communication. Les solutions proposées dans la littérature donnent à l'utilisateur final le choix de minimiser le nombre de rondes de communication ou d'équilibrer la charge des calculs entre les processeurs car ces objectifs sont contradictoires. De plus, leur défaut majeur est de favoriser le temps de latence des processeurs, qui représente la plus grande partie du temps de communication global.

Dans ce travail, nous proposons une technique de partitionnement irrégulier du graphe de dépendances pour apporter une solution à ces objectifs contradictoires. Elle consiste à diviser le graphe de dépendances en sous-graphes (ou blocs) de tailles variables. Elle assure que les blocs des premières étapes (ou diagonales) sont de grandes tailles pour minimiser le nombre de rondes de communication. Ensuite, elle réduit leurs tailles le long des diagonales pour augmenter le nombre de blocs dans ces diagonales et permettre aux processeurs de rester actifs le plus longtemps possible. Ces blocs sont distribués de manière équitable sur les processeurs pour minimiser leur temps d'inactivité et équilibrer leurs charges de calcul. Cependant, cette stratégie induit un temps de latence élevé des processeurs. En effet, la variation de la taille des blocs ne leur permet pas de commencer l'évaluation de certains blocs aussitôt que les données dont ils ont besoin sont disponibles. Pour résoudre ce problème, nous proposons des stratégies consistant à évaluer un bloc comme une succession d'étapes de calcul et de communication d'un ensemble de blocs de petite taille. Les résultats expérimentaux obtenus montrent un gain de performance significatif comparé aux meilleures solutions proposées dans la littérature.

**Mots clés:** Programmation Dynamique, Algorithme Parallèle, Modèle CGM, Graphe de Dépendances, Partitionnement Irrégulier.
List of Acronyms

<table><tr><td>API</td><td>Application Programming Interface;</td></tr><tr><td>BSP</td><td>Bulk Synchronous Parallel;</td></tr><tr><td>CGM</td><td>Coarse-Grained Multicomputer;</td></tr><tr><td>CPU</td><td>Central Processing Unit;</td></tr><tr><td>DAG</td><td>Directed Acyclic Graph;</td></tr><tr><td>DP</td><td>Dynamic Programming;</td></tr><tr><td>DSP</td><td>Digital Signal Processor;</td></tr><tr><td>GPU</td><td>Graphics Processing Unit;</td></tr><tr><td>MCOP</td><td>Matrix Chain Ordering Problem;</td></tr><tr><td>MPI</td><td>Message Passing Interface;</td></tr><tr><td>MPP</td><td>Minimum Cost Parenthesizing Problem;</td></tr><tr><td>OBST</td><td>Optimal Binary Search Tree;</td></tr><tr><td>OpenMP</td><td>Open Multi-Processing;</td></tr><tr><td>PRAM</td><td>Parallel Random Access Memory;</td></tr><tr><td>PU</td><td>Processing Unit;</td></tr><tr><td>TCP</td><td>Triangulation of a Convex Polygon;</td></tr><tr><td>VLSI</td><td>Very-Large-Scale Integration.</td></tr></table>
List of Tables

<table><tr><td>1</td><td>Example of probabilities of three sorted keys</td><td>67</td></tr><tr><td>2</td><td>Total execution time (in seconds) of sequential and CGM-based parallel solutions while solving the TCP problem</td><td>95</td></tr><tr><td>3</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {4096,...,40960}, p ∈ {1,32}, and k ∈ {0,1,2} while solving the MPP with the irregular partitioning technique</td><td>116</td></tr><tr><td>4</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {4096,...,40960}, p ∈ {1,32}, and k ∈ {0,1,2} while solving the OBST problem with the irregular partitioning technique</td><td>116</td></tr><tr><td>5</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {36864,40960}, p ∈ {64,96,128}, and k ∈ {0,1,2} while solving the MPP with the irregular partitioning technique</td><td>117</td></tr><tr><td>6</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {36864,40960}, p ∈ {64,96,128}, and k ∈ {0,1,2} while solving the OBST problem with the irregular partitioning technique</td><td>117</td></tr><tr><td>7</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {4096,...,40960}, p = 32, and k ∈ {1,2} while solving the MPP with the k-block splitting technique</td><td>128</td></tr><tr><td>8</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {36864,40960}, p ∈ {64,96,128}, and k ∈ {1,2} while solving the MPP with the k-block splitting technique</td><td>128</td></tr><tr><td>9</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n = 40960, p ∈ {32,...,128}, and k ∈ {3,4} while solving the MPP with the k-block splitting technique using the k-block by k-block evaluation strategy</td><td>129</td></tr><tr><td>10</td><td>Total execution time (in seconds), speedup, and efficiency (in %) for n ∈ {4096,...,40960}, p = 32, and k ∈ {2,3,4} while solving the MPP with the four-splitting technique on the MatriCS platform</td><td>140</td></tr></table>
11 Total execution time (in seconds), speedup, and efficiency (in %) for $n \in \{36864, 40960\}$, $p \in \{64, 96, 128\}$, and $k \in \{2, 3, 4\}$ while solving the MPP with the four-splitting technique on the MatriCS platform . . . 140

12 Total execution time (in seconds), speedup, and efficiency (in %) for $n \in \{4096, \dots, 40960\}$, $p = 32$, and $k \in \{1, \dots, 5\}$ while solving the OBST problem with the four-splitting technique on the MatriCS platform . . . 141

13 Total execution time (in seconds), speedup, and efficiency (in %) for $n \in \{36864, 40960\}$, $p \in \{64, 96, 128\}$, and $k \in \{1, \dots, 5\}$ while solving the OBST problem with the four-splitting technique on the MatriCS platform 141

14 Total execution time (in seconds), speedup, and efficiency (in %) for $n \in \{4096, \dots, 16384\}$, $p \in \{1, 32\}$, and $k \in \{1, \dots, 5\}$ while solving the OBST problem with the four-splitting technique on our Raspberry Pi cluster141
List of Figures

<table><tr><td>1 Taxonomy of Flynn (1966)</td><td>13</td></tr><tr><td>2 Taxonomy of Raina (1992)</td><td>14</td></tr><tr><td>3 Static and dynamic network topologies</td><td>15</td></tr><tr><td>4 Fork-join paradigm in shared-memory programming model</td><td>20</td></tr><tr><td>5 Network topologies of processors in systolic architectures</td><td>27</td></tr><tr><td>6 Network topologies of processors in hypercube architectures</td><td>28</td></tr><tr><td>7 Description of a BSP superstep</td><td>30</td></tr><tr><td>8 Recursion tree while computing the 4th Catalan number using the divide-and-conquer technique</td><td>35</td></tr><tr><td>9 Recursion tree and the dynamic-programming table while computing the 4th Catalan number using the top-down approach</td><td>37</td></tr><tr><td>10 Dependencies between subproblems while computing the 4th Catalan number using the bottom-up approach</td><td>38</td></tr><tr><td>11 Task graph and dynamic-programming table used to compute Cost[1, 4]</td><td>44</td></tr><tr><td>12 Dynamic-programming and tracking tables filled while computing the product of four matrices with respective dimensions (5 × 10), (10 × 3), (3 × 20), and (20 × 6)</td><td>48</td></tr><tr><td>13 Dynamic graphs D<sub>4</sub> and D'<sub>4</sub> for a problem of size n = 4</td><td>52</td></tr><tr><td>14 Shortest path matrix partitioning strategy proposed by Kechid and Myoupo (2009) for n = 32 and p ∈ {2, 4, 8}</td><td>53</td></tr><tr><td>15 Dependencies of two blocks SM(i, j) and SM(h, l) after applying the partitioning strategy of Kechid and Myoupo (2009)</td><td>55</td></tr><tr><td>16 Distribution schemes of blocks on processors proposed by Kechid and Myoupo (2009)</td><td>56</td></tr><tr><td>17 Alternative bidirectional projection mapping on four and eight processors</td><td>56</td></tr><tr><td>18 Shortest path matrix partitioning strategy proposed by Kengne and Myoupo (2012) for n = 32 and p ∈ {2, 3, 4, 5, 6, 7, 8}</td><td>61</td></tr><tr><td>19 Snake-like mapping on five and eight processors</td><td>61</td></tr><tr><td>20 Example of a binary search tree corresponding to the set of letters (c, e, f, g, h, k, l, n, o, r, s) sorted in alphabetical order</td><td>64</td></tr></table>
<table>
  <tbody>
    <tr>
      <td>21</td>
      <td>Five possible binary search trees obtained from the set of sorted keys &lt;a, b, c&gt;</td>
      <td>64</td>
    </tr>
    <tr>
      <td>22</td>
      <td>Task graph and the dynamic-programming table used to compute Tree[0, 3]</td>
      <td>66</td>
    </tr>
    <tr>
      <td>23</td>
      <td>Dynamic-programming and tracking tables filled while determining the optimal binary search binary tree from probabilities of three sorted keys given in Table 1</td>
      <td>68</td>
    </tr>
    <tr>
      <td>24</td>
      <td>Dependencies and extremities of a block SM(i, j) after applying the partitioning strategy of Kengne et al. (2016)</td>
      <td>72</td>
    </tr>
    <tr>
      <td>25</td>
      <td>Alternative bidirectional projection mapping and snake-like mapping on four processors when g = 8</td>
      <td>73</td>
    </tr>
    <tr>
      <td>26</td>
      <td>A convex polygon and two different triangulations</td>
      <td>76</td>
    </tr>
    <tr>
      <td>27</td>
      <td>Different ways of triangulating the convex polygon P = &lt;5, 10, 3, 20, 6&gt; and the corresponding parenthesis of the product of four matrices M<sub>1</sub>, M<sub>2</sub>, M<sub>3</sub>, and M<sub>4</sub></td>
      <td>77</td>
    </tr>
    <tr>
      <td>28</td>
      <td>A convex polygon P and the corresponding DAG</td>
      <td>79</td>
    </tr>
    <tr>
      <td>29</td>
      <td>Distribution scheme of blocks on three processors proposed by Kechid and Myoupo (2008b)</td>
      <td>84</td>
    </tr>
    <tr>
      <td>30</td>
      <td>Distribution scheme of blocks on four processors proposed by Myoupo and Kengne (2014a)</td>
      <td>89</td>
    </tr>
    <tr>
      <td>31</td>
      <td>Evaluation of the nonleaf node and dependencies between its cones to be computed</td>
      <td>90</td>
    </tr>
    <tr>
      <td>32</td>
      <td>A convex polygon P and the new DAG obtained after building the stack of every node</td>
      <td>91</td>
    </tr>
    <tr>
      <td>33</td>
      <td>DAG in the best and the worst cases of a convex polygon P</td>
      <td>93</td>
    </tr>
    <tr>
      <td>34</td>
      <td>Comparison of the total execution time between our solution versus the sequential solution of Yao (1982) and CGM-based parallel solutions of Kechid and Myoupo (2008b) and Myoupo and Kengne (2014a) while solving the TCP problem</td>
      <td>96</td>
    </tr>
    <tr>
      <td>35</td>
      <td>Comparison of our solution and the best CGM-based parallel solution while solving the TCP problem</td>
      <td>97</td>
    </tr>
    <tr>
      <td>36</td>
      <td>Dynamic graphs D<sub>3</sub> and D'<sub>3</sub> for a problem of size n = 3</td>
      <td>100</td>
    </tr>
    <tr>
      <td>37</td>
      <td>Dynamic graph D'<sub>3</sub> filled while determining the optimal binary search tree from probabilities of three sorted keys given in Table 1</td>
      <td>104</td>
    </tr>
    <tr>
      <td>38</td>
      <td>Fragmentation of a block of size α × β</td>
      <td>104</td>
    </tr>
    <tr>
      <td>39</td>
      <td>Irregular partitioning technique of the shortest path matrix for n = 32, k ∈ {1, 2}, and p ∈ {3, 4, 5, 6, 7, 8}</td>
      <td>106</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>40</td>
      <td>Dependencies of two blocks SM(i, j) and SM(h, l) after applying the ir-regular partitioning technique</td>
      <td>109</td>
    </tr>
    <tr>
      <td>41</td>
      <td>Dependencies and extremities of a block SM(i, j) after applying the ir-regular partitioning technique</td>
      <td>110</td>
    </tr>
    <tr>
      <td>42</td>
      <td>Snake-like mapping on six processors when k = 1</td>
      <td>111</td>
    </tr>
    <tr>
      <td>43</td>
      <td>Global communication time for n ∈ {4096, ..., 40960}, p ∈ {32, ..., 128}, and k ∈ {0, 1, 2} while solving the MPP and the OBST problem with the irregular partitioning technique</td>
      <td>118</td>
    </tr>
    <tr>
      <td>44</td>
      <td>Load imbalance of processors for n ∈ {24576, 32768, 40960}, p = 32, and k ∈ {0, 1, 2} while solving the MPP and the OBST problem with the irregular partitioning technique</td>
      <td>118</td>
    </tr>
    <tr>
      <td>45</td>
      <td>Computation rate versus communication rate for n ∈ {8192, 24576, 40960}, p = 32, and k ∈ {1, 2} while solving the MPP and the OBST problem with the irregular partitioning technique</td>
      <td>119</td>
    </tr>
    <tr>
      <td>46</td>
      <td>Total execution time for n ∈ {4096, ..., 40960}, p ∈ {32, ..., 128}, and k ∈ {0, 1, 2} while solving the MPP and the OBST problem with the irregular partitioning technique</td>
      <td>120</td>
    </tr>
    <tr>
      <td>47</td>
      <td>Drawback of the irregular partitioning technique</td>
      <td>121</td>
    </tr>
    <tr>
      <td>48</td>
      <td>k-block splitting technique of the shortest path matrix for n = 32, k ∈ {1, 2}, and p ∈ {3, 4, 5, 6, 7, 8}</td>
      <td>122</td>
    </tr>
    <tr>
      <td>49</td>
      <td>Dependencies of two k-blocks SM(i, j) and SM(h, l) after applying the k-block splitting technique</td>
      <td>124</td>
    </tr>
    <tr>
      <td>50</td>
      <td>Global communication time for n ∈ {4096, ..., 40960}, p = 32, and k ∈ {0, 1, 2} while solving the MPP with the k-block splitting technique</td>
      <td>129</td>
    </tr>
    <tr>
      <td>51</td>
      <td>Load imbalance of processors for n ∈ {24576, 32768, 40960}, p = 32, and k ∈ {0, 1, 2} while solving the MPP with the k-block splitting technique</td>
      <td>130</td>
    </tr>
    <tr>
      <td>52</td>
      <td>Computation rate versus communication rate for n ∈ {8192, 24576, 40960}, p = 32, and k ∈ {1, 2} while solving the MPP with the k-block splitting technique</td>
      <td>130</td>
    </tr>
    <tr>
      <td>53</td>
      <td>Total execution time for n ∈ {4096, ..., 40960}, p ∈ {32, ..., 128}, and k ∈ {0, 1, 2} while solving the MPP with the k-block splitting technique</td>
      <td>131</td>
    </tr>
    <tr>
      <td>54</td>
      <td>Global communication time and total execution time for n = 40960, p ∈ {32, ..., 128}, and k ∈ {0, ..., 4} while solving the MPP with the k-block splitting technique using the k-block by k-block evaluation strategy</td>
      <td>132</td>
    </tr>
    <tr>
      <td>55</td>
      <td>Four-splitting technique of the shortest path matrix for n = 32, k ∈ {1, 2}, and p ∈ {3, 4, 5, 6, 7, 8}</td>
      <td>134</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>56</td>
      <td>Steps to evaluate the four subblocks of a given block while solving the OBST problem with the four-splitting technique</td>
      <td>139</td>
    </tr>
    <tr>
      <td>57</td>
      <td>Global communication time for n ∈ {4096, ..., 40960}, p ∈ {32, ..., 128}, and k ∈ {0, 1, 2} while solving the MPP and the OBST problem with the four-splitting technique on the MatriCS platform</td>
      <td>142</td>
    </tr>
    <tr>
      <td>58</td>
      <td>Comparison of the overall computation time and the global communication time for n ∈ {24576, 32768, 40960}, p ∈ {32, ..., 128}, and k ∈ {1, ..., 5} while solving the MPP and the OBST problem with the four-splitting technique on the MatriCS platform</td>
      <td>144</td>
    </tr>
    <tr>
      <td>59</td>
      <td>Load imbalance of processors for n ∈ {24576, 32768, 40960}, p = 32, and k ∈ {0, 1, 2} while solving the OBST problem with the four-splitting technique on the MatriCS platform</td>
      <td>146</td>
    </tr>
    <tr>
      <td>60</td>
      <td>Total execution time for n ∈ {4096, ..., 40960}, p ∈ {32, ..., 128}, and k ∈ {0, 1, 2} while solving the MPP and the OBST problem with the four-splitting technique on the MatriCS platform</td>
      <td>146</td>
    </tr>
    <tr>
      <td>61</td>
      <td>Comparison of the overall computation time and the global communication time for n ∈ {8192, 12288, 16384}, p = 32, and k ∈ {1, ..., 5} while solving the OBST problem with the four-splitting technique on the MatriCS platform and our Raspberry Pi cluster</td>
      <td>149</td>
    </tr>
  </tbody>
</table>
List of Algorithms

<table>
  <tbody>
    <tr>
      <td>1</td>
      <td>Divide-and-conquer algorithm to compute the <i>n</i>th Catalan number</td>
      <td>34</td>
    </tr>
    <tr>
      <td>2</td>
      <td>Dynamic-programming algorithm to compute the <i>n</i>th Catalan number using the top-down approach</td>
      <td>36</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Dynamic-programming algorithm to compute the <i>n</i>th Catalan number using the bottom-up approach</td>
      <td>37</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Generic sequential algorithm to solve the studied problems</td>
      <td>43</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Sequential algorithm of Godbole (1973) to solve the MCOP</td>
      <td>48</td>
    </tr>
    <tr>
      <td>6</td>
      <td>CGM-based parallel algorithm of Kechid and Myoupo (2009) to solve the MPP</td>
      <td>58</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Finalization phase using in Algorithm 6 to evaluate blocks SM(i, j)</td>
      <td>59</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Updating phase using in Algorithm 6 to refresh the shortest path costs to nodes of blocks SM(i, j)</td>
      <td>59</td>
    </tr>
    <tr>
      <td>9</td>
      <td>CGM-based parallel algorithm of Kengne and Myoupo (2012) to solve the MPP</td>
      <td>62</td>
    </tr>
    <tr>
      <td>10</td>
      <td>Sequential algorithm of Godbole (1973) to solve the OBST problem</td>
      <td>66</td>
    </tr>
    <tr>
      <td>11</td>
      <td>Sequential algorithm of Knuth (1971) to solve the OBST problem</td>
      <td>67</td>
    </tr>
    <tr>
      <td>12</td>
      <td>CGM-based parallel algorithm of Kengne et al. (2016) to solve the OBST problem</td>
      <td>74</td>
    </tr>
    <tr>
      <td>13</td>
      <td>Recursif algorithm of Yao (1982) to solve the TCP problem</td>
      <td>78</td>
    </tr>
    <tr>
      <td>14</td>
      <td>Dynamic-programming algorithm of Yao (1982) to solve the TCP problem</td>
      <td>80</td>
    </tr>
    <tr>
      <td>15</td>
      <td>CGM-based parallel solution of Kechid and Myoupo (2008b) to solve the TCP problem</td>
      <td>82</td>
    </tr>
    <tr>
      <td>16</td>
      <td>Dependency graph partitioning algorithm of Kechid and Myoupo (2008b)</td>
      <td>83</td>
    </tr>
    <tr>
      <td>17</td>
      <td>Dependency graph partitioning algorithm of Myoupo and Kengne (2014a) balancing the load of processors</td>
      <td>86</td>
    </tr>
    <tr>
      <td>18</td>
      <td>Dependency graph partitioning algorithm of Myoupo and Kengne (2014a) minimizing the number of communication rounds</td>
      <td>88</td>
    </tr>
    <tr>
      <td>19</td>
      <td>Building the stack of every node in the DAG</td>
      <td>91</td>
    </tr>
  </tbody>
</table>
<table>
  <tbody>
    <tr>
      <td>20</td>
      <td>Our dynamic-programming algorithm to solve the TCP problem</td>
      <td>92</td>
    </tr>
    <tr>
      <td>21</td>
      <td>Our CGM-based parallel algorithm based on the irregular partitioning technique to solve the MPP</td>
      <td>112</td>
    </tr>
    <tr>
      <td>22</td>
      <td>Finalization phase using in Algorithm 21 to evaluate blocks SM(i, j) belonging to the lth level of fragmentation</td>
      <td>112</td>
    </tr>
    <tr>
      <td>23</td>
      <td>Updating phase using in Algorithm 21 to refresh the shortest path costs to nodes of blocks SM(i, j) belonging to the lth level of fragmentation</td>
      <td>112</td>
    </tr>
    <tr>
      <td>24</td>
      <td>Our CGM-based parallel algorithm based on the irregular partitioning technique to solve the OBST problem</td>
      <td>113</td>
    </tr>
    <tr>
      <td>25</td>
      <td>Our CGM-based parallel algorithm based on the k-block splitting technique to solve the MPP using the diagonal by diagonal evaluation approach</td>
      <td>125</td>
    </tr>
    <tr>
      <td>26</td>
      <td>Our CGM-based parallel algorithm based on the k-block splitting technique to solve the MPP using the k-block by k-block evaluation approach</td>
      <td>126</td>
    </tr>
    <tr>
      <td>27</td>
      <td>Our CGM-based parallel algorithm based on the four-splitting technique to solve the MPP</td>
      <td>136</td>
    </tr>
    <tr>
      <td>28</td>
      <td>Our CGM-based parallel algorithm based on the four-splitting technique to solve the OBST problem</td>
      <td>138</td>
    </tr>
  </tbody>
</table>
# General Introduction

## CONTENTS

<table><tr><td>Background</td><td>1</td></tr><tr><td>Research problem</td><td>5</td></tr><tr><td>Research aim</td><td>6</td></tr><tr><td>Our contributions</td><td>7</td></tr><tr><td>Thesis outline</td><td>8</td></tr></table>

## Background

### Dynamic programming

Over the last decades, the technological evolution led engineers to confront increasingly complex computational problems in several fields, including artificial intelligence, software engineering, image processing, operations research, economics, and electronics (Collette and Siarry, 2002). Optimization problems, consisting in optimizing (minimizing or maximizing) an objective function from input data, are frequently encountered by engineers. The single-pair shortest path problem is a well-known example. It consists in finding the shortest path between a single pair of vertices in a graph such that the sum of the weights of this path's edges is the smallest. Finding the best solution among a finite and discrete set of feasible solutions is called *combinatorial optimization*. Solving such problems is straightforward because we just need to enumerate all the possible combinations, test them one by one, and choose the best one. Nevertheless, the computational complexity theory shows that this kind of solution is not practical because its time and space complexity is exponential (Cormen et al., 2009). Several techniques have been designed to solve combinatorial optimization problems in polynomial time; one of them is dynamic programming, which is widely studied and used in the literature.
Developed by Bellman Ford in the 1950s, dynamic programming (DP) consists in subdividing a problem into a set of subproblems, organizing their evaluations in such a way that each of them is evaluated just once, and combining their solutions to obtain the optimal solution (Bellman, 1957; Cormen et al., 2009; Lew and Mauch, 2007). The main goal of DP is to avoid recomputing the solutions of subproblems during the solving process of a given problem. Dependencies between subproblems are commonly illustrated by a multi-level directed acyclic graph (DAG), called *task graph* or *dependency graph*. Wah and Li (1988) presented four groups of dynamic-programming problems<fnref n="1" /> based on the nature of the dependency of subproblems and the number of recurrence terms of the objective function: serial monadic, serial polyadic, non-serial monadic, and non-serial polyadic. The latter is characterized by a strong dependency between subproblems, which can lead to an irregular computational load between them. In this thesis, we are interested in a class of non-serial polyadic dynamic-programming problems, in particular the minimum cost parenthesizing problem (MPP), the matrix chain ordering problem (MCOP), the triangulation of a convex polygon (TCP) problem, and the optimal binary search tree (OBST) problem.

These problems can be solved by the sequential algorithm of Godbole (1973). It is sufficient to change the definition of the objective function and the input data according to the semantics of the studied problem. This algorithm runs in $O(n^3)$ time and $O(n^2)$ space by evaluating each of $\Theta(n^2)$ subproblems in $O(n)$ time. However, for the MCOP, the TCP problem, and the OBST problem, there are some speedup sequential algorithms that significantly reduce the execution time and require $O(n^2)$ time and space (Knuth, 1971; Yao, 1982). On the one hand, to solve the OBST problem, the algorithm of Knuth (1971) leverages the table used to store the indices of optimal subproblems to evaluate the $\Theta(n^2)$ subproblems in constant time. On the other hand, the algorithm of Yao (1982) solves the MCOP and the TCP problem<fnref n="2" /> by narrowing down the number of subproblems to $\Theta(n)$ through the quadrangle inequality and evaluating each of them in $O(n)$ time.

**Parallel computing : the master idea**

A few decades ago, general-purpose computers were fitted with single-core processor chips. Nowadays, they are equipped with one or more multi-core processor chips and a large amount of memory. A core is a processing unit capable of per-

<footnote n="1">1. A dynamic-programming problem is a problem that can be solved through the dynamic-programming technique (Kengne, 2014).</footnote>

<footnote n="2">2. There exists a one-to-one correspondence between the MCOP and the TCP problem (Hu and Shing, 1982, 1984). Thus, a solution to one of them remains relevant to the other (Yao, 1982).</footnote>
forming computation operations concurrently (Pacheco and Malensek, 2021). In
spite of this technological leap in hardware, the sequential algorithms of Godbole
(1973), Knuth (1971), and Yao (1982) remain inefficient as the input data size in-
creases. The reason is simple : a sequential algorithm only runs on one processing
unit at a time on a computer, even if that computer has thirty-two, sixty-four, or
one hundred twenty-eight processing units. The performance of such an algorithm
will depend on the power of a processing unit, that is, the number of operations it
is able to perform in one second.

The master idea of parallel computing is to simultaneously use several pro-
cessing units to solve a problem faster by means of a parallel algorithm. Such an
algorithm enables solving a problem in such a way that each processing unit works
on a part of the problem at the same time and cooperates with others to solve the
problem more quickly. It was possible to reduce the execution time of sequential
algorithms that took days, weeks, or months to complete in a few hours, minutes,
or seconds. Parallel computing has emerged at a time when humanity is faced with
increasingly complex problems that need to be solved in a short time, for example
to simulate the reaction of new drugs on the body, to process billions of requests,
to manage millions of banking transactions, and so forth (Sterling et al., 2017).
The need to exploit other computer resources, such as memory to handle problems
with large amounts of data (which could not be handled by a single computer) and
graphics processing unit (GPU) chips to solve problems even faster, has rapidly
become apparent over the years.

A parallel computer is a computer that can execute a parallel algorithm. A
parallel computer architecture consists of one or more parallel computers whose
processing units work together concurrently (Kengne, 2014). An architecture is
said to have a shared memory when the processing units share a common memory
to exchange data, and it is said to have a distributed memory when the processing
units have their own memories and exchange data through a network (Raina, 1992).
These parallel computer architectures are commonly called *supercomputers*.

**Parallel computing model**

A parallel computing model is an abstraction of the architecture's parallel com-
puters intended to be used to formalize in a few parameters the way parallel algo-
rithms will behave in a given architecture (Kengne, 2014). Several models have
been proposed to increase the performance of parallel algorithms. The oldest and
best known of them is the PRAM (Parallel Random Access Memory) model, which
allows the design of parallel algorithms on shared-memory architectures (Czech,
2017; Ferreira and Morvan, 1997; JáJá, 1992; Pacheco and Malensek, 2021). Al-
though this model is simple, it does not reflect the characteristics of real parallel
computers it describes because it assumes that access to information in memory is
achieved in constant time; and thus it leads to designing parallel algorithms that are
efficient in theory but inefficient in practice (Kengne, 2014). Other models, such as
the systolic model and the hypercube model better fit distributed-memory architec-
tures (Czech, 2017; Kung, 1982; Kung and Leiserson, 1978; Leighton, 1992). The
main drawback of these models is that they depend on a specific network topol-
ogy<fnref n="3" />. Thus, parallel algorithms designed in these models are not portable, that is,
the algorithms usually need to be modified when changing architectures to fit the
new architecture (Kengne, 2014).

The bulk synchronous parallel (BSP) model provided a solution to the short-
comings of the previous models (Valiant, 1990). Indeed, it allows better reflecting
the characteristics of real parallel computers by proposing the BSP parallel com-
puter, which is an architecture composed of a set of processors interconnected by a
network and in which two processors can exchange data without going through an
intermediate processor. It is thus easier to design such an architecture being that we
can interconnect general-purpose computers; and thus reduce the cost of design-
ing the architecture compared to systolic and hypercube architectures. Moreover,
through the BSP model, the performance of a parallel algorithm can be estimated
from a few parameters (Kielmann and Gorlatch, 2011; Valiant, 1990).

The coarse-grained multicomputer (CGM) model is a simplified version of the BSP model (Dehne et al., 1993). It allows formalizing in just two parameters the performance of a parallel algorithm : the input data size *n* and the number of processors *p*. A CGM-based parallel algorithm consists in repeating successively two phases until the problem is solved:

* a computation round where processors perform local computations on their data using the best sequential algorithm, and;

* a communication round where processors exchange data through the network.

Because of its simplicity and its compatibility with most recent supercomputers,
this model is suitable for designing parallel algorithms for the class of non-serial
polyadic dynamic-programming problems studied in this thesis.

<footnote n="3">3. A network topology refers to the way in which a set of nodes are connected to each other (Wu and Feng, 1984).</footnote>
Of the problems we study, the MPP is the most classical and generic because it
is widely used to represent a large class of non-serial polyadic DP problems. The
sequential algorithm of Godbole (1973) is besides called *generic sequential algorithm*
in the literature (Kechid and Myoupo, 2009; Kengne, 2014). This algorithm
organizes the subproblems level by level to obtain a DAG having the form of an
upper triangular matrix (see Figure 11a). Subproblems of the same level *i* are evaluated at step (or diagonal) *i*. The parallelization constraint of this problem is the
unevenness of the computational load between different diagonals. In fact, due to
the strong dependency between subproblems, it will have more operations to evaluate the subproblems of diagonal (*i* + 1) than those of diagonal *i*. This constraint
will have heavy consequences on the load balancing.

Moreover, the sequential algorithms of Knuth (1971) and Yao (1982) add new
additional constraints, making the parallelization of the MCOP, the TCP problem,
and the OBST problem a bit more difficult. On the one hand, the sequential al-
gorithm of Knuth (1971) does not change the form of the DAG but neither does it
guarantee that subproblems of the same diagonal will have the same computational
load, compared to the sequential algorithm of Godbole (1973) where subproblems
of the same diagonal have the same computational load (the number of operations
is the same to evaluate these subproblems). On the other hand, according to the
input data, the sequential algorithm of Yao (1982) changes the form of the DAG
into a forest consisting of two binary trees (see Figure 28). Indeed, two input data
of the same size but different values will have two different forms of the DAG; as
a result, designers will be compelled to design a CGM-based parallel solution that
fits all possible forms of the DAGs corresponding to different input data.

The standard methodology for designing CGM-based parallel solutions to solve
these problems is to partition the DAG into subgraphs (or blocks) of same size, then
distribute fairly these blocks among processors, and finally compute them in a suit-
able evaluation order. Some research has investigated the parallelization of these
sequential algorithms on the CGM model in the literature (Cáceres et al., 2010;
Fotso et al., 2010; Higa and Stefanes, 2012; Kechid and Myoupo, 2008a, 2008b,
2009; Kengne and Myoupo, 2012; Kengne et al., 2016; Myoupo and Kengne,
2014a, 2014b). Kengne et al. (2016) highlighted the existence of a relationship
between the total execution time, the load balancing, and the number of commu-
nication rounds from the performance of previous CGM-based parallel solutions.
Indeed, according to the partitioning strategy and the distribution scheme strategy
used when designing these solutions, they showed that minimizing the number of communication rounds and balancing the load of processors are two conflicting objectives when the DAG is partitioned into blocks of the same size :

1. To promote load balancing, the blocks must be small (Kechid and Myoupo, 2008a, 2008b, 2009). Thus, if each processor has one more block to evaluate than another, their load difference will also be low. However, the number of communication rounds will be high.

2. To minimize the number of communication rounds, the blocks must be large (Kengne and Myoupo, 2012; Myoupo and Kengne, 2014a, 2014b). Thus, since there will be few blocks, the number of communication rounds will be reduced. However, the load of processors will be unbalanced.

By generalizing the ideas of the DAG partitioning and distribution scheme introduced in (Kechid and Myoupo, 2008a, 2009; Kengne and Myoupo, 2012; Myoupo and Kengne, 2014b), Kengne et al. (2016) proposed a CGM-based parallel solution that gives the end-user the choice to optimize one criterion according to their own goal, coming for example from the characteristics of parallel computers used.

The main drawback of this solution is the conflicting optimization criteria owing to the fact that the end-user cannot optimize more than one criterion. Moreover, these criteria have a significant impact on the latency time of processors, which in turn has an impact on the global communication time. Recall that the global communication time is obtained by adding up the latency time of processors and the effective transfer time of data. Indeed, when the number of communication rounds is high, excessive communication will lead to communication overhead, which will deteriorate the global communication time. On the other hand, when the load of processors is unbalanced, the evaluation of a block will take a longer time because of its larger size. So, a processor that is waiting for this block to start or continue its computation will wait longer to receive it. As a result, no matter the criterion chosen, the global communication time will decrease the performance of CGM-based parallel solutions.

## Research aim

This thesis aims to design CGM-based parallel solutions that reconcile the conflicting objectives (minimizing the number of communication rounds and balancing the load of processors) to solve a class of non-serial polyadic dynamic-programming
problems. We are confident that solutions where processors evaluate both large and small blocks will reconcile this trade-off. These solutions should minimize the number of communication rounds while balancing the load of processors and reducing the global communication time. These solutions should also minimize the idle time of processors. Indeed, in prior solutions, several processors cannot be active at the same time half in the evaluation of blocks. From this point in time, there are more and more idle processors (one more after each step). However, the closer one gets to the final steps, the higher the load of blocks becomes, and therefore their evaluation requires more computation time.

## Our contributions

We have proposed three contributions to achieve our goal:

* In (Kengne and Lacmou, 2019; Lacmou and Kengne, 2018), we have proposed an irregular partitioning technique to tackle the conflicting objectives. It consists in subdividing the dependency graph into blocks of variable size. It ensures that the blocks of the first diagonals are of large sizes to minimize the number of communication rounds. Thereafter, it decreases these sizes along the diagonals to increase the number of blocks in these diagonals and allow processors to stay active as long as possible. These blocks are fairly distributed among processors to minimize their idle time and balance the load between them. Our CGM-based parallel solutions using this technique minimize also the overall computation time and the latency time of processors; and thus reduce the total execution time. They require $O(n^3/p)$ execution time for the MPP and $O(n^2/\sqrt{p})$ for the OBST problem, each with $O(k\sqrt{p})$ communication rounds. Here, $n$ is the input data size, $p$ is the number of processors, and $k$ is the number of times the size of blocks is subdivided.

* In (Lacmou et al., 2021), we have proposed a fast sequential algorithm for the MCOP and the TCP problem. It consists in organizing the evaluation of the subproblems according to their dependencies, instead of their precedence order as in previous solutions (Myoupo and Kengne, 2014a; Yao, 1982), to solve them fastly by avoiding some unnecessary computations. It requires $O(n)$ time in many cases. Experimental results performed on the MatriCS platform showed that this sequential algorithm is ×18.93 faster than the sequential algorithm of Yao (1982) and ×5.07 faster than the CGM-based parallel solution of Myoupo and Kengne (2014a) on thirty-two processors.
* In (Lacmou et al., 2022a, 2022b, 2022c), we have proposed strategies to reduce the latency time of processors by allowing them to start the evaluation of blocks as soon as possible. Indeed, varying the blocks' sizes does not enable processors to start evaluating small-size blocks as soon as the data they need are available, although these data are available before the end of the evaluation of large-size blocks. We have first proposed the *k*-block splitting technique consisting in splitting the large-size blocks into a set of smaller-size blocks called *k-blocks*. Thus, evaluating a block by a single processor will consist of computing and communicating each *k*-block contained in this block. Experimental results showed that this solution is better than previous ones but leads to communication overhead when *k* increases. That is why we have proposed the four-splitting technique consisting in splitting the large-size blocks into four small-size blocks. It avoids communication overhead and significantly reduces the latency time of processors while preserving the complexity of our solutions using the irregular partitioning technique.

# Thesis outline

The remainder of this document is organized as follows:

## Chapter 1: Parallel Computing and Dynamic Programming

This chapter presents the basic concepts of parallel computing and dynamic programming (DP), and describes the general DP formulation of the studied problems.

## Chapter 2: Parallelization of the Studied Problems: State of the Art

This chapter reviews the state of the art on parallelizing the class of non-serial polyadic DP problems we study, and also presents our fast sequential algorithm which solves the MCOP and the TCP problem with experimental results obtained.

## Chapter 3: Reconciliation of the Minimization of the Number of Communication Rounds and the Load Balancing of Processors

This chapter outlines our CGM-based parallel solutions that reconcile the conflicting objectives to solve the MPP and the OBST problem.

## General Conclusion

This chapter concludes this work and provides some perspectives for future works.
CHAPTER 1

# Parallel Computing and Dynamic Programming

## CONTENTS

<table><tr><td>1.1 - Introduction</td><td>9</td></tr><tr><td>1.2 - High-performance computing</td><td>9</td></tr><tr><td>1.3 - Taxonomies of parallel computer architectures</td><td>11</td></tr><tr><td>1.4 - Designing parallel algorithms</td><td>16</td></tr><tr><td>1.5 - Parallel computing models</td><td>24</td></tr><tr><td>1.6 - Dynamic programming</td><td>32</td></tr><tr><td>1.7 - Taxonomy of dynamic-programming formulations</td><td>39</td></tr><tr><td>1.8 - General dynamic-programming formulation of the studied problems</td><td>41</td></tr><tr><td>1.9 - Summary</td><td>43</td></tr></table>

## 1.1 - Introduction

This chapter outlines the basic concepts of parallel computing and dynamic programming. Section 1.2 introduces the main idea of high-performance computing. Section 1.3 describes the taxonomies of parallel computer architectures. Section 1.4 then explains how to design a parallel algorithm. Parallel computing models are presented in section 1.5 as well as the motivation for choosing the CGM model in this work. Sections 1.6 and 1.7 present dynamic programming and the classification of dynamic-programming formulations, respectively. Finally, Section 1.8 describes the class of dynamic-programming problems studied in this thesis.

## 1.2 - High-performance computing

Toward the end of the 1980s, numerous innovations including the World Wide Web, Google, Microsoft, and Linux led to the growth and creation of many fields
in aeronautical research, space research, economics, medicine, military, and so on. Information and communication technologies have been growing exponentially ever since, owing to the development of diverse materials and computer systems (Ozdamli and Ozdal, 2015). Indeed, they allow accessing, processing, analyzing, storing, or transfering information in different forms to meet a user's need.

In 2020, a study conducted on government-imposed lockdowns of European citizens due to the COVID-19 pandemic showed that the volume of Internet traffic increased by 15-20% almost within a week because the Internet traffic demands of residential users, especially for remote working, entertainment, commerce, and education, increased dramatically (Feldmann et al., 2020). In addition, another study carried out by Li et al. (2019) on the blockchain revealed that fluctuations in the price of Bitcoin have a significant impact on the global computing power of the system. Indeed, a high Bitcoin price will attract more computing power to join the system, while a low Bitcoin price will force computing power to leave the system.

It is therefore virtually impossible to process this massive amount of data in real time with a conventional computer, although nowadays computers are equipped with multi-core processor chips. A multi-core processor is a processor made up of two or more independent processing units, called *cores*, that perform tasks concurrently (Pacheco and Malensek, 2021). Pooling hundreds, thousands, or millions of processing units and making them cooperate to process massive amounts of data and solve complex problems is the master idea of high-performance computing (Sterling et al., 2017). This kind of computer is typically called a *supercomputer*. It has been at the heart of the greatest scientific revolutions and discoveries in the last decades. There are three major types of supercomputers :

* The massive parallel processor, which is a single computer composed of a large number of processors interconnected through a high-speed network. These processors are independent, that is, they do not share memory, and usually each processor runs its own instance of an operating system (Loshin, 2013). Despite the fact that 7.8% of the world's 500 most powerful supercomputers are massive parallel processors in November 2021 according to the TOP500 ranking<fnref n="1" />, Fugaku supercomputer, the first on the list, is a massive parallel processor. It achieved a benchmark score of 442 Pflop/s with 7630848 cores.

* The cluster, which is a set of standalone computers interconnected by a network (Kengne, 2014). According to the TOP500 ranking, 92.2% of super-

<footnote n="1">1. https://www.top500.org</footnote>
computers are clusters. Made up of 4356 computers (for a total of 2414592 cores) interconnected by a dual-rail Mellanox EDR Infiniband, Summit is the second most powerful supercomputer in the world with a performance of 148.8 Pflop/s, which is three times less powerful than the first.

* The quantum computer, which uses properties of quantum physics to perform computations. They aim to solve extremely complex problems the world's most powerful supercomputers cannot and will never solve (Steffen et al., 2011). This technology gained momentum since the first quantum computer was commercialized by D-Wave in 2011, because since then, Google, Intel, IBM, and Microsoft are on a hunt for quantum supremacy<fnref n="2" />.

## 1.3 - Taxonomies of parallel computer architectures

Parallel computing consists in solving a problem by using several processing units simultaneously. Each processing unit works on one part of the problem at the same time and cooperates with others to solve a problem more quickly. A parallel computer is a computer that can perform computation in parallel. A parallel computer architecture consists of one or more parallel computers whose processing units work together concurrently. Nowadays, general-purpose computers are parallel computers since they usually have one or more multi-core processor chips and sophisticated graphics processing unit (GPU) chips. The aspects that distinguish parallel computers have a strong consequence on the design of the parallel algorithms that will be executed on them. A classification according to some criteria, such as the number of instruction streams and data streams, the memory of parallel computers, the network topology, and the granularity, allows highlighting these different aspects.

### 1.3.1 - Classification according to the number of instruction streams and data streams

The first classification relates to the computation activity of the architecture's processing units according to the number of instruction streams and data streams being handled. Recall that the instruction stream is the sequence of instructions executed by a processing unit and the data stream is the sequence of data requested by the

<footnote n="2">². The quantum supremacy is the goal of demonstrating that a quantum computer can solve a problem with such overwhelming speedup that no conventional supercomputer can solve the same problem in a reasonable amount of time (Zhong et al., 2020).</footnote>
instruction stream (Schmidt et al., 2017). Flynn (1966) distinguishes four classes
of architecture :

**SISD (Single Instruction stream, Single Data stream)** where a single-core processor executes a single instruction stream that operates on a single data stream (see Figure 1a). This architecture refers to the traditional von Neumann architecture.

**SIMD (Single Instruction stream, Multiple Data stream)** where many processing units, being part of one or more processors, execute a single instruction stream that operates on different data streams simultaneously (see Figure 1b). This architecture is suitable for image processing and matrix or vector operations, and includes some processor chips like vector processors, array processors, and GPUs.

**MISD (Multiple Instruction stream, Single Data stream)** where many processing units execute different instruction streams that operate on a single data stream (see Figure 1c). This architecture is not very common and is sometimes used to provide fault tolerance with heterogeneous systems running on the same data to provide independent results that are compared with each other (Sitaram and Manjunath, 2012).

**MIMD (Multiple Instruction stream, Multiple Data stream)** where many processing units execute different instruction streams that operates on different data streams (see Figure 1d). Each processing unit has a separate instruction stream and data stream. This architecture is the most common since most modern supercomputers fall into this class.

### 1.3.2 - Classification according to the memory

Based on the address space and physical memory organization of the architecture's processing units, which can be shared or distributed, Raina (1992) proposes three broad classes of architecture:

**SASM (Single Address space, Shared Memory)** which refers to shared-memory architectures. The processing units of a multi-core computer share a common memory through an interconnection network (for example a high-speed memory bus or a crossbar switch, see Figure 2a). A single operating system ensures communication between tasks running on different processing units by writing to and reading from the global memory. There are three
<figure><img src="image_3.png" /><figcaption>**Figure 1** – Taxonomy of Flynn (1966). IS and DS stand for instruction stream and data stream respectively</figcaption></figure>

subcategories of shared memory architectures depending on memory access (Kengne, 2014):

* uniform memory access (UMA): where the shared memory is accessed by all processing units through an interconnection network in the same way that a single processor accesses its memory. Thus, they have equal access time to any memory location.

* non-uniform memory access (NUMA) : each processing unit has a portion of shared memory that has a unique address space. Any processing unit can directly access any memory location using its real address. Therefore, memory access time depends on the distance between the processing units.

* cache-only memory architecture (COMA): which is similar to NUMA except that the shared memory is composed of a cache, so that a data has no proprietary processing unit or specific location in the memory.

**DADM (Distributed Address space, Distributed Memory)** where several computers are interconnected by a network (see Figure 2b). Each computer has a single-core processor chip, its own private memory, and its own operating system. Each processor performs operations on its local memory and if it needs data from other processors, a communication must be established between them to exchange data. This architecture corresponds to distributed-memory architectures.
<figure><img src="image_4.png" /><figcaption>**Figure 2** – *Taxonomy of Raina (1992)*</figcaption></figure>

**SADM (Single Address space, Distributed Memory)** which combines the features of both previous classes (see Figure 2c). This architecture corresponds to distributed shared memory architectures. It is the most used by supercomputers because of its scalability.

### 1.3.3 - Classification according to the network topology

As seen earlier, the parallel computers of an architecture must be interconnected by a fast and reliable network to ensure efficient communication between them. The way they are connected by links (wire or fiber) and switches have a significant impact on the cost, applicability, scalability and performance of the architecture. These networks can be divided into two categories :

**Static networks** where computers are connected by point-to-point links (see Figures 3a, 3b, 3c, and 3d). The topology defined during the design of the architecture does not change. Thus, each computer knows its neighbors as in Figure 3a, where each computer has two neighbors. If a computer sends a data to a non-neighboring computer, it must absolutely transit through all the computers separating them along a given path (Kengne, 2014).

**Dynamic networks** where each computer is connected to a network of switches over one or more links (see Figures 3e, 3f, 3g, and 3h). Thus, the topology can be dynamically updated according to the application requirements running on the architecture. When a data wants to be sent, switches find optimal paths between computers.

### 1.3.4 - Classification according to the granularity

In parallel computation, granularity is the amount of computation related to a task that has been assigned to a processing unit between communication and/or synchronization phases (Grama et al., 2003). It is also defined as the measure of the
<figure><img src="image_5.png" /><figcaption>**Figure 3** – Static and dynamic network topologies. The circle nodes are parallel computers and the square nodes are switches. A link connects either two computers, two switches, or a computer and a switch</figcaption></figure>

ratio of computation to communication (Hwang, 2008). The grain size of a task depends on the type of parallelism that must be performed on the architecture. It has consequences on the performance of a parallel algorithm such as the total number of computation and communication steps, and on the load balancing. Two major types of parallelism exist depending on granularity:

**Fine-grained parallelism** where the computation to communication ratio is low because processing units perform small tasks. It implies high communication overhead as the number of tasks is high. However, this parallelism is suitable for architectures where communications are fast, such as shared-memory architectures, and eases the load balancing (Hwang, 2008). Indeed, the communication overhead can cause a significant drop in performance if the network is not fast.

**Coarse-grained parallelism** where the computation to communication ratio is high because processing units perform large tasks. Since the number of tasks is low, it entails low communication overhead (Hwang, 2008). However, this parallelism makes the load balancing difficult, especially for problems where the amount of computation of tasks is not the same although they may have the same size. In addition, it exacerbates the idleness of processing units because some can perform tasks while others are idle.
## 1.4 - Designing parallel algorithms

Designers of sequential algorithms do not generally care about the characteristics of the computer where these algorithms will be executed because these algorithms exploit only a single processing unit of a multi-core or multi-processor computer. In contrast, a parallel algorithm simultaneously exploits the resources of an architecture's parallel computers (CPU chips, GPU chips, memory, etc.) to solve a problem on a large amount of data more quickly (and thus potentially speed up the production time of a product, save money, etc.). This is why designers of parallel algorithms must know the architecture where these algorithms will be executed. Indeed, the architecture guides designers on the type of parallelization, the source of parallelism, the parallel programming model, and the parallel computing model (described in Section 1.5) to adopt.

### 1.4.1 - Manual parallelization versus automatic parallelization

There are two approaches to parallelizing a sequential algorithm: manual parallelization and automatic parallelization.

#### Manual parallelization

The purpose of manual parallelization is to give designers of parallel algorithms full control over the implementation and execution of the parallel algorithm in a given architecture. More specifically, the designers must analyze the problem to be solved, divide the problem into subproblems (or tasks), analyze the dependencies between tasks, assign the tasks to processing units, define a communication and/or synchronization scheme to allow processing units to exchange data, organize the evaluation of the different tasks until the problem is solved completely, combine the results of different processing units, etc. Of the existing methodologies for manually designing a parallel algorithm, the best known is that of Foster (1995). It is made up of four steps:

1. Partitioning: this process decomposes into tasks the computation that will be performed to solve the whole problem and the data operated on by this computation. The tasks obtained after partitioning are intended to be executed concurrently.

2. Communication: this process determines how tasks will communicate with each other, since the computation to be performed in one task sometimes requires data associated with another task.
3. Agglomeration: this process goes back to the decision made in the two previous steps to group tasks into larger tasks according to the architecture's parallel computers. This step aims to improve performance and reduce the cost of developing the parallel algorithm.

4. Mapping: this process assigns tasks to processing units in an attempt to maximize the use of processing units and minimize communication costs. The tasks are assigned either statically (i.e. determined during the implementation of the algorithm) or dynamically (i.e. determined during the execution of the algorithm) to improve load balancing.

Manual parallelization is usually a time-consuming and complex process. Designers typically encounter many errors and when they correct them, they test their algorithms on large amounts of data (these tests can take days or even weeks) to observe the performance before adopting or redesigning the proposed algorithm. In addition, designers need to have a strong background in parallel computing before designing a parallel algorithm.

## Automatic parallelization

Automatic parallelization is the technique of automatically transforming a sequential algorithm into an equivalent parallel algorithm using the operating system's compiler or specific tools (Mabrouk, 2016). Indeed, it is very common to use iterative structures of nested loops in algorithms; however, these structures consume most of the execution time. It is thus on this part that the compilers ensuring the automatic parallelization will operate by enabling the parallel execution of the loop iterations on processing units. This parallelization consists in making the sequential algorithm undergo a series of legal transformations, in the sense that they preserve the semantics of the initial algorithm, by first carrying out an analysis of dependencies within this sequential algorithm, then detecting the inherent parallelism in the algorithm, and finally determining a series of transformations of this algorithm in respect with its dependencies and making it possible to extract the inherent parallelism (Mabrouk, 2016). There are several loop transformations such as loop distribution, loop interchange, loop skewing, loop fusion, and loop tiling.

Automatic parallelization aims to simplify and reduce the development time of parallel algorithms, which are significantly more difficult to design than sequential algorithms, without necessarily having knowledge of parallel computing. Although it is suitable for shared-memory architectures, it also makes parallel algorithms more portable. However, when faced with a complex sequential algorithm, some semantic information can be lost during the analysis of dependencies and can lead
to the production of an incorrect parallel algorithm or a parallel algorithm that is
less efficient than the initial algorithm (Midkiff, 2012). In some cases, for example
when the loops do not have the same number of iterations, automatic parallelization
still requires the intervention of a designer who has a strong competence in parallel
computing.

1.4.2 - Control parallelism versus data parallelism

Consider a sequential algorithm that executes *m* instructions on each element of an
array of size *n*. To parallelize this algorithm, the designer must spread the tasks
over the resources of the architecture's parallel computers. They have a choice
between two sources of parallelism.

The first is to spread the *m* instructions over processing units. Thus, each pro-
cessing unit will execute a part of instructions on each element of the array (i.e. on
the same data or on different data). This form of parallelism is called *control par-
allelism*. It is widely used in shared-memory architectures. The efficiency of such
a solution depends on the number of instructions to be performed on each element
of the array. The most common control parallelism is called *pipelining*.

The second is to spread the *n* elements of the array over processing units. Each
processing unit will thus execute the *m* instructions on each element of the subset
of the array that has been assigned to it. This form of parallelism is called *data par-
allelism*. It is both used on shared-memory and distributed-memory architectures.
The efficiency of this solution depends on the size of the array (i.e. the input data
size). The advent of GPU chips has expanded the usage of this form of parallelism
by designers of parallel algorithms.

It is obviously possible to combine these two forms of parallelism. In that case,
each processing unit will execute part of the instructions on each element of the
subset of the array that has been assigned to it. This type of solution is used in
distributed shared-memory architectures (see Figure 2c). Its objective is to use the
maximum resources of the architecture's parallel computers to run faster than the
two basic forms.

**1.4.3 - Parallel programming models**

A parallel programming model is a set of program abstractions for fitting parallel
activities from the application to the underlying architecture's parallel computers
(Vitorović et al., 2014). It acts as an interface between the parallel algorithm and
the architecture through its basic operations such as arithmetic operations, spawn-
ing of tasks, reading and writing in shared-memory architectures, or sending and
receiving messages in distributed-memory architectures (Kessler and Keller, 2007). Of the parallel programming models, shared-memory and message-passing models are the best known. The terminology commonly used in literature is as follows: tasks are carried out by threads and processes in the shared-memory programming model and the message-passing programming model, respectively (Czech, 2017; Pacheco and Malensek, 2021; Quinn, 2003).

## Shared-memory programming model

The shared-memory programming model, as the name indicates, is suitable for shared-memory architectures where the processing units share a common memory space (see Figure 2a). In this model, the program corresponding to a parallel algorithm is executed by one or more threads, which are assigned equally to each processing unit by the operating system (Katagiri, 2019b). Indeed, the operating system allocates some resources to the process, including memory space to store instructions and data, the set of registers, etc. All the threads running the process share the address space of a designated area of memory allocated to this process; but they have also a few resources for their exclusive use, such as the stack, program counter register, memory to store private variables (Czech, 2017). So, they implicitly communicate by storing computation results and messages that can be read by other threads asynchronously (Czech, 2017; Vitorović et al., 2014). Mechanisms such as locks, semaphores, and monitors are used to ensure mutual exclusion.

The shared-memory programming model uses the fork-join paradigm, which corresponds to an on-demand parallelism (Czech, 2017; Quinn, 2003). At the beginning of the program execution, only one thread, called *master thread*, executes the sequential portions of the parallel algorithm. In the parts where parallel operations are required (also called *parallel region*), the master thread implicitly spawns or wakes up additional threads; this corresponds to a *fork*. These operations are spread over the master thread and the created threads, and they work simultaneously to accomplish their tasks. At the end of the execution of the parallel region, there is a synchronization barrier to allow the created threads to die or be suspended; this corresponds to a *join*. Thereafter, the master thread resumes the execution that was suspended at the time of the encounter with the parallel region (Czech, 2017; Quinn, 2003). Figure 4 gives an overview.

Nowadays, OpenMP (Open Multi-Processing)<fnref n="3" /> is the mainstream standard for shared-memory programming. It is an application programming interface (API)

<footnote n="3">3. www.openmp.org</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 4** – Fork-join paradigm in shared-memory programming model. α refers to the master thread and β a created thread</figcaption></figure>

that allows writing parallel algorithms in C, C++, and Fortran languages through
compiler directives, library functions, and environment variables. OpenMP does
not provide automatic parallelization (Czech, 2017; Katagiri, 2019b). Indeed, to
design a parallel algorithm with OpenMP, designers usually have to begin with
writing the sequential algorithm before turning it into a parallel algorithm. So
in the first step, designers must identify the parallelism contained in a sequential
algorithm by detecting parallel regions. Then, by applying directives, designers tell
the compiler which sections to run in parallel, and provide information about how
the computation work contained in those sections is spread over the threads. Every
directive begins by the prefix #pragma omp followed by the directive name and the
number of clauses that can be separated by commas or white spaces (Czech, 2017).
The number of threads created by the master thread can be defined by the designer
using the omp_set_num_threads function or the OMP_NUM_THREADS environment
variable.

**Message-passing programming model**

In message-passing programming model, a set of processes (equally distributed on
processing units) executes the same program on different data, and exchanges data
explicitly by sending and receiving messages (Foster, 1995; Hager and Wellein,
2010; Levesque and Wagenbreth, 2010; Quinn, 2003). Each process has its own ad-
dress space. This model targets distributed-memory architectures (see Figure 2b);
nevertheless it can also be used by shared-memory architectures and distributed
shared-memory architectures (see Figure 2c) (Hager and Wellein, 2010). Unlike
the shared-memory programming model, where threads are created when the pro-
gram first enters a parallel region and are blocked while the program executes a
non-parallel region, in this model, processes are created at the beginning of the
execution and are active until the program terminates (Vitorović et al., 2014). This
model gives designers full responsibility for mapping data to processes, coordinat-
ing data communication and synchronization between processes (Nielsen, 2016).
MPI (Message Passing Interface) is the de facto standard for message-passing programming. It is a library that provides functions, beginning by the prefix `MPI_`, for various kinds of communications and for managing processes (Pacheco and Malensek, 2021). MPI is implemented by many free softwares, including MPICH<fnref n="4" /> and OpenMPI<fnref n="5" />, that are supported by most common programming languages such as C, C++, Java, Python, and Fortran.

Two functions are fundamental to designing a parallel algorithm using MPI (Czech, 2017; Foster, 1995). The first function, `MPI_Init`, must be called to perform some setup before any other MPI functions and must be called exactly once per process. The second function is `MPI_Finalize`. It must be called before the end of the parallel algorithm to allow the system to cleanup resources that have been allocated by MPI. No further MPI functions can be called after `MPI_Finalize`.

When a user runs this parallel algorithm, he must specify how many processes should be performed in parallel. Functions `MPI_Comms_size` and `MPI_Comms_rank` respectively determine the number of processes entered by the user and the integer identifier assigned to the current process. Indeed, processes usually identify each other by ranks in the range 0, 1, ..., $p-1$, where $p$ is the number of processes (Pacheco and Malensek, 2021). These ranks are used by the communication functions to send or receive a message. Two modes of communication exist according to the number of processes involved (Czarnul, 2018):

**Point-to-point communication** where two processes are involved in a communication routine. If a process A sends a message to a process B, then the process A must specify the rank of the process B. When the process B needs to receive the message from the process A, the process B can specify the rank of the process A. In case the process B does not know who sent the message, it is possible to allow receiving from any process through the wildcard `MPI_ANY_SOURCE`. Point-to-point communication can be characterized depending on how synchronization is performed (Czarnul, 2018; Katagiri, 2019a):

* The communication is said to be *blocking* when the receiving function `MPI_Recv` does not return to the program after calling the corresponding sending function `MPI Send` until the data are received and the data are entirely copied into a receiving buffer.

<footnote n="4">4. www.mpich.org</footnote>

<footnote n="5">5. www.open-mpi.org</footnote>
* The communication is said to be *non-blocking* when the receiving function `MPI_Irecv` returns to the program immediately even if the corresponding sending function `MPI Send` is not called. There is also the sending function `MPI_Isend` that only starts the sending operation a message without blocking the process. The process may later use on the function `MPI Test` to query the status of the sending operation or the function `MPI_Wait` to wait for its completion (Czech, 2017).

**Collective communication** where several processes are involved in a communication routine. In contrast to the point-to-point communication, all processes involved in a collective call invoke the same function. Some arguments, nevertheless, can be different depending on the rank of a process (Czarnul, 2018). The following are the characteristics of some of the functions:

* the function `MPI Barrier` is used to synchronize processes, i.e. a process is only allowed to pass the barrier once each of the other processes has reached the barrier;

* the function `MPI_Bcast` is used in the case where one process broadcasts data to all other processes;

* the function `MPI_Alltoall` is used in the case where all processes send data to all other processes.

### 1.4.4 - Performance of a parallel algorithm

One of the most important tasks for designers of parallel algorithms is to measure the performance of their algorithms and compare them with those of sequential algorithms that solve the same problem. The goal is to show how the simultaneous use of several resources is more efficient than using only one to solve a problem.

#### Total execution time

The total execution time of a parallel algorithm is usually measured by the difference between the time the first process starts computation and the time the last process finishes the computation (Foster, 1995). This definition is not fully adapted to the shared-memory programming model since the master thread starts and ends the execution of a parallel algorithm by spawning or waking up additional threads if needed. In any case, each process (or thread) computes, communicates, or idles while the algorithm is running. Thus, the execution time is mainly composed of:

1. The computation time, which is the time where a process or a thread performs the computations to carry out the tasks that have been assigned to it.
It depends on several factors such as the problem size, the number of processes, the number of tasks to be performed by a process, the grain size of a task, the characteristics of processing units and memory, and so on. An important challenge for designers is usually to allow processes to perform approximately the same computational load to achieve approximately the same computation time.

2. The communication time, which is the time spent by processes or threads to respectively send and receive messages to each other, or to read and write in the shared memory. In the message-passing programming model, the communication time is often obtained by adding up the effective transfer time of the data and the latency time of processes (Kengne, 2014). Recall that the latency time is the time during which a process waits for a data to continue its computations. Commonly, intra-processor communication, where two depending tasks are located on the same processor chip, is faster than inter-processor communication, where two depending tasks are located on different processor chips. This is because inter-processor communication depends on the physical bandwidth of the communication channel linking the source and destination processors (Foster, 1995).

3. The idle time, which is the time where a process does neither computation nor communication while running the parallel algorithm. It is usually due to a lack of computation in cases where one process performs all the tasks assigned to it (and thus completes its execution) while the others still have tasks to perform. It is also due to the operations of spawning or waking up threads, deleting or suspending threads in the shared-memory programming model. The designers of parallel algorithms seek in most cases to reduce this idle time by allowing processes to remain active, i.e. to perform computation or communication, as long as possible.

Denoting respectively by $T_{comp}^i$, $T_{comm}^i$, and $T_{idle}^i$ the time spent on computation, communication, and idling, on the ith process, Foster (1995) expresses the total execution time $T_{par}$ by Equation (1.1):

$$T_{par} = T_{comp}^i + T_{comm}^i + T_{idle}^i \quad (1.1)$$

## Speedup and efficiency

Speedup is the most commonly used metric for evaluating the performance gain from parallel computing. It shows how many times the execution time of a sequen-
tial algorithm, denoted $T_{seq}$, can be reduced by using a parallel algorithm (Czech, 2017). The speedup is defined as :

$$S = \frac{T_{seq}}{T_{par}} \qquad (1.2)$$

The number of processes (or threads) *p* is the maximum value of the speedup, i.e. a parallel algorithm must be at most *p* times faster than the best sequential algorithm using *p* processes. However, the speedup achieved is typically less than *p*. The reason may be due to the fact that some portions of a sequential algorithm cannot be parallelized, to the unbalanced computational loads of processes, to the overhead associated with the communication time between processes and synchronization of their operations (Czech, 2017; Quinn, 2003).

The efficiency of a parallel algorithm is a metric that measures process utilization. It highlights the efficiency with which a parallel algorithm uses the resources of the architecture's parallel computers (Foster, 1995). It is usually expressed in percentage and defined as :

$$E = \frac{T_{seq}}{p \times T_{par}} = \frac{S}{p} \qquad (1.3)$$

## Scalability

It is important in analyzing the performance of a parallel algorithm to study how the performance of the algorithm varies with parameters such as the problem size and the number of processes (or threads). Let's refer to a parallel algorithm executing on the architecture's parallel computers as a parallel system. The scalability of a parallel system is a measure of its capacity to increase speedup or maintain constant efficiency with a growing number of processes used for computation (Czech, 2017; Foster, 1995; Grama et al., 2003; Quinn, 2003). In general, the efficiency of a parallel system decreases as the number of processes increases. On the other hand, when the problem size is increased while keeping the number of processing units constant, the efficiency of a parallel system increases. A good parallel system is one in which the efficiency can be kept constant as the number of processes and the problem size increase (Grama et al., 2003).

## 1.5 - Parallel computing models

A parallel computing model is an abstraction of an architecture's parallel computers intended to be used by designers of parallel algorithms (Kengne, 2014). It
formalizes in a few parameters the way parallel algorithms will behave in a given
architecture. From the characteristics of a model, it is possible to determine the
time or space complexity of a parallel algorithm without having to implement or
execute it; and thus predict whether or not it will provide good performance, or
compare it with other algorithms solving the same problem to determine which is
the best algorithm for this problem (Kengne, 2014).

A good parallel computing model should meet the following expectations (Ferreira and Morvan, 1997; Foster, 1995; JáJá, 1992):

* a model should be simple and should not have too many parameters to facilitate its usage;

* a model must be well-defined to serve as a common platform for different designers;

* a model must reflect the characteristics of most existing architectures;

* parallel algorithms designed according to a model should not be dependent on a specific architecture;

* theoretical predictions of the performance of parallel algorithms must be in accordance with their execution on the architecture's parallel computers;

* a model must survive the evolution of parallel programming models used by parallel algorithms;

* a model must be deterministic, i.e. when running a parallel algorithm several times with a particular input on a model, it must always lead to the same output.

Unfortunately, it is difficult to design a model that meets all of the above criteria (JáJá, 1992). This section presents the parallel computing models that are most commonly used to design parallel algorithms.

### 1.5.1 - PRAM model

The PRAM (Parallel Random Access Memory) is undoubtedly the most popular shared-memory model in parallel computing (Czech, 2017; Ferreira and Morvan, 1997; JáJá, 1992; Pacheco and Malensek, 2021). This model assumes that a parallel computer is composed of a set of processing units and of a memory to which each processing unit has a constant access time to read and write a data. They operate synchronously under the control of a common clock during computation
steps; hence, access conflicts to the same shared-memory location may occur. Ex-
clusive Read (ER) means that only one processor can read from a memory location
at a given time, and on the other hand, Concurrent Read (CR) means that several
processors can read the contents of the same memory location. A similar classi-
fication has been proposed for write accesses. The writing operation is exclusive
(EW for Exclusive Write) if only one processor can write in a memory cell at a
given time, and concurrent (CW for Concurrent Write) in the case where several
processors can modify the contents of the same memory cell. By combining the
various possibilities of reading and writing, there are four versions of the PRAM
model: ERCW, CRCW, EREW, and CREW (Czech, 2017; JáJá, 1992).

The PRAM model is simple. It allows, in most cases, to know if a problem
can be parallelized or not, to design parallel algorithms without worrying about
communications between processing units, and to easily analyze them (Kengne,
2014). However, this model is too abstract, i.e. it hides the low-level details of
the machines it describes for ease and simplicity, to the point of hiding from de-
signers some that are essential to the design of efficient parallel algorithms. It is
far from real parallel computers because maintaining a constant access time to the
shared memory for a large number of processing units is physically unfeasible and
technically impossible (Czech, 2017). So, a PRAM-based parallel algorithm can-
not be implemented in a real parallel computer and must be readapted to the cho-
sen architecture's parallel computers (Ferreira and Morvan, 1997; Kengne, 2014).
Moreover, a communication operation is much more expensive than a computation
operation and therefore much more influential on the execution time of parallel al-
gorithms in real parallel computers. Thus, hiding this detail from designers leads
to the design of inefficient parallel algorithms (Kengne, 2014).

1.5.2 - Systolic model

Introduced by Kung and Leiserson (1978), the systolic model has proven to be a
powerful tool for the design of specialized embedded processors. A systolic ar-
chitecture is a network of processors that perform simple operations and exchange
data regularly (Kung, 1982). Indeed, systolic-based parallel algorithms operate
synchronously, that is, at each time unit, a processor receives data from some
neighbors, then performs local computations, and finally sends data to some of
its neighbors (JáJá, 1992).

To describe a systolic architecture, it is necessary to specify many details in-
cluding the network topology of processors, the description of processors, the pro-
gram that processors will execute, and the data stream consumed by the network to
<figure><img src="image_7.png" /><figcaption>**Figure 5** – Network topologies of processors in systolic architectures</figcaption></figure>

produce a solution. Figures 5a, 5b, 5c, and 5d illustrate some network topologies of processors in systolic architectures.

The systolic model owes its success to its achievement toward the demand for extremely fast, low-cost supercomputers, and the significant reduction in execution time for more sequential algorithms such as matrix-vector multiplication, matrix-matrix multiplication, matrix inversion, fast Fourier transform, and LU decomposition. This has led to the realization of specialized embedded processors dedicated to many applications like image and signal processing (Kung, 1982). However, this model is too realistic, i.e. it provides designers with the tiniest details (network topologies, routing techniques, etc.) of the architecture's parallel computers it describes in an attempt to increase the efficiency of parallel algorithms, to the point where these details become intrinsic to the validity and efficiency of the parallel algorithms developed. Thus, running these algorithms on computers other than those described by this model will greatly degrade their efficiency because they depend heavily on the chosen network topology; consequently, it is often required to change the parallel algorithm when the network topology evolves (Kengne, 2014).

### 1.5.3 - Hypercube model

A hypercube architecture is composed of $p = 2^d$ processors, indexed from 0 to $p-1$, that are interconnected within a $d$-dimensional cube (Czech, 2017; Ferreira and Morvan, 1997; JáJá, 1992; Leighton, 1992). The latter is a connected graph consisting of $2^d$ nodes and $d \times 2^{d-1}$ edges. Each node corresponds to a $d$-bit binary
<figure><img src="image_8.png" /><figcaption>**Figure 6 – Network topologies of processors in hypercube architectures**</figcaption></figure>

string, and two nodes are linked with an edge if and only if their binary strings differ only at one position (Czech, 2017; Leighton, 1992). Figures 6a, 6b, 6c, and 6d respectively depict a one-dimensional, two-dimensional, three-dimensional, and four-dimensional cube. It can be noticed that $2^d$ processors are needed to extend a $d$-dimensional cube into a $(d+1)$-dimensional cube.

As in the case of the PRAM model and the systolic model, the processors work synchronously. At each step, each processor can send data to one of its neighbors, receive data from one of its neighbors, and perform local computation on that data. A message communicated between two non-neighboring processors is routed through intermediate nodes according to a specific routing algorithm (Mabrouk, 2016). The hypercube model is popular because of its regularity, its versatility, its many interesting graph-theoretic properties, and its ability to handle many com-
putations quickly and simply. However, the logarithmic growth of the number of
connections of each processor with the size of the network is a main drawback
of this model (JáJá, 1992; Leighton, 1992). Moreover, like in the systolic model,
hypercube-based parallel algorithms are not portable (Kengne, 2014).

1.5.4 - BSP model

Valiant (1990) proposed the bulk synchronous parallel (BSP) model as a solution to the efficiency-portability trade-off caused on the one hand by models that are too abstract like the PRAM model, and on the other hand by models that are too realistic like the systolic and hypercube models. Indeed, the idea of this model comes from the observation made on the von Neumann (1945) model, which built a *bridge* between computer designers and sequential algorithm designers. Von Neumann's model allows each of them to focus on their areas of expertise to evolve their different technologies (processors, storage memory, networks, programming languages, compilers, and so on) without worrying about compatibility problems. In the same vein, Valiant (1990) proposed the BSP parallel computer, which is an architecture of parallel computers composed of a set of processors equipped with memory, an interconnection network allowing point-to-point communications between processors, and a synchronization mechanism allowing a barrier synchronization across all processors. The abstraction made on the characteristics of the interconnection network and on the relationship of the memories with processors allows the BSP model to encompass all shared-memory architectures and all distributed-memory architectures; and thus to reflect almost all supercomputers (Kengne, 2014). This is the reason why this model, like the von Neumann model, is widely used till nowadays. It served as a core model for the development of other models such as the LogP model (Culler et al., 1993), the CGM model (Dehne et al., 1993), and the Multi-BSP model (Valiant, 2011); and the development of libraries such as the BSPlib library (Hill et al., 1998), the Green BSP library (Goudreau et al., 1999), the MulticoreBSP library (Yzelman and Roosevelt, 2014).

A BSP-based parallel algorithm consists of a succession of supersteps. As shown in Figure 7, a superstep has three phases (Kielmann and Gorlatch, 2011; Valiant, 1990):

1. A local computation phase where each processor asynchronously performs computation operations by using the best sequential algorithm that solves the problem sequentially, on data in the processor's local memory. These data are put into the local memory either at beginning of the execution of the parallel algorithm or by communication operations of previous supersteps.
<figure><img src="image_9.png" /><figcaption>Figure 7 – Description of a BSP superstep</figcaption></figure>

2 - A global communication phase to exchange data between processors. A communication operation does not become effective before the next super-step begins, i.e. the receiver cannot use received data until the current super-step is finished.

3 - A synchronization phase to finalize communication operations and enable access to received data by the receiving processors before any operation of the next superstep begins.

The performance of a BSP-based parallel algorithm is expressed using three parameters *p*, *l*, and *g*: where *p* is the number of processors, *l* is the time required for a barrier synchronization, and *g* is the time needed for transporting a memory word between two processors (Kielmann and Gorlatch, 2011; McColl, 1995; Valiant, 1990). All communication operations made during a superstep are considered together in what is called an *h*-relation. It is a global communication scheme defined on a set of processors where each one can send or receive at most *h* words (Kengne, 2014). The execution time of a BSP-based parallel algorithm is the sum of the execution times of all supersteps. The execution time of a superstep is obtained by adding up the maximum local computation times of all processors (denoted by *w*), the global communication time, and the synchronization time:

$$T_{\text{superstep}} = w + g \times h + l \quad (1.4)$$
### 1.5.5 - CGM model and motivation behind the choice of this model

Introduced by Dehne et al. (1993) to parallelize sequential algorithms for geometric problems, the coarse-grained multicomputer (CGM) model is a simplified version of the BSP model. It gets rid of the parameters *l*, *g*, *w*, and *h* of the BSP model and retains only the parameters *n* and *p*, which are the input data size and the number of processors, respectively. To represent recent supercomputers composed of several thousand of processors that can process millions or billions of data, this model considers that *p* must be significantly smaller than *n* (*p* ≪ *n*). Unlike the BSP model, this model captures the intrinsic computation and communication features of real parallel computers at two levels:

**At the processors' memory level:** parallel computers of today's architectures are equipped with *state-of-the-art* processors connected to an interconnection network. The CGM model considers that each of the *p* processors is endowed with a sizable amount of local memory so that they can store *O*(*n*/*p*) input data (Olariu, 2008). This is the reason why this model is coarse-grained (Dehne et al., 1993; Ferreira and Morvan, 1997).

**At the communications level:** conventional networks allow the transfer of large amounts of data. The CGM model exploits this feature by allowing processors to send and receive at most *O*(*n*/*p*) data. Indeed, all the information sent from a given processor to another processor during a communication cycle is gathered in a single long message, thereby minimizing the overall message overhead (Dehne et al., 2002). Moreover, during communication operations, synchronization is done implicitly, unlike the BSP model where synchronization is explicitly part of a superstep (Lassous et al., 2000).

The CGM model inherits the portable feature of BSP-based parallel algorithms. Since most personal computers are parallel computers nowadays, the cost of designing a CGM-compatible architecture is relatively low and easy to build. All of these make this model the ideal choice for designing parallel algorithms in this thesis.

A CGM-based parallel algorithm consists of alternating local computation and global communication rounds, a pair of which corresponding to a BSP superstep (Dehne et al., 2002). In each computation round, a processor typically processes its data locally by using the best sequential algorithm to minimize the local computation time per processor. In each communication round, the total data exchanged by each processor are limited to *O*(*n*/*p*). Designers of parallel algorithms will no longer seek to minimize the overall amount of data exchanged, as in the BSP
model, but rather design algorithms with a small number of communication rounds that only depends on *p* and not on *n* (Ferreira, 2001; Lassous et al., 2000). Indeed, for most parallel algorithms based on previous models, the number of communication rounds and the resulting message overhead dramatically increase when the input data size and the number of processors grow, resulting in a considerable performance loss. CGM-based parallel algorithms avoid these problems because they have only a small number of communication rounds whose message size grows with the input data size, which greatly improves performance and scalability (Chan et al., 2008; Dehne, 2006). In a nutshell, an efficient CGM-based parallel algorithm should reduce the number of communication rounds (which should ideally be constant) and the overall computation time.

## 1.6 - Dynamic programming

In mathematics, combinatorics studies the configurations of a set of discrete and finite elements to respond to a given problem (Lovász, 2007). Combinatorial problems occur in many areas of mathematics, including algebra, probability theory, topology, and geometry. A well-known combinatorics problem is to determine the number of possible configurations of a given type (graphs, points, arrays). Another combinatorics problem consists in assigning, from an objective function *f* based on the problem semantics, a numerical value *f*(x) to any configuration x to choose the best configuration that optimizes (maximizes or minimizes) the function *f*. This kind of problem is called *combinatorial optimization problems*. Many algorithm design paradigms have been developed to address these problems.

In the 1950s, greedy algorithms were used to solve some problems in data compression (Huffman, 1952) and graph theory (Dijkstra, 1959; Prim, 1957). This strategy attempts to progressively build the best global solution by choosing at each step the best local solution. Later on, Land and Doig (1960) proposed the Branch and Bound strategy. It consists of recursively splitting the search space into subsets and exploring each of them to find feasible solutions. This strategy eliminates candidate solutions that will not lead to an optimal solution based on the properties of the problem. Tabu search is a heuristic method introduced by Glover (1989, 1990). It consists in determining in a neighborhood the best solution that optimizes the objective function at each iteration. One of its strengths is to maintain a tabu list that prevents visiting the solutions already explored. As for the divide-and-conquer strategy, it recursively decomposes a problem into two or more subproblems until they are sufficiently small to be solved. Then, it solves
these subproblems and gradually builds the original solution from their results. This technique was invented by Gauss (1866) when he proposed an algorithm to solve the discrete Fourier transform for interpolating the trajectories of the asteroids Pallas and Juno (Heideman et al., 1984). Several decades later, many problems have been solved using this technique (such as binary search for search problems or quicksort for sorting problems). Nevertheless, this technique has a drawback when used to solve some problems. Indeed, during the problem-solving process, some subproblems are computed as many times as they are encountered. For this reason, Bellman (1957) proposed a strategy that he called *dynamic programming* to solve this kind of problem.

This section describes the dynamic-programming technique. Section 1.6.1 first recalls the divide-and-conquer technique and presents its main drawback based on a counting problem. Then, Section 1.6.2 shows how the dynamic programming overcomes this drawback by defining a solution. Finally, the principles of dynamic programming are outlined in Section 1.6.3.

### 1.6.1 - Recalling the divide-and-conquer technique

A divide-and-conquer algorithm breaks down a problem into several smaller subproblems that are similar to the original problem, solve them recursively, and then builds the solution of the original problem by combining the solutions of the subproblems (Cormen et al., 2009). There are three main steps in the problem-solving process that are repeated over and over again until the original problem is solved:

**Divide** the problem into subproblems until the base case of the recursive condition is reached. From then on, the problem is no longer divided.

**Conquer** the subproblems by solving them.

**Combine** the solved subproblems together to solve the original problem.

Let us apply this technique to determine the $n$th Catalan number. In combinatorics, the catalan numbers are a sequence of positive integers that appears in many counting problems (Stanley, 2015), such as:

* counting the number of possible binary search trees with $n$ keys;

* counting the number of ways to parenthesize $n$ objects;

* counting the number of ways to connect the 2$n$ points on a circle to form $n$ disjoint chords;
**Algorithm 1** Divide-and-conquer algorithm to compute the *n*th Catalan number

1: **function** CATALAN(*n*)

2: if *n* = 0 or *n* = 1 **then**

3: **return** 1;

4: **else**

5: $c \leftarrow 0$;

6: **for** *k* = 0 to *n* − 1 **do**

7: $c \leftarrow c + \text{CATALAN}(*k*) \times \text{CATALAN}(*n* - *k* - 1);

8: **return** *c*;

* counting the number of ways to triangulate a convex polygon with *n* sides;

* counting the number of non-crossing partitions of a set of *n* elements, or;

* counting the number of semiorders on *n* unlabeled items.

These problems satisfy the recurrence equation (1.5). It is straightforward to deduce Algorithm 1 from this equation.

$$C_n = \begin{cases} 1 & \text{if } n = 0 \text{ or } n = 1, \\ \sum_{k=0}^{n-1} C_k C_{n-k-1} & \text{if } n > 1. \end{cases} \quad (1.5)$$

Figure 8 shows the recursion tree produced by the call `Catalan(4)` to compute $C_4$, the 4th Catalan number:

1. $C_4$ is first divided into two subproblems ($C_0$ and $C_3$);

2. $C_0$ can be conquered easily but not $C_3$;

2.1. $C_3$ is divided into $C_0$ and $C_2$;

2.2. Since $C_2$ can't be conquered;

2.2.1. $C_2$ is divided into $C_0$ and $C_1$;

2.2.2. Since $C_1$ can be conquered, $C_2$ is computed by adding up $C_0C_1$ and $C_1C_0$, the combinations of $C_0$ and $C_1$;

2.3. $C_0$ and $C_2$ is combined to compute $C_0C_2$;

2.4. $C_3$ is reduced to $C_1$, which is conquered and combined to obtain $C_1C_1$;

2.5. $C_3$ is divided into $C_2$ and $C_0$;

2.6. Since $C_2$ can't be conquered;

2.6.1. $C_2$ is divided into $C_0$ and $C_1$;
<figure><img src="image_10.png" /><figcaption>**Figure 8** – Recursion tree while computing the 4th Catalan number using the divide-and-conquer technique</figcaption></figure>

2.6.2 - Since $C_1$ can be conquered, $C_2$ is computed by adding up $C_0C_1$ and $C_1C_0$, the combinations of $C_0$ and $C_1$;

2.7 - $C_2$ and $C_0$ is combined to compute $C_2C_0$;

2.8 - $C_3$ is finally computed by adding up $C_0C_2$, $C_1C_1$, and $C_2C_0$;

3 - $C_0$ and $C_3$ is combined to obtain $C_0C_3$.

It is easy to see that steps 2.2 and 2.6 are similar. Indeed, since the result of $C_2$ is not saved after computing in step 2.2, it must be computed again in step 2.6. In Figure 8, the shaded subtrees represent the subproblems that must be performed more than once when computing $C_4$. $C_3$ is computed two times and $C_2$ is computed six times. This drawback leads to a poor time complexity because the time to compute $C_n$ by this recursive function is exponential in $n$ (Cormen et al., 2009).

## 1.6.2 - Building a dynamic-programming solution

Dynamic programming (DP) is an extension of the divide-and-conquer technique. Unlike a divide-and-conquer algorithm, a DP algorithm stores the results of subproblems when they are first evaluated. Their results are looked up when these subproblems are encountered the next time during the problem solving process. This allows bounding the number of subproblems to be evaluated and the number of operations to be performed to compute a given subproblem. Thus, it is simpler to determine the time and space complexity of a DP algorithm compared to a divide-and-conquer algorithm, which is usually in polynomial time (Cormen et al., 2009).

There are two approaches to obtaining the solution of the original problem from the results of its subproblems: the top-down approach and the bottom-up approach.
**Algorithm 2** Dynamic-programming algorithm to compute the *n*th Catalan number using the top-down approach

1: **function** TOPDOWN-CATALAN(*n*)
2:     let *C*[0..*n*] be a table;
3:     **for** *k* = 0 to *n* **do**
4:         *C*[*k*] ← -∞;
5:     **return** LOOKUP(*C*, *n*);
6: **function** LOOKUP(*C*, *k*)
7:     **if** *C*[*k*] > -∞ **then**
8:         **return** *C*[*k*];
9:     **else**
10:         **if** *k* = 0 or *k* = 1 **then**
11:             *C*[*k*] ← 1;
12:         **else**
13:             *C*[*k*] ← 0;
14:         **for** *i* = 0 to *k* − 1 **do**
15:             *C*[*k*] ← *C*[*k*] + LOOKUP(*C*, *i*) × LOOKUP(*C*, *k* − *i* − 1);
16:     **return** *C*[*k*];

**Top-down approach**

This approach is similar to the divide-and-conquer method in that it decomposes
the original problem until it becomes straightforward to solve. However, it stores
the results of subproblems in a table called *dynamic-programming table*. Algo-
rithm 2 computes the *n*th Catalan number using this approach.

Figure 9a depicts the recursion tree when computing the 4th Catalan number. The shaded subtrees represent the subproblems that are avoided when computing $C_4$ since their results are looked up in the DP table shown in Figure 9b. Indeed, this table is initialized with a value specifying that the results of subproblems are not yet stored (see lines 3-4 in Algorithm 2). Thereafter, the *Lookup* function is called to retrieve the *k*th Catalan number from this table. If the result of $C_k$ is present in the DP table then it is returned (see lines 7-8 in Algorithm 2). Otherwise, $C_k$ is recursively computed by calling *Lookup* on its subproblems, stored in the *k*th entry of the DP table, and returned (see lines 10-16 in Algorithm 2). Compared to Algorithm 1, Algorithm 2 runs in $O(n^2)$ time since it has $\Theta(n)$ entries in the DP table and each of them makes at most $O(n)$ recursive calls.

**Bottom-up approach**

This approach progressively solves the smallest subproblems to build the solution
of the original problem. It is based on the choice of the right order for solving
<figure><img src="image_11.png" /><figcaption>**Figure 9** – **Recursion tree** and the dynamic-programming table while computing the 4th Catalan number using the top-down approach</figcaption></figure>

**Algorithm 3** Dynamic-programming algorithm to compute the *n*th Catalan number using the bottom-up approach

1: **function** BOTTOMUP-CATALAN(*n*)

2: let $C[0..n]$ be a table;

3: $C[0] \leftarrow 1$;

4: $C[1] \leftarrow 1$;

5: **for** $k = 2$ to $n$ **do**

6: $C[k] \leftarrow 0$;

7: **for** $i = 0$ to $k-1$ **do**

8: $C[k] \leftarrow C[k] + C[i] \times C[k-i-1]$;

9: **return** $C[n]$;

the subproblems. This order must guarantee that each subproblem is treated only once, and that it is only processed after all the subproblems on which it depends have been computed. For this purpose, the problem is divided into steps. Each step is composed of several subproblems and has a solving strategy. The solution of each of the subproblems of a given step depends only on those of the subproblems belonging to the previous steps, which are stored in the DP table. The subproblems are solved from the first to the last step. The last step solves usually the original problem. Algorithm 3 illustrates the computation of the 4th Catalan number using the bottom-up approach. Since there are $\Theta(n)$ steps and each of them takes at most $O(n)$ time to be solved, this algorithm runs in $O(n^2)$ time and $O(n)$ space.

Dependencies between the subproblems can be represented by a multi-level, directed and acyclic graph or multi-level DAG (Directed Acyclic Graph) called
<figure><img src="image_12.png" /><figcaption>**Figure 10** – *Dependencies between subproblems while computing the 4th Catalan number using the bottom-up approach*</figcaption></figure>

*task graph* or *dependency graph*. Each node of the DAG represents a subproblem. A level (or a step) *d* is a set of nodes corresponding to the set of subproblems that must be evaluated at step *d*. A node *N*ᵢ is necessary for evaluating the node *N*ⱼ if there is an arc going from *N*ᵢ to *N*ⱼ. A node without outgoing (respectively incoming) arcs corresponds to the original (respectively initial) problem. Figure 10 shows the task graph depicting the dependencies between subproblems while computing the 4th Catalan number. The node *N*ᵢ represents the subproblem *C*ᵢ and a single node is evaluated at each step. Also, the resolution of the subproblem *C*₃ depends to already computed values of *C*₂, *C*₁, and *C*₀.

**1.6.3 - Principles of dynamic programming**

We have applied the dynamic-programming technique to determine the *n*th Catalan number in Section 1.6.2. In this section, two properties required to apply the DP technique to optimization problems are studied : the optimal substructure and the overlapping subproblems.

**Optimal substructure**

The optimal substructure property, also called *principle of optimality* (Lew and Mauch, 2007), is the first criterion that a dynamic-programming problem must satisfy. This principle is stated in (Bellman, 1957, pp. 83) as follows :

> An optimal policy has the property that whatever the initial state and initial decision are, the remaining decisions must constitute an optimal policy with regard to the state resulting from the first decision.

More concretely, this principle states that optimal policies have optimal subpoli-
cies (Lew and Mauch, 2007). This means that each substructure must be optimal
within an optimal structure. Hence, whatever the initial state of a problem, the de-
cisions made in the next states must be optimal and consistent with the decisions
made in the previous states. Therefore, the optimal solution of a problem is made
up of optimal solutions of its subproblems. A DP formulation seeks to define a problem in this way (Grama et al., 2003).

Some combinatorial problems may have the optimal substructure property, but may use too much memory or time to be efficient. Since the DP technique stores the intermediate results of a problem in a DP table for reuse, computing a solution from them can be time intensive when there are too many intermediate results. For example, a set of *n* elements has 2ⁿ subsets and there are *n!* possible permutations of characters of a string of length *n*. Therefore, the DP technique cannot be efficient to find the best solution for each subpermutation/subset since the DP table would be gigantic. However, it can be efficient only when there are not too many results to compute. For example, there are only *n*(*n* − 1)/2 substrings of a string of length *n* and only *n*(*n* − 1)/2 possible subtrees of a binary search tree with *n* keys. In a nutshell, dynamic programming typically produces good results on objects that can be ordered linearly and that cannot be reordered (Kengne, 2014).

## Overlapping subproblems

The property of overlapping subproblems is the second criterion that a DP problem must have. This property means that subproblems have to share common subproblem; that is, a subsubproblem can belong to several subproblems. In this kind of problem, after evaluating a subsubproblem *C* generated by a subproblem *A*, it still needs to be evaluated if it is also generated by a subproblem *B*. Indeed, *C* will be evaluated as many times as it has to be generated by a subproblem. Dynamic programming leverages this property to store in a table the results of subproblems that have already been computed. If this subproblem is encountered once again, then it is looked up in the table instead of being recomputed. The fewer subsubproblems there are in common, the faster the DP algorithm will be. Nevertheless, if subproblems generate new subproblems at each step, then it will be more suitable to drop the DP technique (Cormen et al., 2009).

## 1.7 - Taxonomy of dynamic-programming formulations

The objective function, which characterizes optimization problems, is formulated by a recursive equation called *functional equation* or *optimization equation* (Grama et al., 2003; Kengne, 2014; Lew and Mauch, 2007). The left-hand side represents the unknown value of the optimal solution of a subproblem, and the right-hand side contains the minimization and/or maximization operations on a set of elements (each written in terms of solutions of subproblems belonging to the previous
steps). Defining this equation depends on the semantics of the problem to be addressed. Examples 1 and 2 denote respectively the DP formulations of the all-pairs shortest path problem by Floyd (1962) and the string-to-string correction problem by Wagner and Fischer (1974).

**Example 1** Consider a weighted graph G, which consists of a set of nodes N and a set of edges E. An edge from node i to node j in E has a weight $c_{i,j}$. The cost of a path is the sum of the edges' weights in this path. The algorithm of Floyd (1962) determines the cost $s_{i,j}^k$ of the shortest path between each pair of nodes (i, j) in all intermediate nodes {1, 2, ..., k} through Equation (1.6) defined as follows:

$$s_{i,j}^k = \begin{cases} c_{i,j} & \text{if } k = 0, \\ \min \left\{ s_{i,j}^{k-1}, (s_{i,k}^{k-1} + s_{k,j}^{k-1}) \right\} & \text{if } 1 \le k \le |N| - 1. \end{cases} \quad (1.6)$$

**Example 2** Given two strings a = a₁a₂...aₘ and b = b₁b₂...bₙ on an alphabet Σ, the string-to-string correction problem aims to find the distance between a and b as measured by the minimum number of editing operations (insertions, deletions, or substitutions) required to transform a into b. The edit distance between a and b is given by $d_{m,n}$, and defined by Equation (1.7):

$$d_{i,j} = \begin{cases} d_{i-1,j-1} & \text{if } a_i = b_j, \\ \min \left\{ \begin{array}{l} d_{i-1,j} + \text{deletion}(a_i), \\ d_{i,j-1} + \text{insertion}(b_j), \\ d_{i-1,j-1} + \text{substitution}(a_i, b_j) \end{array} \right. & \text{if } a_i \neq b_j. \end{cases} \quad (1.7)$$

A DP formulation is said to be *monadic* if the functional equation has only one recurrence term, and it is said to be *polyadic* otherwise (Wah and Li, 1988). The DP formulation in Equation (1.7) is monadic since the functional equation involves one recursive term only. Indeed, the evaluation of the subproblem $d_{i,j}$ depends on only one recursive term (subproblem) at a time: either $d_{i-1,j}$, or $d_{i,j-1}$, or $d_{i-1,j-1}$. In contrast, The DP formulation in Equation (1.6) is polyadic since the evaluation of the subproblem $s_{i,j}^k$ involves two subproblems: $s_{i,k}^{k-1}$ and $s_{k,j}^{k-1}$.

A DP formulation can be also classified according to the nature of the dependency of subproblems. If the solution of a subproblem of a given step depends exclusively on the solutions of subproblems of the immediately preceding step, then the DP formulation is said to be *serial*; otherwise, it is said to be *nonserial* (Wah and Li, 1988). The DP formulation in Equation (1.6) is serial since the evaluation of the subproblem $s_{i,j}^k$ belonging to step k required only the solutions of
subproblems $s_{i,j}^{k-1}$, $s_{i,k}^{k-1}$, and $s_{k,j}^{k-1}$, belonging the previous step $k-1$. In contrast, the DP formulation in Equation (1.7) is nonserial since the evaluation of the subproblem $d_{i,j}$ belonging to step $(j-i-1)$ required the solution of the subproblem $d_{i-1,j-1}$ belonging to step $(j-i-3)$, and the solutions of subproblems $d_{i-1,j}$ and $d_{i,j-1}$ belonging step $(j-i-2)$.

Wah and Li (1988) classified DP formulations into four categories based on the above classification criteria: *serial monadic*, *serial polyadic*, *non-serial monadic*, and *non-serial polyadic*. So, the DP formulation of Floyd (1962) for the all-pairs shortest path problem is part of serial polyadic DP formulations, and the DP formulation of Wagner and Fischer (1974) for the string-to-string correction problem is part of non-serial monadic DP formulations. In this thesis, we are interested in non-serial polyadic dynamic-programming problems, especially in a class of problems whose general DP formulation is described in Section 1.8.

## 1.8 - General dynamic-programming formulation of the studied problems

**Definition 1** A finite semigroupoid $(S, R, \bullet)$ is a nonempty finite set $S$ of elements, a binary relation $R \subseteq S \times S$, and an associative binary operator $\bullet$ satisfying the following conditions (Marvins, 1978):

$$1 - \text{if } (a, b) \in R, \text{ then } a \bullet b \in S;$$

$$2 - (a \bullet b, c) \in R \text{ iff } (a, b \bullet c) \in R \text{ and } (a \bullet b) \bullet c = a \bullet (b \bullet c);$$

$$3 - \text{if } (a, b) \in R \text{ and } (b, c) \in R, \text{ then } (a \bullet b, c) \in R.$$

**Definition 2** An associative product is any product of the form $a_1 \bullet a_2 \bullet \cdots \bullet a_n$ such that $(a_i, a_{i+1}) \in R$, $1 \le i \le n$. A linear product is a product of the form $((\cdots (a_1 \bullet a_2) \bullet \cdots) \bullet a_n)$ or $(a_1 \bullet (\cdots \bullet (a_{n-1} \bullet a_n)) \cdots)$.

**Definition 3** A weighted semigroupoid $(S, R, \bullet, pc)$ is a semigroupoid $(S, R, \bullet)$ with a non-negative product cost function $pc$. If $(a_i, a_j) \in R$ then $pc(a_i, a_j)$ is the cost of evaluating $a_i \bullet a_j$. The optimal cost of evaluating an association product $a_i \bullet a_{i+1} \bullet \cdots \bullet a_j$ is denoted by $Cost[i, j]$ and defined by:

$$Cost[i, j] = \begin{cases} \text{Init}(i) & \text{if } 1 \le i = j \le n, \\ \text{Opt}_{i \le k < j} \{Cost[i, k] + Cost[k+1, j] + F(i, k, j)\} & \text{if } 1 \le i < j \le n, \end{cases} \quad (1.8)$$
where n is the problem size and $F(i, k, j) = pc(a_i \bullet \cdots \bullet a_k, a_{k+1} \bullet \cdots \bullet a_j)$.

Equation (1.8) shows that the problems studied in this thesis feature the properties of a DP problem:

* Optimal substructure property: the optimal costs of the subproblems $(k, l)$ are required for computing the optimal cost of a subproblem $(i, j)$, where $i \le k < l \le j$.

* Overlapping subproblems property: consider two subproblems $(i, j)$ and $(i+1, j)$. The evaluation of these subproblems depends on the common values of the evaluation of subproblems $(k, j)$, where $i+1 < k \le j$. It is the same if two subproblems $(i, j)$ and $(i, j+1)$ are considered.

In Equation (1.8), $Cost[i, j]$ corresponds to the value of the optimal solution of the subproblem $(i, j)$. This value is obtained according to the optimum function $Opt$ (for a maximization problem, $Opt = \max$; otherwise $Opt = \min$) among the $(j-i)$ possible combinations of the subsubproblems on which the subproblem $(i, j)$ depends. The value of a combination of two subproblems $(i, k)$ and $(k+1, j)$, where $i \le k < j$, is computed by adding the values of their optimal solutions and the cost corresponding to the combination of these subproblems given by the function $F(i, k, j)$ called union function. The basic subproblems are initialized by the function $Init(i)$.

Several problems can be formulated from Equation (1.8) such as the context-free grammar parsing problem (Kasami, 1965; Younger, 1967) and the Nussinov RNA folding problem (Nussinov and Jacobson, 1980). In this thesis, we are interested in the minimum cost parenthesizing problem (MPP), the matrix chain ordering problem (MCOP), the triangulation of a convex polygon (TCP) problem, and the optimal binary search tree (OBST) problem. The common feature of these problems is the minimization of the optimization equation (this means that $Opt = \min$ in Equation (1.8)). However, they differ in the definition of the union function $F$ and initialization function $Init$ because the definition of these functions depends on the problem to solve.

A solution based on the exhaustive search of all possible combinations would be poor because, for a problem of size $n$, the number of combinations is exponential in $n$ (see Section 1.6.1). To solve the MCOP, Godbole (1973) proposed the first polynomial-time sequential algorithm running in $O(n^3)$ time and $O(n^2)$ space. It became the standard algorithm for solving all problems that can be formulated by Equation (1.8) because the structure and the complexity of this algorithm are
<caption>**Algorithm 4** Generic sequential algorithm to solve the studied problems</caption>

1: **for** *i* = 1 to *n* **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;*Cost*[*i*, *i*] ← *Init*(*i*);
3: **for** *d* = 2 to *n* **do**
4: &nbsp;&nbsp;&nbsp;&nbsp;**for** *i* = 1 to *n* − *d* + 1 **do**
5: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*j* ← *n* − *d* + 1;
6: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*Cost*[*i*, *j*] ← ∞;
7: **for** *k* = *i* to *j* − 1 **do**
8: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*c* ← *Cost*[*i*, *k*] + *Cost*[*k* + 1, *j*] + *F*(*i*, *k*, *j*);
9: **if** *c* < *Cost*[*i*, *j*] **then**
10: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*Cost*[*i*, *j*] ← *c*;
11: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;*Track*[*i*, *j*] ← *k*;

independent of the functions *F* and *Init*. This algorithm is often called *generic sequential algorithm* in the literature (Kengne, 2014). Algorithm 4 draws the big picture.

Computing *Cost*[1,*n*] involves the solution of all subproblems (*i*, *j*), such that 1 ≤ *i* ≤ *j* ≤ *n*. Dependencies between these subproblems can be organized like a DAG, as shown in Figure 11a for a problem of size *n* = 4. This figure reveals, for example, that the value of *Cost*[1,4] is computed from the optimal solutions of the pairs of subproblems (((1,1),(2,4))), (((1,2),(3,4))), and (((1,3),(4,4))). Indeed, Algorithm 4 uses a bottom-up approach to compute the optimal solutions of sub-problems. The value of the optimal solution of each solved subproblem is stored in the DP table depicted in Figure 11b. The index *k* that led to this value among the (*j* − *i*) possible solutions for a subproblem (*i*, *j*) is also saved in a table called *tracking table*. This table, named *Track* in Algorithm 4, is similar to the DP table and is used for the construction of the optimal solution.

## 1.9 - Summary

This chapter was dedicated to the definition of the basic concepts of parallel computing and dynamic programming. In the first part, we started by giving the principle of high-performance computing, which is the *raison d'être* of parallel computing. Then, after presenting the classifications of parallel computer architectures according to several criteria, such as the number of instruction streams and data streams, the memory of parallel computers, the network topology, and the granularity, we studied different techniques to parallelize a sequential algorithm and showed how to measure the performance of a parallel algorithm. Finally, we studied
<figure><img src="image_13.png" /><figcaption>**Figure 11** – *Task graph and dynamic-programming table used to compute Cost* [1, 4]</figcaption></figure>

ied parallel computing models. Different problems of too abstract models like the PRAM model and too realistic models like the systolic model and the hypercube model have been highlighted. The BSP model solved the problems of previous models and the CGM model simplified the BSP model.

In the second part, we first recalled the divide-and-conquer technique and showed its main drawback before presenting the top-down approach and the bottom-up approach of dynamic programming, which overcome this drawback. We then gave the principles of dynamic programming. After classifying the DP formulations, we presented the formulation of the class of non-serial polyadic dynamic-programming problems that interest us.

Chapter 2 will study in detail the parallelization of the sequential algorithms that solve each of these problems.
CHAPTER 2

# Parallelization of the Studied Problems
: State of the Art

## CONTENTS

<table><tr><td>2.1 - Introduction</td><td>45</td></tr><tr><td>2.2 - Minimum cost parenthesizing problem</td><td>45</td></tr><tr><td>2.3 - Optimal binary search tree problem</td><td>62</td></tr><tr><td>2.4 - Triangulation of a convex polygon problem</td><td>76</td></tr><tr><td>2.5 - Summary</td><td>96</td></tr></table>

## 2.1 - Introduction

This chapter reviews the state of the art on parallelization of sequential algorithms for the minimum cost parenthesizing problem, the matrix chain ordering problem, the triangulation of a convex polygon problem, and the optimal binary search tree problem. For each of them, the following points will be presented and studied successively: the description of the problem and its dynamic-programming (DP) formulation, the best sequential algorithm that solves this problem, the parallelization constraints and the literature review on the parallelization of this sequential algorithm on parallel computing models, and finally, the best CGM-based parallel solutions. For the special case of the matrix chain ordering problem, we will present our fast sequential algorithm and experimental results obtained.

## 2.2 - Minimum cost parenthesizing problem

### 2.2.1 - Overview

A parenthesizing process is any process that inserts opening and closing parentheses between a chain of symbols (characters, numbers, matrices, or objects) in
order to define the best sequential processing order according to an intended goal (Kengne, 2014). The minimum cost parenthesizing problem (MPP) of a chain of symbols consists in finding the parenthesizing that will minimize the cost of the computations involved on this chain. Depending on the kind of entities in the chain to parenthesize and the treatment to perform (see the general formulation in Section 1.8), this problem appears in the literature under several variants. A well-known variant is the parenthesizing problem of the lexical ordering of the computation blocks for embedded DSP (Digital Signal Processor) applications (Bhattacharyya and Murthy, 1995). The goal is to find the parenthesizing that minimizes the memory requirement on the DSP. The most popular variant is the matrix chain ordering problem (MCOP), which consists in finding the parenthesizing that minimizes the cost of the product of a chain of matrices (Cormen et al., 2009; Kengne, 2014).

First, consider the cost of multiplying two matrices A and B with respective dimensions ($a_0 \times a_1$) and ($b_0 \times b_1$). They can be multiplied only if $a_1 = b_0$. This operation requires $a_0 \times c \times b_1$ scalar multiplications, where $c = a_1 = b_0$, and produces in output a matrix C with dimensions ($a_0 \times b_1$). The way in which a chain of matrices is bracketed can have a significant impact on the evaluation of the product cost. To illustrate this, consider four matrices $M_1, M_2, M_3$, and $M_4$ with respective dimensions ($5 \times 10$), ($10 \times 3$), ($3 \times 20$), and ($20 \times 6$). Evaluating their product can be done in different ways:

$$1 - (M_1 \times (M_2 \times M_3)) \times M_4, \text{ which requires } 10 \times 3 \times 20 + 5 \times 10 \times 20 + 5 \times 20 \times 6 = 2200 \text{ scalar multiplications;}$$

$$2 - M_1 \times (M_2 \times (M_3 \times M_4)), \text{ which requires } 3 \times 20 \times 6 + 10 \times 3 \times 6 + 5 \times 10 \times 6 = 2460 \text{ scalar multiplications;}$$

$$3 - (M_1 \times M_2) \times (M_3 \times M_4), \text{ which requires } 5 \times 10 \times 3 + 3 \times 20 \times 6 + 5 \times 3 \times 6 = 600 \text{ scalar multiplications.}$$

Thus, the product cost depends on the chosen order although all bracketings yield the same output matrix M with dimensions ($5 \times 6$)$^1$. Therefore, the MCOP consists in finding the optimal order having the lowest cost for multiplying matrices. The input is a list of ($n+1$) natural numbers $d_0, d_1, \dots, d_n$, representing the dimensions of the matrices (a matrix $M_i$ has dimensions ($d_{i-1} \times d_i$)). The output is the minimum cost to compute the product of the $n$ matrices.

Note that the MCOP does not multiply the matrices. It simply looks for an order of multiplication that minimizes the cost. The time spent to determine this optimal

<footnote n="1">1. The matrix multiplication is associative.</footnote>
order is commonly offset by the time saved when performing matrix multiplications using this optimal order; for example, performing 600 scalar multiplications will theoretically take four less time than 2460 (Cormen et al., 2009). Traditionally, the matrix multiplication was introduced in linear algebra to facilitate and clarify computations (Marvins, 1978; Schreier and Sperner, 2011). This strong relationship between matrix multiplication and linear algebra remains fundamental in mathematics, as well as in physics, engineering, and computer science, more precisely in scientific computing (Lee et al., 2003; Lin, 1994; Yau and Lu, 1993).

Back to the MCOP, the minimum cost for evaluating $M_{i,j} = M_i \times M_{i+1} \times \cdots \times M_j$ is denoted by $Cost[i, j]$ and defined by:

$$ Cost[i, j] = \begin{cases} 0 & \text{if } 1 \le i = j \le n, \\ \min_{i \le k < j} \{Cost[i, k] + Cost[k+1, j] + d_{i-1} \times d_k \times d_j\} & \text{if } 1 \le i < j \le n. \end{cases} \tag{2.1} $$

Solving $M = M_1 \times M_2 \times \cdots \times M_i \times \cdots \times M_n$ reduces to computing $Cost[1, n]$. Equation (2.1) is equivalent to Equation (1.8) as $Opt = \min, Init(i) = 0$, and $F(i, k, j) = d_{i-1} \times d_k \times d_j$, where $1 \le i \le k \le j \le n$.

### 2.2.2 - Sequential algorithm of Godbole (1973)

Algorithm 5 gives an overview of the sequential algorithm of Godbole (1973) to solve the MCOP. Since the optimal product cost of $M_{i,j}$ never depends on a product of length greater than $(j-i)$, the subproducts can be organized and computed in $n$ levels (steps or diagonals) according to their length. Thus, products of length $d$, which correspond to $M_{i,i+d-1}$ such that $1 \le i \le i+d-1$, are evaluated at diagonal $d$ in Algorithm 5. The solution of Godbole (1973) aims at filling, diagonal by diagonal, in cells of the DP table (named Cost in Algorithm 5) and the tracking table (named Track in Algorithm 5) respectively the value of the optimal product cost of $M_{i,j}$ (see line 10 in Algorithm 5) and the value of the index $k$ corresponding to the pair of subproducts $(M_{i,k}, M_{k+1,j})$, which minimizes $Cost[i, j]$ (see line 11 in Algorithm 5). The first diagonal of the DP table is initialized to 0 since there is no product of length 1 (see lines 1-2 in Algorithm 5).

Figures 12a and 12b illustrate, respectively, the DP table and the tracking table that have been filled while computing the product of four matrices with respective dimensions $(5 \times 10)$, $(10 \times 3)$, $(3 \times 20)$, and $(20 \times 6)$. Each table is a square matrix of order $n=4$. However, only the part of the table strictly above the main diagonal (the first diagonal) is used. It is obvious that one can use a single table in which the part of the table above the main diagonal will be used to store the values of
<caption>**Algorithm 5** Sequential algorithm of Godbole (1973) to solve the MCOP</caption>

1: **for** i = 1 to *n* **do**
2:     *Cost*[i, i] ← 0;
3:     **for** d = 2 to *n* **do**
4:         **for** i = 1 to *n* - *d* + 1 **do**
5:             *j* ← *n* - *d* + 1;
6:             *Cost*[i, *j*] ← ∞;
7:             **for** *k* = *i* to *j* − 1 **do**
8:                 *c* ← *Cost*[i, *k*] + *Cost*[*k* + 1, *j*] + *d*ᵢ₋₁ × *d*ₖ × *d*ⱼ;
9:                 **if** *c* < *Cost*[i, *j*] **then**
10:                 *Cost*[i, *j*] ← *c*;
11:             **Track*[i, *j*] ← *k*;

<figure><img src="image_14.png" /><figcaption>**Figure 12** – Dynamic-programming and tracking tables filled while computing the product of four matrices with respective dimensions (5 × 10), (10 × 3), (3 × 20), and (20 × 6)</figcaption></figure>

*Cost*[i, j], and the part of the table below the main diagonal will be used to store
the values of *Track*[j, i]. The shaded entries in Figure 12a depict the computation
of *Cost*[1, 4]. The pairs that have the same shading are taken together in line 8 of
Algorithm 5 when computing

$$
\begin{align*}
\text{Cost}[1,4] &= \min \left\{ \begin{array}{l} \text{Cost}[1,1] + \text{Cost}[2,4] + d_0 \times d_1 \times d_4 = 0 + 540 + 300 = 840, \\ \text{Cost}[1,2] + \text{Cost}[3,4] + d_0 \times d_2 \times d_4 = 150 + 360 + 90 = 600, \\ \text{Cost}[1,3] + \text{Cost}[4,4] + d_0 \times d_3 \times d_4 = 450 + 0 + 600 = 1050 \end{array} \right. \\
&= 600.
\end{align*}
$$
To solve the MCOP, the sequential algorithm of Godbole (1973) requires $O(n^3)$ time. This is deducible since the loops are nested at three deep, and each of them takes at most $(n-1)$ values. It requires $O(n^2)$ memory space to store the DP table and the tracking table.

## Parallelization constraint of the sequential algorithm of Godbole (1973)

The main parallelization constraint of this algorithm comes from the irregular computational load for the evaluation of subproducts of a diagonal $d$ since the evaluation of subproducts of length $d$ is smaller than subproducts of length $(d+1)$ (it has more operations to evaluate the subproblems of length $(d+1)$ than those of length $d$). Thus, in the parallelization process, this irregular computational load between the diagonals of the sequential algorithm of Godbole (1973) must be considered. It can cause the load imbalance between processors and lead to poor efficiency of the parallel algorithm. Some research has investigated the parallelization of this sequential solution on different parallel computing models in the literature.

## Literature review on the parallelization of the sequential algorithm of Godbole (1973)

On the PRAM model, Valiant et al. (1983) proposed a general method for parallelizing non-serial polyadic dynamic-programming problems. It requires $Θ(lg^2 n)$ computation time on $O(n^9)$ processors. Rytter (1988) reduced the number of processors down to $O(n^6/lnn)$ by using the specific features of this class of problems. Later on, Huang et al. (1994) improved this solution by reducing the number of processors to $O(n^6/ log^5 n)$. On CREW-PRAM machines, Bradford (1994) and Czumaj (1993) designed solutions running in $O(log^3 n)$ time using respectively $O(n^3/ log n)$ and $O(n^2/ log^3 n)$ processors. Tang and Gupta (1995) presented an algorithm requiring $Θ(n)$ computation time with $Θ(n^2)$ processors. Ramanan (1996) is the first to design an algorithm using only $O(n)$ processors with $O(log^3 n)$ time. Based on the row minima of totally monotone matrices, Bradford et al. (1998) proposed the first parallel polylogarithmic algorithm. It requires $O(n log^{1.5} n)$ time and processors. On CRCW-PRAM machines, Bradford (1994) proposed an $O(lg^4 n)$-time algorithm on $O(n/lnn)$ processors. Later on, Bradford et al. (1998) improved it to $O(lg^{1.5} n \sqrt{lg lg n})$ time using $O(n \sqrt{lg n})$ processors. On the systolic model, Guibas et al. (1979) proposed a bidirectional systolic algorithm on $O(n^2)$ processors with $O(n)$ running time. Karypis and Kumar (1993) mapped a systolic table onto a mesh-connected array of size $n^2$. Myoupo (1992, 1993) proposed a technique for mapping the MCOP to a linear systolic array. On
the hypercube model, Ibarra et al. (1991) presented a solution in $O(n^2)$ time on $O(n)$ processors.

On realistic models of parallel machines, Nishida et al. (2011) presented an efficient parallel implementation of the sequential algorithm of Godbole (1973) on GPU architectures. Ito and Nakano (2013) accelerated this solution. Their idea consists in partitioning the sequential algorithm of Godbole (1973) into many sequential kernel calls, selecting the best values for the size and the number of blocks for each kernel call, and minimizing the memory access overhead. On a NVIDIA GeForce GTX 680 chip, their solution executes in 5.57 seconds for a problem of size 8192. Shyamala et al. (2017) accelerated the computation time through C++ high-performance accelerated massive parallel code. More recently, Diwan and Tembhurne (2019) designed an adaptive generalized mapping method to parallelize non-serial polyadic dynamic-programming problems that utilize GPUs, for efficient mapping of subproblems onto processing threads in each phase. Biswas and Mukherjee (2021) proposed a new memory optimized technique and a versatile technique of utilizing shared memory in blocks of threads to minimize time for accessing dimensions of matrices on GPU architectures. On shared-memory architectures, Tan et al. (2007) proposed a parallel algorithm based on a pipeline to fill the DP table by decomposing the computation operators. Mabrouk (2016) designed solutions based on loop transformations. She showed that the associative expression evaluation technique with the loop interchange transformation provides efficient parallel solutions. Many researchers proposed parallel solutions on the CGM model on distributed-memory architectures. For designing these solutions, the standard methodology consists of subdividing the dependency graph into subgraphs of same size, then distributing these subgraphs fairly among processors, and finally computing them in a suitable evaluation order.

Based on the sequential algorithm of Godbole (1973), Cáceres et al. (2010) proposed a CGM-based parallel solution running in $O(n^3/p)$ execution time with $O(p)$ communication rounds, where $n$ is the input data size and $p$ the number of processors. Higa and Stefanes (2012) used several concepts from the PRAM-based parallel algorithm of Bradford (1994) and designed an $O(n^3/p^3)$ execution time with $O(1)$ communication round. They replaced the traditional DP table used in the sequential algorithm of Godbole (1973) with a weighted DAG. The processors used in the local computation of the algorithm of Dijkstra (1959) and the algorithm of Floyd (1962). Kechid and Myoupo (2009) proposed a CGM-based parallel solution, which runs in $O(n^3/p)$ execution time with $O(p)$ communication rounds. Likewise, the CGM-based parallel solution proposed by Kengne and My-
oupo (2012) runs in $O(n^3/p)$ execution time with $\lceil\sqrt{2p}\rceil$ communication rounds. The common point between these last two solutions is that they use the graph model proposed by Bradford (1994). Before detailing the latter two solutions in Sections 2.2.4 and 2.2.5 respectively, the graph model proposed by Bradford (1994), which enables to transform the MCOP as a shortest path problem, is presented in Section 2.2.3.

### 2.2.3 - Dynamic graph model of Bradford (1994)

Dynamic-programming problems are often solved through the shortest path problems on weighted DAGs. Indeed, Bradford (1994) proposed a graph model called *dynamic graph* for the problem that can be formulated by Equation (1.8). For a problem of size *n*, this graph is denoted by $D_n$ and defined as follows:

**Definition 4** A dynamic graph $D_n = (V, E \cup E')$ is defined as a set of vertices,

$$V = \{(i, j) : 1 \le i \le j \le n\} \cup \{(0, 0)\}$$

a set of unit edges,

$$E = \{(i, j) \to (i, j+1) : 1 \le i \le j < n\} \cup \{(i, j) \uparrow (i-1, j) : 1 < i \le j \le n\} \cup \\ \{(0, 0) \uparrow (i, i) : 1 \le i \le n\}$$

a set of jumpers,

$$E' = \{(i, j) \Rightarrow (i, t) : 1 \le i < j < t \le n\} \cup \{(s, t) \uparrow (i, t) : 1 \le i < s < t \le n\}$$

a weight function $W$ such that

$$W((i, j) \to (i, j+1)) = d_{i-1} \times d_j \times d_{j+1} \qquad 1 \le i \le j < n$$

$$W((i, j) \uparrow (i-1, j)) = d_{i-2} \times d_{i-1} \times d_j \qquad 1 < i \le j \le n$$

$$W((0, 0) \uparrow (i, i)) = 0 \qquad 1 \le i \le n$$

$$W((i, k) \Rightarrow (i, j)) = SP[k+1, j] + d_{i-1} \times d_k \times d_j \quad 1 \le i < k < j \le n$$

$$W((k+1, j) \uparrow (i, j)) = SP[i, k] + d_{i-1} \times d_k \times d_j \quad 1 \le i \le k < j \le n$$

A square matrix of size *n* called *shortest path matrix*, and denoted by *SP*, is used to store in the cell *SP[i, j]* the shortest path from node (0,0) to (i, j). Bradford (1994) showed that the computation of *Cost[i, j]* is equivalent to search in $D_n$ the shortest path from node (0,0) to (i, j). Therefore, it is straightforward to prove that the algorithm of Godbole (1973) presented in Section 2.2.2 (Algorithm 5) is equivalent to compute the shortest paths from (0,0) to the other vertices in a dynamic graph $D_n$, incrementally, diagonal after diagonal, from left to right (Kechid and
<figure><img src="image_15.png" /><figcaption>Figure 13 – Dynamic graphs $D_4$ and $D'_4$ for a problem of size $n = 4$</figcaption></figure>

Myoupo, 2009; Kengne and Myoupo, 2012). Indeed, each path from the root to an edge node corresponds to one of the possible parenthesizes in a dynamic graph given that :

* $(i, j) \to (i, j+1)$ represents the product $(M_i \times M_{i+1} \times \dots \times M_j) \times M_{j+1}$;

* $(i, j) \uparrow (i-1, j)$ represents the product $M_{i-1} \times (M_i \times M_{i+1} \times \dots \times M_j)$;

* $(i, k) \Rightarrow (i, j)$ and $(k+1, j) \uparrow (i, j)$ represent the product $(M_i \times \dots \times M_k) \times (M_{k+1} \times \dots \times M_j)$.

Thus, the shortest path corresponds to the optimal parenthesizing. Figure 13a shows a dynamic graph $D_n$ for $n = 4$. It has the same DAG form representing the dependency between subproblems depicted in Figure 11a, with an additional node $(0,0)$.

Given a problem of size $n$ and its corresponding dynamic graph $D_n$, Theorem 1 can be stated.

**Theorem 1 (Duality theorem)** If the shortest path from node $(0,0)$ to $(i, j)$ needs the edge from node $(i, k)$ to $(i, j)$, then there exists a dual shortest path with the same cost needing the edge from node $(k+1, j)$ to $(i, j)$.

This is a fundamental element of CGM-based parallel algorithms proposed by (Kechid and Myoupo, 2009; Kengne and Myoupo, 2012) because it avoids computation redundancy of the shortest path costs in $D_n$. The input graph of these algorithms is therefore a subgraph of $D_n$ denoted by $D'_n$, in which the set of edges
<figure><img src="image_16.png" /><figcaption>Figure 14 – Shortest path matrix partitioning strategy proposed by Kechid and Myoupo (2009) for $n = 32$ and $p \in \{2, 4, 8\}$</figcaption></figure>

from node $(k+1, j)$ to $(i, j)$, such that $1 \le i < k < j \le n$, is removed. Figure 13b shows the dynamic graph $D'_4$.

### 2.2.4 - CGM-based parallel solution of Kechid and Myoupo (2009)

#### Dynamic graph partitioning

Kechid and Myoupo (2009) partition the shortest path matrix $SP$ into $p(p+1)/2$ submatrices called *blocks*. A block, denoted by $SM(i, j)$, is a matrix of size $\theta(n, p) \times \theta(n, p)$, where $\theta(n, p) = \lceil n/p \rceil$. Equation (2.2) shows entries of the shortest path matrix $SP$ delimiting a block $SM(i, j)$, such that $1 \le i \le j \le n$.

$$SM(i, j) = \begin{pmatrix} SP[i, j - \theta(n, p) + 1] & \cdots & SP[i, j] \\ \vdots & \ddots & \vdots \\ SP[i + \theta(n, p) - 1, j - \theta(n, p) + 1] & \cdots & SP[i + \theta(n, p) - 1, j] \end{pmatrix} \quad (2.2)$$

Figures 14a, 14b, and 14c depict three scenarios of this partitioning for $n = 32$ and $p \in \{2, 4, 8\}$ (the node (0,0) and unit edges of the dynamic graph $D'_{32}$ have not been schematized for readability reasons). In each of them, $SP$ is subdivided into $p$ rows and $p$ columns of blocks. The number of blocks in row $i$ (respectively in column $j$) is $p - i + 1$ (respectively $j$). Usually, all blocks are not full when $n \mod p \neq 0$. Recall that a block is said *full* when the number of rows is equal to the number of columns. In these figures, the number in each block represents the diagonal in which it belongs. The number of blocks in the diagonal $d$ is $p - d + 1$. Thus, there are $p$ blocks at the first diagonal and one block at the last. The blocks of the first diagonal are upper triangular matrices of sizes $\theta(n, p) \times \theta(n, p)$.
**Blocks' dependency analysis**

Kechid and Myoupo (2009) showed that the evaluation of shortest paths for different node blocks of the same diagonal can be carried out in parallel. Indeed, the dependency relationship between blocks, given by Theorem 2, proved that those on the same diagonal are independent.

**Lemma 1 (Nodes' dependency)** *To find the shortest path cost to a node (i, j) in graph $D'_n$, it is necessary to know the shortest path cost to nodes (i, i), ..., (i, j-1) and (i+1, j), ..., (j, j).*

**Proof.** To compute the shortest path to a vertex (i, j), it is necessary, for each of its incoming edges, to have the value of the shortest path to the starting vertex and the weight of this edge. In the dynamic graph $D'_n$, the edges that reach a vertex (i, j) come from the set of vertices : $S'_{i,j} = (i,i),..., (i,j-1) \cup (i+1,j)$. The weight of a jump is $W((i,k) \Rightarrow (i,j)) = SP[k+1,j] + d_{i-1} \times d_k \times d_j$. The computation of weights of edges coming from the vertices of $S'_{i,j}$ requires the values of the shortest paths to the vertices $S''_{i,j} = (i+1,j),..., (j,j)$ of the column j. Thus, the computation of $SP[i,j]$ involves the values of the shortest paths to the vertices of the groups $S'_{i,j}$ and $S''_{i,j}$.

**Lemma 2 (Weights of jumps to a block)** *Computing the weights of jumps from nodes of block SM(i,h) to nodes of block SM(i,j), such that $1 \le i \le h \le j \le n$, only involves the shortest path costs to nodes of block $SM(h - \Theta(n,p) + 2, j)$.*

**Proof.** Computing the weights of jumps from a set of successive nodes in the row i, {(i,h), (i,h+1), ..., (i,l)}, to a node (i,m) such that $i \le h < l < m$, requires only the shortest path costs to nodes {(h+1,m), ..., (l+1,m)}, which is a set of successive nodes of column m. Thus, computing the weight of any row i' of block SM(i,h), i.e. (i',h - Θ(n,p) + 1), (i',h - Θ(n,p) + 2), ..., (i',h) to each node (i',j') of row i' in block SM(i,j), i.e. (i',j - Θ(n,p) + 1), (i',j - Θ(n,p) + 2), ..., (i',j'), ..., (i',j), requires only the shortest path costs to the set of nodes: {(h - Θ(n,p) + 2, j'), (h - Θ(n,p) + 3, j'), ..., (h+1, j')}.

But this set is exactly the set of vertices of column j' in block $SM(h - \Theta(n,p) + 2, j)$ with $j' \in \{(j - \Theta(n,p) + 1), (j - \Theta(n,p) + 2), \dots, (j-1), j\}$. Thus, computing the weights of jumps of each row in $SM(i,h)$ requires only the shortest path costs of nodes in block $SM(h - \Theta(n,p) + 2, j)$.

**Theorem 2 (Blocks' dependency)** *The shortest path costs to every node in blocks $SM(i, j - \Theta(n, p))$, $SM(i, j - 2 \times \Theta(n, p))$, ..., $SM(i, j - (u-1) \times \Theta(n, p))$, and $SM(i$
<figure><img src="image_17.png" /><figcaption>**Figure 15** – Dependencies of two blocks SM(i, j) and SM(h, l) after applying the partitioning strategy of Kechid and Myoupo (2009)</figcaption></figure>

+θ(n, p), j), SM(i + 2 × θ(n, p), j), ..., SM(i + (u − 1) × θ(n, p), j) are mandatory
to computing the shortest paths to nodes of the block SM(i, j), where u = ⌈(j −
i)/θ(n, p)⌋.

Figure 15 illustrates an example of dependencies of two blocks SM(i, j) and SM(h, l). The most shaded blocks are required to evaluate them.

Mapping blocks onto processors

Several distribution schemes have been proposed by Kechid and Myoupo (2009).
The first is a straightforward mapping that consists in assigning, horizontally (re-
spectively vertically), the blocks of the line (respectively column) *i* to processor
*P*<sub>*i*−1</sub>. The corresponding communication scheme, depicted in Figure 16a (respec-
tively in Figure 16b), is easily implementable. However, these mappings make
processors idle since after each step, one processor will not have assigned block.
For example, Figure 16a (respectively Figure 16b) shows that the processor *P*<sub>*p*−1</sub>
(respectively *P*<sub>0</sub>) will be idle after the first step, then it will be the turn of the pro-
cessor *P*<sub>*p*−2</sub> (respectively *P*<sub>1</sub>) after the second step, and so forth. Moreover, these
mappings yield an unbalanced load between them since the processor having the
highest load will evaluate *p* blocks and the one having the lowest load will evaluate
one block.

The second mapping, called *alternative bidirectional projection mapping*, fa-
cilitates the cost prediction and implementation of the algorithm while trying to
<figure><img src="image_18.png" /><figcaption>**Figure 16** – Distribution schemes of blocks on processors proposed by Kechid and Myoupo (2009)</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 17** – Alternative bidirectional projection mapping on four and eight processors</figcaption></figure>

balance the loads between processors. Its principle is illustrated in Figure 16c. An arrow from a block points to the processor to which it is assigned. To distribute the $p(p+1)/2$ blocks of the dynamic graph $D'_n$ on $p$ processors, one of the $p$ blocks of the first diagonal is assigned to a processor. The blocks on other diagonals are projected onto the first diagonal, alternatively, once horizontally and once vertically. This approach allows obtaining horizontal and vertical arrows. Figures 17a and 17b give an example of this mapping on four and eight processors.

The communication scheme corresponding to this mapping is simple and easy to implement like the latter two mappings. However, this mapping halves the number of processors to which a block is sent after being computed. Indeed, a processor owns half of the blocks on the row to which it has been assigned its first block of
the first diagonal, moving from left to right. It also owns half the blocks of the
blocks on the column by moving from bottom to top. In addition, it balances the
loads between processors better than the last two mappings for the reasons stated
above.

The alternative bidirectional projection mapping, by contrast, promotes the idle-
ness of processors from step ( ⌈ *p* / 2 ⌉ + 1). Indeed, starting from this step, processor
*P*_{ ⌈ *p* / 2 ⌉ } becomes idle; for example in Figure 17b, processor *P*₄ after step 4 becomes
idle. Subsequently, after each step, two new processors become idle; for example
in Figure 17b, processors *P*₂ and *P*₃ become idle after step 5. In addition, some
processors from step ( ⌈ *p* / 2 ⌉ + 1) have a large computational load while others are
idle; for example processors *P*₀, *P*₁, and *P*₇ in Figure 17b, which have two blocks
to be performed after step 4.

CGM-based parallel algorithm

Evaluating blocks of the dynamic graph $D'_n$ must be done according to a well-adapted order following dependencies between them. Indeed, according to Theorem 2, the shortest path costs to nodes of a given block belonging to the diagonal $d$ cannot be performed before those contained in each of blocks on which they depend along the preceding diagonals (the diagonals 1, 2, ..., $d-1$). Hence, the evaluation of blocks of a diagonal $d$ can start after computing blocks of the last diagonal of blocks on which they depend, i.e. those of diagonal ($d-1$).

However, Kechid and Myoupo (2009) have shown that the manner of comput-
ing the shortest path costs to nodes of a block allows starting before the end of
the evaluation of blocks belonging to diagonal (d - 1). Thus, the evaluation of the
blocks is done in a progressive fashion and starts as soon as possible. Theorem 3
indicates when computations of the shortest path costs to nodes of a block SM(i, j)
can be started.

**Theorem 3** After computing solutions of each diagonal *h*, ⌈(*j* − *i*)/2⌋ + 1 ≤ *h* ≤ *j* − *i* + 1, at least two possible values of the shortest path from each node in block *SM*(*i*, *j*) can be evaluated.

**Proof.** To evaluate the shortest paths to nodes of a block *SM*(*i*, *j*) of the diagonal *k* = *j* − *i* + 1, for each node, the cost of different edges and the value of the shortest path to the source node of each of these edges are required. *SM*(*i*, *j*) receives edges from blocks *SM*(*i*, *j*'), such that *j'* ∈ [*i*, *j* − 1], and the costs of edges from a block *SM*(*i*, *j*') can be evaluated as soon as the values of the shortest path of nodes of *SM*(*j*', *j*) are computed (from Lemma 2). This is only possible from
1: **for** *d* = 1 to *p* **do**

2: **Finalization** of computations of the shortest path costs to nodes in blocks belonging to the diagonal *d*;

3: **Communication** of block *SM(i, j)* of current diagonal to processors that hold upper and right blocks;

4: **Update** the shortest path costs to each block belonging to diagonals (*d* + 1, *d* + 2, ..., min{2 × (*d* − 1), *p*});

the diagonal *t*, such that (*t* ≥ *j*' - *i* + 1) and (*t* ≥ *j* - *j*' + 1). This implies that: (*t* ≥ ⌊*j* - *i* / 2⌋ + 1. Thus, after evaluating the diagonal *h*, such that *t* ≤ *h* ≤ *j* - *i* + 1, the blocks *SM(i, h + *i* - 1)*, *SM(h + *i* - 1, *j*)*, *SM(i, j* - *h* + 1)*, and *SM(j* - *h* + 1, *j*)* are evaluated. And thus, two possible values of the shortest path to each of nodes of *SM(i, j)* are computable. Those whose incoming edges come from *SM(i, h + *i* - 1)* and *SM(i, j* - *h* + 1)*.

The CGM-based parallel algorithm of Kechid and Myoupo (2009) is a succession of *p* similar steps. The blocks are evaluated in parallel from the first diagonal of blocks, followed by the second, and so on till the last. The overall structure is given by Algorithm 6.

Evaluating the shortest path cost to a node of a block belonging to the diagonal *d* starts at the diagonal ⌊*d*/2⌋. After the computation of blocks on diagonal *d* (line 2 in Algorithm 6), each block is forwarded (line 3 in Algorithm 6) to processors that need these blocks for updating (line 4 in Algorithm 6) or for finalizing (line 2 in Algorithm 6) the computations of values in next steps. In fact, two tasks have to be done at step (*d* + 1):

1. for each block belonging to the diagonal *m*, (*d* + 2 ≤ *m* ≤ 2*d*), some new values may be computed and an update of its temporary values is done (path relaxation principle (Cormen et al., 2009));

2. some new values may be computed and a final update of its temporary values is done for each block of the diagonal (*d* + 1).

These processes are repeated at each step until step *p* where only the finalization phase is carried out. At step *d* = 1, only finalization and communication phases are executed. It is not difficult to notice that computing the shortest path costs (in which jumps from blocks *SM(i, k)* are involved) to blocks *SM(i, j)* is equivalent
<caption>**Algorithm 7** Finalization phase using in Algorithm 6 to evaluate blocks $SM(i, j)$</caption>

1: **for** $d = (j - i - \theta(n, p))$ to $(j - i)$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;**for each** node $(a, b)$ of diagonal $d$ belonging to $SM(i, j)$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;$SP[a, b] \leftarrow \min\{SP[a, b]$, weight of paths whose final edge are jumps coming from block $SM(i, i + \theta(n, p))$, weights of paths whose final edges are internal jumps, weights of paths whose final edges are unit edges}\;

<caption>**Algorithm 8** Updating phase using in Algorithm 6 to refresh the shortest path costs to nodes of blocks $SM(i, j)$</caption>

1: **for** $d = (j - i - \theta(n, p))$ to $(j - i)$ **do**
2: &nbsp;&nbsp;&nbsp;&nbsp;$M_1 \leftarrow$ matrix-multiplication(+, min) $(SM(i, h + i - 1), SM(h + i - 1, j))$;
3: &nbsp;&nbsp;&nbsp;&nbsp;$M_2 \leftarrow$ matrix-multiplication(+, min) $(SM(i, j - h + 1), SM(j - h + 1, j))$;
4: &nbsp;&nbsp;&nbsp;&nbsp;$SP[i, j] \leftarrow \min\{SP[i, j], M_1, M_2\}$;

to the sequential matrix-multiplication $(+, \min)^2$ of the two matrices $SM(i, k)$ and $SM(k - \theta(n, p) + 2, j)$. Pseudocodes of finalization and updating phases in Algorithm 6 are given by Algorithms 7 and 8, respectively.

After the $(j-i)$th iteration of the above algorithms, the only paths which remain to evaluate for nodes in $SM(i, j)$ are those whose last edge is:

1 - either an unit edge (vertical or horizontal);

2 - or a horizontal jump, which comes from an internal node in $SM(i, j)$;

3 - or a horizontal jump, which comes from a node in $SM(i, i)$.

In any case, computing the weight induced by each of these paths (due to these edges) to a node $(i', j')$ of $SM(i, j)$ needs the shortest path cost from a node $(e', f')$ of $SM(i, j)$ such that $f' - e' < j' - i'$. In cases (1) and (2), this value is necessary to compute the shortest path cost of the start node of the last edge. In case (3), this value is necessary for computing the weight of the last edge. Therefore, Algorithm 7 is in fact a classical algorithm of the shortest path in $SM(i, j)$, in which each node can receive a simple edge or a jump from an internal node in the dynamic graph $D'_n$. In summary, Theorem 4 yields the complexity of the CGM-based parallel solution of Kechid and Myoupo (2009).

**Theorem 4** *The CGM-based parallel solution of Kechid and Myoupo (2009) runs in $O(n^3/p)$ execution time with $O(p)$ communication rounds in the worst case.*

<footnote n="2">2. Matrix-multiplication $(+, \min)$ is a matrix multiplication in which the multiplication and summation operations are replaced by addition and the minimum, respectively.</footnote>
### 2.2.5 - CGM-based parallel solution of Kengne and Myoupo (2012)

This section presents the CGM-based parallel solution of Kengne and Myoupo (2012). They proposed a partitioning strategy that reduces the number of sub-graphs (or blocks) of the dynamic graph $D'_n$ as the number of processors increases to minimize the number of communication rounds. The dependencies between blocks are similar to those of Kechid and Myoupo (2009) since it is only the size of blocks that changes. Then, they propose a distribution scheme that fairly assigns blocks onto processors. This mapping balances better the loads between the processors than the alternative bidirectional projection mapping proposed by Kechid and Myoupo (2009) (described in Section 2.2.4). Moreover, it limits to at most two the number of blocks that a processor must evaluate. Nonetheless, the overall structure of their CGM-based parallel algorithm is similar to that of Kechid and Myoupo (2009).

#### Dynamic graph partitioning

Kengne and Myoupo (2012) subdivide the shortest path matrix $SP$ into $f(p)$ rows and $f(p)$ columns of blocks, where $f(p) = ⌊\sqrt{2p}⌋$. A block $SM(i, j)$ is a $\theta(n, p) \times \theta(n, p)$ matrix, where $\theta(n, p) = ⌊n/f(p)⌋$. The size of this block is larger than those obtained after applying the partitioning strategy of Kechid and Myoupo (2009) since $f(p) < p$ when $p > 2$. This means that entries of a block, illustrated by Equation (2.2), are more extensive and therefore, evaluating a block will require more computing time.

Figures 18a, 18b, and 18c show three scenarios of this partitioning for $n = 32$ and $p \in \{2, 3, 4, 5, 6, 7, 8\}$. $SP$ is partitioned into $f(p)(f(p) + 1)/2$ blocks. For example, in Figure 18b (respectively Figure 18c), $SP$ is partitioned into six blocks (respectively ten blocks) when $p \in \{3, 4\}$ (respectively $p \in \{5, 6, 7, 8\}$). Usually, all blocks are not full when $n \bmod f(p) \neq 0$. For example in Figure 18b, while blocks in the first two columns are full, those in the third column are not. There are $p$ blocks at the first diagonal and one block at the last. Recall that in these figures, the number in each block represents the diagonal in which it belongs.

#### Mapping blocks onto processors

The mapping proposed by Kengne and Myoupo (2012) consists in assigning the blocks of a given diagonal from the upper leftmost corner to the lower rightmost corner. First, the blocks of the first diagonal are assigned from processor $P_0$ to processor $P_{f(p)-1}$. Then, the process is repeated on the next diagonal starting with processor $P_{f(p)}$, and so on until a block has been assigned to each processor. The mapping starts again with processor $P_0$, and continues along the diagonals with a
<figure><img src="image_20.png" /><figcaption>**Figure 18** – Shortest path matrix partitioning strategy proposed by Kengne and Myoupo (2012) for n = 32 and p ∈ {2, 3, 4, 5, 6, 7, 8}</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 19** – Snake-like mapping on five and eight processors</figcaption></figure>

snake-like path until the last diagonal $f(p)$. Figures 19a and 19b give an example of this mapping on five and eight processors.

The communication scheme corresponding to this mapping is simple and easy to implement like the alternative bidirectional projection mapping proposed by Kechid and Myoupo (2009). Moreover, this mapping allows processors to remain active as soon as possible and minimizes the number of communication rounds. It also ensures load balancing because it allows some processors to evaluate at most one block more than others. For example, Figure 19b shows that only processors $P_0$ and $P_1$ evaluate one more block than others. However, this mapping does not optimize communications since a processor does not usually hold the upper blocks that are on the same row and column as the block that has been assigned to that processor.
**Algorithm 9** CGM-based parallel algorithm of Kengne and Myoupo (2012) to solve the MPP

1: **for** $d = 1$ to $f(p)$ **do**

2: **Finalization** of computations of the shortest path costs to nodes in blocks belonging to the diagonal $d$;

3: **Communication** of block $SM(i, j)$ of current diagonal to processors that hold upper and right blocks;

4: **Update** the shortest path costs to each block belonging to diagonals $(d+1, d+2, \dots, \min\{2 \times (d-1), f(p)\})$;

**Lemma 3** After partitioning the shortest path matrix into $b$ submatrices (blocks) such that $p < b \le 2p$ and applying the snake-like mapping, each processor evaluates at most two blocks.

*Proof.* The shortest path matrix is partitioned into $b = f(p)(f(p) + 1)/2$ blocks, where $f(p) = \lceil \sqrt{2p} \rceil$. It is deducible that $p < b \le 2p$. Given that the snake-like mapping allows processors to evaluate at most one block more than others, each processor evaluates at most two blocks. ■

## CGM-based parallel algorithm

The CGM-based parallel algorithm of Kengne and Myoupo (2012) is given by Algorithm 9. It consists of $f(p)$ similar steps as the one in Kechid and Myoupo (2009). As seen in Section 2.2.4, and according to Theorem 3, the updates for a block $SM(i, j)$ (which concern jumps from block $SM(i, k)$) are equivalent to a matrix multiplication $(+, \min)$ of matrices $SM(i, k)$ and $SM(k - \theta(n, p) + 2, j)$. Evaluating a block of diagonal $d$ starts as soon as the diagonal $\lceil d/2 \rceil$ is evaluated, i.e. at step $(\lceil d/2 \rceil + 1)$.

**Theorem 5** The CGM-based parallel solution of Kengne and Myoupo (2012) runs in $O(n^3/p)$ execution time with $\lceil \sqrt{2p} \rceil$ communication rounds in the worst case.

## 2.3 - Optimal binary search tree problem

### 2.3.1 - Overview

A tree is a data structure that hierarchically represents a set of elements (often called keys) (Deo, 1974; Knuth, 1997). It consists of a set of vertices and edges. A node stores a key and an edge is a link between two vertices. A path is a sequence of
vertices in which two consecutive vertices are connected by an edge (Deo, 1974).
Any two vertices in a tree are connected by a single path. In fact, a tree is an
undirected acyclic connected graph (Knuth, 1997).

A rooted tree is a tree in which a vertex has been designated as the root (Knuth, 1997). Consider that there is a path from the root vertex to a vertex x. If there is a path from x to a vertex y, then x is an ancestor of y, and y is a descendant of x (Deo, 1974). If there are no intermediate vertices on the path between x and y, then x is a father of y, and y is a child of x (Deo, 1974). The root vertex is the only vertex that has no parent in a rooted tree. If there is no path between vertex x and another vertex, then x is a leaf (Deo, 1974).

A binary tree is a rooted tree where each vertex has at most two children: a left child and a right child (Deo, 1974; Knuth, 1997). The subtree of a vertex x in a tree is the part of that tree containing x, as well as its descendants. Thus, the left subtree (respectively the right subtree) of a vertex x is the subtree of the left child (respectively the right child) of x (Deo, 1974; Knuth, 1997).

A binary search tree is a binary tree that maintains a set of sorted keys according to the following rule: for each vertex x of the tree, the key of x is larger than all the keys contained in the left subtree of x and it is smaller than all the keys contained in the right subtree of x (Nagaraj, 1997). Figure 20 shows an example of a binary search tree corresponding to the set of keys {c, e, f, g, h, k, l, n, o, r, s} sorted in alphabetical order. A binary search tree is typically used to efficiently search for a particular key among a set of sorted keys. The cost of searching for this key is equal to length of the path from the root to this key plus one, which is in fact, the depth of this key (Cormen et al., 2009). For example, in Figure 20, the cost of searching for the key "n" (the root vertex) is equal to 1 and the cost of searching for the key "h" is equal to 5.

The overall cost of searching, which is equal to the sum of the cost of search-
ing of each key, must be usually as less as possible in many applications. Self-
balancing binary search trees<fnref n="3" /> minimize this cost by ensuring that all keys are as
near the root vertex as possible. However, when a query frequency is associated
with each key, this tree will not be the most efficient in some cases. Indeed, it
would be better to place the most frequently queried keys closer to the root vertex,
and to place the less frequently queried keys further from the root vertex. For ex-
ample, consider the set of sorted keys {a, b, c} with respective frequencies of 3, 1,
and 7. The cost of searching for a key x is equal to the depth of x multiplied by

<footnote n="3">3. A self-balancing binary search tree is a binary search tree that automatically attempts to keep its height (maximal number of levels below the root vertex) as small as possible at all times (Knuth, 1998).</footnote>
<figure><img src="image_22.png" /><figcaption>**Figure 20** – Example of a binary search tree corresponding to the set of letters `c, e, f, g, h, k, l, n, o, r, s` sorted in alphabetical order</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21** – Five possible binary search trees obtained from the set of sorted keys `{a, b, c}`</figcaption></figure>

the frequency of `x`. Of the five possible binary search trees, depicted in Figure 21, that can be obtained from these keys, the tree that minimizes the overall cost of searching is determined as follows:

* for the first tree, it is equal to $(1 \times 3) + (2 \times 7) + (3 \times 1) = 20$;

* for the second tree, it is equal to $(1 \times 3) + (2 \times 1) + (3 \times 7) = 26$;

* for the third tree, it is equal to $(1 \times 1) + (2 \times 3) + (2 \times 7) = 21$;

* for the fourth tree, it is equal to $(1 \times 7) + (2 \times 1) + (3 \times 3) = 18$;

* for the fifth tree, it is equal to $(1 \times 7) + (2 \times 3) + (3 \times 1) = 16$.

As shown in Figure 21, the third binary search tree is balanced but the fifth is not. Nevertheless, the overall search cost of the latter is lower than the self-balanced binary search tree (the third tree in Figure 21) and is the lowest of all; thus, it is the optimal binary search tree for the above example.
More formally, the optimal binary search tree (OBST) problem can be defined as follows: consider a set of $n$ sorted keys $K = \langle k_1, k_2, \dots, k_n \rangle$, such that $k_1 < k_2 < \dots < k_n$. The probability of finding a key $k_i$ is denoted by $p_i$. In the case where the searched element is not in $K$, consider a set of $(n+1)$ dummy keys $D = \langle d_0, d_1, \dots, d_n \rangle$. Indeed, $d_0$ represents the set of values that are smaller than $k_1$ and $d_n$ the set of values that are larger than $k_n$. A dummy key $d_i$, such that $1 \le i < n$, represents the set of values between $k_i$ and $k_{i+1}$. The probability of finding a dummy key $d_i$ is denoted by $q_i$. In summary, either the search ends in success (i.e. finding the searched key $k_i$) or in failure (i.e. finding a dummy key $d_i$); thus the sum of probabilities of success and failure is equal to 1. A binary search tree $T$ constructed from these keys is composed of $n$ internal vertices in the set $K$ and $(n+1)$ leaves in the set $D$. Equation (2.3) gives the overall cost of searching of $T$:

$$ Cost(T) = \sum_{i=1}^{n} (depth(k_i) \times p_i) + \sum_{i=0}^{n} (depth(d_i) \times q_i) \quad (2.3) $$

where $depth(x)$ is the depth of the key $x$. The OBST problem consists in finding the binary search tree that will have the lowest overall cost of searching.

Let $w(i, j) = p_{i+1} + \dots + p_j + q_i + q_{i+1} + \dots + q_j$. The minimum cost of searching of a tree $T_{i,j}$ from the set of keys $K_{i,j} = \langle k_{i+1}, k_{i+2}, \dots, k_j \rangle$ and the set of dummy keys $D_{i,j} = \langle d_i, d_{i+1}, \dots, d_j \rangle$ is denoted by $Tree[i, j]$ and defined by:

$$ Tree[i, j] = \begin{cases} q_i & \text{if } 0 \le i = j \le n, \\ \min_{i \le k < j} \{Tree[i, k] + Tree[k+1, j] + w(i, j)\} & \text{if } 0 \le i < j \le n. \end{cases} \quad (2.4) $$

Finding the optimal binary search tree from the set of keys $K$ and the set of dummy keys $D$ is reduced to compute $Tree[0, n]$. Equation (2.4) is equivalent to Equation (1.8) as $Opt = \min, Init(i) = q_i$, and $F(i, k, j) = w(i, j)$, where $0 \le i \le k \le j \le n$.

The straightforward sequential algorithm of Godbole (1973), running in $O(n^3)$ time and $O(n^2)$ space, is given by Algorithm 10. The DP table (named $Tree$ in Algorithm 10) stores the value of the optimal cost of searching of $T_{i,j}$ (see line 10 in Algorithm 10). The tracking table (named $Cut$ in Algorithm 10) stores the value of the index $k$, which minimizes $Tree[i, j]$ (see line 11 in Algorithm 10). It is the optimal decomposition value of $T_{i,j}$ into two subtrees. For a problem of size $n=3$, the DAG and the DP table are illustrated in Figures 22a and 22b, respectively.

Applications of optimal search binary trees are numerous. One of the most obvious is the search for a word in a dictionary. Indeed, from the words of a dictionary and the access frequency of each word, an optimal binary search tree can be built to quickly and efficiently answer to a query. An application derived from the latter
<caption>**Algorithm 10** Sequential algorithm of Godbole (1973) to solve the OBST problem</caption>

1: **for** i = 0 to n **do**
2:     Tree[i, i] ← qi;
3:     **for** d = 1 to n **do**
4:         **for** i = 0 to n - d + 1 **do**
5:             j ← n - d + 1;
6:             Tree[i, j] ← ∞;
7:             **for** k = i to j − 1 **do**
8:                 c ← Tree[i, k] + Tree[k + 1, j] + w(i, j);
9:                 **if** c < Tree[i, j] **then**
10:                 Tree[i, j] ← c;
11:             Cut[i, j] ← k;

<figure><img src="image_24.png" /><figcaption>**Figure 22** – Task graph and the dynamic-programming table used to compute Tree[0,3]</figcaption></figure>

is the translation of a word from one language to another (Cormen et al., 2009). El-
Qawasmeh (2004) used an optimal binary search tree to solve the word prediction
problem. This problem attempts to guess and update the next word in a sentence
as it is typed. In view of these different applications, an optimal binary search tree
should be constructed from the sequential algorithm of Knuth (1971), which is an
improved version of the sequential algorithm of Godbole (1973), running in O(n²)
time and space.

**2.3.2 - Sequential algorithm of Knuth (1971)**

Knuth (1971) noticed entries of tracking table *Cut* satisfy the *monotonicity property*
such that *Cut*[i, j − 1] ≤ *Cut*[i, j] ≤ *Cut*[i + 1, j] for all 0 ≤ i < j ≤ n. This property
means that computing *Cut*[i, j] consists in finding all indices between *Cut*[i, j − 1]
<caption>**Algorithm 11** Sequential algorithm of Knuth (1971) to solve the OBST problem</caption>

1: **for** *i* = 0 to *n* **do**
2:     *Tree[i,i]* ← *q*<sub>*i*</sub>;
3:     **for** *d* = 1 to *n* **do**
4:         **for** *i* = 0 to *n* − *d* + 1 **do**
5:             *j* ← *n* − *d* + 1;
6:             *Tree[i,j]* ← ∞;
7:             **for** *k* = *Cut[i,j* − 1] to *Cut[i* + 1, *j*] **do**
8:                 *c* ← *Tree[i,k]* + *Tree[k* + 1, *j*] + *w(i,j)*;
9:                 **if** *c* < *Tree[i,j]* **then**
10:                 *Tree[i,j]* ← *c*;
11:             *Cut[i,j]* ← *k*;

<table><caption>**Table 1** – Example of probabilities of three sorted keys</caption><thead><tr><th><i>i</i></th><th>0</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th><i>p</i><sub><i>i</i></sub></th><td></td><td>0.25</td><td>0.4</td><td>0.15</td></tr><tr><th><i>q</i><sub><i>i</i></sub></th><td>0.03</td><td>0.06</td><td>0.07</td><td>0.04</td></tr></tbody></table>

and *Cut[i+1, j]*, instead of between *i* and *j* as it is done in line 7 of Algorithm 10. Knuth (1971) thus achieves to evaluate the Θ(*n*²) subproblems in constant time. Algorithm 11 draws a big picture.

## Parallelization constraint of the sequential algorithm of Knuth (1971)

One downside to the speedup of Knuth (1971), however, is that the number of comparison operations for entries in the same diagonal varies from one entry to another compared to the classical version. Consider for example the set of three sorted keys with probabilities given in Table 1. Figures 23a and 23b depict, respectively, the DP table and the tracking table that have been filled while determining the optimal binary search tree from these probabilities. Figure 23a shows that the value of *Tree*[0,3] is computed only from the optimal solutions of the pair of subproblems ((0,1), (2,3)), since *Cut*[0,2] and *Cut*[1,3] are equal to 1. Yet in the classical version, the value of *Tree*[0,3] should be computed from the optimal solutions of the pairs of subproblems ((0,0), (1,3)), ((0,1), (2,3)), and ((0,2), (3,3)). Therefore, when designing parallel algorithms based on the sequential algorithm of Knuth (1971), there is no guarantee that the processors will have the same load if entries on a diagonal are fairly distributed among them, as shown in Figures 19a and 19b for example. It is the main parallelization constraint of this algorithm.
<figure><img src="image_25.png" /><figcaption>**Figure 23** – Dynamic-programming and tracking tables filled while determining the optimal binary search binary tree from probabilities of three sorted keys given in Table 1</figcaption></figure>

## Literature review on the parallelization of the sequential algorithm of Knuth (1971)

Although parallelization of the sequential algorithm of Godbole (1973) has been widely studied in the literature on different parallel computing models (see Section 2.2.2), little work has been done on parallelization of the sequential algorithm of Knuth (1971). On the PRAM model, Karpinski and Rytter (1994) designed a sublinear time parallel algorithm to construct optimal binary search trees. It requires $O(n^{1-\epsilon}\log n)$ computation time with the total work $O(n^{2+2\epsilon})$ for an arbitrarily small constant $0 < \epsilon \le 0.5$. Recall that the total work is equal to the computation time multiplied by the number of processors (Karpinski and Rytter, 1994). On CREW-PRAM machines, Karpinski et al. (1996) presented an algorithm running in $O(n^{0.6})$ computation time with $n$ processors. On realistic models of parallel machines, Craus (2002) proposed an $O(n)$-time parallel algorithm using $O(n^2)$ processors on very-large-scale integration (VLSI) architectures. Wani and Ahmad (2019) proposed a parallel implementation of the sequential algorithm of Knuth (1971) on GPU architectures. For a problem of size 16384, they achieved a speedup factor of 409 on a NVIDIA GeForce GTX 570 chip and a speedup factor of 745 on a NVIDIA GeForce GTX 1060 chip. On shared-memory architectures, after demonstrating that dependencies of subproblems available in the code implementing the sequential algorithm of Knuth (1971) allow generating only 2D tiled code, using the polyhedral model, Bielecki et al. (2021) proposed a way of trans-
forming this algorithm to a modified one exposing dependencies to generate 3D
parallel tiled code. Experimentations conducted using OpenMP demonstrated that
the 3D tiled code considerably outperforms the 2D tiled code.

On distributed-memory architectures, Kechid and Myoupo (2008a) proposed the first CGM-based parallel solution to solve the OBST problem. It requires $O(n^2/p)$ execution time with $O(p)$ communication rounds. They improved this solution in their works done in (Kechid and Myoupo, 2009) to solve the MPP. Myoupo and Kengne (2014b) built a CGM-based parallel solution, from their proposed solution in (Kengne and Myoupo, 2012) to solve the MPP, running in $O(n^2/\sqrt{p})$ execution time with $\lceil\sqrt{2p}\rceil$ communication rounds. This solution uses the same partitioning strategy and the same distribution scheme strategy; except that the evaluation of blocks is not done in a progressive manner. Indeed, Myoupo and Kengne (2014b) showed that the progressive evaluation is not suitable because of the speedup of Knuth (1971). Later on, Kengne et al. (2016) highlighted the existence of a relationship between the execution time, the load balancing, and the number of communication rounds in a family of CGM-based parallel solutions solving the OBST problem. They proposed a general methodology for deriving a CGM-based parallel algorithm in accordance with the user's parameters. This solution is described more succinctly in Section 2.3.3.

2.3.3 - CGM-based parallel solution of Kengne et al. (2016)

Looking at the performance of CGM-based parallel solutions proposed by Kechid
and Myoupo (2008a) and Myoupo and Kengne (2014b) to solve the OBST prob-
lem (and the solutions proposed by Kechid and Myoupo (2009) and Kengne and
Myoupo (2012) to solve the MPP), Kengne et al. (2016) noticed that the execution
time was related to the load balancing and the number of communication rounds.
These criteria depend on the partitioning strategy and the distribution scheme strat-
egy used when designing the CGM-based parallel solutions:

1. When the dependency graph is subdivided into small-size blocks (Kechid and Myoupo, 2008a, 2009), the load difference between processors is small if one processor has one more block than another. However, the number of communication rounds will be high.

2. When the dependency graph is subdivided into large-size blocks (Kengne and Myoupo, 2012; Myoupo and Kengne, 2014b), the number of communication rounds of the corresponding algorithm is reduced since there are few blocks. However, the load of processors will be unbalanced.
These criteria in fact have a significant impact on the global communication time, which in turn has an impact on the total execution time. Indeed, when the blocks are small, the number of communication rounds is high since there are many blocks. Thus, excessive communication will lead to communication overhead, which will deteriorate the global communication time. On the other hand, when the blocks are large, the evaluation of a block will take a long time because of the load imbalance. So a processor that is waiting for this block to start or continue its computation will wait longer to receive this block. As a result, Kengne et al. (2016) showed that these criteria are contradictory because it is difficult to simultaneously optimize them in a CGM-based parallel solution.

Kengne et al. (2016) proposed a CGM-based parallel solution that gives the end-user the choice to optimize one criterion according to their own goal. They generalized the ideas of dependency graph partitioning and distribution scheme introduced in (Kechid and Myoupo, 2008a, 2009; Kengne and Myoupo, 2012; Myoupo and Kengne, 2014b). The end-user derives a CGM-based parallel algorithm from these parameters:

* $g$ (the granularity) : it is used in the partitioning strategy to define the size and number of blocks;

* $p$ (the number of processors) : it is used in the distribution scheme strategy to assign blocks onto processors.

The complexity of the derived CGM-based parallel algorithm depends on these parameters. Kengne et al. (2016) also proposed a simulator based on the work of Fotso et al. (2010), which allows users to modify these parameters and observe the performance (number of blocks, number of communication rounds, load balancing, and efficiency) of the derived CGM-based parallel algorithm.

## Dependency graph partitioning

Kengne et al. (2016) partition the DAG into $g(g+1)$ blocks. A block $SM(i, j)$ is a matrix of size $\theta(n,g) \times \theta(n,g)$, where $\theta(n,g) = \lceil(n+1)/g\rceil$. This strategy encompasses all the previous ones because:

* it is equivalent to Kechid and Myoupo (2008a) when $g = 2p$;

* it is equivalent to Kechid and Myoupo (2009) when $g = p$;

* it is equivalent to Myoupo and Kengne (2014b) when $g = \lceil\sqrt{2p}\rceil$.
It is straightforward to notice that g can be larger, smaller or equal to p. Nevertheless, the blocks obtained after the partitioning must be distributed fairly to processors whatever the value of g.

## Blocks' dependency analysis

Figures 24a and 24b respectively show an example of dependencies and extremities of a block $SM(i, j)$ after applying the partitioning strategy of Kengne et al. (2016). The extremities of $SM(i, j)$ are defined by:

* the leftmost upper entry $LUE_{ij} = (i, j - \theta(n, g) + 1)$;

* the rightmost upper entry $RUE_{ij} = (i, j)$;

* the leftmost lower entry $LLE_{ij} = (i + \theta(n, g) - 1, j - \theta(n, g) + 1)$;

* the rightmost lower entry $RLE_{ij} = (i + \theta(n, g) - 1, j)$.

Figure 24a illustrates eight points (A, B, C, D, E, F, G, and H) that identify blocks on which the block $SM(i, j)$ depends:

* $A = \text{Tree}[LUE_{i,j-3 \times \theta(n,g)-1}]$

* $B = \text{Tree}[LUE_{i,j-\theta(n,g)-2}]$

* $C = \text{Tree}[LLE_{i,j-\theta(n,g)-2}]$

* $D = \text{Tree}[LLE_{i,j-3 \times \theta(n,g)-1}]$

* $E = \text{Tree}[LLE_{i+\theta(n,g)+1,j}]$

* $F = \text{Tree}[RLE_{i+\theta(n,g)+1,j}]$

* $G = \text{Tree}[RLE_{i+3 \times \theta(n,g),j}]$

* $H = \text{Tree}[LLE_{i+3 \times \theta(n,g),j}]$

All lower blocks that are in the same row and column as the block $SM(i, j)$ are no longer absolutely required to evaluate $SM(i, j)$ as a consequence of the speedup of Knuth (1971) (Kengne, 2014). In Figure 24a for example, only the five most shaded blocks will be needed. In other cases, it could have been just four, three or two blocks. Moreover, it may happen that some nodes of blocks located in extremities are needed instead of the whole blocks. It depends on the input data. Therefore, the speedup of Knuth (1971) does not allow estimating the exact load of a block before its evaluation (Kengne, 2014).
<figure><img src="image_26.png" /><figcaption>**Figure 24** – Dependencies and extremities of a block SM(i, j) after applying the partitioning strategy of Kengne et al. (2016)</figcaption></figure>

However, the blocks on which the block SM(i, j) depends are not located on the same diagonal as SM(i, j). Consequently, they can be carried out in parallel.

## Mapping blocks onto processors

Kengne et al. (2016) implemented the alternative bidirectional projection mapping (described in Section 2.2.4, page 55) and snake-like mapping (described in Section 2.2.5, page 60). They extended these strategies and let the end-user choose one of them according to their own goal.

The alternative bidirectional projection mapping is not practical when $g < p$ because no blocks will be assigned to processors $P_g, P_{g+1}, \dots, P_{p-1}$; these processors will remain idle throughout the problem-solving process. It is only practical when $g = k \times p$, where k is a positive integer. k blocks of the first diagonal are assigned to a processor. The blocks on other diagonals are then projected onto the first diagonal. Figure 25a shows an example of this mapping when $g = 8$ and $p = 4$. This mapping minimizes communications and makes the processors less idle compared to Figure 17b when $g = p = 8$. Indeed, a processor has more than half of the blocks in the row and in the column to which it has been assigned its first blocks of the first diagonal. Furthermore, in Figure 25a, the first processor (processor $P_2$)
<figure><img src="image_27.png" /><figcaption>**Figure 25** – Alternative bidirectional projection mapping and snake-like mapping on four processors when $g = 8$</figcaption></figure>

becomes idle after step 6; yet in Figure 17b, it is after step 4 that the first processor ($P_4$) becomes idle.

The snake-like mapping is convenient whatever the value of $g$. When $g \ge p$, the blocks of the first diagonal are assigned from processor $P_0$ to processor $P_{p-1}$. Then, the process is repeated starting again with processor $P_0$ until all blocks of the first diagonal are assigned (when $g > p$), continues to the next diagonal and so on until a block has been assigned to each processor. Figure 25b gives an example of this mapping when $g = 8$ and $p = 4$. This mapping does not minimize communications enough compared to alternative bidirectional projection mapping. However, Kengne et al. (2016) have shown that the snake-like mapping is more efficient and balances the load of processors better than alternative bidirectional projection mapping.

## CGM-based parallel algorithm

The analysis of the dependency of blocks showed that the evaluation of a block does not always depend on all lower blocks that are in the same row and column as this block. Evaluating blocks in a progressive fashion, like the MPP (see Section 2.2.4, page 57), would not be ideal because some unnecessary computations could be performed (Kengne, 2014). Kengne et al. (2016) evaluated blocks in a non-progressive fashion, i.e. the evaluation of blocks of the diagonal $d$ starts after computing blocks of the diagonal ($d-1$).
**Algorithm 12** CGM-based parallel algorithm of Kengne et al. (2016) to solve the OBST problem

1: **for** d = 1 to *R* **do**

2: **Computation** of blocks belonging to the round *d* using Algorithm 11;

3: **Communication** of entries (*Tree* and *Cut* tables) required for computing each block of rounds {*d* + 1, *d* + 2, ..., *R*};

Lemma 4 states the number of communication rounds of the CGM-based par-
allel algorithm of Kengne et al. (2016). It is equal to *g* when *g* ≤ *p* because each
diagonal of blocks contains at most *p* blocks; this means that each processor will
evaluate at most one block per diagonal. This is not the case when *g* > *p*. In fact, a
processor will evaluate at most *k* blocks per diagonal when *g* = *k* × *p*. For example
in Figures 25a and 25b, a processor evaluates at most two blocks from diagonal 1
to diagonal 4; and there is 12 communication rounds. Considering *R* as the number
of communication rounds regardless of the values of *g* and *p*, Algorithm 12 gives
an overview of the CGM-based parallel algorithm of Kengne et al. (2016).

**Lemma 4** *The number of communication rounds of the CGM-based parallel algorithm of Kengne et al. (2016) is defined by:*

$$
R = \begin{cases}
g & \text{if } g \le p, \\
\left\lfloor \frac{g(g+1)-p(p-1)}{2p} \right\rfloor + p & \text{if } g > p.
\end{cases}
$$

*Proof.* When *g* ≤ *p*, it is clear that there is *g* communication rounds. When *g* > *p*, diagonals that contain at least *p* blocks are evaluated in $\left\lfloor \frac{g(g+1)-p(p-1)}{2p} \right\rfloor$ communication rounds. The last (*p* − 1) diagonal of blocks are evaluated in (*p* − 1) rounds. The potential blocks located in the diagonal that contains *p* blocks and has not yet been considered can be evaluated in at most one round.

■

**Theorem 6** *The CGM-based parallel solution of Kengne et al. (2016) runs in O(n² × R/g²) execution time with R communication rounds in the worst case.*

**Simulator of Kengne et al. (2016)**

Kengne et al. (2016) proposed a simulator to allow users to modify the parameters *g* and *p* and observe the performance of the derived CGM-based parallel algorithm. From the partitioning strategy and the distribution scheme strategy, they easily deduce the size and number of blocks, and the number of communication rounds. They relied on the work of Fotso et al. (2010), which computes the load of blocks