After almost 40 years of active research on parallel scientific computing, researchers are still struggling with the implementation of Gaussian elimination. The problem centers on how to implement pivoting, that is, row and/or column exchanges for stability. If concurrency and data movement were not issues, the alternatives would be 1) complete pivoting, or Gaussian elimination with complete pivoting (GECP), which involves relocating the largest element into the pivot position with row and column exchanges, or 2) Gaussian elimination with partial pivoting (GEPP), which means moving the largest element in the current column into the pivot position. GECP requires O(n³) comparisons, thus approximately doubling the effort, but is stable. GEPP requires O(n²) comparisons and can be unstable, but only if you are very unlucky. For reasons no one completely understands, examples of matrices where GEPP is unstable have been found using the imagination of computational scientists, but never by accident in computational experience. For that reason, GEPP is considered stable “in practice.” Unfortunately, in distributed memory computing environments, the data movement involved for even the modest number of comparisons and row exchanges required by GEPP can be burdensome. This paper plots the pivoting overhead in a central processing unit/graphics processing unit (CPU/GPU) architecture as ranging from 40 percent of factorization time on small matrices to 17 percent on matrices of size 10,000. Although this is an improvement for large matrices, it still includes a much greater overhead than suggested by the O(n²) comparisons for pivoting versus the O(n³) operations for factorization.

As an alternative, the authors propose a heuristic called the partial random butterfly transformation (PRBT), which comes from an idea by Parker [1] and Parker and Pierce [2]. The techniques in constructing PRBT are based on those used to construct fast Fourier transforms, thus allowing the use of similar parallel algorithms and data structures. The authors propose transforming a matrix A to be factored by Gaussian elimination into A_r = U^TAV where U and V are PRBTs. Then A_r is factored by Gaussian elimination without pivoting. Care is taken to make certain the condition numbers of U and V are bounded by modest constants, thereby insuring that A_r has a condition number not too different from that of A. The authors perform exhaustive tests on a number of matrices from the literature, including three well-known ones that are known to perform poorly with GEPP. Iterative refinement is performed when the residual is unsatisfactory. In no case was more than one step of iterative refinement needed, and in all cases PRBT produced solutions with small residuals. The operational complexity of PRBT is slightly greater than GEPP, theoretically O(n² log n), but in practice O(n²) with a larger constant. Using the CPU/GPU architecture described in the paper, PRBT was much faster for matrices of sizes 1,000 to 3,000. However, for matrices larger than 9,000, the improvement over GEPP was only around 10 percent. The authors believe that asymptotically, PRBT and GEPP will be very close.

From this paper, we can conclude that PRBT works well as a proposed alternative to GEPP, but its advantage fades as the matrices get larger. As a method for keeping Gaussian elimination stable, there is no theory explaining why PRBT works at all. In fact, there is no known theoretical reason why this method should be better or worse than Gaussian elimination without pivoting. This paper proposes an interesting alternative to partial pivoting, but significant theoretical work is needed to understand it properly.