#ifndef _BGC_MATRICES_H_
#define _BGC_MATRICES_H_

// ================== Matrix2x2 ================= //

typedef struct {
    float row1_col1, row1_col2;
    float row2_col1, row2_col2;
} BGC_FP32_Matrix2x2;

typedef struct {
    double row1_col1, row1_col2;
    double row2_col1, row2_col2;
} BGC_FP64_Matrix2x2;

// ================== Matrix2x3 ================= //

typedef struct {
    float row1_col1, row1_col2;
    float row2_col1, row2_col2;
    float row3_col1, row3_col2;
} BGC_FP32_Matrix2x3;

typedef struct {
    double row1_col1, row1_col2;
    double row2_col1, row2_col2;
    double row3_col1, row3_col2;
} BGC_FP64_Matrix2x3;

// ================== Matrix3x2 ================= //

typedef struct {
    float row1_col1, row1_col2, row1_col3;
    float row2_col1, row2_col2, row2_col3;
} BGC_FP32_Matrix3x2;

typedef struct {
    double row1_col1, row1_col2, row1_col3;
    double row2_col1, row2_col2, row2_col3;
} BGC_FP64_Matrix3x2;

// ================== Matrix3x3 ================= //

typedef struct {
    float row1_col1, row1_col2, row1_col3;
    float row2_col1, row2_col2, row2_col3;
    float row3_col1, row3_col2, row3_col3;
} BGC_FP32_Matrix3x3;

typedef struct {
    double row1_col1, row1_col2, row1_col3;
    double row2_col1, row2_col2, row2_col3;
    double row3_col1, row3_col2, row3_col3;
} BGC_FP64_Matrix3x3;

// ========== Matrix Product 2x2 at 2x2 ========= //

inline void bgc_fp32_multiply_matrix2x2_by_matrix2x2(BGC_FP32_Matrix2x2* product, const BGC_FP32_Matrix2x2* matrix1, const BGC_FP32_Matrix2x2* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
}

inline void bgc_fp64_multiply_matrix2x2_by_matrix2x2(BGC_FP64_Matrix2x2* product, const BGC_FP64_Matrix2x2* matrix1, const BGC_FP64_Matrix2x2* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
}

// ========== Matrix Product 2x2 at 3x2 ========= //

inline void bgc_fp32_multiply_matrix2x2_by_matrix3x2(BGC_FP32_Matrix3x2* product, const BGC_FP32_Matrix2x2* matrix1, const BGC_FP32_Matrix3x2* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;
    const float row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;
    const float row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;
}

inline void bgc_fp64_multiply_matrix2x2_by_matrix3x2(BGC_FP64_Matrix3x2* product, const BGC_FP64_Matrix2x2* matrix1, const BGC_FP64_Matrix3x2* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;
    const double row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;
    const double row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;
}

// ========== Matrix Product 2x3 at 2x2 ========= //

inline void bgc_fp32_multiply_matrix2x3_by_matrix2x2(BGC_FP32_Matrix2x3* product, const BGC_FP32_Matrix2x3* matrix1, const BGC_FP32_Matrix2x2* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;

    const float row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1;
    const float row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
}

inline void bgc_fp64_multiply_matrix2x3_by_matrix2x2(BGC_FP64_Matrix2x3* product, const BGC_FP64_Matrix2x3* matrix1, const BGC_FP64_Matrix2x2* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;

    const double row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1;
    const double row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
}

// ========== Matrix Product 2x3 at 3x2 ========= //

inline void bgc_fp32_multiply_matrix2x3_by_matrix3x2(BGC_FP32_Matrix3x3* product, const BGC_FP32_Matrix2x3* matrix1, const BGC_FP32_Matrix3x2* matrix2)
{
    product->row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    product->row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;
    product->row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3;

    product->row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    product->row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;
    product->row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3;

    product->row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1;
    product->row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2;
    product->row3_col3 = matrix1->row3_col1 * matrix2->row1_col3 + matrix1->row3_col2 * matrix2->row2_col3;
}

inline void bgc_fp64_multiply_matrix2x3_by_matrix3x2(BGC_FP64_Matrix3x3* product, const BGC_FP64_Matrix2x3* matrix1, const BGC_FP64_Matrix3x2* matrix2)
{
    product->row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1;
    product->row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2;
    product->row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3;

    product->row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1;
    product->row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2;
    product->row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3;

    product->row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1;
    product->row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2;
    product->row3_col3 = matrix1->row3_col1 * matrix2->row1_col3 + matrix1->row3_col2 * matrix2->row2_col3;
}

// ========== Matrix Product 3x2 at 2x3 ========= //

inline void bgc_fp32_multiply_matrix3x2_by_matrix2x3(BGC_FP32_Matrix2x2* product, const BGC_FP32_Matrix3x2* matrix1, const BGC_FP32_Matrix2x3* matrix2)
{
    product->row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    product->row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;

    product->row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    product->row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
}

inline void bgc_fp64_multiply_matrix3x2_by_matrix2x3(BGC_FP64_Matrix2x2* product, const BGC_FP64_Matrix3x2* matrix1, const BGC_FP64_Matrix2x3* matrix2)
{
    product->row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    product->row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;

    product->row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    product->row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
}

