How to calculate Matrix exponential with Tailor series PARALLEL using MPI c++

Question

Need to calculate Matrix exponential with Tailor series with MPI matrix is small 3 x 3 for example

Meanwhile

vector<vector<double>> matrixMult(const vector<vector<double>>& A, const vector<vector<double>>& B) {
    int n = A.size();
    vector<vector<double>> C(n, vector<double>(n, 0));
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            for (int k = 0; k < n; k++)
                C[i][j] += A[i][k] * B[k][j];
    return C;
}
vector<vector<double>> matrixExp(const vector<vector<double>>& A) {
    int n = A.size();
    vector<vector<double>> E(n, vector<double>(n, 0));
    vector<vector<double>> T(n, vector<double>(n, 0));
    vector<vector<double>> localE(n, vector<double>(n, 0));

    int rank, size;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    for (int i = 0; i < n; i++)
        E[i][i] = 1;

    for (int i = 0; i < n; i++)
        localE[i][i] = 0;

    T = E;
    for (int j = 1; j <= rank; j++) 
    {
        T = matrixMult(T, A);
        T = matrixDiv(T, j);
    }
    localE = T;
    for (int i = rank + 1; i <= N; i += size) {
        for (int j = i; j < i + size; j++) {
            T = matrixMult(T, A);
            T = matrixDiv(T, j);
        }
        localE = matrixSum(localE, T);
    }
    MPI_Reduce(localE[0].data(), E[0].data(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    MPI_Reduce(localE[1].data(), E[1].data(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    MPI_Reduce(localE[2].data(), E[2].data(), n, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);
    return E;
}

But i dont know how to optimize this

for (int j = i; j < i + size; j++) {
            T = matrixMult(T, A);
            T = matrixDiv(T, j);
        }

Maybe it's impossible with this implementation scheme