Сложение транспонированных матриц (последовательный алгоритм)

/// <summary>
/// Перегрузка операции сложения матриц с присваиванием
/// </summary>
/// <typeparam name="T"> тип элементов матрицы </typeparam>
/// <param name="A"> слагаемое-результат </param>
/// <param name="B"> второе слагаемое </param>
template <class T>
void operator += (CSR_Matrix<T> &A, CSR_Matrix<T> B)
{
    assert(A.n() == B.n() && A.m() == B.m());

    int _N = A.n(), _M = A.m(), _NNZ = 0;

    vector <T> values;
    vector <int> cols, rowptr(1);

    vector <T> values_A = A.values(), values_B = B.values();
    vector <int> cols_A = A.cols(), cols_B = B.cols();
    vector <int> rowptr_A = A.rowptr(), rowptr_B = B.rowptr();

    for (int i = 0; i < _N; ++i)
    {
        int left_A = rowptr_A[i], right_A = rowptr_A[i + 1];
        int left_B = rowptr_B[i], right_B = rowptr_B[i + 1];

        int ptr_A = left_A, ptr_B = left_B;

        while (ptr_A < right_A && ptr_B < right_B)
        {
            int col_A = cols_A[ptr_A], col_B = cols_B[ptr_B];
            T value_A = values_A[ptr_A], value_B = values_B[ptr_B];

            if (col_A == col_B)
            {
                if (value_A + value_B)
                {
                    values.push_back(value_A + value_B);
                    cols.push_back(col_A);  ++_NNZ;
                }

                ++ptr_A; ++ptr_B;
            }
            else if (col_A < col_B)
            {
                values.push_back(value_A);
                cols.push_back(col_A);

                ++ptr_A; ++_NNZ;
            }
            else
            {
                values.push_back(value_B);
                cols.push_back(col_B);

                ++ptr_B; ++_NNZ;
            }
        }

        while (ptr_A < right_A)
        {
            values.push_back(values_A[ptr_A]);
            cols.push_back(cols_A[ptr_A]);

            ++ptr_A; ++_NNZ;
        }

        while (ptr_B < right_B)
        {
            values.push_back(values_B[ptr_B]);
            cols.push_back(cols_B[ptr_B]);

            ++ptr_B; ++_NNZ;
        }

        rowptr.push_back(_NNZ);
    }

    A = CSR_Matrix<T>(_N, _M, _NNZ, values, cols, rowptr);
}