Repeated single-subscript contraction is a (rather inefficient) way to sum the components of a matrix. For example, this program:

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

struct add_asgn_double {
    double operator() (double & z, double const y) { return z += y; }
};

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

    form const frm_a {blt{1,4},blt{1,4},blt{1,4}};
    con_matrix<double> mat_a {sub_matrix_random (frm_a,rand_double(0.0, 1.0))};
    cout << "\n matrix a: " << mat_a;

    {
        cout << "\n\n three contractions on one subscript each: \n";

        vec_bool const which_b1 {false,false,true};
        con_matrix<double> mat_b1 (contraction<add_asgn_double> (mat_a, which_b1));
        cout << "\n matrix b1: " << mat_b1;

        vec_bool const which_b2 {false,true};
        con_matrix<double> mat_b2 (contraction<add_asgn_double> (mat_b1, which_b2));
        cout << "\n matrix b2: " << mat_b2;
        
        vec_bool const which_b3 {true};
        con_matrix<double> mat_b3 (contraction<add_asgn_double> (mat_b2, which_b3));
        cout << "\n matrix b3: " << mat_b3;
    } {
        cout << "\n\n one contraction on two subscripts, ";
        cout << "\n then one contraction on one subscript: \n";

        vec_bool const which_c1 {false,true,true};
        con_matrix<double> mat_c1 (contraction<add_asgn_double> (mat_a, which_c1));
        cout << "\n matrix c1: " << mat_c1;

        vec_bool const which_c2 {true};
        con_matrix<double> mat_c2 (contraction<add_asgn_double> (mat_c1, which_c2));
        cout << "\n matrix c2: " << mat_c2;
    } {
        cout << "\n\n one contraction on one subscript, ";
        cout << "\n then one contraction on two subscripts: \n";

        vec_bool const which_d1 {false,false,true};
        con_matrix<double> mat_d1 (contraction<add_asgn_double> (mat_a, which_d1));
        cout << "\n matrix d1: " << mat_d1;

        vec_bool const which_d2 {true,true};
        con_matrix<double> mat_d2 (contraction<add_asgn_double> (mat_d1, which_d2));
        cout << "\n matrix d2: " << mat_d2;
    } {
        cout << "\n\n one contraction on all three subscripts: \n";

        vec_bool const which_e {true,true,true};
        con_matrix<double> mat_e (contraction<add_asgn_double> (mat_a, which_e));
        cout << "\n matrix e: " << mat_e;
    }

    return 0;
}

 matrix a: 
    ( 1 1 1 ) = 0.0850    ( 1 1 2 ) = 0.8916    ( 1 1 3 ) = 0.1897
    ( 1 2 1 ) = 0.3980    ( 1 2 2 ) = 0.7435    ( 1 2 3 ) = 0.5604
    ( 1 3 1 ) = 0.8096    ( 1 3 2 ) = 0.5117    ( 1 3 3 ) = 0.9951
    ( 2 1 1 ) = 0.9666    ( 2 1 2 ) = 0.4261    ( 2 1 3 ) = 0.6530
    ( 2 2 1 ) = 0.9615    ( 2 2 2 ) = 0.8580    ( 2 2 3 ) = 0.2940
    ( 2 3 1 ) = 0.4146    ( 2 3 2 ) = 0.5149    ( 2 3 3 ) = 0.7898
    ( 3 1 1 ) = 0.5443    ( 3 1 2 ) = 0.0936    ( 3 1 3 ) = 0.4323
    ( 3 2 1 ) = 0.8449    ( 3 2 2 ) = 0.7728    ( 3 2 3 ) = 0.1919
    ( 3 3 1 ) = 0.7804    ( 3 3 2 ) = 0.1813    ( 3 3 3 ) = 0.5791

 three contractions on one subscript each: 

 matrix b1: 
    ( 1 1 ) = 1.1663    ( 1 2 ) = 1.7019    ( 1 3 ) = 2.3164
    ( 2 1 ) = 2.0457    ( 2 2 ) = 2.1135    ( 2 3 ) = 1.7193
    ( 3 1 ) = 1.0702    ( 3 2 ) = 1.8096    ( 3 3 ) = 1.5409
 matrix b2: 
    ( 1 ) = 5.1846    ( 2 ) = 5.8785    ( 3 ) = 4.4206
 matrix b3: 
    ( ) = 15.4838

 one contraction on two subscripts, 
 then one contraction on one subscript: 

 matrix c1: 
    ( 1 ) = 1.8236    ( 2 ) = 2.6144    ( 3 ) = 1.8963
 matrix c2: 
    ( ) = 6.3343

 one contraction on one subscript, 
 then one contraction on two subscripts: 

 matrix d1: 
    ( 1 1 ) = 1.1663    ( 1 2 ) = 1.7019    ( 1 3 ) = 2.3164
    ( 2 1 ) = 2.0457    ( 2 2 ) = 2.1135    ( 2 3 ) = 1.7193
    ( 3 1 ) = 1.0702    ( 3 2 ) = 1.8096    ( 3 3 ) = 1.5409
 matrix d2: 
    ( ) = 4.8207

 one contraction on all three subscripts: 

 matrix e: 
    ( ) = 1.5222