// ========== Matrix Product 3x2 at 3x3 ========= //

inline void bgc_fp32_multiply_matrix3x2_by_matrix3x3(BGC_FP32_Matrix3x2* product, const BGC_FP32_Matrix3x2* matrix1, const BGC_FP32_Matrix3x3* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;
    const float row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3 + matrix1->row1_col3 * matrix2->row3_col3;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
    const float row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3 + matrix1->row2_col3 * matrix2->row3_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;
}

inline void bgc_fp64_multiply_matrix3x2_by_matrix3x3(BGC_FP64_Matrix3x2* product, const BGC_FP64_Matrix3x2* matrix1, const BGC_FP64_Matrix3x3* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;
    const double row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3 + matrix1->row1_col3 * matrix2->row3_col3;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
    const double row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3 + matrix1->row2_col3 * matrix2->row3_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;
}

// ========== Matrix Product 3x3 at 2x3 ========= //

inline void bgc_fp32_multiply_matrix3x3_by_matrix2x3(BGC_FP32_Matrix2x3* product, const BGC_FP32_Matrix3x3* matrix1, const BGC_FP32_Matrix2x3* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;

    const float row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1 + matrix1->row3_col3 * matrix2->row3_col1;
    const float row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2 + matrix1->row3_col3 * matrix2->row3_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
}

inline void bgc_fp64_multiply_matrix3x3_by_matrix2x3(BGC_FP64_Matrix2x3* product, const BGC_FP64_Matrix3x3* matrix1, const BGC_FP64_Matrix2x3* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;

    const double row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1 + matrix1->row3_col3 * matrix2->row3_col1;
    const double row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2 + matrix1->row3_col3 * matrix2->row3_col2;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
}

// ========== Matrix Product 3x3 at 3x3 ========= //

inline void bgc_fp32_multiply_matrix3x3_by_matrix3x3(BGC_FP32_Matrix3x3* product, const BGC_FP32_Matrix3x3* matrix1, const BGC_FP32_Matrix3x3* matrix2)
{
    const float row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const float row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;
    const float row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3 + matrix1->row1_col3 * matrix2->row3_col3;

    const float row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const float row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
    const float row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3 + matrix1->row2_col3 * matrix2->row3_col3;

    const float row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1 + matrix1->row3_col3 * matrix2->row3_col1;
    const float row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2 + matrix1->row3_col3 * matrix2->row3_col2;
    const float row3_col3 = matrix1->row3_col1 * matrix2->row1_col3 + matrix1->row3_col2 * matrix2->row2_col3 + matrix1->row3_col3 * matrix2->row3_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
    product->row3_col3 = row3_col3;
}

inline void bgc_fp64_multiply_matrix3x3_by_matrix3x3(BGC_FP64_Matrix3x3* product, const BGC_FP64_Matrix3x3* matrix1, const BGC_FP64_Matrix3x3* matrix2)
{
    const double row1_col1 = matrix1->row1_col1 * matrix2->row1_col1 + matrix1->row1_col2 * matrix2->row2_col1 + matrix1->row1_col3 * matrix2->row3_col1;
    const double row1_col2 = matrix1->row1_col1 * matrix2->row1_col2 + matrix1->row1_col2 * matrix2->row2_col2 + matrix1->row1_col3 * matrix2->row3_col2;
    const double row1_col3 = matrix1->row1_col1 * matrix2->row1_col3 + matrix1->row1_col2 * matrix2->row2_col3 + matrix1->row1_col3 * matrix2->row3_col3;

    const double row2_col1 = matrix1->row2_col1 * matrix2->row1_col1 + matrix1->row2_col2 * matrix2->row2_col1 + matrix1->row2_col3 * matrix2->row3_col1;
    const double row2_col2 = matrix1->row2_col1 * matrix2->row1_col2 + matrix1->row2_col2 * matrix2->row2_col2 + matrix1->row2_col3 * matrix2->row3_col2;
    const double row2_col3 = matrix1->row2_col1 * matrix2->row1_col3 + matrix1->row2_col2 * matrix2->row2_col3 + matrix1->row2_col3 * matrix2->row3_col3;

    const double row3_col1 = matrix1->row3_col1 * matrix2->row1_col1 + matrix1->row3_col2 * matrix2->row2_col1 + matrix1->row3_col3 * matrix2->row3_col1;
    const double row3_col2 = matrix1->row3_col1 * matrix2->row1_col2 + matrix1->row3_col2 * matrix2->row2_col2 + matrix1->row3_col3 * matrix2->row3_col2;
    const double row3_col3 = matrix1->row3_col1 * matrix2->row1_col3 + matrix1->row3_col2 * matrix2->row2_col3 + matrix1->row3_col3 * matrix2->row3_col3;

    product->row1_col1 = row1_col1;
    product->row1_col2 = row1_col2;
    product->row1_col3 = row1_col3;

    product->row2_col1 = row2_col1;
    product->row2_col2 = row2_col2;
    product->row2_col3 = row2_col3;

    product->row3_col1 = row3_col1;
    product->row3_col2 = row3_col2;
    product->row3_col3 = row3_col3;
}

#endif // _BGC_MATRIX_TYPES_H_