matrix_utils.h

//------------------------------------------------------------------------------
//  
//      MATRIX UTILS
//      Utils to work with Eigen::MatrixXd
//
//------------------------------------------------------------------------------
 
 
#ifndef MATRIX_UTILS_H
#define MATRIX_UTILS_H
 
//Eigen libs
#include"Eigen/Dense"
#include"Eigen/Core"
 
//OpenMP libs
#include"omp.h"
 
//------------------------------------------------------------------------------
//  Functions to work with reduced form of the jacobians.
//  The jacobians of the layers map are sparse matrices
//  with the same number of (consecutive) non-zero entries 
//  in each row (See the documentation).
//  The following functions allow to convert 
//  a standard matrix from/to reduced form and implement 
//  mutliplications between different forms of matrices
//------------------------------------------------------------------------------
 
//Multiplication between a reduced matrix and a standard one
inline Eigen::MatrixXd reduced_standard_mul(Eigen::MatrixXd &red, Eigen::MatrixXd &std)
{
    int rows_red = red.rows();
    int non_null_per_row = rows_red*(red.cols());
    int rows_std = std.rows();
    int cols_std = std.cols();
    Eigen::MatrixXd result(rows_red,cols_std);
    #pragma omp parallel for collapse(2) shared(red,std,result,rows_red,non_null_per_row,rows_std,cols_std)
    for (int i = 0; i < rows_red; i++)
    {
        for (int j = 0; j < cols_std; j++)
        {
            #pragma omp simd
            for (int k = 0; k < rows_std; k++)
            {
                result(i,j) = red(i,k)*std(k*non_null_per_row,j);
            }
        }
    }
    return result;
}
//------------------------------------------------------------------------------
//Multiplication between a standard matrix and a reduced one
inline Eigen::MatrixXd standard_reduced_mul(Eigen::MatrixXd &std, Eigen::MatrixXd &red)
{
    int rows_std = std.rows();
    int cols_std = std.cols();
    int rows_red = red.rows();
    int cols_red = red.cols();
    int non_null_per_row = red.cols();
    Eigen::MatrixXd result(rows_std,cols_red*rows_red);
    result = move(Eigen::MatrixXd::Zero(rows_std,cols_red*rows_red));
    for (int i = 0; i < rows_std; i++)
    {
        for (int j = 0; j < cols_std; j++)
        {
            #pragma omp simd
            for (int k = 0; k < non_null_per_row; k++)
            {
                result(i,j*non_null_per_row+k) = std(i,j)*red(j,k);
            }
        }
    }
    return result;
}
//------------------------------------------------------------------------------
//Conversion from reduced form to standard form
inline Eigen::MatrixXd reduced_to_standard(Eigen::MatrixXd &a)
{
    
    int rows_a = a.rows();
    int cols_a = a.cols();
    int a_non_null_per_row = cols_a;
    int cols_res = rows_a*a_non_null_per_row;
    Eigen::MatrixXd result;
    result = move(Eigen::MatrixXd::Zero(rows_a,cols_res));
    #pragma omp parallel for shared(a,result,rows_a,cols_a,a_non_null_per_row) 
    for (int i = 0; i < rows_a; i++)
    {
        #pragma omp simd
        for (int j = 0; j < cols_a; j++)
        {
            result(i,i*a_non_null_per_row+j) = a(i,j);
        }
        
    }
    //cout << "OK!" << endl;
    return result;
}
//------------------------------------------------------------------------------
//Conversion from standard form to reduced form
inline Eigen::MatrixXd standard_to_reduced(Eigen::MatrixXd &a)
{
    int rows_a = a.rows();
    int cols_a = a.cols();
    int a_non_null_per_row = cols_a/rows_a;
    int cols_res = a_non_null_per_row;
    Eigen::MatrixXd result(rows_a,cols_res);
    #pragma omp parallel for shared(a,result,rows_a,cols_a,a_non_null_per_row,cols_res)
    for (int i = 0; i < rows_a; i++)
    {
        #pragma omp simd
        for (int j = 0; j < cols_res; j++)
        {
            result(i,j) = a(i,i*a_non_null_per_row+j);
        }
    }
    return result;
}
//------------------------------------------------------------------------------
 
 
//------------------------------------------------------------------------------
//      Operations on generic matrices
//------------------------------------------------------------------------------
 
//Compute the Frobenius norm of a matrix
inline double Frobenius_norm(Eigen::MatrixXd &mat)
{
    const int rows = mat.rows();
    const int cols = mat.cols();
    double result = 0.;
    #pragma omp parallel for collapse(2) shared(mat) reduction(+: result)
    for (int i = 0; i < rows; i++)
    {
        for (int j = 0; j < cols; j++)
        {
 
            result += mat(i,j)*mat(i,j);
        }
    }
    result = sqrt(result);
    return result;
}
//------------------------------------------------------------------------------
 
//------------------------------------------------------------------------------
//      Write a matrix to file
//------------------------------------------------------------------------------
 
//Write a matrix to file
inline void write_eigen_matrixxf_to_file(string filename, char separating_char, Eigen::MatrixXd &mat)
{
    std::ofstream file;
    file.open(filename);
    const int nrows = mat.rows();
    const int ncols = mat.cols();
    for (int i = 0; i < nrows; i++)
    {
        for (int j = 0; j < ncols; j++)
        {
            file << mat(i,j);
            file << separating_char;
        }
        file << endl;
    }
    file.close();
}
#endif /* MATRIX_UTILS_H */