Algebra.h

Go to the documentation of this file.

/*
 *  CepGen: a central exclusive processes event generator
 *  Copyright (C) 2020-2024  Laurent Forthomme
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
#ifndef CepGen_Utils_Algebra_h
#define CepGen_Utils_Algebra_h
 
#include <gsl/gsl_matrix_double.h>
 
#include <initializer_list>
#include <memory>
 
namespace cepgen {
  class Vector;
  class VectorRef;
  class Matrix {
  public:
    explicit Matrix(size_t num_rows, size_t num_cols);
    Matrix(const std::initializer_list<Vector>& vecs);
 
    Matrix(const Matrix&);             
    Matrix& operator=(const Matrix&);  
    operator Vector() const;           
 
    static Matrix zero(size_t num_rows, size_t num_cols = 0ull);
    static Matrix uniform(size_t num_rows, size_t num_cols, double value = 1.);
    static Matrix identity(size_t);         
    static Matrix diagonal(const Vector&);  
 
    size_t numColumns() const;  
    size_t numRows() const;     
 
    bool operator==(const Matrix&) const;                                           
    inline bool operator!=(const Matrix& oth) const { return !(operator==(oth)); }  
 
    Matrix subset(size_t min_y, size_t min_x, size_t max_y = 0ull, size_t max_x = 0ull) const;
 
    Matrix& operator*=(double);               
    Matrix& operator*=(const Vector&);        
    Matrix& operator*=(const Matrix&);        
    Matrix& operator/=(double);               
    Matrix operator-() const;                 
    Matrix& operator+=(const Matrix&);        
    Matrix& operator-=(const Matrix&);        
    double& operator()(size_t, size_t);       
    double operator()(size_t, size_t) const;  
 
    Vector operator%(const Vector&) const;  
 
    friend Matrix operator*(double, const Matrix&);         
    friend Matrix operator*(const Matrix&, double);         
    friend Vector operator*(const Matrix&, const Vector&);  
    friend Matrix operator*(const Matrix&, const Matrix&);  
    friend Matrix operator/(const Matrix&, double);         
    friend Matrix operator+(const Matrix&, const Matrix&);  
    friend Matrix operator-(const Matrix&, const Matrix&);  
 
    typedef std::pair<size_t, size_t> Indices;
    Indices imin() const;  
    Indices imax() const;  
 
    double min() const;  
    double max() const;  
 
    bool null() const;         
    bool positive() const;     
    bool negative() const;     
    bool nonNegative() const;  
 
    friend std::ostream& operator<<(std::ostream&, const Matrix&);  
 
    Matrix& truncate(double min = 1.e-14);  
 
    Matrix& transpose();        
    Matrix transposed() const;  
    Matrix& invert();           
    Matrix inverted() const;    
 
    VectorRef column(size_t);     
    Vector column(size_t) const;  
    VectorRef row(size_t);        
    Vector row(size_t) const;     
    VectorRef diagonal();         
    Vector diagonal() const;      
 
  private:
    std::unique_ptr<gsl_matrix, void (*)(gsl_matrix*)> gsl_mat_;
  };
 
  class VectorRef : public gsl_vector_view {
  public:
    VectorRef(const gsl_vector_view&);
 
    VectorRef& operator=(const Vector&);                                          
    operator Vector() const;                                                      
    bool operator==(const Vector&) const;                                         
    inline bool operator!=(const Vector& vec) const { return !operator==(vec); }  
    double& operator()(size_t);                                                   
    double operator()(size_t) const;                                              
    friend std::ostream& operator<<(std::ostream&, const VectorRef&);             
 
  private:
    friend class Vector;
  };
 
  class Vector : public Matrix {
  public:
    explicit Vector(size_t num_coord, double def = 0.);
    Vector(const std::initializer_list<double>&);  
    Vector(const VectorRef&);                      
    Vector(const gsl_vector_const_view&);          
 
    size_t size() const;  
 
    Vector subset(size_t min, size_t max = 0ull) const;
 
    double& operator()(size_t);         
    double operator()(size_t) const;    
    double dot(const Vector&) const;    
    Vector cross(const Vector&) const;  
 
    friend std::ostream& operator<<(std::ostream&, const Vector&);  
  };
}  // namespace cepgen
 
#endif