Iterative linear algebra methods are the important parts of the overall computing time of applications in various fields since decades. Recent research related to social networking, big data, machine learning and artificial intelligence has increased the necessity for non-hermitian solvers associated with much larger sparse matrices and graphs. The analysis of the iterative method behaviors for such problems is complex, and it is necessary to evaluate their convergence to solve extremely large non-Hermitian eigenvalue and linear problems on parallel and/or distributed machines. This convergence depends on the properties of spectra. Then, it is necessary to generate large matrices with known spectra to benchmark the methods. These matrices should be non-Hermitian and nontrivial, with very high dimension. A scalable parallel matrix generator SMG2S that uses the user-defined spectrum to construct large-scale sparse matrices and ensures their eigenvalues as the given ones with high accuracy is implemented based on MPI and C++11. This report gives the manual of SMG2S.

Chapter 1 Introduction 1.1 Getting Started SMG2S (Scalable Matrix Generator with Given Spectrum) [START_REF] Galicher | Generate very large sparse matrices starting from a given spectrum[END_REF][START_REF] Wu | A parallel generator of non-hermitian matrices computed from given spectra[END_REF] is a software which provides for generating the non-Hermitian Matrices with user-customized eigenvalues. SMG2S is implemented in parallel based on MPI (Message Passing Interface) and C++11 to support efficiently the generation of test matrices in parallel on distributed memory platforms.

Iterative linear algebra methods are essential for the applications in various fields. The analysis of the iterative method behaviors is complex, and it is necessary to evaluate their convergence to solve extremely large non-Hermitian eigenvalue and linear problems on parallel and/or distributed machines. This convergence depends on the properties of spectra. Thus, we propose SMG2S to generate large matrices with known spectra to benchmark these methods. The generated matrices are non-Hermitian and non-trivial, with very high dimension.

The functionality proposed inside SMG2S can verify the ability of SMG2S to keep the accuracy of a given spectrum. This function is based on the shift inverse power method. SMG2S also gives a graphic user interface to compare the given and final spectral distribution for the verification.

We will describe the following subset of the SMG2S.

• Parallel Vector and Matrix: this part presents the functions implemented in SMG2S to establish parallel vector and matrix over distributed memory platforms.

• Nilpotent Matrix Object: this part presents a special nilpotent matrix object for the matrix generation procedure in SMG2S.

• Generating Matrix with prescribed eigenvalues: this part gives the way to use SMG2S to generate required test matrices.

• Interface to Other Languages/Libraries: this part introduces the interface of SMG2S to other languages and existing scientific computational libraries such as PETSc and Trilinos.

• Verification of Eigenvalues of Generated Matrix: this part gives the way to verify the accuracy of eigenvalues of generated matrices comparing with given spectrum. A graphic user interface is also provided to facilitate the comparison.

Installation

To obtain SMG2S, please follow the instructions at the SMG2S download page: https://smg2s.github.io/download.html. Prerequisites:

• C++ Compiler with c++11 support;

• Cmake (version minimum 3.6);

• (Optional) PETSc and SLPEc are necessary for the verification of the ability to keep the given spectrum.

In the main directory:

Programming Language in SMG2S

SMG2S is a collection of templated header files written in C++. The wrappers to C and Python codes are provided. The users of PETSc or Trilinos can directly use SMG2S with the interfaces implemented inside.

Referencing SMG2S

Place cite these papers if you want to reference SMG2S. The distributed vector inside SMG2S is implemented using the C++ programming language and MPI. The parallel vector implementation is composed of two main classes: parV ectorM ap and parV ector. parV ectorM ap class is a vector index map which controls the partitioning and distribution over the processes, and parV ector is the parallel vector itself, which contains the actual distributed data and the corresponding functions controlling the data.

Vector Map

The distribution of a set of integer labels (or elements) across the processes is here called a parV ectorM ap. In SMG2S, it is implemented with the help of std :: map, which maps the proc number and related integers. Basically, parV ectorM ap handles the definition of global and local indices for the mapping across the processes. Here we give some methods of parV ectorM ap < S >: S G e t G l o b a l S i z e ( ) ; // Get t o t a l i n t e r g e r number f o r a l l p r o c s .

Creating a Distributed Vector

A parallel vector parV ector across the processes can be created with the help of parV ectorM ap. The entries of parV ector are distributed over different processes referring to the parV ectorM ap. parV ector < T, S > is implemented with the C++ template, with T the data type of entry, and S the data type of index. This is the constructor of parV ector, with M P I Comm ncomm the working MPI communicator, Slbound and Subound respectively the lower and upper bound of vector global indices on each proc. S G e t G l o b a l S i z e ( ) ; // Get t o t a l i n t e r g e r number f o r a l l p r o c s ;

T * GetArray ( ) ; // Get t h e a r r a y c o n t a i n g t h e e n t r i e s on each pro c ;

void S e t V a l u e L o c a l ( S row , T v a l u e ) ; // i n s e r t v a l u e i n t h e l o c a l i n d e x named row ;

void S e t V a l u e s L o c a l ( S nindex , S * rows , T * v a l u e s ) ; // i n s e r t a r r a y i n t h e l o c a l i n d i c e s named rows ;

void SetValueGocal ( S row , T v a l u e ) ; // i n s e r t v a l u e i n t h e g o c a l i n d e x named row ;

void S e t V a l u e s G o c a l ( S nindex , S * rows , T * v a l u e s ) ; // i n s e r t a r r a y i n t h e g o c a l i n d i c e s named rows ;

void SetTovalue (T v a l u e ) ; // S e t t h e e n t r i e s o f p a r V e c t o r a l l t o same g i v e n v a l u e ; void MA( N i l p o t e n c y <S> n i l p , p a r M a t r i x S p a r s e <T, S> * prod ) ; // m a t r i x m u l t i p l e a s p e c i a l n i l p o t e n t m a t r i x ;

void AM( N i l p o t e n c y <S> n i l p , p a r M a t r i x S p a r s e <T, S> * prod ) ; // s p e c i a l n i l p o t e n t m a t r i x m u l t i p l e a n o t h e r m a t r i x ;

Creating a Distributed Matrix

The is an example of creating a distributed matrix with the mapping of parallel vector:

/

Introduction

The nilpotent matrix is very important for the generation of test matrices with given spectrum. Il can be defined by several parameters, the explicit implementation is not necessary. The three parameters defined a nilpotent matrix is listed as:

• dIagPostion: the distance of the off-diagonal to the diagonal, refering to p in Fig. 3.1;

• length: the continuous one on the off-diagonal of nilpotent matrix, refering to d in Fig. 3.1;

• probSize: the number of row/column of nilpotent matrix, refering to n in Fig. 3.1. 14

Different Types of Nilpotent Matrix

The different nilpotent matrix will influence the sparsity pattern of the final generated matrix.

• NilpType1: diagPostion = 2

• NilpType2: diagPostion = 3

• NilpType3: diagPostion > 3 3.3

Parameter Validation for Nilpotent Matrix

• NilpType1: parameter length can be any integer value > 0;

• NilpType2: parameter length sould be even;

• NilpType3: validation of parameter length is complex. length should be divisible by p.

Chapter 4

Generating Matrix with SMG2S

SMG2S Class

The header file ./smg2s/smg2s.h implement the matrix generation method. It is defined as:

template<typename T, typename S> p a r M a t r i x S p a r s e <T, S> * smg2s ( S p r o b S i z e , N i l p o t e n c y <S> n i l p , S lbandwidth , s t d : : s t r i n g spectrum , MPI Comm comm )

Inside the definition, typename T is to define the size of matrix, and typename S is to define the scalar types of entries of matrix. We give the meaning of the input parameter as below:

• S ProbSize: the size of matrix to generate;

• Nilpotency<S> nilp: the nilpotent matirx object for generation;

• S lbandwidth: the bandwidth of lower-diagonal band of initial matrix;

• std::string spectrum: the file path of spectra file;

• MPI Comm comm: the working MPI communicator.

Generation Workflow

Customize the Low Band of Initial Matrix

We know that the low band bandwidth of initial matrix can be set by the parameter lbandwidth of smg2s. Additionaly, the distribution of entries of initial matrix can also be customized by the function matInit provided by the file ./verification/tests/specGen.h. En default, these entries are filled in random. The different mechanism to fill them will influence the sparsity of final generated sparse matrix. template<typename T, typename S> void m a t I n i t ( p a r M a t r i x S p a r s e <T, S> * Am, p a r M a t r i x S p a r s e <T, S> * matAop , S p r o b S i z e , S lbandwidth )

In this function, distributed matrix Am and matAop should be filled with the same way. And these entries of matrix can be filled by the method Loc SetValue implemented in parMatrixSparse. Loc SetValue uses the global indices of matrix to set values.

Chapter 5

Interface to Other Languages/Libraries

Until now, SMG2S provides interfaces to C, Python, PETSc and Trilinos.

Interface to C

SMG2S install command will generate a shared library libsmg2s.so (libsmg2s2c.dylib on OS X platform) into ${INSTALL DIRECTORY}/lib. It can be used to profit the C wrapper of SMG2S.

The way to use: smg2sRealDoubleInt (m, 1 0 , n , 3 , " " ,MPI COMM WORLD) ; 6. Release Nilpotency Object and parMatrixSparse Object : R e l e a s e N i l p o t e n c y I n t (&n ) ; R e l e a s e P a r M a t r i x S p a r s e R e a l D o u b l e I n t (&m) ; 19 SMG2S provides the C interface to different data types. For the data type of matrix size, it can be either int or longint; for the data type of matrix entries, it can be either complex or real with single or double precision.

The Nilpotent Matrix object is implemented for both int and longint as below:

struct N i l p o t e n c y I n t ; struct N i l p o t e n c y L o n g I n t ;

The interface of C for parMatrixSparse Object and smg2s function can be defined as below, SUFFIX can be replaced by one of the selected data types:

• ComplexDoubleLongInt; • ComplexDoubleInt; • ComplexSingleLongInt; • ComplexSingleInt; • RealDoubleLongInt; • RealDoubleInt; • RealSingleLongInt; • RealSignleInt. // l

