matrixmath.hpp

/******************************************************************************
* Copyright (c) ATTX LLC 2024. All Rights Reserved.
*
* This software and associated documentation (the "Software") are the 
* proprietary and confidential information of ATTX, LLC. The Software is 
* furnished under a license agreement between ATTX and the user organization 
* and may be used or copied only in accordance with the terms of the agreement.
* Refer to 'license/attx_license.adoc' for standard license terms.
*
* EXPORT CONTROL NOTICE: THIS SOFTWARE MAY INCLUDE CONTENT CONTROLLED UNDER THE
* INTERNATIONAL TRAFFIC IN ARMS REGULATIONS (ITAR) OR THE EXPORT ADMINISTRATION 
* REGULATIONS (EAR99). No part of the Software may be used, reproduced, or 
* transmitted in any form or by any means, for any purpose, without the express 
* written permission of ATTX, LLC.
******************************************************************************/
/*
Matrix math header file
-----------------------
This file defines basic functions for mathematical operations on
the matrix class
Rather than raising an error on fail cases (i.e. divide by zero),
the following libraries are designed to check for errors and return
corresponding codes.
 
Note: The naming convention observed for matrix math is as follows:
- Capital letters represent matrices
- Lowercase letters represent scalars
- All vectors are 1-dimensional matrices for the purpose of math
- Pass by reference return values are all named "result"
- The letter used in math indicates the order of operations - A/a
  comes first, then B/b. A/a or B/b always represent the "other"
  value in an operation
 
Author: Alex Reynolds
*/
 
#ifndef MATRIXMATH_HPP
#define MATRIXMATH_HPP
 
#include "Matrix.hpp"
#include "clockwerkerrors.h"
 
namespace clockwerk {
 
    /// ---------------------------------------------------------------------------------------------------
    /// Standard functions
    /// ---------------------------------------------------------------------------------------------------
 
    /// @brief   Function to multiply two matrices
    /// @param   A   First matrix in the multiplication operation
    /// @param   B   Second matrix in the multiplication operation
    /// @param   result  PBR return of A*B result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    /// @note    Here we are using the "standard" O^3 matrix multiplication algorithm. Quick lookup
    ///          indicates that more highly optimized matrix multiplication algorithms, such as strassen,
    ///          require much larger sizes than we will typically be working with to be valuable.
    template<typename T, unsigned int R1, unsigned int C1R2, unsigned int C2>
    void multiply(Matrix<T, R1, C1R2> A, Matrix<T, C1R2, C2> B, Matrix<T, R1, C2> &result) {
        /// Check and enforce matrix size constraints
        /// Note: A number of checks necessary for matrix multiplication are performed at compile time.
        /// The number of rows and columns in matrix result needs to correcpond to the number of rows
        /// in matrix one and columns in matrix two -- this is enforced by templating size constraints
        /// and will not compile if incorrect
        /// The number of columns in matrix one need to correspond to the number of rows in matrix 2.
        /// This constraint is enforced by the template variable C1R2 and will not compile if incorrect
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j, k;
        for(i = 0; i < R1; i++) {
            for(j = 0; j < C2; j++) {
                result.values[i][j] = T();      /// Note: here we're explicitly calling the default constructor.
                                                /// should be zero for basic types but essentially enforces assumption
                                                /// that object in question has a default constructor.
                for(k = 0; k < C1R2; k++) {
                    /// NOTE: Per the above dimension check 
                    result.values[i][j] += A.values[i][k]*B.values[k][j];
                }
            }
        }
    }
 
