C语言 优化循环展开以通过高斯约当方法找到矩阵的逆

trnvg8h3  于 2023-01-29  发布在  其他
关注(0)|答案(2)|浏览(125)

我尝试应用循环展开来通过Gauss Jorda方法找到矩阵的逆矩阵,当矩阵的大小非常大,无法放入缓存时,减少内存访问的次数(瓶颈)。我让它运行得更快,但我得到的结果是错误的,我不知道为什么。

for(k=0;k<size;k++)                                  
{                                                       
    
    pivot=original[k][k];
    for(j=0;j<size;j++)                             
    {
        original[k][j]/=pivot;                                  
        inverse[k][j]/=pivot;                                   
    }   
            
    for(i=0;i<size;i++)                                 
    {
           
        if(i!=k)
        {
            pivot = original[i][k];                                 
            for(j=0;j<size;j++)                         
            {                                                                                   
                original[i][j] -= original[k][j]*pivot;                     
                inverse[i][j] -= inverse[k][j]*pivot;   
                                            
            }
        }
            
    }     
}

我希望通过接收内存访问次数来加快问题的执行速度。
我的循环展开版本如下:

for(k=0;k<sizeOfMatrix;k += 2)                                   
{                                                           
    pivot=original[k][k];
    for(j=0;j<sizeOfMatrix;j++)                             
    {
        original[k][j]/=pivot;                                  
        inverse[k][j]/=pivot;                                   
    }   
           
    for(i=0;i<sizeOfMatrix;i++)                                 
    {
       
        if(i!=k && i!=k+1)
        {
            pivot = original[i][k];                                 
            for(j=0;j<sizeOfMatrix;j++)                         
            {                                                                                   
                original[i][j] -= original[k][j]*pivot;                     
                inverse[i][j] -= inverse[k][j]*pivot;                              
            }
        }
               
        if(i!=k+1)
        {
            pivot=original[k][k];
            for(j=0;j<sizeOfMatrix;j++)                             
            {
                original[k+1][j]/=pivot;                                    
                inverse[k+1][j]/=pivot;                                 
            }
                              
            pivot = original[i][k+1];                                   
            for(j=0;j<sizeOfMatrix;j++)                         
            {                                                                                   
                original[i][j] -= original[k+1][j]*pivot;                       
                inverse[i][j] -= inverse[k+1][j]*pivot;                              
            }
        }        

    }     
        
}
watbbzwu

watbbzwu1#

现在,让我们先不考虑原始代码的数值适应性问题,以及展开是否是您应该尝试的练习。您展示了通过重写 middle 循环体来展开 outer 循环的尝试。尽管可以沿着与转换结果类似的方式构建与原始代码等效的计算。它不会构成一个展开,因为循环展开不会改变在循环迭代范围内循环体中执行的操作的顺序。2你所做的转换确实会产生这样的改变。
但在您的例子中,转换后的代码就是错误的,至少因为

  • sizeOfMatrix为奇数时,它在k循环迭代期间执行越界数组访问,其中k == sizeOfMatrix - 1
  • 它在if(i!=k+1)块的开始处使用了错误的主元值。
  • 原则上,正确的值将被存储在x1M4 N1 x,但是在x1M5 N1 x的迭代中还没有计算出所需的值。
  • X1 M6 N1 X块对于除了X1 M9 N1 X之外的 * 每个 * X1 M8 N1 X缩放原始矩阵和逆矩阵的行X1 M7 N1 X,但是它们每个应当仅缩放一次。

在不进一步考虑您所尝试的转换是否有助于性能改进的情况下,我观察到,如果您想执行真正的展开,那么尝试展开的循环是错误的。如果真正的展开有任何好处,那么它将来自于展开最内层的循环。例如:

for (k = 0; k < size; k++) {                                                       
    // ...
    for (i = 0; i < size; i++) {
        // ...
            for (j = 0; j < size ^ 1; j += 2) {
                // two rows per iteration                                                                                  
                original[i][j] -= original[k][j] * pivot;
                original[i][j + 1] -= original[k][j + 1] * pivot;
                inverse[i][j] -= inverse[k][j] * pivot;
                inverse[i][j + 1] -= inverse[k][j + 1] * pivot;
            }
            if (size & 1) {
                // a single additional row when needed
                assert(j == size - 1);
                original[i][j] -= original[k][j] * pivot;
                inverse[i][j] -= inverse[k][j] * pivot;
            }
        // ...
    }
}

但是正如我在注解中所说的,如果你的编译器期望这样的展开来加速程序,那么当你在编译时启用优化时,它会自动这样做,而且编译器对这些事情的判断是相当不错的。
实际上,手动优化很容易使性能恶化,而不是提高性能,这可能是固有的,也可能是因为手动优化会干扰编译器自己执行更有效的优化。
如果最佳速度是你的目标,那么你可能应该使用一个调优的第三方实现,比如lapack或者atlas。
如果您坚持自己滚动,那么您可能只需改进算法(甚至坚持使用G-J)就可以获得更高的加速比,但您可能应该给予其中一些返回以获得更高的数值稳定性。
1.在每次外循环迭代时,选择当前消除列中具有最大量值的可用行,并在必要时将其与当前行交换。这将最小化大数除以非常小的数(尤其是零)的发生率。
1.(在上面的(1)之后),仅计算当前透视表右侧列中元素的更新。原始矩阵的其他元素已经对结果产生了所有影响。(但您不能使用逆矩阵来简化这种方法。)
如果你不对(1)做一些改变,那么你的逆矩阵将在一些实际上可逆的矩阵上失败。
沿着如此,例如:

// Temporary space for row swaps
element_type *tmp = malloc(size * sizeof(*tmp));
if (!tmp) {
    abort();  // you probably want something more graceful here
}

for (int pivot_col = 0; pivot_col < size; pivot_col++) {

    // Choose the remaining row that maximizes the magnitude of the pivot element

    int pivot_row = pivot_col;
    element_type pivot = fabs(original[pivot_row][pivot_col]);  // or fabsf() or fabsl() ...
                                                     
    for (row = pivot_row + 1; row < size; row++) {
        element_type test_el = fabs(original[row][pivot_col]);

        if (test_el > pivot) {
            pivot_row = row;
            pivot = test_el;
        }
    }

    if (pivot_row != pivot_col) {

        // This assumes array-of-array structure.
        // There are better options for array-of-pointers structure.

        int current_row = pivot_col;

        // Only the tails of the original rows, starting at pivot_col, need to be swapped
        memcpy(tmp, &original[current_row][pivot_col], (size - pivot_col) * sizeof(*tmp));
        memcpy(&original[current_row][pivot_col], &original[pivot_row][pivot_col], (size - pivot_col) * sizeof(*tmp));
        memcpy(&original[pivot_row][pivot_col], tmp, (size - pivot_col) * sizeof(*tmp));

        // The full inverse rows need to be swapped
        memcpy(tmp, inverse[current_row], sizeof(inverse[current_row]));
        memcpy(inverse[current_row], inverse[pivot_row], sizeof(inverse[current_row]));
        memcpy(inverse[pivot_row], tmp, sizeof(inverse[current_row]));
    }
    pivot_row = pivot_col;
    pivot = original[pivot_row][pivot_col];  // get the correctly-signed pivot value

    // All columns of the inverse row need to be scaled, but
    // only the tail columns of the original row need to be.
    int col = 0;

    for (; col < pivot_col + 1; col++) {
        inverse[pivot_row][col] /= pivot;                                   
    }   
    for (; col < size; col++) {
        original[pivot_row][col] /= pivot;                                  
        inverse[pivot_row][col] /= pivot;                                   
    }   
            
    // Compute the effect on the other rows of canceling the pivot column
    for (int row = 0; row < size; row++) {
        if (row == pivot_row) continue;

        element_type cancel_val = original[row][pivot_col];

        // where col < pivot_col, original[pivot_row][col] is (logically) zero
        // where col == pivot_col, we are (logically) setting original[row][col] to zero
        for (col = 0; col <= pivot_col; col++) {
            inverse[row][col] -= inverse[pivot_row][col] * cancel_val;   
        }
        // we need actually to write updates to the remaining columns of the original
        for(; col < size; col++) {
            original[row][col] -= original[pivot_row][col] * cancel_val;
            inverse[row][col] -= inverse[pivot_row][col] * cancel_val;
        }
    }     
}

free(tmp);

这是从您的第一个代码中得到的灵感,改进了名称,并实现了算法改进(1)和(2)。

pbpqsu0x

pbpqsu0x2#

你已经找到了问题的解决方案。因为目标是学术性的。目标是减少执行时间,减少内存访问的次数,在我的情况下减少一半。
1.我修改了第k行,我用它的主元划分了该行的元素
1.用新的行k修改矩阵的行k+1
1.修改行k+1,将行的元素除以其透视。
1.我修改了值为k和k + 1的其他行

for (k = 0; k < size; k += 2)
    {
        pivot = original[k][k];
        for (j = 0; j < size; j++)
        {
            original[k][j] /= pivot;
            inverse[k][j] /= pivot;

        }
        pivot = original[k + 1][k];
        for (i = 0; i < size; i++)
        {
            original[k + 1][i] -= original[k][i] * pivot;
            inverse[k + 1][i] -= inverse[k][i] * pivot;
        }
        pivot = original[k+1][k+1];

        for (j = 0; j < size; j++)
        {
            original[k+1][j] /= pivot;
            inverse[k+1][j] /= pivot;
        }
        for (i = 0; i < size; i++)
        {
            if (i != k && i != k + 1)
            {
                pivot = original[i][k];
                
                    for (j = 0; j < size; j++)
                    {
                        original[i][j] -= original[k][j] * pivot;
                        inverse[i][j] -= inverse[k][j] * pivot;
                    }
            }

            if (i != k + 1)
            {
                pivot = original[i][k+1];
                for (j = 0; j < size; j++)
                {
                    original[i][j] -= original[k + 1][j] * pivot;
                    inverse[i][j] -= inverse[k + 1][j] * pivot;

                }
            }
        }
    }

相关问题