Create Your Inferface

On each process, the submatrix is stored by the std::map<T,S> provided by C++, which can be gotten through the function GetDynM atLoc implemented in the sparse matrix implementation. The column index and related entry value can be gotten by the C++ iterator. Verification of Eigenvalues SMG2S provides the functionality to verify the abilibity to keep given spectrum. In the directory of ./verification/. The implementation of the functionality is powerInverse.cpp.

Prerequisites

The verification method is implemented based on the shifted inverse method proposed by PETSc/SLEPc. Before the starting of verification, it is necessary to have the two on the platforms.

If not, the download and installation of PETSc can be found: LENGTH=$ (awk 'NR==2{ p r i n t $1 } ' v e c t o r . t x t ) echo " Test E i g e n v a l u e s number = " ${LENGTH} f o r ( ( i =3; i<=$ {LENGTH}+2; i ++)) do r e a l=$ (awk 'NR=='$ i ' { p r i n t $2 } ' v e c t o r . t x t ) imag=$ (awk 'NR=='$ i ' { p r i n t $3 } ' v e c t o r . t x t ) s r u n -n 1 $ {EXEC} -n $ {N} -l $ {L} -e p s m o n i t o r c o n v \ -e p s p o w e r s h i f t t y p e c o n s t a n t -s t t y p e s i n v e r t \ -e x a c t v a l u e ${ r e a l }+$ { imag } i -t e s t t o l $ {TEST TOL} \ -d e g r e e $ {DEGREE} done Here we list the meaning of the critical parameters in the script above:

• N: the size of matrix to generate, which should be equal to the number of given eigenvalues;

• L: the bandwidth of low part diagonal of matrix to generate;

• TEST TOL: the tolerance to check if the accuracy of one eigenvalue can be accepted or not;

• DEGREE: the continuous one for the nilpotency matrix.

Script for result cleaning

The result file generated during the verification can be cleaned into the pseudo-Matrix Market Vector by the script below:

#! / b i n / b a s h g r e p "@> The e i g e n v a l u e " $1 > tmp . t x t awk ' { p r i n t $5 " " $7 } ' tmp . t x t > tmp2 . t x t awk ' { p r i n t s u b s t r ( $0 , 1 , l e n g t h ( $0 ) -1)} ' tmp2 . . / t r a i t e m e n t . sh r e s u l t s . t x t r e s u l t s c l e a n . t x t

In this command, the 1st and 2nd arguments for the execution are separately the initial results file and the final cleaned and formatted file.

  cmake . -DCMAKE INSTALL PREFIX=${INSTALL DIRECTORY} The main.cpp will generate an exectutable smg2s.exe to demonstrate a minimum sample : make The main part of SMG2S is a collection of header files. Install the header files into ${INSTALL DIRECTORY} make i n s t a l l For testing the software in your platforms: make t e s t The output of the test should be like: Running t e s t s . . . Test p r o j e c t / User /name/SMG2S S t a r t 1 : T e s t S i z e 1 0 0 0 0 w p r o c 1 1/4 Test #1: T e s t S i z e 1 0 0 0 0 w p r o c 1 . . Passed 1 . 2 0 s e c S t a r t 2 : T e s t S i z e 2 0 0 0 0 w p r o c 2 2/4 Test #2: T e s t S i z e 2 0 0 0 0 w p r o c 2 . . Passed 1 . 2 2 s e c S t a r t 3 : T e s t S i z e 1 0 0 0 0 s p r o c 1 3/4 Test #3: T e s t S i z e 1 0 0 0 0 s p r o c 1 . . Passed 1 . 2 0 s e c S t a r t 4 : T e s t S i z e 1 0 0 0 0 s p r o c 2 4/4 Test #4: T e s t S i z e 1 0 0 0 0 s p r o c 2 . . Passed 0 . 6 6 s e c 100% t e s t s passed , 0 t e s t s f a i l e d out o f 4 T o t a l Test time ( r e a l ) = 4 . 2 9 s e c

S

  Loc2Glob ( S l o c a l i n d e x ) ; // c o n v e r t l o c a l i n d e x t o g l o b a l i n d e x ; S Glob2Loc ( S g l o b a l i n d e x ) ; // c o n v e r t l o c a l i n d e x t o g l o b a l i n d e x ; int GetRank ( ) ; // Get pro c i n d e x ; S GetLowerBound ( ) ; // Get l o w e r bound i n d e x f o r each proc ; S GetUpperBound ( ) ; // Get upper bound i n d e x f o r each pro c ; S G e t L o c a l S i z e ( ) ; 9 // Get i n d e x number f o r each pr oc ;

/

  * c o n s t r u c t o r * / p a r V e c t o r (MPI Comm ncomm , S lbound , S ubound ) ; Here we give some methods of parV ector: parVectorMap<S> * GetVecMap ( ) ; // r e t u r n t h e r e l a t e d parVectorMap o f p a r V e c t o r ; S Loc2Glob ( S l o c a l i n d e x ) ; // c o n v e r t l o c a l i n d e x t o g l o b a l i n d e x ; S Glob2Loc ( S g l o b a l i n d e x ) ; // c o n v e r t l o c a l i n d e x t o g l o b a l i n d e x ; int GetRank ( ) ; // Get pro c i n d e x ; S GetLowerBound ( ) ; // Get l o w e r bound i n d e x f o r each proc ; S GetUpperBound ( ) ; // Get upper bound i n d e x f o r each pro c ; S G e t L o c a l S i z e ( ) ; // Get i n d e x number f o r each pr oc ;

void

  VecAdd ( p a r V e c t o r * v ) ; //Add a n o t h e r v e c t o r v w i t h same mapping ; void V e c S c a l e (T s c a l e ) ; // S c a l i n g t h e v e c t o r ; T VecDot ( p a r V e c t o r * v ) ; // v e c t o r d o t p r o d u c t o p e r a t i o n ; void ReadExtVec ( s t d : : s t r i n g spectrum ) ; // read v e c t o r from l o c a l f i l e ; void VecView ( ) ; // d i s p l a y t h e v e c t o r ; void specGen ( s t d : : s t r i n g spectrum ) ; // Generate / l o a l a s p e c i a l v e c t o r c o n t a i n i n g g i v e n spectrum Here we give an example to generate a parV ector by SMG2S: int w o r l d s i z e ; int w o r l d r a n k ; int span , l o w e r b , upper b ; MPI Comm size (comm, &w o r l d s i z e ) ; MPI Comm rank (comm, &w o r l d r a n k ) ; span = int ( c e i l ( double ( p r o b S i z e ) / double ( w o r l d s i z e ) ) ) ; i f ( w o r l d r a n k == w o r l d s i z e -1 ) { l o w e r b = w o r l d r a n k * span ; upper b = p r o b S i z e -1 + 1 ; } e l s e { l o w e r b = w o r l d r a n k * span ; upper b = ( w o r l d r a n k + 1 ) * span -1 + 1 ; } parVector<T, S> * vec = new parVector<T, S>( \ comm, l o w e r b , upper b ) ; vec->SetTovalue (T v a l ) ;2.1.3 Parallel MatrixThe one-dimensional row-major parallel matrix parM atrixSparse in SMG2S is distributed with the same parV ectorM ap of a given parV ector. On each process, the columns indices and the entries values are stored by a std :: map < T, S > with S data type of indices and T data type of entries. This is the constructor of parM atrixSparse:/ * c o n s t r u c t o r * / p a r M a t r i x S p a r s e ( parVector<T, S> * Vec , parVector<T, S> * Vec ) ;Here we given some methods of the parM atrixSparse object: parVectorMap<S> * GetYMap ( ) ; // r e t u r n t h e m a t r i x mapping o f t h e columns ;MPI Comm GetComm ( ) ; // Get c u r r e n t w o r k i n g MPI communicator ; s t d : : map<S , T> * GetDynMatLoc ( ) ; // Get t h e map s t o r i n g c o l s and e n t r i e s on each pr oc ; void LOC MatView ( ) ; // d i s p l a y t h e p a r a l l e l m a t r i x void L o c S e t V a l u e L o c a l ( S row , S c o l , T v a l u e ) ; // i n s e r t v a l u e i n t h e l o c a l i n d e x row ; void L o c S e t V a l u e s L o c a l ( S nindex , S * rows , S * c o l s , \ T * v a l u e s ) ; // i n s e r t a r r a y i n t h e l o c a l i n d e x rows ; void SetValueGocal ( S row , T v a l u e ) ; // i n s e r t v a l u e i n t h e g o c a l i n d e x named row ; void S e t V a l u e s G o c a l ( S nindex , S * rows , T * v a l u e s ) ; // i n s e r t a r r a y i n t h e g o c a l i n d i c e s named rows ; void L o c S e t V a l u e ( S row , S c o l , T v a l u e ) ; // S e t t h e e n t r i y ( row , c o l ) o f p a r M a t r i x w i t h v a l u e ; void L o c S e t D i a g o n a l ( parVector<T, S> * d i a g ) ; // s e t t h e d i a g o n a l o f m a t r i x t o a g i v e n v e c t o r ; void Loc MatScale (T s c a l e ) ; // S c a l i n g t h e m a t r i x ; void Loc MatAXPY ( p a r M a t r i x S p a r s e <T, S> * X, T s c a l e ) ; //AXPY o p e r a t i o n ; void Loc MatAYPX ( p a r M a t r i x S p a r s e <T, S> * X, T s c a l e ) ; //AYPX o p e r a t i o n ; void L o c Z e r o E n t r i e s ( ) ; // Zeros a l l e n t r i e s and k e e p t h e p r e v i o u s m a t r i x p a t t e r n ;

  * l o w e r b and u p p e r b o f each proc i s g i v e n * / parVector<T, S> * vec = new parVector<T, S>( \ comm, l o w e r b , upper b ) ;p a r M a t r i x S p a r s e <T, S> * A=new p a r M a t r i x S p a r s e <T, S>(vec , vec ) ;

Figure 3 . 1 :

 31 Figure 3.1: Nilpotent Matrix

