Selectively multiplying entries in a block matrix

Question

I know that the question I am about to ask is a tall order, but for the time being I have no solution ahead and I am hopping someone can shed some light into this.

I am working on a code that has a block matrix structure. I need to do some vector-matrix multiplication with sub-matrices and sub-vectors. However I have no idea where to even start with this.

I have a mock problem illustrating the issue:

#include <iostream>
#include <vector>
#include <array>

typedef unsigned int label;

// Matrix class
template <class Type, int length>
class matrixN
{
    public: 
       matrixN():
       v_(std::array<double, 16>{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16})
       {}

       const Type& operator[](const int& i ) const {return v_[i];}

        enum
        {
            size = length
        };

    private:
        std::array<Type, 16> v_;
};

// Vector class
template <class Type, int length>
class VectorN
{
    public: 
       VectorN():v_(std::array<double, 4>{1,2,3,4}){}
        enum
        {
            size = length
        };
        
       const Type& operator[](const int& i ) const {return v_[i];}

    private:
        std::array<Type, length> v_;
};

// InnerProduct for a 4x4 matrix and a 4x1 vector
template <class Type, int length>
std::array<Type, length>
operator&(const matrixN<Type, length>& t, const VectorN<Type, length>& v)
{
    std::array<Type, length> result{0};

    label i = 0;
    for (label row = 0; row < matrixN<Type, length>::size; row++)
    {
        Type& r = result[row];

        for (label col = 0; col < matrixN<Type, length>::size ; col++)
        {
            r += t[i]*v[col];
            i++;
        }
    }
    std::cout << "I am the original function" << std::endl;

    return result;
}

// Vector of matrices
template <class Type, int length>
class matrixFieldN
{
    public:
        matrixFieldN(const int& size):tf_(size){}

        label size() const{return tf_.size();}

        const matrixN<Type,length>& operator[](const int& i ) const {return tf_[i];}

    private:
        std::vector<matrixN<Type,length> > tf_;
};

// Vector of vectors
template <class Type, int length>
class vectorFieldN 
{
    public:
       vectorFieldN (const int& size):vf_(size){}

       label size()const {return vf_.size();}

        const VectorN<Type,length>& operator[](const int& i ) const {return vf_[i];}

    private:
        std::vector<VectorN <Type, length> > vf_;
};

// Inner product function for vectors of stuff (calls the individual &)
template <class Type, int length>
std::vector<std::array<Type, length>>
operator&(const matrixFieldN<Type, length>& tf, const vectorFieldN <Type, length>& vf)
{
    std::vector<std::array<Type, length>> result (tf.size());

    for(label i=0; i<tf.size(); i++)
    {
        result[i] = tf[i] & vf[i];
    }
    return result;
}

// Someway to slice the matrix

int main(int argCount, char *args[])
{
    matrixN<double, 4> t1;
    VectorN<double, 4> v1;

    auto test1 =t1 & v1;

    matrixFieldN<double,4> tf1 (5);
    vectorFieldN<double,4> vf1 (5);

    std::vector<std::array<double,4>> test2 = tf1 & vf1;

    return 0;
}

I have created a base matrix class which always initializes to 1...16, representing a 4x4 matrix. This class has a retrieval parameter to index each element and a length enum. A similar approach was made for the vector class but with elements from 1...4.

Afterwards, I have defined the inner product operation templated and defined for matrix-vector multiplication.

The last two classes are vectors of the previous two: a vector of matrices and a vector of vectors. The inner product function was defined to call the & operator for each element in the vector.

How can I now write something that will only perform sub-block multiplication: e.g., 3x3 multiplication {elements: 1...9} of the matrix with elements {1..3} of the vector. Or more generally, by selecting which entries should be multiplied:

Does anyone know a code doing something similar so that I can take a look?