nummatrix.h

Go to the documentation of this file.

/*--------------------------------*- VM2D -*-----------------*---------------*\
| ##  ## ##   ##  ####  #####   |                            | Version 1.14   |
| ##  ## ### ### ##  ## ##  ##  |  VM2D: Vortex Method       | 2026/03/06     |
| ##  ## ## # ##    ##  ##  ##  |  for 2D Flow Simulation    *----------------*
|  ####  ##   ##   ##   ##  ##  |  Open Source Code                           |
|   ##   ##   ## ###### #####   |  https://www.github.com/vortexmethods/VM2D  |
|                                                                             |
| Copyright (C) 2017-2026 I. Marchevsky, K. Sokol, E. Ryatina, A. Kolganova   |
*-----------------------------------------------------------------------------*
| File name: nummatrix.h                                                      |
| Info: Source code of VM2D                                                   |
|                                                                             |
| This file is part of VM2D.                                                  |
| VM2D 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.                                         |
|                                                                             |
| VM2D 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 VM2D.  If not, see <http://www.gnu.org/licenses/>.               |
\*---------------------------------------------------------------------------*/
 
 
#ifndef NUMMATRIX_H_
#define NUMMATRIX_H_
 
#include "numvector.h"
 
namespace VMlib
{
 
    template<typename T, size_t n, size_t m>
    class nummatrix : public std::array<T, n*m>
    {
 
    public:
 
    public:
 
        numvector<T, m>& operator[](size_t i)
        {
            return *(reinterpret_cast<numvector<T, m>*>(this->data() + (i*m)));
        }
 
        const numvector<T, m>& operator[](size_t i) const
        {
            return *(reinterpret_cast<const numvector<T, m>*>(this->data() + (i*m)));
        }
 
 
 
        template <typename P>
        nummatrix<T, n, m>& operator*=(const P c)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] *= c;
            return *this;
        }//operator*=
 
 
        template <typename P>
        nummatrix<T, n, m>& operator/=(const P c)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] /= c;
            return *this;
        }//operator/=
 
 
        template <typename P>
        nummatrix<T, n, m>& operator+=(const nummatrix<P, n, m>& y)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] += y.data()[i];
            return *this;
        }//operator+=
 
 
        template <typename P>
        nummatrix<T, n, m>& operator-=(const nummatrix<P, n, m>& y)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] -= y.data()[i];
            return *this;
        }//operator-=
 
 
        template <typename P>
        auto operator+(const nummatrix<P, n, m>& y) const -> nummatrix<typename std::remove_const<decltype(this->data()[0] + y.data()[0])>::type, n, m>
        {
            nummatrix<typename std::remove_const<decltype(this->data()[0] + y.data()[0])>::type, n, m> res;
            for (size_t i = 0; i < n*m; ++i)
                res.data()[i] = this->data()[i] + y.data()[i];
            return res;
        }//operator+
 
 
        template <typename P>
        auto operator-(const nummatrix<P, n, m>& y) const -> nummatrix<typename std::remove_const<decltype(this->data()[0] - y.data()[0])>::type, n, m>
        {
            nummatrix<typename std::remove_const<decltype(this->data()[0] - y.data()[0])>::type, n, m> res;
            for (size_t i = 0; i < n*m; ++i)
                res.data()[i] = this->data()[i] - y.data()[i];
            return res;
        }//operator-
 
 
        template <typename P>
        auto operator*(const P c) const -> nummatrix<typename std::remove_const<decltype(this->data()[0] * c)>::type, n, m>
        {
            nummatrix<typename std::remove_const<decltype(this->data()[0] * c)>::type, n, m> res;
            for (size_t i = 0; i < n*m; ++i)
                res.data()[i] = c * this->data()[i];
            return res;
        }//operator*
 
 
 
        nummatrix<T, n, m> operator-() const
        {
            nummatrix<T, n, m> res;
            for (size_t i = 0; i < n*m; ++i)
                res.data()[i] = -this->data()[i];
            return res;
        }//operator-
 
 
        template <typename P>
        bool operator==(const nummatrix<P, n, m>& y) const
        {
            for (size_t i = 0; i < n*m; ++i)
                if (this->data()[i] != y.data()[i])
                    return false;
            return true;
        }//operator==
 
 
        template <typename P>
        bool operator!=(const nummatrix<P, n, m>& y) const
        {
            return !(*this == y);
        }//operator!=
 
 
        std::pair<size_t, size_t> size() const { return{ n, m }; }
 
 
        auto norm1() const -> typename std::remove_const<typename std::remove_reference<decltype(this->data()[0])>::type>::type
        {
            std::vector<typename std::remove_const<typename std::remove_reference<decltype(this->data()[0])>::type>::type> ms(m, 0);
            for (size_t i = 0; i < m; ++i)
                for (size_t j = 0; j < n; ++j)
                    ms[i] += abs(this->data()[j*m + i]);
            return std::max_element(ms.begin(), ms.end());
        }
 
        auto norminf() const -> typename std::remove_const<typename std::remove_reference<decltype(this->data()[0])>::type>::type
        {
            std::vector<typename std::remove_const<typename std::remove_reference<decltype(this->data()[0])>::type>::type> ms(n, 0);
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    ms[i] += abs(this->data()[i*m + j]);
            return std::max_element(ms.begin(), ms.end());
        }
 
 
        nummatrix<T, n, m>& toZero(T val = 0.0)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] = val;
            return *this;
        }
 
 
        nummatrix<T, n, n>& toIdentity()
        {
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    this->data()[i*m + j] = (i == j) ? 1 : 0;
            return *this;
        }
 
        nummatrix<T, m, n> t() const
        {
            nummatrix<T, m, n> res;
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    res[j][i] = (*this)[i][j];
            return res;
        }
 
        nummatrix<T, n, n> sym() const
        {
            return (*this + this->t())*0.5;
        }
 
        nummatrix() { };
 
 
        template <typename P>
        explicit nummatrix(const P z)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] = z;
        }//nummatrix(...)
 
 
        template <typename P>
        nummatrix(const nummatrix<P, n, m>& mtr)
        {
            for (size_t i = 0; i < n*m; ++i)
                this->data()[i] = mtr.data()[i];
        }//nummatrix(...)
 
 
        template <typename P>
        nummatrix(const std::vector<std::vector<P>>& vec)
        {
            if ((vec.size() != n) || (vec[0].size() != m))
                throw;
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    this->data()[i*m + j] = vec[i][j];
        }//nummatrix(...)
 
 
        template <typename P>
        nummatrix(const numvector<numvector<P, m>, n>& mtr)
        {
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    this->data()[i*m + j] = mtr[i][j];
        }//nummatrix(...)
 
 