1 .

 1 Include the head file #i n c l u d e <smg2s / smg2s . h> 16 template<typename T, typename S> void parVector<T, S > : : specGen ( s t d : : s t r i n g spectrum ) In this function, the eigenvalues are stored by the distributed vector text-colorblueparVector. And the filling of values on this parVector can be done by the method SetValueGlobal implemented in parVector, which takes the global indices to set values.

1 .

 1 Add this shared library to LD LIBRARY PATH: export LD LIBRARY PATH=${INSTALL DIRECTORY}/ l i b 2. Include the header file: #i n c l u d e < i n t e r f a c e /C/ c w r a p p e r . h> 3. create Nilpotency object : struct N i l p o t e n c y I n t * n ; n = n e w N i l p o t e n c y I n t ( ) ; NilpType1 ( n , 2 , 1 0 ) ; 4. After that, you need to create the parallel Sparse Matrix Object Mt like this : struct p a r M a t r i x S p a r s e R e a l D o u b l e I n t * m; m = newParMatrixSparseRealDoubleInt ( ) ; 5. Generate by SMG2S :

  p a r M a t r i x S p a r s e <S , T > * M T c o l ; S v a l ; / * On each p r o c e s s * / parVectorMap<T> * pmap = M->GetYMap ( ) ; / * Get row number on each pro c * / T l s i z e = pmap->G e t L o c a l S i z e ( ) ; s t d : : map<T, S> * d y n l o c ; s t d : : map<T, S > : : i t e r a t o r i t ; d y n l o c = M ->GetDynMatLoc ( ) ; / * Get c o l i n d i c e s and v a l u e s * / f o r (T i = 0 ; i < l s i z e ; i ++){ f o r ( i t = d y n l o c [ i ] . b e g i n ( ) ; i t != d y n l o c [ i ] . end ( ) ; ++i t ) { c o l = i t -> f i r s t ; v a l = i t ->s e c o n d ; } } Chapter 6

  [PETSc Download] and [SLPEc Installation]. The download and installation of PETSc can be found: [SLEPc Download] and [SLEPc Installation].

6. 2

 2 Verifcation by Shifted Inverse Power Method 1. compile the file powerInverse.cpp by the command make This will generate an executable powerInverse.exe. 2. Suppose the given eigenvalues are stored in the file vector.txt by the pseudo-Matrix Matrix Vector format, run the verification script as below: #! / b i n / b a s h EXEC=./ p o w e r I n v e r s e . exe N=100 L=10 TEST TOL=0.00001 DEGREE=4

Figure 6 . 1 :

 61 Figure 6.1: Home Screen Capture

  SMG2S is an open source software published under the GNU Lesser General Public License v3.0. SMG2S can be redistributed and modified under the terms of this license.SMG2S is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. SMG2S is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MER-CHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. You should have received a copy of the GNU Lesser General Public License along with SMG2S. If not, see ¡http://www.gnu.org/licenses/¿.

	1.3 CMake Options
	We use CMake to build, test and package SMG2S. If you do not have PETSc
	and SLEPc in your platform, please make sure the option below is OFF in
	CMakeLists.txt.
	o p t i o n (INSTALL TO USE " I n s t a l l SMG2S i n c l u d e f i l e s ? " OFF)
	1.4 Copyright and Licensing of SMG2S

  o n g i n t c a s e

	5.2 Interface to Python
	SMG2S uses SWIG to generate the wrapper of SMG2S to Python. Generate
	the shared library and install the python module of smg2s.
	#i n s t a l l o n l i n e from p y p i
	CC=mpicxx p i p i n s t a l l smg2s
	#b u l i d i n l o c a l
	cd . / i n t e r f a c e / python ;
	CC=mpicxx python s e t u p . py b u i l d e x t --i n p l a c e
	#or
	CC=mpicxx python s e t u p . py b u i l d
	#or
	CC=mpicxx python s e t u p . py i n s t a l l
	#run
	mpirun -np 2 python g e n e r a t e . py
	Before the utilization, make sure that mpi4py is installed.
	This is a little example of usage :
	struct parMatrixSparseSUFFIX ;
	/ * parVectorMap C wrapper * / from mpi4py import MPI
	import smg2s struct parVectorMapLongInt * newparVectorMapLongInt ( void ) ; import s y s
	/ * complex d o u b l e l o n g i n t * /
		s i z e = MPI .COMMWORLD. G e t s i z e ( )
	struct parMatrixSparseComplexSUFFIX * newPar\ rank = MPI .COMMWORLD. Get rank ( ) MatrixSparseSUFFIX ( void ) ; name = MPI . G e t p r o c e s s o r n a m e ( )
	void ReleaseParMatrixSparseSUFFIX ( struct \
	parMatrixSparseSUFFIX * * p p I n s t a n c e ) ; s y s . s t d o u t . w r i t e (
	" H e l l o , World ! I am p r o c e s s %d o f %d on %s . \ n" void LOC MatViewSUFFIX ( struct p a r M a t r i x \ SparseSUFFIX * m) ; % ( rank , s i z e , name ) )
	void GetLocalSizeSUFFIX ( struct p a r M a t r i x \ i f rank == 0 :
	SparseSUFFIX * m, print ( 'INFO ]> S t a r t i n g . . . ' ) i n t 6 4 t * r s , i n t 6 4 t * c s ) ;
	void Loc ConvertToCSRSUFFIX ( struct p a r M a t r i x \ print ( "INFO ]> The MPI World S i z e i s %d" %s i z e )
	SparseComplexDoubleLongInt * m) ;
	#b a n d w i d t h f o r t h e l o w e r band o f i n i t i a l m a t r i x void Loc CSRGetRowsArraySizesSUFFIX ( struct p a r M a t r i \ xSparseSUFFIX * m, i n t 6 4 t * s i z e , i n t 6 4 t * s i z e 2 ) ; lbandwidth = 3
	void Loc CSRGetRowsArraysSUFFIX ( struct par \ #c r e a t e t h e n i l p o t e n t m a t r i x
	MatrixSparseSUFFIX * m, i n t 6 4 t s i z e 2 , i n t * * c o l s , double * * r e a l , double * * imag ) ; i n t 6 4 t s i z e , i n t * * rows , \ n i l p = smg2s . N i l p o t e n c y I n t ( )
	void smg2sSUFFIX ( struct parMatrixSparseSUFFIX \ #s e t u p t h e n i l p o t e n t m a t r i x :
	* m,	i n t 6 4 t p r o b S i z e , struct N i l p o t e n c y L o n g I n t \ n i l p . NilpType1 ( 2 , 1 0 )
	* n i l p ,	i n t 6 4 t lbandwidth , char * spectrum , MPI Comm comm ) ;
	Mt = smg2s . p a r M a t r i x S p a r s e D o u b l e I n t ( )
	#Generate Mt by SMG2S

#v e c t o r . t x t i s t h e f i l e t h a t s t o r e s t h e g i v e n \

Here, the probsize parameter represent the matrix size, nilp is the nilpotency matrix object that we have declared previously, lbandwidth is the bandwidth of lower-diagonal band. spectrum is the file path of spectra file, if spectrum is set as " ", SMG2S will use the mechanism inside to generate the spectral distribution. comm is the basic object used by MPI to determine which processes are involved in a communication.

The given spectra file is in pseudo-Matrix Market Vector format. For the complex eigenvalues, the given spectrum is stored in three columns, the first column is the coordinates, the second column is the real part of complex values, and the third column is the imaginary part of complex values. %%MatrixMarket matrix c o o r d i n a t e complex g e n e r a l 3 3 3 1 10 6 . 5 1 5 4 2 1 0 . 6 2 8 8 3 . 4 7 9 0 3 1 0 . 7 6 2 1 5 . 0 5 4 0 For the eigenvalues values, the given spectrum is stored in two columns, the first column is the coordinates, the second column is related values. %%MatrixMarket matrix c o o r d i n a t e complex g e n e r a l 3 3 1 10 2 1 0 . 6 2 8 8 3 1 0 . 7 6 2 1

Creation of Given Spectrum

In the directory ./verification/tests, we give an example to generate the file of given spectrum, which can be reused by the users to create their own values. This is a small C++ file named as vecgen.cpp.

If the users want to generate the eigenvalues in time without loading from local file, they can customize their eigenvalues generation by the function spec-Gen in the file ./verification/tests/specGen.h, and set the parameter spectrum of smg2s to be " ".

s p e c t r a l d i s t r i b u t i o n in l o c a l f i l e s y s t e m . Mt=smg2s . smg2sDoubleInt ( 1 0 , n i l p , lbandwidth , \ " v e c t o r . t x t " , MPI .COMMWORLD) Here are the example of Arnoldi, GMRES, and another Krylov method. You need to have Python2.X or Python3.X to run it. Moreover, UI uses some libraries to support a dynamic and intuitive graphical user interface, you can see the list of libraries. Normally, some of them are included in Python distribution. You can find below the list of necessary libraries of the UI.

Interface to PETSc

Interface to Trilinos/Teptra

• Modules which are bundled in the Python installation: Tkinter, re, sys, decimal ;

• Modules which need to be installed in addition to Python: NumPy & SciPy, Matplotlib, Pillow(PIL)

Install modules to Python 2.X: apt-g e t i n s t a l l python-tk python-imaging-tk p i p -m i n s t a l l P i l l o w python -mpip i n s t a l l -U p i p python -mpip i n s t a l l -U m a t p l o t l i b p i p i n s t a l l -U numpy s c i p y Install modules to Python 3.X: sudo apt-g e t i n s t a l l python3-tk python-imaging-tk p i p -m i n s t a l l P i l l o w python -mpip i n s t a l l -U p i p python -mpip i n s t a l l -U m a t p l o t l i b p i p i n s t a l l -U numpy s c i p y It will be generate with automatic lens scaling, but your can generate it with your own scales as Fig. 6.4.2:

How to use the GUI


Mathis Fleury

Pierre Maurel

Marsel Mano

Elise Bannier

Christian Barillot

Automatic Electrodes Detection during simultaneous EEG/fMRI acquisition

Synopsis

