matdatastructures.hpp

Go to the documentation of this file.

/*****************************************************************************
    NumeRe: Framework fuer Numerische Rechnungen
    Copyright (C) 2021  Erik Haenel et al.
 
    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
    (at your option) 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 MATDATASTRUCTURES_HPP
#define MATDATASTRUCTURES_HPP
 
#include <vector>
#include <string>
#include <complex>
#include <utility>
 
#include "../ui/error.hpp"
#include "../utils/stringtools.hpp"
 
class Matrix
{
    private:
        std::vector<mu::value_type> m_storage;       
        size_t m_rows;                               
        size_t m_cols;                               
        bool m_transpose;                            
 
        mu::value_type& get(size_t i, size_t j)
        {
            return m_transpose ? m_storage[i*m_cols + j] : m_storage[i + j*m_rows];
        }
 
        const mu::value_type& get(size_t i, size_t j) const
        {
            return m_transpose ? m_storage[i*m_cols + j] : m_storage[i + j*m_rows];
        }
 
        void bakeTransposition()
        {
            if (m_transpose)
            {
                std::vector<mu::value_type> buffer(m_cols*m_rows);
 
                for (size_t i = 0; i < m_rows; i++)
                {
                    for (size_t j = 0; j < m_cols; j++)
                    {
                        buffer[i + j*m_rows] = m_storage[i*m_cols + j];
                    }
                }
 
                std::swap(buffer, m_storage);
                m_transpose = false;
            }
        }
 
    public:
        Matrix() : m_rows(0u), m_cols(0u), m_transpose(false) {}
 
        Matrix(size_t r, size_t c = 1u, const mu::value_type& init = NAN) : m_rows(r), m_cols(c), m_transpose(false)
        {
            m_storage.resize(r*c, init);
        }
 
        Matrix(const Matrix& mat)
        {
            assign(mat);
        }
 
        // We do not declare the move constructor by ourselves
        Matrix(Matrix&& mat) = default;
 
        Matrix(const std::vector<std::vector<mu::value_type>>& vectMatrix) : Matrix()
        {
            m_rows = vectMatrix.size();
            m_cols = vectMatrix[0].size();
 
            for (size_t i = 1; i < m_rows; i++)
            {
                m_cols = std::max(m_cols, vectMatrix[i].size());
            }
 
            m_storage.resize(m_rows*m_cols, NAN);
 
            for (size_t i = 0; i < m_rows; i++)
            {
                for (size_t j = 0; j < vectMatrix[i].size(); j++)
                    get(i, j) = vectMatrix[i][j];
            }
        }
 
        Matrix(size_t r, size_t c, mu::value_type* vVals, int nVals)
        {
            assign(r, c, vVals, nVals);
        }
 
        Matrix& operator=(const Matrix& mat)
        {
            assign(mat);
            return *this;
        }
 
        void assign(const Matrix& mat)
        {
            m_rows = mat.m_rows;
            m_cols = mat.m_cols;
            m_storage = mat.m_storage;
            m_transpose = mat.m_transpose;
        }
 
        void assign(size_t r, size_t c, mu::value_type* vVals, int nVals)
        {
            m_storage.assign(vVals, vVals+std::min(r*c, (size_t)nVals));
            m_rows = r;
            m_cols = c;
            m_storage.resize(m_rows*m_cols);
            m_transpose = false;
        }
 
        void assign(size_t r, size_t c, const std::vector<mu::value_type>& vVals)
        {
            m_storage.assign(vVals.begin(), vVals.begin()+std::min(r*c, vVals.size()));
            m_rows = r;
            m_cols = c;
            m_storage.resize(m_rows*m_cols);
            m_transpose = false;
        }
 
        std::vector<mu::value_type>& data()
        {
            return m_storage;
        }
 
        const std::vector<mu::value_type>& data() const
        {
            return m_storage;
        }
 
        void transpose()
        {
            std::swap(m_rows, m_cols);
            m_transpose = !m_transpose;
        }
 
        size_t rows() const
        {
            return m_rows;
        }
 
        size_t cols() const
        {
            return m_cols;
        }
 
        std::pair<size_t,size_t> dims() const
        {
            return std::make_pair(m_rows, m_cols);
        }
 
        bool isScalar() const
        {
            return m_rows == 1 && m_cols == 1;
        }
 
        bool isSquare() const
        {
            return m_rows == m_cols;
        }
 
        bool isEmpty() const
        {
            return !m_rows || !m_cols;
        }
 
        bool containsInvalidValues() const
        {
            for (auto& v : m_storage)
            {
                if (mu::isnan(v) || mu::isinf(v))
                    return true;
            }
 
            return false;
        }
 
        std::string printDims() const
        {
            return toString(m_rows) + "x" + toString(m_cols);
        }
 
        mu::value_type& operator()(size_t i, size_t j = 0u)
        {
            if (i >= m_rows || j >= m_cols)
                throw SyntaxError(SyntaxError::INVALID_INDEX, "INTERNAL INDEXING ERROR",
                                  SyntaxError::invalid_position,
                                  "MAT(" + toString(i) + "," + toString(j) + ") vs. size = " + printDims());
 
            return get(i, j);
        }
 
        const mu::value_type& operator()(size_t i, size_t j = 0u) const
        {
            if (i >= m_rows || j >= m_cols)
                throw SyntaxError(SyntaxError::INVALID_INDEX, "INTERNAL INDEXING ERROR",
                                  SyntaxError::invalid_position,
                                  "MAT(" + toString(i) + "," + toString(j) + ") vs. size = " + printDims());
 
            return get(i, j);
        }
 
        void resize(size_t r, size_t c)
        {
            if (r == m_rows && c == m_cols)
                return;
 
            // Use a buffer
            std::vector<mu::value_type> buffer(r*c);
 
            // Copy the remaining elements
            #pragma omp parallel for
            for (size_t i = 0; i < std::min(r, m_rows); i++)
            {
                for (size_t j = 0; j < std::min(c, m_cols); j++)
                {
                    buffer[i + j*r] = m_storage[i + j*m_rows];
                }
            }
 
            // Swap storage with the buffer
            std::swap(m_storage, buffer);
 
            // Update the sizes
            m_rows = r;
            m_cols = c;
        }
 
        void extend(size_t r, size_t c)
        {
            if (isScalar())
                return;
 
            // We need to bake the transposition
            bakeTransposition();
 
            if (r != m_rows || c != m_cols)
            {
                // Keep the relevant scalar dimensions
                bool rowScalar = m_rows == 1 && r > 1;
                bool colScalar = m_cols == 1 && c > 1;
 
                // Resize the matrix
                resize(r, c);
 
                if (rowScalar)
                {
                    for (size_t i = 1; i < m_rows; i++)
                    {
                        for (size_t j = 0; j < m_cols; j++)
                        {
                            get(i, j) = get(0, j);
                        }
                    }
                }
                else if (colScalar)
                {
                    for (size_t j = 1; j < m_cols; j++)
                    {
                        for (size_t i = 0; i < m_rows; i++)
                        {
                            get(i, j) = get(i, 0);
                        }
                    }
                }
            }
        }
 
        Matrix& operator*=(const Matrix& mat)
        {
            assign(operator*(mat));
            return *this;
        }
 
        Matrix operator*(const Matrix& mat) const
        {
            if (mat.m_rows != m_cols)
                throw SyntaxError(SyntaxError::WRONG_MATRIX_DIMENSIONS_FOR_MATOP, "INTERNAL INDEXING ERROR",
                                  SyntaxError::invalid_position,
                                  printDims() + " vs. " + mat.printDims());
 
            Matrix ret(m_rows, mat.m_cols, 0.0);
 
            #pragma omp parallel for
            for (size_t i = 0; i < ret.m_rows; i++)
            {
                for (size_t j = 0; j < ret.m_cols; j++)
                {
                    for (size_t k = 0; k < m_cols; k++)
                    {
                        // Using the private method avoids addressing
                        // issues occuring with transposed matrices
                        ret.get(i, j) += get(i, k) * mat.get(k, j);
                    }
                }
            }
 
            return ret;
        }
 
        Matrix hcat(const Matrix& mat) const
        {
            size_t cols = m_cols + mat.m_cols;
            size_t rows = std::max(m_rows, mat.m_rows);
 
            Matrix mRet(rows, cols);
 
            // Copy this matrix
            for (size_t i = 0; i < m_rows; i++)
            {
                for (size_t j = 0; j < m_cols; j++)
                {
                    mRet.get(i, j) = get(i, j);
                }
            }
 
            // Copy other matrix
            for (size_t i = 0; i < mat.m_rows; i++)
            {
                for (size_t j = 0; j < mat.m_cols; j++)
                {
                    mRet.get(i, j+m_cols) = mat.get(i, j);
                }
            }
 
            return mRet;
        }
 
        Matrix vcat(const Matrix& mat) const
        {
            size_t cols = std::max(m_cols, mat.m_cols);
            size_t rows = m_rows + mat.m_rows;
 
            Matrix mRet(rows, cols);
 
            // Copy this matrix
            for (size_t i = 0; i < m_rows; i++)
            {
                for (size_t j = 0; j < m_cols; j++)
                {
                    mRet.get(i, j) = get(i, j);
                }
            }
 
            // Copy other matrix
            for (size_t i = 0; i < mat.m_rows; i++)
            {
                for (size_t j = 0; j < mat.m_cols; j++)
                {
                    mRet.get(i+m_rows, j) = mat.get(i, j);
                }
            }
 
            return mRet;
        }
};
 
 
 
// Erster Index: No. of Line; zweiter Index: No. of Col (push_back verwendet dazu stets zeilen!)
//typedef std::vector<std::vector<std::complex<double>> > Matrix;
 
 
struct MatFuncData
{
    const Matrix& mat1;
    const Matrix& mat2;
    const Matrix& mat3;
    mu::value_type fVal;
    int nVal;
    int mVal;
 
    MatFuncData() : mat1(Matrix()), mat2(Matrix()), mat3(Matrix()), fVal(NAN), nVal(0), mVal(0) {}
    MatFuncData(const Matrix& _mat, int n = 0) : mat1(_mat), mat2(Matrix()), mat3(Matrix()), fVal(NAN), nVal(n), mVal(0) {}
    MatFuncData(const Matrix& _mat1, const Matrix& _mat2, int n) : mat1(_mat1), mat2(_mat2), mat3(Matrix()), fVal(NAN), nVal(n), mVal(0) {}
    MatFuncData(const Matrix& _mat1, const Matrix& _mat2) : mat1(_mat1), mat2(_mat2), mat3(Matrix()), fVal(NAN), nVal(0), mVal(0) {}
    MatFuncData(const Matrix& _mat1, const Matrix& _mat2, const Matrix& _mat3) : mat1(_mat1), mat2(_mat2), mat3(_mat3), fVal(NAN), nVal(0), mVal(0) {}
    MatFuncData(const Matrix& _mat, mu::value_type f) : mat1(_mat), mat2(Matrix()), mat3(Matrix()), fVal(f), nVal(0), mVal(0) {}
    MatFuncData(const Matrix& _mat, mu::value_type f, int n) : mat1(_mat), mat2(Matrix()), mat3(Matrix()), fVal(f), nVal(n), mVal(0) {}
    MatFuncData(const Matrix& _mat, int n, int m) : mat1(_mat), mat2(Matrix()), mat3(Matrix()), fVal(NAN), nVal(n), mVal(m) {}
    MatFuncData(int n, int m = 0) : mat1(Matrix()), mat2(Matrix()), mat3(Matrix()), fVal(NAN), nVal(n), mVal(m) {}
};
 
 
struct MatFuncErrorInfo
{
    const std::string& command;
    const std::string& expression;
    size_t position;
 
    MatFuncErrorInfo(const std::string& _sCmd, const std::string& _sExpr, size_t pos) : command(_sCmd), expression(_sExpr), position(pos) {}
};
 
 
typedef Matrix (*MatFunc)(const MatFuncData&, const MatFuncErrorInfo&);
 
 
enum MatFuncSignature
{
    MATSIG_INVALID,
    MATSIG_MAT,
    MATSIG_MAT_NOPT,
    MATSIG_MAT_MAT,
    MATSIG_MAT_MAT_MAT,
    MATSIG_MAT_MAT_N,
    MATSIG_MAT_F,
    MATSIG_MAT_F_N,
    MATSIG_MAT_N_MOPT,
    MATSIG_N_MOPT
};
 
 
struct MatFuncDef
{
    MatFuncSignature signature;
    MatFunc func;
    bool isPureMatFunc;
 
    MatFuncDef() : signature(MATSIG_INVALID), func(nullptr), isPureMatFunc(false) {}
    MatFuncDef(MatFuncSignature sig, MatFunc f, bool isMatFunc = true) : signature(sig), func(f), isPureMatFunc(isMatFunc) {}
};
 
#endif // MATDATASTRUCTURES_HPP