#if !defined(__CUDACC__)
        nummatrix(const std::initializer_list<numvector<T, m>>& z)
        {
            if (z.size() != n)
                throw;
            for (size_t i = 0; i < n; ++i)
            {
                const numvector<T, m>& numv = *(z.begin() + i);
                for (size_t j = 0; j < m; ++j)
                    this->data()[i*m + j] = numv[j];
            }
        }//nummatrix(...)
 
 
        template <typename P>
        nummatrix(const std::initializer_list<numvector<P, m>>& z)
        {
            if (z.size() != n)
                throw;
            for (size_t i = 0; i < n; ++i)
            {
                const numvector<T, m>& numv = *(z.begin() + i);
                for (size_t j = 0; j < m; ++j)
                    this->data()[i*m + j] = numv[j];
            }
        }//nummatrix(...)
#endif 
 
 
        operator std::vector<std::vector<T>>() const
        {
            std::vector<std::vector<T>> vec;
            vec.reserve(n);
            for (size_t i = 0; i < n; ++i)
            {
                vec[i].reserve(m);
                for (size_t j = 0; j < m; ++j)
                    vec[i].push_back(this->data()[i*m + j]);
            }
            return vec;
        }
 
 
        inline nummatrix<T, m, n> transpose() const
        {
            nummatrix<T, m, n> res;
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < m; ++j)
                    res[j][i] = this->data()[i*m + j];
            return res;
        }//Transpose(...)
 
 
        template<typename P>
        auto operator&(const numvector<P, m>& x) -> \
            numvector<typename std::remove_const<decltype(this->data()[0] * x[0])>::type, n>
        {
            numvector<typename std::remove_const<decltype(this->data()[0] * x[0])>::type, n> res;
 
            for (size_t i = 0; i < n; ++i)
                res[i] = (*this)[i] & x;
 
            return res;
        }//operator&
 
        template<typename P, size_t p>
        auto operator&(const nummatrix<P, m, p>& B) const -> \
            nummatrix<typename std::remove_const<decltype(this->data()[0] * B.data()[0])>::type, n, p>
        {
            nummatrix<typename std::remove_const<decltype(this->data()[0] * B.data()[0])>::type, n, p> res;
 
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < p; ++j)
                    res[i][j] = 0;
 
            for (size_t i = 0; i < n; ++i)
                for (size_t j = 0; j < p; ++j)
                    for (size_t k = 0; k < m; ++k)
                        res[i][j] += (*this)[i][k] * B[k][j];
            return res;
        }//operator&
 
    };
 
 
    template<typename T, typename P, size_t n, size_t m>
    auto operator&(const numvector<P, n>& x, const nummatrix<T, n, m>& A) -> \
        numvector<typename std::remove_const<decltype(x[0] * A.data()[0])>::type, m>
    {
        numvector<typename std::remove_const<decltype(x[0] * A.data()[0])>::type, m> res;
        for (size_t i = 0; i < m; ++i)
        {
            res[i] = 0;
            for (size_t j = 0; j < n; ++j)
                res[i] += x[j] * A[j][i];
        }
        return res;
    }//operator&(...)
 
 
    template<typename T, typename P, size_t n, size_t m>
    inline auto operator| (const numvector<T, n>& x, const numvector<P, m>& y) -> \
        nummatrix<typename std::remove_const<decltype(x[0] * y[0])>::type, n, m>
    {
        nummatrix<typename std::remove_const<decltype(x[0] * y[0])>::type, n, m> res;
 
        for (size_t i = 0; i < n; ++i)
            for (size_t j = 0; j < m; ++j)
                res[i][j] = x[i] * y[j];
 
        return res;
    }//operator|(...)
 
 
 
    template<typename T, size_t n, size_t m>
    std::ostream& operator<< (std::ostream& str, const nummatrix<T, n, m>& x)
    {
        str << "{ ";
        for (size_t j = 0; j < n - 1; ++j)
            str << x[j] << ", ";
        str << x[n - 1];
        str << " }";
        return str;
    }//operator<<
 
}//namespace VMlib
 
using VMlib::nummatrix;
 
#endif