Simultaneous EEG/fMRI acquisition allows to measure brain activity at high spatial-temporal resolution. The localisation of EEG sources depends on several parameters including the position of the electrodes on the scalp. The position of the MR electrodes during its acquisitions is obtained with the use of the UTE sequence allowing their visualisation. The retrieval of the electrodes consists in obtaining the volume where the electrodes are located by applying a sphere detection algorithm. We detect around 90% of electrodes for each subject, and our UTE-based electrode detection showed an average position error of 3.7mm for all subjects.

Purpose

The coupling of EEG/fMRI allows to measure brain activity at high spatial-temporal resolution. The localisation of EEG sources depends on several parameters including the position of the electrodes on the scalp. An accurate knowledge about these informations is important for source reconstruction. Currently, when acquiring EEG/fMRI together, the position of the electrodes is generally estimated according to fiducial points by using a template. In the context of simultaneous EEG/fMRI acquisition, a natural idea is to use MR images to localise EEG electrodes. However MR compatible electrodes are built to be almost invisible on MRI Images. Taking advantage of a recently proposed Ultra-Short Echo Time (UTE) sequence , we propose a fully automatic method to detect those electrodes in MR images.

Materials &methods

The UTE sequence was acquired in five healthy volunteers and all underwent a simultaneous EEG/fMRI examination. A 64channels MR compatible cap (Brainproduct, Gilching, Germany) was used for all the subjects. The retrieval of the electrodes consisted in two part, in the one hand we provided a mask that consists in obtaining the volume where the electrodes are located. In the second hand we performed the electrode detection inside the volume of interest (VOI).

First, the anatomical T1 image is segmented using FSL ,and a mask of the scalp is computed from the segmentation. This mask is then dilated to isolate the volume where the electrodes are located. Figure 2 shows the result of this first step.

The second step consists in obtaining a list of potential candidates for electrodes position. This list is given by a spherical Hough algorithm that will detect spherical shapes in our VOI. Detection of spherical shapes are presented in figure 3.a.

Third, we will now filter these potential electrodes to keep only the real electrodes by applying a criterion based on their theoretical positions. To do this, we use a spherical 64-template EEG net and the Iterative Closest Point algorithm (ICP) to register the template to the detected points. For each of the 64 electrodes of the registered template, the closest point detected by the spherical Hough is selected. Unselected points are discarded. Once the filtering step has been performed the number of electrodes will be equal to 64. Detected electrodes for one subject are shown in figure 3.b.

Finally the points that are far from the registered template are discarded and replaced by the local maxima of the UTE image around the theoretical position given by the registered template. The 64 obtained detection are then labelled using the template.

Results

The performance of our method was evaluated by comparing our detections to a ground truth (GT). Following , this GT was manually obtained by picking up the position of each 64 electrodes for each subject on the pancake view, which is roughly a 2D projection of the scalp.

We also compared our method to a more classical semi automatic one : five fiducial points were obtained manually and the spherical template was adjusted to these points.

The figures 4 and 5 illustrate that we are able to detect almost 90% of electrodes for each subject. Our UTE-based electrode detection showed an average position error of 3.7mm for all subjects. Paired t-test was computed between the two groups with p <0.0001. Today, the time for a UTE sequence is almost five minutes. Lowering the sample of this sequence makes the acquisition time can going down to 1 min 30. We expect to evaluate the impact of this down-sampling in the near future.

Conclusion

We present a method to detect and label automatically EEG electrodes (64-channels) during an EEG/fMRI acquisition. We have used a UTE sequence to obtain the electrode position on a MR-volume. This method does not require additional cost and can be used by just adding a UTE sequence to the MR protocol.

For future research, this technique will be used to extract the position of the electrodes in real time. Indeed, this technique is interesting for applications requiring immediate knowledge of the position of the electrodes. We believe this method will be useful in the improvement of the fusion of EEG-fMRI data-sets.

Figures

Figure 1: Visualisation of the electrodes around the scalp, the T1 image (on the right) show that the MR compatible electrodes does not appear however they emerge on the UTE sequence.

Figure 2 :

 2 Figure 2: Extraction of the VOI on the T1 image (right) and UTE image masked with the VOI.

Figure 3 :

 3 Figure 3: Example of electrode detection before (on the right) appliance of the filtering and after filtering.

Figure 4 :

 4 Figure4: Number of false positive (FP) and false negative (FN) for every subject. The detection accuracy is the percentage of electrodes that have been detected. We consider that an electrode is detected when the position error (PE) is bellow 1cm[START_REF] Kavanagk | Evaluation of methods forthree-dimensional localization of electrical sources in the human brain[END_REF] ).

Figure 5 :

 5 Figure 5: Position error for UTE-based detection and semi-automatic detection for all subjects.

4

  https://submissions.mirasmart.com/ISMRM2018/... 3 of 3 11/8/17, 6:29 PM

Acknowledgements

This work has received a French government support granted to the CominLabs excellence laboratory and managed by the National Research Agency in the "Investing for the Future" program under reference ANR-10-LABX-07-01. It was also financed by Brittany region under HEMISFER project, MRI data acquisition was supported by the Neurinfo MRI research facility from the University of Rennes I. Neurinfo is granted by the European Union (FEDER),the French State, the Brittany Council, Rennes.


Sophie Kaleba
email: sophie.kaleba@etudiant.univ-lille1.fr

Clément Béra
email: clement.bera@vub.ac.be

Stéphane Ducasse
email: stepahne.ducasse@inria.fr

Assessing primitives performance on multi-stage execution

Keywords: Time compiler, Primitive, Virtual machine, Managed runtime

Virtual machines, besides the interpreter and just-in-time compiler optimization facilities, also include a set of primitive operations that the client language can use. Some of these are essential and cannot be performed in any other way. Others are optional: they can be expressed in the client language but are often implemented in the virtual machine to improve performance when the just-in-time compiler is unable to do so (start-up performance, speculative optimizations not implemented or not mature enough, etc.).

In a hybrid runtime, where code is executed by an interpreter and a just-in-time compiler, the implementor can choose to implement optional primitives in the client language, in the virtual machine implementation language (typically C or C++), or on top of the just-in-time compiler backend. This raises the question of the maintenance and performance trade-offs of the different alternatives. As a case study, we implemented the String comparison optional primitive in each case. This paper describes the different implementations, discusses the maintenance cost of each of them and evaluates for different string sizes the execution time in Cog, a Smalltalk virtual machine.

Introduction

High-level object-oriented programming languages are often implemented on top of a virtual machine (VM). As VMs have become more popular, different techniques have been set up using Just-In-Time (JIT) compilation to improve the overall runtime performance: method JITs [START_REF] Hölzle | Adaptive optimization for Self: reconciling high performance with exploratory programming[END_REF] which usually compile methods to native code and tracing [Copyright notice will appear here once 'preprint' option is removed.] JITs [START_REF] Bala | Dynamo: A Transparent Dynamic Optimization System[END_REF][START_REF] Bolz | Tracing the meta-level: PyPy's tracing JIT compiler[END_REF]) which usually compile linear traces of execution into native code.

On top of its virtual machine, the client language can use a set of primitives 1 , which are performed directly by the interpreter rather than by evaluating expressions in a method. We distinguish two kinds of primitives:

• Essential primitives cannot be performed in any other way. A high-level object-oriented language without primitives can move values from one variable to another, but cannot add two integers together. Many arithmetic and comparison operations between numbers are primitives. Some primitives allow one to communicate with I/O devices such as the disk, the display and the keyboard.

• Optional primitives exist only to make the system run faster [START_REF] Goldberg | Smalltalk-80: The Language and Its Implementation[END_REF]. They can be implemented in the client language directly, making its implementation as a primitive optional, but they are often implemented directly in the VM to improve performance.

Optional primitives as client code. The implementation of an optional primitive in the VM rather than the client language improves performance but also often increases the implementation engineering cost. For example, accessing objects in the VM usually requires understanding of the memory layout or implementation details, which are not needed in the client language. For this reason, some virtual machine implementors attempted to remove most or all optional primitives. Their goal was to use a JIT performing speculative optimizations to be able to implement the optional primitives in the client language without performance loss. Ideally they wanted to get the performance of the same primitives written in the VM. Self [START_REF] Hölzle | Adaptive optimization for Self: reconciling high performance with exploratory programming[END_REF] and Strongtalk (Sun Microsystems 2006) were the first to try this approach. Early versions of the Javascript engine V8 (Google 2008) also implemented most of the primitives, including for example Array operations in Javascript itself. However, without optional primitives implemented in the VM, several problems rise:

• High performance relies on a JIT with speculative optimizations, which is hard to implement, maintain and evolve.

• There is a big performance gap between the baseline performance and the peak performance.

• Even mature JITs fail to optimize some narrow cases, where the performance is drastically slower [START_REF] Barrett | Virtual machine warmup blows hot and cold[END_REF].

Overall, this approach requires high engineering time to get a mature optimizing JIT with speculative optimizations and even then, it leads in practice to unreliable performance.

World border cost. VMs are traditionally implemented in a low-level language such as C or C++. To balance between memory footprint, start-up performance and peak performance, the execution of code is usually done through multiple execution tiers: the first few executions of a code snippet are done through an interpreter and the JIT compiler optimizes at runtime the frequently used code snippets.

This leads to annoying concerns: for example, let's say we implement a primitive in C. Frequently used portion of code calling that primitive are going to be compiled to native code by the JIT. This is problematic since the native code generated by the JIT is not directly compatible with native code generated by the C compiler: a call from one to another may require to edit the stack pointer and the frame pointer from the client stack to the C stack and to spill or move multiple registers. Switching between both usually wastes up to around a dozen native instructions. If the primitive is going to execute many instructions, the switch overhead might be negligible, but if the primitive executes only a few instructions, the overhead can be noticeable.

Measurements.

In this paper, we want to measure the complexities and gain of different implementations. For this purpose, we take the example of a string comparison optional primitive. The primitive compares two strings and answers if one string is greater, equal or lower than the other string. The exact specification is detailed in Section 2.4. The primitive is convenient for measurements since it can be executed with strings of small sizes or large sizes, leading the primitive to execute many processor instructions or only a few of them.

