This function calculates the outer product of two matrices. If mat_a has p dimensions and mat_b has q dimensions, then their outer product mat_y has p + q dimensions. In detail:

• if the subscript of some component of mat_a is ( a₀, a₁, … a_p−1 ),
• and the subscript of some component of mat_b is ( b₀, b₁, … b_q−1 ),
• then the subscript of some component of mat_y will be ( a₀, a₁, … a_p−1, b₀, b₁, … b_q−1 )
• and that component's value will be mat_a ( a₀, a₁, … a_p−1 ) times mat_b ( b₀, b₁, … b_q−1 )

template <
    typename mul_oper,
    typename compo_a,
    typename compo_b = compo_a,
    typename compo_y = compo_a,
    typename alloc_a = std::allocator <compo_a>,
    typename alloc_b = std::allocator <compo_b>,
    typename alloc_y = std::allocator <compo_y>
> var_matrix<compo_y,alloc_y> outer_product (
    sub_matrix_base<compo_a,alloc_a> const & mat_a,
    sub_matrix_base<compo_b,alloc_b> const & mat_b
);

#include "SOME_DIRECTORY/mat_gen_dim.h"
using namespace mat_gen_dim;

form const frm_a {blt{3,5},blt{3,6},blt{5,7}};
form const frm_b {blt{2,4},blt{1,3}};

struct mul_double {
    double operator() (double const a, double const b) { return a * b; }
};

int main() {
    cout.setf (std::ios::fixed, std::ios::floatfield);

    cout.precision (3);
    var_matrix<double> mat_a {sub_matrix_linear (frm_a,1.0)};
    cout << "\n matrix a: " << mat_a;
    con_matrix<double> mat_acon {mat_a};

    cout.precision (2);
    var_matrix<double> mat_b {sub_matrix_linear (frm_b,2.0)};
    cout << "\n matrix b: " << mat_b;
    con_matrix<double> mat_bcon {mat_b};

    cout.precision (5);
    var_matrix<double> const mat_y { outer_product<mul_double> (mat_a, mat_bcon) };
    cout << "\n outer product of a and b: " << mat_y;
    con_matrix<double> mat_ycon {mat_y};

    return 0;
}

 matrix a: 
    ( 3 3 5 ) = 1.335    ( 3 3 6 ) = 1.336    ( 3 4 5 ) = 1.345    ( 3 4 6 ) = 1.346
    ( 3 5 5 ) = 1.355    ( 3 5 6 ) = 1.356    ( 4 3 5 ) = 1.435    ( 4 3 6 ) = 1.436
    ( 4 4 5 ) = 1.445    ( 4 4 6 ) = 1.446    ( 4 5 5 ) = 1.455    ( 4 5 6 ) = 1.456
 matrix b: 
    ( 2 1 ) = 2.21    ( 2 2 ) = 2.22    ( 3 1 ) = 2.31    ( 3 2 ) = 2.32
 outer product of a and b: 
    ( 3 3 5 2 1 ) = 2.95035    ( 3 3 5 2 2 ) = 2.96370    ( 3 3 5 3 1 ) = 3.08385    ( 3 3 5 3 2 ) = 3.09720
    ( 3 3 6 2 1 ) = 2.95256    ( 3 3 6 2 2 ) = 2.96592    ( 3 3 6 3 1 ) = 3.08616    ( 3 3 6 3 2 ) = 3.09952
    ( 3 4 5 2 1 ) = 2.97245    ( 3 4 5 2 2 ) = 2.98590    ( 3 4 5 3 1 ) = 3.10695    ( 3 4 5 3 2 ) = 3.12040
    ( 3 4 6 2 1 ) = 2.97466    ( 3 4 6 2 2 ) = 2.98812    ( 3 4 6 3 1 ) = 3.10926    ( 3 4 6 3 2 ) = 3.12272
    ( 3 5 5 2 1 ) = 2.99455    ( 3 5 5 2 2 ) = 3.00810    ( 3 5 5 3 1 ) = 3.13005    ( 3 5 5 3 2 ) = 3.14360
    ( 3 5 6 2 1 ) = 2.99676    ( 3 5 6 2 2 ) = 3.01032    ( 3 5 6 3 1 ) = 3.13236    ( 3 5 6 3 2 ) = 3.14592
    ( 4 3 5 2 1 ) = 3.17135    ( 4 3 5 2 2 ) = 3.18570    ( 4 3 5 3 1 ) = 3.31485    ( 4 3 5 3 2 ) = 3.32920
    ( 4 3 6 2 1 ) = 3.17356    ( 4 3 6 2 2 ) = 3.18792    ( 4 3 6 3 1 ) = 3.31716    ( 4 3 6 3 2 ) = 3.33152
    ( 4 4 5 2 1 ) = 3.19345    ( 4 4 5 2 2 ) = 3.20790    ( 4 4 5 3 1 ) = 3.33795    ( 4 4 5 3 2 ) = 3.35240
    ( 4 4 6 2 1 ) = 3.19566    ( 4 4 6 2 2 ) = 3.21012    ( 4 4 6 3 1 ) = 3.34026    ( 4 4 6 3 2 ) = 3.35472
    ( 4 5 5 2 1 ) = 3.21555    ( 4 5 5 2 2 ) = 3.23010    ( 4 5 5 3 1 ) = 3.36105    ( 4 5 5 3 2 ) = 3.37560
    ( 4 5 6 2 1 ) = 3.21776    ( 4 5 6 2 2 ) = 3.23232    ( 4 5 6 3 1 ) = 3.36336    ( 4 5 6 3 2 ) = 3.37792