std_library_matrix.cmajor

//
//     ,ad888ba,                              88
//    d8"'    "8b
//   d8            88,dba,,adba,   ,aPP8A.A8  88     The Cmajor Standard Library
//   Y8,           88    88    88  88     88  88
//    Y8a.   .a8P  88    88    88  88,   ,88  88     (C)2024 Cmajor Software Ltd
//     '"Y888Y"'   88    88    88  '"8bbP"Y8  88     https://cmajor.dev
//                                           ,88
//                                        888P"
//
//  The Cmajor standard library may be used under the terms of the ISC license:
//
//  Permission to use, copy, modify, and/or distribute this software for any purpose with or
//  without fee is hereby granted, provided that the above copyright notice and this permission
//  notice appear in all copies. THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL
//  WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
//  AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, INDIRECT, OR
//  CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS,
//  WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN
//  CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

/// std.matrix

/**
    This namespace contains various matrix manipulation functions.
*/
namespace std::matrix
{
    //==============================================================================
    /// Returns the matrix product of two 2-dimensional arrays.
    ElementType[n, m] multiply<ElementType, n, m, k> (ElementType[n, k] a,
                                                      ElementType[k, m] b)
    {
        ElementType[n, m] result;

        for (wrap<n> i)
            for (wrap<m> j)
                for (wrap<k> o)
                    result[i, j] += a[i, o] * b[o, j];

        return result;
    }

    /// Returns the dot-product of two 1-dimensional arrays.
    ElementType dot<ElementType, n> (ElementType[n] a,
                                     ElementType[n] b)
    {
        ElementType product;

        for (wrap<n> i)
            product += a[i] * b[i];

        return product;
    }

    /// Returns the dot-product of two 2-dimensional arrays.
    /// (The dot product for 2D arrays is the same as a matrix multiplication).
    ElementType[n, m] dot<ElementType, n, m, k> (ElementType[n, k] a,
                                                 ElementType[k, m] b)
    {
        return multiply (a, b);
    }

    /// Returns the inverse of the input 2-dimensional array.
    /// If no inverse is possible, it returns an empty array.
    /// This algorithm uses the Gaussian elimination method to calculate the inverse.
    ElementType[n, n] inverse<ElementType, n> (ElementType[n, n] matrix)
    {
        ElementType<n * 2>[n] m;

        for (wrap<n> i)
        {
            for (wrap<n> j)
                m[i][j] = matrix [i, j];

            m[i][n + i] = 1;
        }

        for (wrap<n> c)
        {
            // Partial pivot to ensure the largest unprocessed value is in [c][c]
            for (wrap<n> r = n - 1; r > c; r--)
                if (m.at (r-1)[0] < m[r][0])
                    swap (m.at (r-1), m[r]);

            // Inverse not possible, return empty matrix
            if (m[c][c] == 0)
                return ();

            m[c] = m[c] / m[c][c];

            for (wrap<n> r)
            {
                if (r != c)
                    m[r] = m[r] - m[r][c] * m[c];
            }
        }

        ElementType[n, n] result;

        for (wrap<n> i)
            for (wrap<n> j)
                result[i, j] = m[i][n + j];

        return result;
    }

    /// Returns the sum of two 2-dimensional arrays.
    ElementType[n, m] add<ElementType, n, m> (ElementType[n, m] a,
                                              ElementType[n, m] b)
    {
        ElementType[n, m] result;

        for (wrap<n> i)
            for (wrap<m> j)
                result[i, j] = a[i, j] + b[i, j];

        return result;
    }

    /// Returns the difference between two 2-dimensional arrays.
    ElementType[n, m] subtract<ElementType, n, m> (ElementType[n, m] a,
                                                   ElementType[n, m] b)
    {
        ElementType[n, m] result;

        for (wrap<n> i)
            for (wrap<m> j)
                result[i, j] = a[i, j] - b[i, j];

        return result;
    }

} // namespace std::matrix