We implemented a first version of the optional primitive in Pharo. We then implemented two versions compiling through C compilers to native code, using different parts of our VM infrastructure. Lastly, we implemented a version in the back-end of the JIT. Each time, we tried to write the code in the most efficient way possible. Then, we compared the execution time of all versions on strings of different sizes, showing that overall the most complex implementation is the fastest to execute.

Section 2 describes our implementation context, i.e., the execution model of the VM we used for our evaluation, how it executes primitives and the specification of our string comparison primitive. Section 3 details the various implementations of the primitive we used for the evaluation. In Section 4, we evaluate the performance of the different implementations with different C compilers and string sizes. We also show the native code of the performance critical part of the primitive generated by each version and discuss it. Further sections compare our work to related works, discuss future work and conclude.

Implementation context

Our evaluation is based on the Cog VM: a Smalltalk virtual machine. The Cog VM is based on Dan Ingalls' Smalltalk VM [START_REF] Ingalls | Back to the Future: The Story of Squeak, a Practical Smalltalk Written in Itself[END_REF]) and has been enhanced with a JIT [START_REF] Miranda | The cog smalltalk virtual machine writing a jit in a high-level dynamic language[END_REF]) compiling one method or closure at a time to native code. Historically, the VM implemented a 32-bits adaptation of the Smalltalk-80 specifications [START_REF] Goldberg | Smalltalk-80: The Language and Its Implementation[END_REF], but the Cog VM is now the default VM for multiple programming languages such as Pharo [START_REF] Black | Pharo by Example[END_REF], Squeak [START_REF] Black | Squeak by Example[END_REF]) and Newspeak [START_REF] Bracha | Modules As Objects in Newspeak[END_REF].

Section 2.1 describes briefly the compilation process of the Cog VM. Section 2.2 explains briefly the hybrid execution model with an interpreter and a JIT. In Section 2.3, we detail how primitives are executed. Section 2.4 discusses the string comparison primitive specifications used for our evaluation.

VM compilation

Core VM code

Slang

C

Generated execution engine code Platform code

Executable VM

Native code

Slang-to-C compilation

C compilation

Step 1

Step 2

Figure 1: Cog VM compilation

Most of the Cog VM code base is written in Slang, a subset of Smalltalk. Slang is compiled to C and then to native code through standard C compilers. The execution engine (the memory manager, the interpreter and the baseline JIT) is entirely written in Slang. The two main purposes of using Slang code over plain C are:

• to specify with annotations what function needs to be inlined or duplicated with constant parameters, and

• to be able to simulate VM execution, executing the Slang code as Smalltalk code on top of a compiled VM for debugging purposes.

The executable is generated in two steps as shown on Figure 1, similarly to the RPython toolchain [START_REF] Rigo | Pypy's approach to virtual machine construction[END_REF]. The first step is to generate the two C files representing the whole execution engine written in Slang using the Slang-to-C compiler. During the second step, a C compiler is called to compile the execution engine and the platformspecific code written directly in C to the executable VM.

Plugins. The Cog VM can be extended using plugins [START_REF] Guzdial | Squeak -Open Personal Computing and Multimedia[END_REF]. VM plugins enable the addition of new features to the VM. The main benefit of writing a plugin over extending the VM is modularity, plugins are in separate code bases and can be easily added or removed from the VM. Plugins can be compiled as internal or external:

• External plugins are compiled as a dynamic libraries distributed with the VM. They can be added or removed from a compiled VM: each dynamic library can be removed, or recompiled separately and modified.

• Internal plugins are compiled with the VM executable. They can be added or removed at VM compilation time.

There is a performance overhead for the plugins compared to normal VM code. For example, accesses to object headers are dependent on the memory representation and therefore need to be done indirectly through calls to main VM executable, since the plugin cannot know ahead of time which memory representation is used. Such calls cannot be inlined by the C compiler2 .

VM execution

The Cog VM executes the client's code using an hybrid runtime with an interpreter and a JIT. The interpreter uses a global look-up cache to speed-up virtual calls. On a cache hit, it requests the JIT to compile that method and the native code version is used instead. In most case, the JIT is used at the second activation of the method. Closures have similar heuristics.

The VM code, including the interpreter code, is compiled by a C compiler and uses the C stack at runtime. During its execution, the interpreter modifies the execution stack of the client language, which is different from its own C execution stack as shown on Figure 2. This means for example that the frame pointer register in the processor refers to a stack frame on the C stack and that a processor push instruction pushes a value on the C stack. When executing the native code generated by JIT, only the client language stack is used. This means for example that the frame pointer register in the processor refers to a stack frame on the client stack and that a processor push instruction pushes a value on the client stack. Because of this design, switching the execution between the interpreter and the native code generated by the JIT has a cost. Effectively, the VM needs to switch the frame pointer, the stack pointer and the link register (the latter only on some architectures) from one stack to the other one. It also needs to save those pointers to be able to resume execution from one stack to the other one. In addition, caller-saved registers on the C conventions or the JIT convention need to be saved since the VM cannot tell ahead what code is going to be used in the other execution model and what registers are required. We estimate the overhead to up to around a dozen of processor instructions, the exact number of numbers varies a lot (the number of registers to edit depends on the processor used and on the code executed).

Most of the primitives are implemented either in Slang or directly in C. They benefit from the C compiler optimizations to be efficient. When called from the interpreter, activating a primitive is just a C function call (the primitive function pointer is cached in the look-up cache next to the method to activate). When called from the native code generated by the JIT, activating a primitive requires to switch the execution from the client stack to the C stack, perform the primitive, and switch back to resume execution. If the primitive takes a significant amount of instructions to execute, for example, the primitive is copying 100kb from one array to another one, the stack switch dance overhead is negligible. If the primitive is very quick to perform, for example, the primitive is adding two integers which do not overflow, the overhead is massive: the execution cost of the addition goes from a few processor instructions to a few dozen instructions.

To solve this performance bottleneck, the Cog VM allows one to redefine primitives in Cog's Register Transfer Language (RTL), an abstract assembly which is compiled to native code through Cog's JIT back-end. Such primitives are then generated to native code at runtime when a method annotated with the primitive is compiled to native code. This allows the primitive to be performed in the client stack with the client calling convention, saving register moves and the stack switch dance. Primitives redefined in such way can be entirely or partially redefined. In the latter case, only the common cases are re-implemented on top of the JIT backend. If an unimplemented case happens at runtime, the generated native code calls the C primitive instead, adding overhead only in uncommon cases.

Primitive execution

In the Cog VM, primitives are always associated with compiled methods. A compiled method has information in its header to inform the virtual machine if it has a primitive operation or not. When activating a method with a primitive, the primitive function is executed before the method's bytecode. If the primitive succeeds, the primitive returns a result, as if it were the result of the virtual call. If the primitive fails, it does so without side-effects and execution continues executing the method's bytecode, as if the primitive had not been present.

To fail without side-effects a primitive must validate any parameters and any state fetched from them, before changing execution state by performing its operations. Validation involves any of testing for a specific class, testing for bit vs pointer objects, bounds checks, and recursively applying these tests to substructure of the parameters. For example, the primitive that installs a cursor examines the first parameter to check that it represents a valid cursor object, comprised of two bitmaps, one for the image and one for the shape, plus a point to specify the cursor's hotspot.

String comparison primitive

In this section we describe first the string representations then the string comparison primitive specification.

String representations. Strings are represented on top of the Cog VM in two possible forms: ByteStrings and WideStrings. ByteStrings encode each character of the string in a single byte. The encoding usually follows the Latin-1 (ISO 8859-1) standard. WideStrings encode each character in 32 bits. The encoding usually follows Unicode specifications. In this paper, we will discuss only the implementation of the string comparison primitive for ByteStrings. WideStrings use an implementation entirely implemented in the client language, and so far, it has not been reported as a performance bottleneck for any program of any users. To sort arrays of strings, strings can be compared. They can be compared in the default order (Latin-1) or using different orders, for example, the case insensitive order.

Primitive specification. The String comparison optional primitive takes two or three parameters. The two first operands are the two ByteStrings to compare and the third operand, optional, is the order table (a ByteArray of size 256 encoding the order of characters). If no order is specified, the VM assumes the byte ordering (Latin-1 order) is the order to compare against, and compares directly the bytes of the ByteStrings. If an order is specified, for example a case insensitive order (the entries for the 41 and 61 in the order byte array, respectively characters A and a in the ASCII table, both answer the same value), the primitive compares the bytes through the indirection order. Having the order optional allows to have quicker string comparison in the common case where the Latin-1 order is used. The primitive answers a negative integer, zero, or a positive integer as the first parameter is less than, equal to, or greater than the second parameter.

In the paper we will focus on the comparison of two strings with the default Latin-1 ordering and measure the performance only of this case. The goal of this work was to improve the performance of common string operations, typically string equality which internally uses the primitive, and the cases where the order is non-Latin-1 are considered uncommon.

Conceptually, the primitive with two parameters first checks if the operands are ByteStrings and fails if they are not. It then extract the ByteString sizes from their object header and computes the minimum size. Lastly, it iterates over the two ByteStrings until the minimum size is reached. If a difference is found, if answers the difference between the two different characters. If no difference is found, it answers the difference between the two string sizes.

Different primitive implementations

This section describes the different implementations of the string comparison optional primitive and then compare them.

We then explain how each of these versions is executed by the VM and the pros and cons of each and every.

Different implementations

We distinguish four implementations: the Baseline in pure Smalltalk, the Plugin version, the Slang version and the Slang+RTL version.

Baseline. To evaluate the primitive performance, we use as baseline an implementation in pure Smalltalk. There are multiple ways of writing the comparison of two strings in Smalltalk, following Smalltalk coding conventions, a standard developer would likely use high-level constructs which apply a closure on each element of the strings. We did not use high-level construct and chose to implement the Smalltalk version the most optimized way Smalltalk allows, i.e., we used the constructs to:do: and ifFalse:, which are both recognized by the Smalltalk to bytecode compiler and compiled respectively to loops and branches at the bytecode level3 . The implementation follows the new specification, hence two methods are available, baselineCompareWith: and baselineCompareWith:order:. Since the focus is only on Latin-1 ordered comparison, we show the code only of baselineCompareWith:. Plugin. Written in Slang, the VM plugin implementation is compiled to C independently to the VM and is then compiled to native code through the C compiler as an internal plugin, without enabling linking time optimizations. Accesses to the ByteString object header (to check if they are actually ByteString and not other objects as well as extracting their size), and accesses to the stack (to read the operands) are done indirectly through the main VM API since they're dependent on the memory representation of objects and stack layout used. These calls generate some overhead since they cannot be inlined. To avoid such calls in the main comparison loop, the primitive requests the main VM to provide a pointer to the first byte of each string and then assumes all bytes are then contiguous.

