我应该如何提高这段代码的缓存友好性？(3D复数组到2D矩阵的一维数组)

Question

我正在开发一个处理多声道音频处理的Android应用程序。STFT函数产生三维复数组形式的复频率-时间表示，其维数为nFrames by nFreq的nChannels。

然而，在下一步，我需要执行盲源分离，它的运行时受益于我将每个频率段的通道和帧移动到矩阵中。目前，在读取STFTin条目时，该代码对缓存不友好。有没有办法让它对缓存更友好呢？

    Complex[][] temp = new Complex[nFrames][nChannels];
    Complex[][] tempConj = new Complex[nFrames][nChannels];

    X = new Array2DRowFieldMatrix[nFreqs];
    Xcopy = new Array2DRowFieldMatrix[nFreqs];
    Xconj = new Array2DRowFieldMatrix[nFreqs];
    Y = new Array2DRowFieldMatrix[nFreqs];
    for (int f = 0; f < nFreqs; f++) {
        for (int t = 0; t < this.nFrames; t++) {
            for (int c = 0; c < this.nChannels; c++) {
                temp[t][c] = STFTin[c][t][f];
                tempConj[t][c] = STFTin[c][t][f].conjugate();
                //STFTin is nChannels by nFrames by nFreq
        }
        X[f] = new Array2DRowFieldMatrix<>(temp);
        Xconj[f] = new Array2DRowFieldMatrix<>(tempConj);
        Xcopy[f] = (Array2DRowFieldMatrix<Complex>) X[f].copy();
        Y[f] = (Array2DRowFieldMatrix<Complex>) X[f].copy();
    }

Answer 1

正如我在注释中提到的，您可能希望采用现有的缓存友好的矩阵转置算法来实现您的逻辑(因为这是您正在做的大部分工作-矩阵转置)。其他人已经充实了关于最优块大小、循环排序和如何考虑边缘情况(例如不规则形状的矩阵)的实现细节，因此调整现有代码将导致代码更快、错误更少。

然而，如果你真的需要推出自己的解决方案，以下是我会尝试的。首先，重新排序循环，使频率和通道循环位于frames for-loop内。您直接从频率子数组访问元素，并将它们存储在通道子数组中，因此在将它们加载到缓存中时，充分利用这两个元素非常重要；这就是为什么我们需要在外部保持帧循环的原因。

接下来，根据缓存大小的一定比例对数组访问进行分块。这样，我们操作的所有频率和通道子数组都可以同时存在于缓存中-我们不会通过一次从频率数组中读取太多值来清除通道子数组，反之亦然。有一些方法可以计算这个大小，但老实说，运行它并计时会更可靠、更快-当事情停止变得更快时，停止增加块大小。

粗略的代码概要如下：

Complex[][][] temp = new Complex[nFreqs][nFrames][nChannels];
Complex[][][] tempConj = new Complex[nFreqs][nFrames][nChannels];

int blockSizeF = 1 << 2;  // Increase these until you see no speedup
int blockSizeC = 1 << 3;

X = new Array2DRowFieldMatrix[nFreqs];
Xcopy = new Array2DRowFieldMatrix[nFreqs];
Xconj = new Array2DRowFieldMatrix[nFreqs];
Y = new Array2DRowFieldMatrix[nFreqs];
for (int t = 0; t < this.nFrames; t++) {
    for (int fBlock = 0; fBlock < nFreqs; fBlock += blockSizeF) {
        for (int cBlock = 0; cBlock < this.nChannels; cBlock += blockSizeC) {
            for (int f = fBlock; f < fBlock + blockSizeF; f++) {
                for (int c = cBlock; c < cBlock + blockSizeC; c++) {
                    temp[f][t][c] = STFTin[c][t][f];
                    tempConj[f][t][c] = STFTin[c][t][f].conjugate();
                    //STFTin is nChannels by nFrames by nFreq
                }
            }
        }
    }
}
for (int f = 0; f < nFreqs; f++) {
    X[f] = new Array2DRowFieldMatrix<>(temp[f]);
    Xconj[f] = new Array2DRowFieldMatrix<>(tempConj[f]);
    Xcopy[f] = (Array2DRowFieldMatrix<Complex>) X[f].copy();
    Y[f] = (Array2DRowFieldMatrix<Complex>) X[f].copy();
}

通常情况下，blockSizeF和blockSizeC是相同的，但在这种情况下，每次读取STFTin的频率子数组时，您都要对temp和tempConj这两个单独的通道子数组执行两次写入。这意味着你想要一个比频率更大的通道块大小--可能是2的一个因素，也可能是sqrt(2)的一个因素--我不太确定。我会想象它会是这两个中的一个，所以我会进行实验，找出最好的。无论采用哪种方式，您都可能希望将块大小四舍五入到最接近的2的幂(或者至少是2的大幂的倍数)，以帮助与缓存线或页面边界对齐。

但是请注意，blockSizeF和blockSizeC 必须分别是nFreqs和nChannels的因子。有一些方法可以绕过这一规定，但它复杂、速度慢且容易出错。通常，只需填充矩阵，并在转换后剥离多余的部分，通常会更容易。