The following program first obtains the conventional matrix product of two ordinary matrices (con_matrix<double>). After that, it constructs equivalent block matrices (var_matrix<var_matrix<double>>) and multiples those, obtaining the equivalent product.

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

// the following will always be the initialization of "which" for conventional matrix multiplication
vec_bool const which {false,true,true,false};

struct mul_dou {
    double operator() (double const x, double const y) { return x * y; }
};

struct add_asgn_dou {
    double operator() (double & x, double const y) { return x += y; }
};

struct mul_mat_dou {
    var_matrix<double> operator() (
        con_matrix<double> const & x,
        con_matrix<double> const & y)
    {
        return inner_product<mul_dou, add_asgn_dou> (x, y, which);
    }
};

struct add_asgn_mat_dou {
    var_matrix<double> operator() (
        var_matrix<double> & x,
        con_matrix<double> const & y)
    {
        return x += y;
    }
};

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

    // you can change these numbers:
    sign_int const A (3); // non-negative
    sign_int const B (5); // non-negative
    sign_int const C (6); // positive
    sign_int const D (4); // positive
    sign_int const E (2); // non-negative
    sign_int const F (1); // non-negative

    form const f_AB_CD { blt (0, A + B), blt (0, C + D) };
    form const f_A_C   { blt (0, A    ), blt (0, C    ) };
    form const f_A_D   { blt (0, A    ), blt (C, C + D) };
    form const f_B_C   { blt (A, A + B), blt (0, C    ) };
    form const f_B_D   { blt (A, A + B), blt (C, C + D) };

    form const f_CD_EF { blt (0, C + D), blt (0, E + F) };
    form const f_C_E   { blt (0, C    ), blt (0, E    ) };
    form const f_C_F   { blt (0, C    ), blt (E, E + F) };
    form const f_D_E   { blt (C, C + D), blt (0, E    ) };
    form const f_D_F   { blt (C, C + D), blt (E, E + F) };

    cout << "\n\n the unblocked matrix: "; 

    con_matrix<double> m_AB_CD {sub_matrix_linear<double> (f_AB_CD,1.0)};
    con_matrix<double> m_CD_EF {sub_matrix_linear<double> (f_CD_EF,2.0)};
    con_matrix<double> m_AB_EF {inner_product
        <mul_dou, add_asgn_dou> (m_AB_CD, m_CD_EF, which)};

    cout << "\n m_AB_CD: " << m_AB_CD;
    cout << "\n m_CD_EF: " << m_CD_EF;
    cout << "\n m_AB_EF: " << m_AB_EF;

    cout << "\n\n extraction of blocks: ";

    con_matrix<double> m_A_C (con_lens_less_range<double> (m_AB_CD, f_A_C));
    con_matrix<double> m_A_D (con_lens_less_range<double> (m_AB_CD, f_A_D));
    con_matrix<double> m_B_C (con_lens_less_range<double> (m_AB_CD, f_B_C));
    con_matrix<double> m_B_D (con_lens_less_range<double> (m_AB_CD, f_B_D));

    cout << "\n m_A_C: " << m_A_C;
    cout << "\n m_A_D: " << m_A_D;
    cout << "\n m_B_C: " << m_B_C;
    cout << "\n m_B_D: " << m_B_D;

    con_matrix<double> m_C_E (con_lens_less_range<double> (m_CD_EF, f_C_E));
    con_matrix<double> m_C_F (con_lens_less_range<double> (m_CD_EF, f_C_F));
    con_matrix<double> m_D_E (con_lens_less_range<double> (m_CD_EF, f_D_E));
    con_matrix<double> m_D_F (con_lens_less_range<double> (m_CD_EF, f_D_F));

    cout << "\n m_C_E: " << m_C_E;
    cout << "\n m_C_F: " << m_C_F;
    cout << "\n m_D_E: " << m_D_E;
    cout << "\n m_D_F: " << m_D_F;

    cout << "\n\n assembly of blocks into compound matrix: ";

    form const f_22 { blt (0, 2), blt (0, 2) };

    var_matrix<var_matrix<double>> m_m_AB_CD (f_22, {m_A_C, m_A_D, m_B_C, m_B_D});
    var_matrix<var_matrix<double>> m_m_CD_EF (f_22, {m_C_E, m_C_F, m_D_E, m_D_F});
    var_matrix<var_matrix<double>> m_m_AB_EF {inner_product
        <mul_mat_dou, add_asgn_mat_dou> (m_m_AB_CD, m_m_CD_EF, which)};

    cout << "\n m_m_AB_CD: " << m_m_AB_CD;
    cout << "\n m_m_CD_EF: " << m_m_CD_EF;
    cout << "\n m_m_AB_EF: " << m_m_AB_EF;

    return 0;
}

 m_AB_EF: 
    ( 0 0 ) = { 25.6850 }  ( 0 1 ) = { 25.7895 }  ( 0 2 ) = { 25.8940 }
    ( 1 0 ) = { 28.1350 }  ( 1 1 ) = { 28.2495 }  ( 1 2 ) = { 28.3640 }
    ( 2 0 ) = { 30.5850 }  ( 2 1 ) = { 30.7095 }  ( 2 2 ) = { 30.8340 }
    ( 3 0 ) = { 33.0350 }  ( 3 1 ) = { 33.1695 }  ( 3 2 ) = { 33.3040 }
    ( 4 0 ) = { 35.4850 }  ( 4 1 ) = { 35.6295 }  ( 4 2 ) = { 35.7740 }
    ( 5 0 ) = { 37.9350 }  ( 5 1 ) = { 38.0895 }  ( 5 2 ) = { 38.2440 }
    ( 6 0 ) = { 40.3850 }  ( 6 1 ) = { 40.5495 }  ( 6 2 ) = { 40.7140 }
    ( 7 0 ) = { 42.8350 }  ( 7 1 ) = { 43.0095 }  ( 7 2 ) = { 43.1840 }

 m_m_AB_EF: 
    ( 0 0 ) = { 
                ( 0 0 ) = { 25.6850 }  ( 0 1 ) = { 25.7895 }
                ( 1 0 ) = { 28.1350 }  ( 1 1 ) = { 28.2495 }
                ( 2 0 ) = { 30.5850 }  ( 2 1 ) = { 30.7095 } 
              }
    ( 0 1 ) = { 
                ( 0 2 ) = { 25.8940 }
                ( 1 2 ) = { 28.3640 }
                ( 2 2 ) = { 30.8340 } 
              }
    ( 1 0 ) = { 
                ( 3 0 ) = { 33.0350 }  ( 3 1 ) = { 33.1695 }
                ( 4 0 ) = { 35.4850 }  ( 4 1 ) = { 35.6295 }
                ( 5 0 ) = { 37.9350 }  ( 5 1 ) = { 38.0895 }
                ( 6 0 ) = { 40.3850 }  ( 6 1 ) = { 40.5495 }
                ( 7 0 ) = { 42.8350 }  ( 7 1 ) = { 43.0095 } 
              }
    ( 1 1 ) = { 
                ( 3 2 ) = { 33.3040 }
                ( 4 2 ) = { 35.7740 }
                ( 5 2 ) = { 38.2440 }
                ( 6 2 ) = { 40.7140 }
                ( 7 2 ) = { 43.1840 } 
              }