Slang. The Slang implementation is written inside the interpreter and has access to the full interpreter APIs. It is compiled to C and then to native code as part of the main VM. There is no overhead in accessing the object memory representation or the stack as in the Plugin version.

Slang+RTL. In addition to the previous implementation, we've re-implemented here the primitive on top of Cog's JIT back-end. The primitive is written in Cog's RTL, an abstract assembly whose instructions compile almost oneto-one to native instructions. This primitive requires a more important amount of work to be written: we had to be careful to generate efficient machine code (there is no C compiler to optimize the code) and debugging such code requires the use of a processor simulator such as Bochs. Since only performance critical implementations matter in this context, we've only re-implemented the case where there are two parameters. For the three parameter primitive, the VM falls back to the Slang version.

Comparison

In this section we compare the four different implementations according to the following criteria:

Recompilation required: Changing the implementation requires only Smalltalk bytecode recompilation (No), requires to recompile a plugin to a dynamic library if compiled as external or the recompilation of the full VM if compiled as internal (Plugin) or requires to recompile the full VM (Yes).

Language familiar to developers: Smalltalk developers are considered to be familiar with Smalltalk only, while Cog VM developers are considered being familiar with Smalltalk, Slang, C and Assembly code. Only Cog's JIT implementors are familiar with its RTL, which is a subset of the VM developers. We use All if all communities are familiar with the language, VM if all VM developers are and VM-JIT if only a subset of VM developers are.

Number of lines of code:

We measured the number of lines of code of each implementation.

Number of implementations to maintain:

In each case, the client language has to implement a fall-back code in Smalltalk to be executed if the primitive fails. This code is maintained by the client language implementors, and not by the VM implementors. Hence we do not count the fallback code in the number of implementation to maintain, nor the pure Smalltalk implementation. However, as a VM implementor, we count all the Slang implementations as well as the Cog's RTL one.

Stack switch overhead: Yes if activating the primitive from the native code generated by the JIT requires to switch between the C and the client stack, inducing execution time overhead.

Optimization potential: C cross-file inlining is not performed by default for internal plugins. External plugins don't have access to the whole VM API and have to go through an interface to do so, leading to an overhead, limiting the performance unless the string size is really high.

Control over the generated code: The implementations written in Slang are compiled by the C compilers into native code. The generated code can therefore differ according the compilers and the optimizations they performed. On the contrary, the JIT generates the same native code whichever C compiler is used.

Table 1 summarizes the main pros and cons related to each implementation: the engineering time tends to increase when more low-level languages are used. VM recompilation slows slightly the development process, but the resulting code is usually faster to run. Likewise, the use of plugins brings modularity at the cost of some performance loss.

Evaluation

The Plugin, Slang and Slang+RTL implementations of this string comparison primitive have been implemented in the Cog VM. In this section, we compare execution speed of these primitives against the other existing implementations. We measure their respective execution times for different string sizes.

Set-up

To evaluate the performance, we had to compile a special VM with all the different primitive implementations included. The production VM normally only features one string comparison primitive implementation, not three of them. We compiled one VM for Linux with GCC, and one for Mac with LLVM (exact version numbers below). We ran all benchmarks on the same run of the same VM on both platforms.

GCC-Linux. The GCC-Linux evaluation was performed on an Asus Zenbook with Ubuntu 16.04.4 LTS, a 2.2 Ghz processor Intel Core 5, 8 Gb 1600MHz DDR3 of RAM. The Linux VM was compiled using GCC version 5.4.0.

LLVM-Mac. The LLVM-Mac evaluation was performed on a MacBook pro with Mac OS 10.11.6, a 2.9 Ghz processor Intel Core 5, 8 Gb 1867MHz DDR3 of RAM. The Mac VM was compiled using Apple LLVM version 8.0.0 (clang-800.0.38).

The following evaluation was performed on 32 bits Intel VMs (x86).

Methodology

We assessed the execution time of the different implementations by running the following benchmark: Each benchmark was measured 30 times. Each time, we run a string comparison in a loop and measure the total execution time.

The number of iterations of the loop depends on the string size: the comparison is repeated between 1414 and 10000 times for a short string (5 characters or less), and between 100 and 316 times for a long string (100 characters or more). The warm-up time of the JIT is in this context not noticeable on the performance results: the primitive is compiled to native code usually at the second call as explained in Section 2.2 and the VM sampling profiler [START_REF] Kaleba | A detailed vm profiler for the cog vm[END_REF]) confirms no significant time is spent in the interpreter loop (Between 0 and 1 sample over dozens of thousands are shown in the interpreter loop). The number of iterations allows the benchmarks to run for at least 1 second in the slower implementations to avoid processor and OS noises.

The loop overhead is measured and removed from the measurements (We made sure the JIT did not remove the dead loop to measure it). The two compared strings are equal as this situation corresponds to the worst scenario for most implementations.

Results

We performed the benchmarks upon five different string sizes: 0, 1, 5, 100 and 1000 characters. In the case of small strings, most of the time is spent activating the primitive and in the code before the main comparison loop, the primitive is performed with a little overall number of native instructions. In the case of large strings, the time spent executing the primitive is dominated by the main comparison loop. Figure 3 and Table 2 show the results in terms of relative speed-up towards the baseline implementation, respectively for GCC-Linux and LLVM-Mac.

These results can be explained based on the last 3 criteria detailed in Section 3.2.

Stack switch overhead. As mentioned in Section 2.2, all implementations written in Slang come with an overhead due to the switch between the C and Smalltalk stacks. This overhead is clearly visible in the bench results in Figure 3 for small strings (size 1 or less): the Plugin implementation is between 58,7 (GCC-Linux) to 66,1% (LLVM-Mac) slower than the baseline's for empty strings.The problem is even more important for the Plugin version, since the beginning of its code, before the main comparison loop, is slower due to the calls to the main VM code. As the string size goes larger, the overhead is less significant because it is absorbed by other performance gains. Regarding this criteria, Slang+RTL version is always the fastest as it is not impacted by the switch overhead. Optimization potential. When it comes to comparing large strings (100 characters and more), the execution time is mainly spent in the comparison loop (Figure 4). The stack switch overhead is still noticeable in the Slang implementation when strings of 100 characters are compared: in GCC-Linux, Slang is 17.2 times faster than Baseline while Slang+RTL is 18,6 times faster than Baseline. Due to the beginning of its code calling the VM code, the Plugin version is still around 0.68 times slower than the Slang and Slang+RTL version for a string of size 100 (Still in GCC-Linux). For strings of size 1000, the performance of all primitives is very similar, since they have almost the same comparison loop, except for the Slang version on LLVM-Mac, as we are about to discuss in the next paragraph.

Control over generated code. One can notice a rather striking performance difference between the 2 set-ups we have. In GCC-Linux, the performance between the Slang and Slang+RTL implementations is really close for long strings (1000 characters wide), with respectively 25.24 and 26.76 relative speed-up factors for a 1000-characters string. Yet performance greatly differs when the LLVM-Mac set-up is used, with respectively 12,94 and 24,88 relative speed-up factors. This difference is mainly due to the optimizations performed or not by the GCC and Clang compilers for the Slang version. As depicted in Figure 4, GCC-Linux generates instructions that are very close to the ones generated by the JIT, hence the execution time is similar. On the contrary, LLVM-Mac generates more instructions: it fails at keeping all the operands live across the comparison loops, leading to two extra memory reads per iteration (bolded in Figure 4). This is an example of one native function in the whole VM and we do not conclude anything about the C compiler capabilites from it. VMs compiled with LLVM and GCC are around equally fast except in narrow cases like this one and LLVM is able to optimize correctly the Plugin version. LLVM likely fails at optimizing the Slang version due to the old GNU C extensions we use in the interpreter code present in the same file to force variables into specific registers, which are better supported by GCC, the only compiler we supported back in the days.

However, we can see that the Slang+RTL version has the most reliable performance since it does not depend on what optimization the C compiler is able to perform or not. This is even more relevant in less common back-ends not shown in this evaluation, such as the MIPSEL back-end we support. On such uncommon processors, C compilers are usually less clever at performing back-end specific optimizations since less engineering time has been invested into it. Cog's RTL is closer to the native instructions and can therefore often generate better code in this case.

Conclusion

In this section, we compared the performance of the different implementations against a baseline (i.e., pure Smalltalk code) implementation. The results showed the Slang+RTL implementation was always the fastest. Otherwise, most of the performance gaps between the primitive versions could be explained by 3 criteria: stack switching, optimization for large strings, and the control over generated code.

Discussion and related work

