nummatrix.h

Go to the documentation of this file.

/*--------------------------------*- VMlib -*----------------*---------------*\
| ##  ## ##   ## ##   ##  ##    |                            | Version 1.12   |
| ##  ## ### ### ##       ##    |  VMlib: VM2D/VM3D Library  | 2024/01/14     |
| ##  ## ## # ## ##   ##  ####  |  Open Source Code          *----------------*
|  ####  ##   ## ##   ##  ## ## |  https://www.github.com/vortexmethods/VM2D  |
|   ##   ##   ## #### ### ####  |  https://www.github.com/vortexmethods/VM3D  |
|                                                                             |
| Copyright (C) 2017-2024 Ilia Marchevsky                                     |
*-----------------------------------------------------------------------------*
| File name: nummatrix.h                                                      |
| Info: Source code of VMlib                                                  |
|                                                                             |
| This file is part of VMlib.                                                 |
| VMLib 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.                                         |
|                                                                             |
| VMlib 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 VMlib.  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];
        }

        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