c++ - 理论和实践矩阵乘法

我的系统:

系统规格:
Intel core2duo E4500 3700g内存L2高速缓存2M x64 fedora 17

我如何测量人字拖鞋

好吧，我使用papi库(读取硬件性能计数器)来测量我的代码的flops和mflops。它返回实时处理时间，flops以及最终的flops /处理时间，等于MFLOPS.library使用硬件计数器来计算浮点数指令或浮点运算以及“总周期”以获得包含触发器和MFLOPS的最终结果。

我的计算内核

我使用了三个循环矩阵矩阵乘法(方阵)和三个嵌套循环，它们在其内部循环中对一维数组进行了一些运算。

第一个内核MM

    float a[size][size];
    float b[size][size];
    float c[size][size];

 start_calculate_MFlops();

for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            for (int k = 0; k < size; **k+=1**) {
                *c[i][j]=c[i][j]+a[i][k] * b[k][j];*
                     }
            }
 }
  stop_calculate_MFlops();

具有一维数组的第二个内核

    float d[size];
    float e[size];
    float f[size];
    float g[size];
    float r = 3.6541;

 start_calculate_MFlops();

for (int i = 0; i < size; ++i) {
    for (int j = 0; j < size; ++j) {
        for (int k = 0; k < size; ++k) {
            d[k]=d[k]+e[k]+f[k]+g[k]+r;
        }
    }
}

stop_calculate_MFlops();

我对拖鞋的了解

矩阵矩阵乘法(MM)在其内部循环中执行2次运算(此处为浮点运算)，并且由于存在3个循环，因此针对大小X进行迭代，因此从理论上讲，MM的总触发器为2 * n ^ 3。

在第二个内核中，我们有3个循环，在最内层的循环中，我们有1d数组进行一些计算。此循环中有4个浮点运算。因此，理论上我们的总触发器为4 * n ^ 3触发器

我知道我们上面计算出的触发器与实际机器中的触发器并不完全相同。在实际的机器中，还有其他操作，例如加载和存储，这将导致理论上的失败。

问题？:

当我在第二个内核理论触发器中使用1d数组时，

是
通过执行代码并进行测量得到的相同或左右的翻牌
它。实际上，当我使用一维数组触发器等于
最内层循环乘以n ^ 3，但是当我使用第一个内核MM时
使用2d数组的理论触发器是2n ^ 3，但是当我运行代码时
，测量值比理论值高太多
大约4+(矩阵乘法的最内层循环中的2个运算)* n ^ 3 + = 6n ^ 3。
我只用下面的代码更改了最内层循环中的矩阵乘法行:

A[i][j]++;

在3个嵌套循环中，此代码的理论触发器为1次操作* n ^ 3 = n ^ 3再次运行代码时，结果比预期的2+(最内部循环的1次操作)* n高^ 3 = 3 * n ^ 3

大小为512X512的矩阵的示例结果:

实时:1.718368处理时间:1.227672总flpops:
807,107,072 MFLOPS:657.429016

实时:3.608078处理时间:3.042272播放总数:
807,024,448 MFLOPS:265.270355

理论上的翻牌:2 * 512 * 512 * 512 = 268,435,456

实测触发器= 6 * 512 ^ 3 = 807,107,072

3个嵌套循环中1d数组操作的样本结果

实时:1.282257进行时间:1.155990总flpops:
536,872,000 MFLOPS:464.426117

理论底数: 4n ^ 3 = 536,870,912

被测翻牌: 4n ^ 3 = 4 * 512 ^ 3 + overheads(other operation？)= 536,872,000

我找不到上述行为的任何原因吗？
我的假设是真的吗？

希望使其比之前的描述更加简单。

实际上，我的意思是通过执行代码来衡量真实的失败。

代码:

 void countFlops() {

    int size = 512;
    int itr = 20;
    float a[size][size];
    float b[size][size];
    float c[size][size];
/*  float d[size];
    float e[size];
    float f[size];
    float g[size];*/
        float r = 3.6541;

    float real_time, proc_time, mflops;
    long long flpops;
    float ireal_time, iproc_time, imflops;
    long long iflpops;
    int retval;

    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            a[j][j] = b[j][j] = c[j][j] = 1.0125;
        }
    }