In some programming languages (Strongtalk, Self, RSqueak), the implementors tried to decrease the overall number of optional primitives by improving the client language performance, reducing the need for such primitives [START_REF] Microsystems | Strongtalk official website[END_REF][START_REF] Hölzle | Adaptive optimization for Self: reconciling high performance with exploratory programming[END_REF][START_REF] Felgentreff | Towards reducing the need for algorithmic primitives in dynamic language vms through a tracing jit[END_REF]. In other work [START_REF] Ballard | Quicktalk: A Smalltalk-80 dialect for defining primitive methods[END_REF][START_REF] Chari | Waterfall: Primitives generation on the fly[END_REF], optional primitives can be written directly in a subset of the client language and compiled at runtime to native code, hence they can be modified in the client language without recompiling the VM. Lastly, several VMs implemented their interpreter on top of an abstract assembly (Google 2008;[START_REF] Wingo | Inside Javascriptcore's Low-Level Interpreter[END_REF], allowing the interpreter to use the same register conventions as the code generated by the JIT and avoiding the stack switch cost.

Decreasing the number of primitives

The philosophy of Self and Strongtalk, both high performance VMs, were to improve the performance of the client language with a mature optimizing JIT featuring advanced optimizations such as speculative optimizations [START_REF] Microsystems | Strongtalk official website[END_REF][START_REF] Hölzle | Adaptive optimization for Self: reconciling high performance with exploratory programming[END_REF]. In this context, the need of primitives is reduced, since the performance gain from using a primitive over the client language is low or inexistant. More recently, Felgentreff et al. [START_REF] Felgentreff | Towards reducing the need for algorithmic primitives in dynamic language vms through a tracing jit[END_REF] analysed the speed difference between primitives implemented in C, RPython, and Smalltalk. They showed that with the RPython toolchain framework, the speed difference between all three implementations is not that big, hence they could choose to implement some optional primitives in the client language over C or RPython without performance drop to reduce the number of primitives to maintain. These approaches have the advantage of decreasing the primitive maintenance cost since fewer primitives have to be maintained. The performance problem is however partially solved: only the peak performance is similar with and without the primitives. Since those systems rely on JITs, startup performance is worse than peak performance, increasing the performance difference at start-up between the client code and the primitive. In addition, the performance can be unreliable with optimizing JITs. The implementors of Morphic [START_REF] Maloney | Directness and liveness in the morphic user interface construction environment[END_REF] on top of the Self VM were complaining that changing a single line of code could lead to important performance drop and they would not understand why. Those performance drops were due to the optimizing JIT taking different optimization decisions. Aside from the start-up performance and the performance unreliability problems described, these approaches also require a mature optimizing JIT, hard to implement and maintain.

Primitives written in the client language

In QuickTalk [START_REF] Ballard | Quicktalk: A Smalltalk-80 dialect for defining primitive methods[END_REF]) and Waterfall [START_REF] Chari | Waterfall: Primitives generation on the fly[END_REF], the implementors are able to write the primitives in a subset of the client language. The primitives are not part of the VM anymore, but are instead generated as part of the compiled code bytecode representation with a special encoding. They are executed by generating at runtime native code for the special primitive encoding. They can be modified without recompiling the VM and can be performed by the client language implementors instead of the VM implementors, lowering the maintenance cost for the VM implementor.

However, this design has some issues that can lead to performance drops. Many VMs support multiple processor back-ends and multiple memory representation of objects. Primitives written in such DSL are required to be independent from the processor used and the memory representation of objects. If they were written in the VM, a different version could exist for each different processor and memory representation and a different version would be picked at compile-time. Having different versions allows to refine carefully the instructions generated for performance critical primitives (Array accesses for example) and can lead to noticeable speed-ups.

No stack switch overhead

To avoid paying the stack switch cost between the client Stack and the C stack, some VMs include an interpreter written on top of an abstract assembly. It allows the interpreter to have the same register convention as the native code generated by the JIT, removing the stack switch overhead. For example, the V8 team recently wrote an interpreter named Ignition. This interpreter is compiled ahead of time through their JIT back-end to native code. With their primitives and their interpreter written in such way, the V8 runtime rarely needs to switch between the client and the C stack (no need to switch to interpret code, no need to switch for primitives). In Javascript Core, a similar approach is taken with the interpreter named LLInt [START_REF] Wingo | Inside Javascriptcore's Low-Level Interpreter[END_REF]. LLInt is written in an abstract assembly code compiled ahead of time to native code compatible with the JIT's calling conventions.

In our context, implementing a similar approach would require us to compile Slang to native code through our JIT compiler back-end or to rewrite the interpreter on top of Cog JIT's RTL. Slang has some abstractions over native code, similarly to C. Compiling Slang ahead of time to native code efficiently requires us to re-implement optimizations implemented in C compilers and we don't want that complexity. Alternatively, we could try to lower the abstractions level of Slang to express native code optimizations, but that would increase the complexity of Slang, which we don't want either. Rewriting the interpreter on top of Cog JIT's RTL is possible but requires a significant amount of work. It would be nice to do so as a future work.

Future work and conclusion

In this section we discuss two future work directions and conclude.

Future work: back-end specific. All the primitive implementations we compared are written in a processor independent way, due to the C compiler abstraction for the Slang versions, the abstract assembly abstraction for Cog's RTL version or the VM abstraction for the baseline implementation in Smalltalk. We could rewrite the primitive differently for each back-end we currently support in production (x86, x64, MIPSEL, ARMv6) and for the back-end yet to support (ARMv8, ...). It could be relevant in some cases performance-wise, for example, on Intel processors, the JVM implements string comparison using the SSE4.2 string comparison instructions when available [START_REF] Davis | How the jvm compares your strings using the craziest x86 instruction you've never heard of[END_REF]. We did not go in this direction because the maintenance cost is too high since we would have to maintain a different implementation per back-end.

Future work: adaptive optimizer integration. For the past few years, an experimental adaptive optimizer (a JIT featuring advanced optimizations such as speculative inlining) has been developed on top of the Cog VM [START_REF] Béra | Sista: Saving optimized code in snapshots for fast start-up[END_REF][START_REF] Bera | Sista: a Metacircular Architecture for Runtime Optimisation Persistence[END_REF]. The optimizer has specific rules to be able to in-Assessing primitives performance on multi-stage execution 9 2018/9/15 line primitive operations without performance loss. We need to carefully measure, in different benchmarks, the performance of the inlined string comparison primitive. We also need to discuss the maintenance cost, for example, do we need a new representation of the primitive for the optimizing JIT to inline it? The micro-benchmarks used in the paper were relevant to show the client stack to C stack switch overhead as well as baseline performance. In the optimizing JIT context, such micro-benchmarks are not really relevant anymore since the optimizing JIT usually entirely optimizes away such benchmarks and real-application performance depends on how well the primitive is inlined and specialized using information from the optimized method (typically, if one operand of the string comparison is a constant string, can the generated code benefit from this knowledge?).

Conclusion.

In this paper we compared four different implementations of the string comparison primitive on top of the Cog VM. We showed the performance difference between each implementation for different string sizes. Choosing the right implementation for one's VM is a trade-off between maintenance cost and execution time, the quickest to execute the primitive is, the harder it is to maintain.