    /// @brief   Function to multiply a scalar by a matrix
    /// @param   a   Scalar by which to multiply A
    /// @param   A   Matrix in the scalar multiplication operation
    /// @param   result  PBR return of a*A result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    template<typename T, unsigned int R, unsigned int C>
    void multiply(const T &a, Matrix<T, R, C> A, Matrix<T, R, C> &result) {
        /// Check and enforce matrix size constraints
        /// Note: The check necessary to ensure safe operations here is that the size of result
        /// and A match -- That is enforced via templating and the compiler at compile time
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = a*A.values[i][j];
            }
        }
    }
 
    /// @brief   Function to multiply two matrices element-wise
    /// @param   A   First matrix in the element-wise multiplication operation
    /// @param   B   Second matrix in the element-wise multiplication operation
    /// @param   result  PBR return of A.*B result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    template<typename T, unsigned int R, unsigned int C>
    void eMultiply(Matrix<T, R, C> A, Matrix<T, R, C> B, Matrix<T, R, C> &result) {
        /// Check and enforce matrix size constraints
        /// Note: The check necessary to ensure safe operations here is that the sizes of result,
        /// A, and B match -- That is enforced via templating and the compiler at compile time
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = A.values[i][j]*B.values[i][j];
            }
        }
    }
 
    /// @brief   Function to add a scalar to a matrix
    /// @param   a   Scalar to add to A
    /// @param   A   Matrix in the scalar addition operation
    /// @param   result  PBR return of a+A result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    template<typename T, unsigned int R, unsigned int C>
    void add(const T &a, Matrix<T, R, C> A, Matrix<T, R, C> &result) {
        /// Check and enforce matrix size constraints
        /// Note: The check necessary to ensure safe operations here is that the size of result
        /// and A match -- That is enforced via templating and the compiler at compile time
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = a + A.values[i][j];
            }
        }
    }
 
    /// @brief   Function to add two matrices element-wise
    /// @param   A   First matrix in the element-wise addition operation
    /// @param   B   Second matrix in the element-wise addition operation
    /// @param   result  PBR return of A+B result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    template<typename T, unsigned int R, unsigned int C>
    void eAdd(Matrix<T, R, C> A, Matrix<T, R, C> B, Matrix<T, R, C> &result) {
        /// Check and enforce matrix size constraints
        /// Note: The check necessary to ensure safe operations here is that the sizes of result,
        /// A, and B match -- That is enforced via templating and the compiler at compile time
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = A.values[i][j] + B.values[i][j];
            }
        }
    }
 
    /// @brief   Function to subtract B from A, element-wise
    /// @param   A   First matrix in the element-wise subtraction operation
    /// @param   B   Second matrix in the element-wise subtraction operation
    /// @param   result  PBR return of A-B result
    /// @return  Error code corresponding to errors in clockwerkerrors.h
    template<typename T, unsigned int R, unsigned int C>
    void eSubtract(Matrix<T, R, C> A, Matrix<T, R, C> B, Matrix<T, R, C> &result) {
        /// Check and enforce matrix size constraints
        /// Note: The check necessary to ensure safe operations here is that the sizes of result,
        /// A, and B match -- That is enforced via templating and the compiler at compile time
        /// Note: Because size checks are performed by templating, this function does not return an error code.
        /// If a size error exists, function will not compile
 
        /// Now we can perform our multiplication operation
        /// Note: Because of templating size checks, we know that no operations
        /// will request an index out of bounds in a matrix. Therefore, we can ignore these checks
        /// for speed in this function
        unsigned int i, j;
        for(i = 0; i < R; i++) {
            for(j = 0; j < C; j++) {
                result.values[i][j] = A.values[i][j] - B.values[i][j];
            }
        }
    }
 
    /// ------------------------------------------------------------------------------------------------------
    /// Operator overloads -- NOTE: May be less efficient than equivalent PBR function in some cases due to 
    ///                             return by value. Where overload is less efficient, it is explicitly noted.
    /// ------------------------------------------------------------------------------------------------------
    /// @brief   Overload of matrix multiplication - Two matrices
    /// @note    Returns value -- less efficient in some cases
    template<typename T, unsigned int R1, unsigned int C1R2, unsigned int C2>
    Matrix<T, R1, C2> operator*(const Matrix<T, R1, C1R2> &A, const Matrix<T, C1R2, C2> &B) {
        /// Call our function to add two matrices and return the result
        Matrix<T, R1, C2> result;
        multiply(A, B, result);
        return result;
    }
 
    /// @brief   Overload of matrix multiplication - Matrix and scalar
    /// @note    Returns value -- less efficient in some cases
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator*(const T &a, const Matrix<T, R, C> &A) {
        /// Call our function to multiply and return the result
        Matrix<T, R, C> result;
        multiply(a, A, result);
        return result;
    }
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator*(const Matrix<T, R, C> &A, const T &a) {
        /// Call our function to multiply and return the result
        Matrix<T, R, C> result;
        multiply(a, A, result);
        return result;
    }
 
    /// @brief   Overload of matrix addition - Two matrices
    /// @note    Returns value -- less efficient in some cases
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator+(const Matrix<T, R, C> &A, const Matrix<T, R, C> &B) {
        /// Call our function to add two matrices and return the result
        Matrix<T, R, C> result;
        eAdd(A, B, result);
        return result;
    }
 
    /// @brief   Overload of matrix addition - Matrix and scalar
    /// @note    Returns value -- less efficient in some cases
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator+(const T &a, const Matrix<T, R, C> &A) {
        /// Call our function to add two matrices and return the result
        Matrix<T, R, C> result;
        add(a, A, result);
        return result;
    }
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator+(const Matrix<T, R, C> &A, const T &a) {
        /// Call our function to add two matrices and return the result
        Matrix<T, R, C> result;
        add(a, A, result);
        return result;
    }
 
    /// Returns value -- less efficient in some cases
    template<typename T, unsigned int R, unsigned int C>
    Matrix<T, R, C> operator-(const Matrix<T, R, C> &A, const Matrix<T, R, C> &B) {
        /// Call our function to add two matrices and return the result
        Matrix<T, R, C> result;
        eSubtract(A, B, result);
        return result;
    }
 
}
 
#endif