Matrix.h

Go to the documentation of this file.

// SPDX-License-Identifier: LGPL-3.0-or-later
/***********************************************************************************/
/*  Copyright 2010 WSL Institute for Snow and Avalanche Research    SLF-DAVOS      */
/***********************************************************************************/
/* This file is part of MeteoIO.
    MeteoIO 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.
 
    MeteoIO 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 Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public License
    along with MeteoIO.  If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef MATRIX_H
#define MATRIX_H
 
#include <vector>
#include <iostream>
 
namespace mio {
 
class Matrix {
        public:
                Matrix() : vecData(), ncols(0), nrows(0) {}
 
                Matrix(const int& rows, const int& cols);
                Matrix(const size_t& rows, const size_t& cols) : vecData(rows*cols), ncols(cols), nrows(rows) {}
 
                Matrix(const size_t& rows, const size_t& cols, const double& init) : vecData(rows*cols, init), ncols(cols), nrows(rows) {}
 
                Matrix(const size_t& n, const double& init);
 
                Matrix(const size_t& rows, const size_t& cols, const std::vector<double> data) : vecData(data), ncols(cols), nrows(rows) {}
 
                void identity(const size_t& n, const double& init = 1.);
 
                void resize(const size_t& rows, const size_t& cols);
                void resize(const size_t& rows, const size_t& cols, const double& init);
                void resize(const size_t& rows, const size_t& cols, const std::vector<double>& data);
 
                void size(size_t& rows, size_t& cols) const;
 
                size_t getNx() const {return ncols;}
                size_t getNy() const {return nrows;}
 
                void clear();
 
                void random(const double& range);
 
                double& operator ()(const size_t& x, const size_t& y);
                double operator ()(const size_t& x, const size_t& y) const;
 
                double scalar() const;
                static double scalar(const Matrix& m);
 
                static double dot(const Matrix& A, const Matrix& B);
 
                Matrix getT() const;
                static Matrix T(const Matrix& m);
                void T();
 
                Matrix getInv() const;
                bool inv();
 
                Matrix getRow(const size_t i) const;
 
                void setRow(const size_t i, const Matrix& row);
 
                Matrix getCol(const size_t j) const;
 
                void setCol(const size_t j, const Matrix& col);
 
                Matrix getDiagonal() const;
 
                Matrix extract(size_t r_low, size_t r_high, size_t c_low, size_t c_high) const;
                
                double maxCoeff(size_t& max_row, size_t& max_col) const;
 
                static Matrix solve(const Matrix& A, const Matrix& B);
 
                static bool solve(const Matrix& A, const Matrix& B, Matrix& X);
 
                static Matrix TDMA_solve(const Matrix& A, const Matrix& B);
 
                static bool TDMA_solve(const Matrix& A, const Matrix& B, Matrix& X);
 
                double det() const;
 
                static void gaussElimination(Matrix& M, std::vector<size_t>& p);
 
                static bool gaussSolve(Matrix& M, Matrix& A, Matrix& X);
 
                static bool gaussSolve(const Matrix& M, const Matrix& A, Matrix& X);
 
                static bool gaussInverse(Matrix& M);
 
                static bool gaussInverse(const Matrix& M, Matrix& Inv);
 
                static double gaussDet(Matrix& M);
 
                static double normEuclid(const Matrix& vv);
 
                bool LU(Matrix& L, Matrix& U) const;
 
                void partialPivoting(std::vector<size_t>& pivot_idx);
                void partialPivoting();
 
                void maximalPivoting();
 
                //void bidiagonalize();
 
                static unsigned int eigenvaluesJacobi(Matrix& AA, Matrix& DD);
                static double jacobiEpsilon(Matrix& AA);
 
                static void svdJacobi(const Matrix& MM, Matrix& UU, Matrix& SS, Matrix& VV);
                static void sortEigenvalues(Matrix& EE, Matrix& VV);
                void swapCols(const size_t &j1, const size_t &j2); //used by svd
 
                const std::string toString(const int& precision=2, const bool& prettify=true) const;
 
                Matrix& operator+=(const Matrix& rhs);
                const Matrix operator+(const Matrix& rhs) const;
                Matrix& operator+=(const double& rhs);
                const Matrix operator+(const double& rhs) const;
 
                Matrix& operator-=(const Matrix& rhs);
                const Matrix operator-(const Matrix& rhs) const;
                Matrix& operator-=(const double& rhs);
                const Matrix operator-(const double& rhs) const;
 
                Matrix& operator*=(const Matrix& rhs);
                const Matrix operator*(const Matrix& rhs) const;
                Matrix& operator*=(const double& rhs);
                const Matrix operator*(const double& rhs) const;
 
                Matrix& operator/=(const double& rhs);
                const Matrix operator/(const double& rhs) const;
 
                bool operator==(const Matrix&) const; 
                bool operator!=(const Matrix&) const; 
 
                bool isIdentity() const;
                static bool isIdentity(const Matrix& A);
 
                static const double epsilon, epsilon_mtr;
 
        protected:
                std::vector<double> vecData;
                size_t ncols;
                size_t nrows;
 
                size_t findMaxInCol(const size_t &col);
                size_t findMaxInRow(const size_t &row);
                void swapRows(const size_t &i1, const size_t &i2);
 
};
 
} //end namespace
 
#endif