Commit b8980663 authored by CHAMONT David's avatar CHAMONT David
Browse files

Transformation de NatationSynchronisee en sous-groupe

parents
# Tutoriel sur Xtensor
## Recette d'exécution du comparatif GSoC, via docker
Ici, le code soumis par les différents candidats a été modifié le moins possible.
Pour se placer dans le conteneur docker, et compiler le banc d'essai :
```shell
> docker/run.sh
>> cd gsoc
>> make bench.exe
```
Pour inverser 1000 fois un tableau aléatoire de
100 matrices :
```shell
>> ./bench.exe 1000 1000
```
## Recette d'exécution du banc d'essai, via docker
Ici, les codes soumis par les différents candidats ont été modifiés
pour travailler aussi bien en float qu'en double, aussi bien en
xt::xarray qu'en xt::xtensor. Et nous avons ajouté une implémentation
à base de xtensor-blas.
Pour se placer dans le conteneur docker, et compiler le banc d'essai :
```shell
> docker/run.sh
>> cd bench
>> make bench.exe
```
Pour inverser 1000 fois un xarray aléatoire de 100 matrices en float :
```shell
>> ./bench.exe 0 1000 1000
```
Pour inverser 1000 fois un xarray aléatoire de 100 matrices en double :
```shell
>> ./bench.exe 1 1000 1000
```
Pour inverser 1000 fois un xtensor aléatoire de 100 matrices en float :
```shell
>> ./bench.exe 2 1000 1000
```
Pour inverser 1000 fois un xtensor aléatoire de 100 matrices en double :
```shell
>> ./bench.exe 3 1000 1000
```
## A propos du répertoires /tests
Ce répertoire contient des tests de très bas niveau pour essayer
xtensor, et pour mettre en évidence certains des bugs rencontrés.
## A propos du portage des codes de xt::xarray à xt::xtensor
Cela s'est avéré particulièrement compliqué pour le code d'Ibrahim :
* on ne connait pas le nombre de dimensions des différents xarray,
ce qui m'a obligé à des impressions de shape pour savoir le
nombre de dimensions des xt::xtensor.
* Plus de possibilités de faire des reshape()
* Probleme de preparation des shapes pour faire les eyes
(il semble impossible de donner n'importe quel type "shape"
au constructeur de eye).
* J'ai souvent oublié la "lazy evaluation" et essayer d'imprimer une
expression on encore évaluée. En cas d'erreur de compilation
du type "could not convert 'index' from ...", un appel
à xt::eval est souvent la solution.
* Les messages d'erreur de compil des templates sont toujours aussi
indigestes.
* Trois bugs de xtensor ont été levés à cette occasion.
## A propos du portage des codes de xt::xtensor à xt::xtensorf
1. Maintenant que la taille des MatricesSet est fixe, il faut un nouveau paramètre pour dire quel sous-ensemble
d'une MatricesSet on veut traiter, notamment pour invert_pair().
2. A chaque fois qu'une implémentation utilise des matrices de taille spéciale, il y a problème.
Pour commencer, Ayoub construit une matrice double pour ajouter une matrice Identité.
En fait, il faut traquer tous les usages de shape()...
3. Utile : la redéfinition en intene des types, et l'utilisation des noms redéfinis,
facilite l'écriture des versions spécialisées.
\ No newline at end of file
*.o
*.exe
\ No newline at end of file
//
// David Chamont Comments
// The adjonction of an identity matrix before "echelon"
// prevent the use of xtensor ?
//
/**
* @author Ayoub Chouak (@ntauth)
*/
#pragma once
#include <climits>
#include <vector>
#include <iostream>
#include <random>
#include <chrono>
#include <string>
#include "bench.hh"
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
class AyoubInverter
{
public:
// types
using RealType = Real ;
using VectorType = Vector ;
using MatrixType = Matrix ;
using MatricesSetType = MatricesSet ;
using size_type = typename MatricesSetType::shape_type::size_type ;
static void invert( MatricesSetType & mats, size_type nbmats ) {
for ( size_type i=0 ; i<nbmats ; ++i ) {
auto matview = xt::view(mats,i,xt::all(),xt::all()) ;
MatrixType mat = matview ;
invert_gauss_jordan(mat) ;
matview = mat ;
}
}
/**
*
* @tparam Ty type of the matrix elements
* @param mat0 the matrix whose row is to be swapped
* @param mat1 the matrix whose row is to be swapped
* @param row0 the index of the row of mat0 that is to be swapped
* @param row1 the index of the row of mat1 that is to be swapped
*/
static void swap_row(MatrixType & mat0, MatrixType& mat1, size_t row0, size_t row1);
/**
* @brief Reduces a matrix to the row echelon form
* @tparam Ty matrix element type
* @param mat the matrix to be reduced
*/
static void reduce_echelon(MatrixType & mat);
/**
* @brief Inverts a matrix using the Gauss-Jordan method
* @tparam Ty matrix element type
* @param mat the matrix to be reduced
*/
static void invert_gauss_jordan(MatrixType & mat);
};
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
void AyoubInverter<Real,Vector,Matrix,MatricesSet>::invert_gauss_jordan(Matrix & mat)
{
MatrixType mat_;
// Get the matrix shape
auto const & shape = mat.shape();
// Augment the matrix by adjoining the identity matrix to the right
typename MatrixType::shape_type am_shape = { shape[0], 2 * shape[0] };
mat_ = MatrixType(am_shape);
// Fill the original matrix
for (size_t i = 0; i < am_shape[0]; i++)
for (size_t j = 0; j < am_shape[0]; j++)
mat_(i, j) = mat(i, j);
// Fill the identity matrix
for (size_t i = 0; i < am_shape[0]; i++)
for (size_t j = am_shape[0]; j < am_shape[1]; j++)
mat_(i, j) = (j - am_shape[0] == i) ? 1 : 0;
// Reduce the augmented matrix to echelon form and extract the inverse
reduce_echelon(mat_);
// Copy the inverse to mat
for (size_t i = 0; i < am_shape[0]; i++)
{
for (size_t j = am_shape[0]; j < am_shape[1]; j++) {
mat(i, j - am_shape[0]) = mat_(i, j);
}
}
}
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
void AyoubInverter<Real,Vector,Matrix,MatricesSet>::reduce_echelon( Matrix & mat )
{
auto const & shape = mat.shape();
size_t last_nz_idx = std::numeric_limits<size_t>::max();
for (size_t col = 0; col < shape[1]; col++)
{
size_t nz_idx = std::numeric_limits<size_t>::max();
bool nz_found = false;
for (size_t row = col; row < shape[0] && !nz_found; row++)
{
if (mat(row, col) != 0)
{
nz_idx = row;
nz_found = true;
}
}
if (nz_found)
{
// Move the row to allow for the zero rows to cascade down
if (nz_idx != col)
swap_row(mat, mat, nz_idx, col);
RealType pivot = mat(col, col);
// Reduce the row to echelon form by dividing each element by pivot
for (size_t i = 0; i < shape[1]; i++)
{
mat(nz_idx, i) /= pivot;
}
// Apply a gauss move to zero out the elements in the rows above and below
for (size_t i = 0; i < shape[0]; i++)
{
if (i != col && mat(i, col) != 0)
{
RealType lambda = mat(i, col); // Reduction factor
for (size_t j = 0; j < shape[1]; j++)
mat(i, j) -= mat(col, j) * lambda;
}
}
}
}
}
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
void AyoubInverter<Real,Vector,Matrix,MatricesSet>::swap_row(Matrix& mat0, Matrix& mat1, size_t row0, size_t row1)
{
VectorType x_view0 = xt::view(mat0, row0);
VectorType x_view1 = xt::view(mat1, row1);
auto const& shape0 = mat0.shape();
auto const& shape1 = mat1.shape();
// Make sure the rows are not out of bound
xt::check_access(shape0, row0);
xt::check_access(shape1, row1);
// Make sure the row dimensions are compatible
if (shape0[0] != shape1[0])
throw std::string("Incompatible row dimension!");
for (size_t i = 0; i < x_view0.size(); i++) {
mat0(row0, i) = x_view1(i);
mat1(row1, i) = x_view0(i);
}
}
// Specialisation for xtensorf
template<typename Real>
class AyoubInverter<Real, xt::xtensorf<Real,xt::xshape<5>>, xt::xtensorf<Real,xt::xshape<5,5>>, xt::xtensorf<Real,xt::xshape<NB_RAND_MATRICES,5,5>>>
{
public:
// types
using RealType = Real ;
using VectorType = xt::xtensorf<Real,xt::xshape<5>> ;
using MatrixType = xt::xtensorf<Real,xt::xshape<5,5>> ;
using MatricesSetType = xt::xtensorf<Real,xt::xshape<NB_RAND_MATRICES,5,5>> ;
using size_type = typename MatricesSetType::shape_type::size_type ;
static void invert( MatricesSetType & mats, size_type nbmats ) {
for ( size_type i=0 ; i<nbmats ; ++i ) {
auto matview = xt::view(mats,i,xt::all(),xt::all()) ;
MatrixType mat = matview ;
invert_gauss_jordan(mat) ;
matview = mat ;
}
}
static void swap_row(xt::xtensorf<Real,xt::xshape<5,10>> & mat0, xt::xtensorf<Real,xt::xshape<5,10>>& mat1, size_t row0, size_t row1)
{
VectorType x_view0 = xt::view(mat0, row0);
VectorType x_view1 = xt::view(mat1, row1);
auto const& shape0 = mat0.shape();
auto const& shape1 = mat1.shape();
// Make sure the rows are not out of bound
xt::check_access(shape0, row0);
xt::check_access(shape1, row1);
for (size_t i = 0; i < x_view0.size(); i++) {
mat0(row0, i) = x_view1(i);
mat1(row1, i) = x_view0(i);
}
}
static void reduce_echelon(xt::xtensorf<Real,xt::xshape<5,10>> & mat)
{
auto const & shape = mat.shape();
size_t last_nz_idx = std::numeric_limits<size_t>::max();
for (size_t col = 0; col < shape[1]; col++)
{
size_t nz_idx = std::numeric_limits<size_t>::max();
bool nz_found = false;
for (size_t row = col; row < shape[0] && !nz_found; row++)
{
if (mat(row, col) != 0)
{
nz_idx = row;
nz_found = true;
}
}
if (nz_found)
{
// Move the row to allow for the zero rows to cascade down
if (nz_idx != col)
swap_row(mat, mat, nz_idx, col);
Real pivot = mat(col, col);
// Reduce the row to echelon form by dividing each element by pivot
for (size_t i = 0; i < shape[1]; i++)
{
mat(nz_idx, i) /= pivot;
}
// Apply a gauss move to zero out the elements in the rows above and below
for (size_t i = 0; i < shape[0]; i++)
{
if (i != col && mat(i, col) != 0)
{
Real lambda = mat(i, col); // Reduction factor
for (size_t j = 0; j < shape[1]; j++)
mat(i, j) -= mat(col, j) * lambda;
}
}
}
}
}
static void invert_gauss_jordan(MatrixType & mat)
{
xt::xtensorf<Real,xt::xshape<5,10>> mat_;
// Fill the original matrix
for (size_t i = 0; i < 5; i++)
for (size_t j = 0; j < 5; j++)
mat_(i, j) = mat(i, j);
// Fill the identity matrix
for (size_t i = 0; i < 5; i++)
for (size_t j = 5; j < 10; j++)
mat_(i, j) = (j - 5 == i) ? 1 : 0;
// Reduce the augmented matrix to echelon form and extract the inverse
reduce_echelon(mat_);
// Copy the inverse to mat
for (size_t i = 0; i < 5; i++)
{
for (size_t j = 5; j < 10; j++) {
mat(i, j - 5) = mat_(i, j);
}
}
}
};
//
// David Chamont Comments
// The use of reshape() prevent the simple replacement of xarray with xtensor.
// And the fact to set =0 for some Matrix ?
//
//
// Created by Ibrahim Radwan on 3/25/18.
//
#pragma once
#include <xtensor/xarray.hpp>
#include <xtensor/xio.hpp>
#include <xtensor/xeval.hpp>
using namespace xt;
using namespace std;
#include "bench.hh"
/**
* Notes:
* [1] This program heavily uses matrices operations to get things done, which means for a better utilization,
* you may invert matrices in a bulky fashion
* [2] Every function here works as if there're many matrices(e.g. depth)
*
*/
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
class IbrahimInverter
{
public :
// types
using RealType = Real ;
using VectorType = Vector ;
using MatrixType = Matrix ;
using MatricesSetType = MatricesSet ;
static void invert( MatricesSet & mats ) {
mats = invertMatrices(mats) ;
}
/**
* This function takes matrices and calculates their inverses
*
* @param matrices This array of shape (#matrices, dim, dim), where dim is the matrix length/width
*
* @return Array of shape (#matrices, dim, dim) contains the inverses
*/
static MatricesSet invertMatrices(const MatricesSet &);
private:
/**
* Decompose the given matrices using Crout decomposition
*
* @param matrices Array of matrices to be decomposed
* @param depth Count of the matrices
* @param dim Size of the matrix
*
* @return pair of 2 arrays {@first contains all matrices L triangles}, {@second contains all matrices U triangles}
*/
static pair<MatricesSet, MatricesSet> matricesLUDecomposition(const MatricesSet &matrices, const size_t depth, const size_t dim);
/**
* Calculate the value to be subtracted from the current L-Col
*
* @param L Matrices Lower triangle
* @param R Matrices Upper triangle
* @param LColIdx Index of the column we are calculating
* @param dim Matrix dimension
*
* @return xarray the values to be subtracted from the original matrices columns
*/
static Matrix getLColSubtractedValue(const MatricesSet &L, const MatricesSet &U, size_t LColIdx, size_t dim);
/**
* Calculate the value to be subtracted from the current L-Col
*
* @param L Matrices Lower triangle
* @param R Matrices Upper triangle
* @param URowIdx Index of the column we are calculating
* @param dim Matrix dimension
*
* @return xarray the values to be subtracted from the original matrices columns
*/
static Matrix getURowSubtractedValue(const MatricesSet &L, const MatricesSet &U, size_t URowIdx, size_t dim);
/**
* Calculate the intermediate solutions of the system
*
* @param L
* @param U
* @param depth Count of the matrices
* @param dim Size of the matrix
*
* @return xarray The intermediate solutions array
*/
static MatricesSet calculateIntermediateSolution(const MatricesSet &L, const MatricesSet &U, const size_t depth,
const size_t dim);
/**
* Calculate the final solutions of the system
*
* @param L Lower triangle
* @param U Upper triangle
* @param D Intermediate solution matrix
* @param depth Count of the matrices
* @param dim Matrix dimension
*
* @return xarray The intermediate solutions array
*/
static MatricesSet calculateFinalSolution(const MatricesSet &L, const MatricesSet &U,
const MatricesSet &D, const size_t depth, const size_t dim);
};
/**
* Calculate matrices inverses
*
* Assumptions:
* [1] The matrices are squares
* [2] The given matrices are invertible (the non invertible ones
* will yeild a nan. matrix)
*
* The general idea is to follow the following steps:
* [1] Decompose the given each matrix into its LU matrices (using Crout Decomposition)
* [2] Calculate the intermediate solution
* [3] Calculate the final solution (matrix inverse)
*
* The steps are being taken from this online article with some modifications:
* https://www.gamedev.net/articles/programming/math-and-physics/matrix-inversion-using-lu-decomposition-r3637/
*
* @param matrices This array of shape (#matrices, dim, dim), where dim is the matrix length/width
* @return xarray of shape (#matrices, dim, dim) contains the inverses
*/
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
MatricesSet IbrahimInverter<Real,Vector,Matrix,MatricesSet>::invertMatrices(const MatricesSet &matrices)
{
// Get matrices count (a.k.a depth), and each matrix dimension
size_t depth = matrices.shape()[0];
size_t dim = matrices.shape()[1];
// Decompose matrices to LU matrices
pair<MatricesSet, MatricesSet> && LU = matricesLUDecomposition(matrices, depth, dim);
// Calculate intermediate solutions
MatricesSet && D = calculateIntermediateSolution(LU.first, LU.second, depth, dim);
// Return the final matrices
return calculateFinalSolution(LU.first, LU.second, D, depth, dim);
}
/**
* Decompose the given matrices using Crout decomposition
*
*
* This function follows these steps:
* [A0] The first L-col and U-row don't need some operations, so we skip those operation in case of i == 0
* [A1] Fill the L, U matrices:
* [A11] First fill the L-col then the U-row (as U-row depends on the previous L-col)
*
* How L-cols are filled (check the document mentioned above for more info):
* [B1] Multiply the previous L-cols with the corresponding U-col (Each L-col gets multiplied by a single value from the
* corresponding U-col)
* [B2] Sum results from step 1
* [B3] Subtract the result from step 2, from the corresponding original matrix column
*
* How U-rows are filled (check the document mentioned above for more info):
* [C1] Multiply the previous R-rows with the corresponding L-row (Each U-row gets multiplied by a single value from the
* corresponding L-row)
* [C2] Sum results from step 1
* [C3] Subtract the result from step 2, from the corresponding original matrix column
* [C4] Divide the result from step 3 by L(i,i) where i = index of the current U-row we fill
*
*
* The formulas:
* l(i1) = a(i1), for i=1,2,⋯,n
* u1j = a(1j)/l(11), for j=2,3,⋯,n
*
* For j=2,3,⋯,n−1
* l(ij) = a(ij)− ∑{k=[1, j−1]} l(ik) * u(kj), for i=j,j+1,⋯,n
*
* u(jk) = a(jk) − (∑{i=[1,j−1]} l(ji)u(ik)) / l(jj), for k=j+1,j+2,⋯,n
*
*
* @param matrices Array of matrices to be decomposed
* @param depth Count of the matrices
* @param dim Size of the matrix
*
* @return pair of 2 xarrays {@first contains all matrices L triangles}, {@second contains all matrices U triangles}
*/
template<typename Real, typename Vector, typename Matrix, typename MatricesSet>
pair<MatricesSet, MatricesSet> IbrahimInverter<Real,Vector,Matrix,MatricesSet>::matricesLUDecomposition(const MatricesSet &matrices, const size_t depth, const size_t dim)
{
// Each matrix will have its own Lower triangle matrix and Upper triangle matrix
// Crout assumes initial U as an eye matrix
std::vector<std::size_t> matrices_shape { matrices.shape()[0], matrices.shape()[1], matrices.shape()[2] } ;
MatricesSet L = xt::zeros<Real>(matrices_shape);
MatricesSet U = xt::eye<Real>(matrices_shape);
// In this loop we fill L, U
// We use same loop to reduce complexity
for (size_t i = 0; i < dim; ++i)
{
// ========================================
// Fill L cols
// ========================================
// [B3] Update the L column
auto && LCol = xt::view(L, xt::all(), xt::range(i, dim), i);
Matrix LColSubtractedValue = getLColSubtractedValue(L, U, i, dim) ; // sans ca, ca plante avec les xt::xtsensor !
Matrix NewLView = xt::view(matrices, xt::all(), xt::range(i, dim), i) ;
Matrix NewLCol = NewLView - LColSubtractedValue ;
LCol = NewLCol ;
// ========================================
// Fill U rows
// ========================================
// U rows finish calculations earlier than L-cols by 1 iteration
if (dim - i - 1 <= 0) continue;
// [C4] Get the corresponding L(i,i) value
typename Matrix::shape_type shape { depth, 1 } ;
Matrix LCorrespondingDivisorValue { shape } ;
xt::view(LCorrespondingDivisorValue, xt::all(), 1)
= xt::view(L, xt::all(), i, i) ;
//Vector LCorrespondingDivisorValue= xt::view(L, xt::all(), i, i);
//LCorrespondingDivisorValue.reshape({LCorrespondingDivisorValue.shape()[0], 1}); // Vector => Matrix
// [C3,4]