Figure 2 :

 2 Figure 2: Client and C stacks.

  min: length2) do: [ :i | (c1 := self basicAt: i) = (c2 := aString basicAt: i) ifFalse: [ ĉ1 -c2 ] ]. ˆlength1 -length2

  1000) collect: [ :size | | iterations time overhead | collection := ByteString new: size. collection2 := ByteString new: size. size < 100 ifTrue: [iterations := (100000000 // (size*10 max: 1) sqrt) floor.] ifFalse: [iterations := (10000000 // size sqrt) floor.]. overhead := [ 1 to: iterations do: [ :i | ] ] timeToRun. time := [ 1 to: iterations do: [ :i | collection primitiveStringCompare: collection2 ] ] timeToRun. stream nextPutAll: ((timeoverhead) asString)].

Figure 3 :

 3 Figure 3: String comparison benchmarks for strings of different sizes on different architectures. Average relative speed-up regarding baseline. The higher the better.

  Figure 4: GCC-Linux and LLVM-Mac: different optimizations on the comparison loop.

Table 1 :

 1 Comparison between implementations.

	Assessing primitives performance on multi-stage execution	5

Table 2 :

 2 Linux 0.413 ± 0.0005 0.746 ± 0.0008 1.804 ± 0.0005 12.05 ± 0.0088 25.36 ± 0.1742 LLVM-Mac 0.339 ± 0.0007 0.624 ± 0.0012 1.581 ± 0.0092 11.41 ± 0.0474 23.21 ± 0.1789 Slang GCC-Linux 1.287 ± 0.0048 2.257 ± 0.0116 4.958 ± 0.0111 17.24 ± 0.0515 25.43 ± 0.1621 LLVM-Mac 0.932 ± 0.0015 1.690 ± 0.0030 3.241 ± 0.0067 10.15 ± 0.0321 12.94 ± 0.0177 Benchmark results with standard errors. Average relative speed-up regarding baseline.

	0	1	5	100	1000
	Plugin GCC-Slang+RTL GCC-Linux 3.343 ± 0.0211 5.171 ± 0.0265 8.715 ± 0.0221 17.66 ± 0.0481 26.76 ± 0.1861 LLVM-Mac 3.255 ± 0.0079 5.216 ± 0.0153 8.971 ± 0.0359 17.64 ± 0.06 24.88 ± 0.1983

In the case of specific internal plugins, C compiler linking time optimizations can be enabled to limit this overhead.

Assessing primitives performance on multi-stage execution 3 2018/9/15

Assessing primitives performance on multi-stage execution 4 2018/9/15

The bytecode compiler moves the upper limit of the loop ahead of the loop, so (length1 min: length2) is not compute at each iteration.

Assessing primitives performance on multi-stage execution 6 2018/9/15

Assessing primitives performance on multi-stage execution 7 2018/9/15

Acknowledgments

We thank Levente Uzonyi for providing the benchmark template used in this article and Eliot Miranda for his advice on the implementation.

This work was supported by Ministry of Higher Education and Research, Nord-Pas de Calais Regional Council, CPER Nord-Pas de Calais/FEDER DATA Advanced data science and technologies 2015-2020.


Y El

M Elomari
email: elomari@uca.ac.ma

J Absi

F Fakhari Tehrani
email: fateh.fakhari-tehrani@unilim.fr

Proceeding of The International Conference on Optimization and Applications (ICOA2017) Numerical Simulation of Fracture at Asphalt Mastic Materials

In this paper, numerical simulations have been conducted to investigate how damage initiates and propagates at mastic materials. Mastic is known as the matrix component of asphalt concrete. The 2D specimen digital model has been created by using a layer of mastic material which is inserted between two coarse aggregate. Cohesive elements have been inserted into mastic to simulate crack initiation and propagation. The effects of loading and stiffness modulus will be investigated. Many important conclusions will be given.

I. INTRODUCTION

Asphalt concretes (AC) are complex biphasic materials composed of asphalt binder, fine and coarse aggregate. Fine aggregate are mixed with asphalt binder to formulate the mastic [START_REF] Fakhari Tehrani | Twodimensional/three-dimensional biphasic modelling of the dynamic modulus of bituminous materials[END_REF], [START_REF] Fakhari Tehrani | Heterogeneous numerical modeling of asphalt concrete through use of a biphasic approach: Porous matrix/inclusions[END_REF], [START_REF] Fakhari Tehrani | Investigation into the impact of the use of 2D/3D digital models on the numerical calculation of the bituminous composites complex modulus[END_REF]. In general, fracture occurs mostly at mastic and at the interface between coarse aggregate and mastic [START_REF] Song | A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material[END_REF], [START_REF] Yin | Multiscale fracture simulation of three-point bending asphalt mixture beam considering material heterogeneity[END_REF]. Recently, many researchers have been concentrating on the simulation of damage and fracture within AC. Song et al [START_REF] Song | A bilinear cohesive zone model tailored for fracture of asphalt concrete considering viscoelastic bulk material[END_REF] and Yin et al [START_REF] Yin | Multiscale fracture simulation of three-point bending asphalt mixture beam considering material heterogeneity[END_REF] investigated about fracture by simulation of the three point bending beam, in this studies, cohesive element approach was used.

In this paper, 2D fracture simulation at mastic will be modelled. The RVE is represented as a layer of mastic which is placed between two coarse aggregate. Crack initiation and propagation within mastic will be shown. The effects of the loading rate and stiffness modulus will be investigated.

The objectives of this work are :

• Numerical simulations of damage initiation and fracture within mastic. • Evaluating the impact of loading on the crack propagation.

• Identifying influence of stiffness modulus on damage initiation and evolution.

A. Numerical model description 1) Cohesive element approach: Cohesive element approach is a numerical method that can be used to simulate cohesive damage at mastic. Cohesive element behaviour is shown in Figure 1, cohesive element behaves elastically until the maximum stress criterion is reached ( Initiation of damage).

The elastic behaviour can be written as described in equation ( 1) [START_REF]Defining the constitutive response of cohesive elements using a traction-separation description[END_REF].

  τ n τ s τ t   =   K nn 0 0 0 K ss 0 0 0 K tt   *   δ n δ s δ t   (1)
Where τ n , τ s and τ t are the normal and two shear traction direction stresses respectively. K nn is the stiffness in the normal compressive direction, K ss is the stiffness in the first shear direction and K tt is the stiffness in the second shear direction.

The strain components are denoted as ε n ,ε s and ε t . The separation values can be defined as [START_REF]Defining the constitutive response of cohesive elements using a traction-separation description[END_REF]:

ε n = δ n T , ε n = δ s T , ε n = δ t T ( 2
)
Where T is the initial thickness of the interface element. δ n , δ s and δ t are the separation values toward each direction.

The criterion of damage is given in equation 3 :

M ax{ τ n τ 0 n ; τ s τ 0 s ; τ t τ 0 t } = 1 (3)
Where τ 0 n τ 0 s and τ 0 t are the strength toward normal, first and second shear directions respectively.

After damage initiates, the stiffness of mastic will be degraded. Stiffness degradation ratio is evaluated by a scalar variable named D. D varies from 0 to 1. When D equals 1, the stiffness of mastic is totally degraded. The stress components at the stiffness degradation are given in equations 4,5 and 6 as :

τ n = (1 -D)τ n (4) τ s = (1 -D)τ s ( 5
)
τ t = (1 -D)τ t ( 6
)
Where τ n τ s and τ t are the stress components predicted by the elastic traction-separation behaviour for the current strains without damage [START_REF]Defining the constitutive response of cohesive elements using a traction-separation description[END_REF].

After degradation of the stiffness of mastic, fracture will occur. The fracture energy can be expressed as below in equation 7 :

φ = 1 2 τ 0 (n,s,t) δ f (n,s,t) (7)
Where δ f (n,s,t) is the separation value at fracture.

2) Cohesive elements insertion: Cohesive elements are inserted between the original solid elements of mastic. The method of insertion can be found in [START_REF] Tarasovs | [END_REF]. Figure 2 illustrates the cohesive element insertion steps. The first step is to mesh the layer of mastic on ABAQUS, triangle finite elements for 2D specimen is used (figure 2a). The second step is to modify the initial mesh, every common node between two or more solid elements is duplicated. Moreover, cohesive elements with zero thickness are inserted between faces of neighbouring solid elements and joining them together (figure 2b).

3) Numerical model: The 2D numerical model chosen for this paper is represented in Figure 3. The model is composed of two coarse aggregate and a layer of mastic which is inserted between them. The 2D specimen with dimensions of 40 mm × 25 mm. The dimensions of each coarse aggregate is 40 mm × 10 mm. The mastic is with dimensions of 40 mm × 5 mm.

In general, aggregate behaves elastically. While mastic is thermal-sensitive material, it is viscous at higher temperature, visco-elastic at medium temperature and elastic at lower temperature respectively [START_REF] Fakhari Tehrani | Twodimensional/three-dimensional biphasic modelling of the dynamic modulus of bituminous materials[END_REF], [START_REF] Fakhari Tehrani | Heterogeneous numerical modeling of asphalt concrete through use of a biphasic approach: Porous matrix/inclusions[END_REF]. Fracture tests of asphalt concrete were conducted at -10 • C by H. Kim et al [START_REF] Kim | Discrete fracture modeling of asphalt concrete[END_REF], the test results are given in table I. Abaqus software is used to carry out the finite element analysis. The quasi-static simulations are considered with a step time of 1 s. The displacement is applied at the top surface of the upper aggregate, it increases linearly from 0 to 0.5 mm The free meshing method is considered. An overall mesh

B. Results and discussion

The maximum nominal stress criterion (Equation 3) is evaluated by MAXSCRT (Output of ABAQUS varies from 0 to 1). If MAXSCRT is inferior than 1, damage will not initiate, whereas damage is assumed to begin when MAXSCRT achieves 1.

The scalar damage variable D is evaluated by SDEG (Output of ABAQUS varies from 0 to 1). For example, if SDEG=0.2, 20% of the stiffness of mastic will degrade. Fracture is assumed to begin when SDEG achieves 1.

1) Effect of loading: In this section, the effect of loading on fracture will be shown. The applied displacement increases linearly from 0 to 0.5 mm during the step time of 1 s. Cracks initiation and propagation are observed in 2D model of AC specimen with upon increment of loading. The detailed crack initiation and evolution are presented in figure 5, in this figure, only mastic bloc is shown without the aggregate and the mastic size scale is increased to visualize the cracking paths. Figs 5a, 5b, 5c and 5d show the status of the 2D digital specimen at step time of 0, 0.2, 0.3 and 1 s respectively. It can be seen that two cracks are initiated and propagated as the loading increases. The first crack is initiated at the top-left region of mastic, it propagates horizontally from the left side to the right one of specimen. However, the second crack is created at the bottom-right region of mastic and propagates horizontally from the right side to the left one for a distance almost equals 10 mm.

MAXSCRT (Damage criterion) and SDEG (Stiffness degradation) for two cohesive elements are plotted in figure 6 . On the one hand, the first cohesive element is where the first crack tip appeared (top-left region of mastic). On the other hand, the second cohesive element is where the crack tip of the first crack ended (top-right region of mastic). For the first cohesive element, damage initiates at t=0.01 s In the meantime, SDEG value exceeds 0 and achieves 1 at the time t=0.1 s. For the second cohesive element, damage starts at t=0.17 s and SDEG value exceeds 0 and achieves 0.9997 at t=1 s. The second cohesive element has not been fractured yet (SDEG is still less than 1). Fracture might occur if the loading rate exceeds 0.5mm. 2) Effect of stiffness modulus: One of the important parameter affecting the fracture behaviour is the stiffness modulus of mastic. The stiffness modulus K used in the previous case is 18.2 GP a/mm, to investigate how K affects the overall fracture behaviour, 50 % of K, K and 150% of K are assigned to the mastic material at the same loading condition (0.5 mm). Figure 7 illustrates the crack paths at mastic. Figures 7a and7b show the crack paths of stiffness of 9100 MPa (50 % of K) and 27300 MPa (150 % of K) respectively. For stiffness of 9100 MPa, it reveals that only one crack appears and the mastic is divided to two parts. For stiffness of 27300 MPa, it reveals that the crack paths are almost identical to the crack path found in the case of stiffness equals K. The difference between the two crack paths is that in the case of the second crack (crack at the bottom-right region), there is one more cohesive element which is separated in the case of stiffness which equals 27300 MPa (see figs 7b and 5d).

II. CONCLUSION

The present work presents numerical simulations of fracture at asphalt materials using cohesive element method. Cohesive elements are inserted into mastic to inquire into cracking initiation and propagation. Based on the above results, the following conclusions are drown :

• Cracks were initiated and propagated at mastic with upon increment of loading. • The stiffness modulus value can influence the cracking paths at mastic.

Fig. 1 :

 1 Fig. 1: Cohesive element behaviour

Fig. 2 :

 2 Fig.2: Original mesh (a) and embedded cohesive element (b)[START_REF] Tarasovs | [END_REF]

Fig. 3 :

 3 Fig. 3: Numerical model, loading and boundary condition

Fig. 4 :

 4 Fig. 4: Mesh configuration of the numerical model

  (a) The status of mastic at 0 s (b) The status of mastic at 0.2 s (c) The status of mastic at 0.3 s (d) The status of mastic at 1 s

Fig. 5 :

 5 Fig. 5: The status of cracks initiation and propagation within mastic at different step times

  MAXSCRT and SDEG) MAXSCRT of the first cohesive element SDED of the first cohesive element MAXSCRT of the second cohesive element SDEG of the second cohesive element

Fig. 6 :

 6 Fig. 6: MAXSCRT and SDEG variation

TABLE I :

 I The input parameters of the numerical simulation

	Aggregate	Mastic
	Elastic modulus	56.8 GP a	Initial stiffness	18.2 GP a/mm
			Fracture energy	270 J/m 2
	Poisson's ratio	0.15	Tensile strength	3.78 MPa