Trung tâm cá nhân
sách gpt4 ăn đã đi

c - 用于大型密集矩阵乘法的循环平铺/阻塞

In lại 作者：太空狗 更新时间：2023-10-29 16:32:45 27 4
mua khóa gpt4 giày nike

我想知道是否有人可以告诉我如何有效地使用循环平铺/循环阻塞进行大型密集矩阵乘法。我正在用 1000x1000 矩阵做C = AB。我按照 Wikipedia 上的循环平铺示例进行操作，但使用平铺得到的结果比不使用时更差。

http://en.wikipedia.org/wiki/Loop_tiling

http://software.intel.com/en-us/articles/how-to-use-loop-blocking-to-optimize-memory-use-on-32-bit-intel-architecture

我在下面提供了一些代码。由于缓存未命中，天真的方法非常慢。 transpose 方法在缓冲区中创建 B 的转置。这种方法给出最快的结果(矩阵乘法为 O(n^3)，转置为 O(n^2)，因此转置至少快 1000 倍)。没有阻塞的 wiki 方法也很快，不需要缓冲区。阻塞方法较慢。阻塞的另一个问题是它必须多次更新 block 。这对线程/OpenMP 来说是一个挑战，因为如果不小心，它可能会导致竞争条件。

我应该指出，使用 AVX 并修改转置方法可以比 Eigen 更快地获得结果。但是，我使用 SSE 的结果比 Eigen 慢一点，所以我认为我可以更好地使用缓存。

编辑:我想我知道我想做什么。它来自于 1969 年的 Cannon 算法。
http://en.wikipedia.org/wiki/Matrix_multiplication#Communication-avoiding_and_distributed_algorithms

我需要进行 block 矩阵乘法，并让每个线程处理 C thay vì MỘT 的子矩阵。 strong> 和 B。例如，如果我将矩阵分成四个 block 。我会这样做:

+-+ +-+ +-+ +-+ +-+ +-+
| | | | | |
| C1 C2 | | A1 A2 | | B1 B2 |
| | = | | x | |
| C3 C4 | | A3 A4 | | B3 B4 |
| | | | | |
+-+ +-+ +-+ +-+ +-+ +-+


thread 1: C1 = A1B1 + A2B3
thread 2: C2 = A1B2 + A2B4
thread 3: C3 = A3B1 + A4B3
thread 4: C4 = A3B2 + A4B4

这消除了竞争条件。我得考虑一下。

void matrix_mult_naive(const float*A , const float* B, float* C, const int N, const int M, const int K) {
    #pragma omp parallel for
    for(int i=0; i
        for(int j=0; j
            float tmp = 0;
            for(int l=0; l
                tmp += A[M*i+l]*B[K*l+j];
            }
            C[K*i + j] = tmp;
        }
    }
}
void matrix_mult_transpose(const float*A , const float* B, float* C, const int N, const int M, const int K) {
    float *B2 = (float*)aligned_malloc(M*K*sizeof(float), 32);
    transpose(B, B2, M, K, 1);
    #pragma omp parallel for
    for(int i=0; i
        for(int j=0; j
            float tmp = 0;
            for(int l=0; l
                tmp += A[M*i+l]*B2[M*j+l];
            }
            C[K*i + j] = tmp;
        }
    }
    aligned_free(B2);
}

void matrix_mult_wiki(const float*A , const float* B, float* C, const int N, const int M, const int K) {
    for(int i=0; i
        for(int j=0; j
            C[K*i + j] = 0;
        }  
    }
    #pragma omp parallel for
    for(int i=0; i
        for(int l=0; l
            float a = A[M*i+l];
            for(int j=0; j
                C[K*i + j] += a*B[K*l+j];
            }
        }
    }
}

void matrix_mult_wiki_block(const float*A , const float* B, float* C, const int N, const int M, const int K) {
   const int block_size = 8; //I have tried several different block sizes
   for(int i=0; i
       for(int j=0; j
           C[K*i + j] = 0;
       }
    }
    for(int l2=0; l2
        for(int j2=0; j2
        #pragma omp parallel for
            for(int i=0; i
                for(int l=l2; l
                    for(int j=j2; j
                        C[K*i + j] += A[M*i+l]*B[K*l+j];
                    }
                }
            }
        }
    }
}

câu trả lời hay nhất

我得到的最佳结果是再添加一个 vì 循环来阻塞 N，然后重新排列循环。我还提升了循环不变代码，但编译器的优化器应该会自动执行此操作。 block 大小应该是缓存行大小除以 sizeof(float)。这比转置方法快 50%。

如果您只能选择 AVX 或阻塞之一，使用 AVX 扩展(vfmadd###ps Và haddps)仍然会快得多。考虑到您已经在测试 block 大小是否是 64/sizeof(float) == 16 个 float == 两个 256 位 AVX 寄存器的倍数，因此最好和直接使用两者。

  • 转置:1,816,522 个刻度
  • 平铺:892,431(快 51%)
  • AVX+平铺:230,512(快 87%)

平铺:

void matrix_mult_wiki_block(const float*A , const float* B, float* C,
                            const int N, const int M, const int K) {
    const int block_size = 64 / sizeof(float); // 64 = common cache line size
    for(int i=0; i
        for(int j=0; j
            C[K*i + j] = 0;
        }
    }
    for (int i0 = 0; i0 < N; i0 += block_size) {
        int imax = i0 + block_size > N ? N : i0 + block_size;

        for (int j0 = 0; j0 < M; j0 += block_size) {
            int jmax = j0 + block_size > M ? M : j0 + block_size;

            for (int k0 = 0; k0 < K; k0 += block_size) {
                int kmax = k0 + block_size > K ? K : k0 + block_size;

                for (int j1 = j0; j1 < jmax; ++j1) {
                    int sj = M * j1;

                    for (int i1 = i0; i1 < imax; ++i1) {
                        int mi = M * i1;
                        int ki = K * i1;
                        int kij = ki + j1;

                        for (int k1 = k0; k1 < kmax; ++k1) {
                            C[kij] += A[mi + k1] * B[sj + k1];
                        }
                    }
                }
            }
        }
    }
}

至于 Cannon 引用，SUMMA algorithm是一个更好的遵循。

万一其他人正在优化高瘦乘法({~1e9 x 50} x {50 x 50}，我是如何在这里结束的)，转置方法在性能上几乎与阻塞方法相同 N=18( float )。 N=18 是一种病态情况(比 17 或 19 更糟)，我不太清楚导致这种情况的缓存访问模式。所有较大的 N 都通过阻塞方法得到改进。

关于c - 用于大型密集矩阵乘法的循环平铺/阻塞，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/15829223/

27 4 0
Đề xuất bài viết: c - 强制 64 位长 double ？
Đề xuất bài viết: c - 无锁队列中的内存屏障使用
Đề xuất bài viết: c - 无法理解格式字符串利用代码
Đề xuất bài viết: c - 如何在 Clion 中自动生成 .h 文件的函数头？
太空狗
Hồ sơ

Tôi là một lập trình viên xuất sắc, rất giỏi!

Nhận phiếu giảm giá taxi Didi miễn phí
Phiếu giảm giá taxi Didi
Các bài viết phổ biến trên toàn bộ trang web
Chứng chỉ ICP Bắc Kinh số 000000
Hợp tác quảng cáo: 1813099741@qq.com 6ren.com
Xem sitemap của VNExpress