/*  for (int i = 0; i < size; ++i) {
                d[i]=e[i]=f[i]=g[i]=10.235;
        }*/

    if ((retval = PAPI_flops(&ireal_time, &iproc_time, &iflpops, &imflops))
            < PAPI_OK) {
        printf("Could not initialise PAPI_flops \n");
        printf("Your platform may not support floating point operation event.\n");
        printf("retval: %d\n", retval);
        exit(1);
    }
    for (int i = 0; i < size; ++i) {
        for (int j = 0; j < size; ++j) {
            for (int k = 0; k < size; k+=16) {
                c[i][j]=c[i][j]+a[i][k] * b[k][j];
            }
        }
    }

/*  for (int i = 0; i < size; ++i) {
    for (int j = 0; j < size; ++j) {
        for (int k = 0; k < size; ++k) {
            d[k]=d[k]+e[k]+f[k]+g[k]+r;
        }
    }
    }*/

    if ((retval = PAPI_flops(&real_time, &proc_time, &flpops, &mflops))
            < PAPI_OK) {
        printf("retval: %d\n", retval);
        exit(1);
    }
    string flpops_tmp;
    flpops_tmp = output_formatted_string(flpops);
    printf(
            "calculation: Real_time: %f Proc_time: %f Total flpops: %s MFLOPS: %f\n",
            real_time, proc_time, flpops_tmp.c_str(), mflops);

}

谢谢

最佳答案

如果您需要计算操作数-您可以制作简单的类，其作用类似于浮点值并收集统计信息。它可以与内置类型互换。

LIVE DEMO:

#include <boost/numeric/ublas/matrix.hpp>
#include <boost/operators.hpp>
#include <iostream>
#include <ostream>
#include <utility>
#include <cstddef>
#include <vector>

using namespace boost;
using namespace std;

class Statistic
{
    size_t ops = 0;
public:
    Statistic &increment()
    {
        ++ops;
        return *this;
    }
    size_t count() const
    {
        return ops;
    }
};

template<typename Domain>
class Profiled: field_operators<Profiled<Domain>>
{
    Domain value;
    static vector<Statistic> stat;
    void stat_increment()
    {
        stat.back().increment();
    }
public:
    struct StatisticScope
    {
        StatisticScope()
        {
            stat.emplace_back();
        }
        Statistic &current()
        {
            return stat.back();
        }
        ~StatisticScope()
        {
            stat.pop_back();
        }
    };
    template<typename ...Args>
    Profiled(Args&& ...args)
        : value{forward<Args>(args)...}
    {}
    Profiled& operator+=(const Profiled& x)
    {
        stat_increment();
        value+=x.value;
        return *this;
    }
    Profiled& operator-=(const Profiled& x)
    {
        stat_increment();
        value-=x.value;
        return *this;
    }
    Profiled& operator*=(const Profiled& x)
    {
        stat_increment();
        value*=x.value;
        return *this;
    }
    Profiled& operator/=(const Profiled& x)
    {
        stat_increment();
        value/=x.value;
        return *this;
    }
};
template<typename Domain>
vector<Statistic> Profiled<Domain>::stat{1};

int main()
{
    typedef Profiled<double> Float;
    {
        Float::StatisticScope s;
        Float x = 1.0, y = 2.0, res = 0.0;
        res = x+y*x+y;
        cout << s.current().count() << endl;
    }
    {
        using namespace numeric::ublas;
        Float::StatisticScope s;
        matrix<Float> x{10, 20},y{20,5},res{10,5};
        res = prod(x,y);
        cout << s.current().count() << endl;
    }
}

输出为:

3
2000

附言您的矩阵循环不适合缓存，结果是very inefficient。

点对点

int size = 512;
float a[size][size];

这不是合法的C++代码。 C++不支持VLA。