Gaussian elimination assume system of equations is represented by an nx(n+1) matrix. assume indices start at 1 (not 0) 1 for i = 1 to n // i is the pivot row 2 parfor j = 1 to n // j is the elimination row 3 r = a[j,i]/a[i,i] 4 parfor k = i to n+1 // k is the elimination column 5 if j != i 6 a[j,k] -= r*a[i,k] 7 for j = 1 to n 8 s[j] = a[j,n+1]/a[j,j] This is simple to write, and appropriate for a GPU or a SIMD machine like the old CM-2, but it does more work than it needs to: multiplications line 3 n^2 line 6 (n-1)\sum_{2<=i<=n+1}{i} = (n-1)(n^2+3n)/2 = (n^3+2n^2-3n)/2 line 8 n --------------------- total (n^3+4n^2-n)/2 additions line 6 (n-1)\sum_{2<=i<=n+1}{i} = (n-1)(n^2+3n)/2 = (n^3+2n^2-3n)/2 The algorithm you learned in linear algebra does less total work 1 for i = 1 to n // i is the pivot row 2 for k = i+1 to n+1 3 a[i,k] /= a[i,i] 4 a[i,i] = 1.0 5 for j = i+1 to n // j is the elimination row 6 p = a[j,i] 7 a[j,i] = 0.0 8 for k = i+1 to n+1 // k is the elimination column 9 a[j,k] -= p*a[i,k] // a is now upper triangular // back-substitute to find solution vector s 10 for j = n downto 1 11 s[j] = a[j,n+1] 12 for i = n downto j+1 13 s[j] -= a[j,i]*s[i] This requires multiplications: line 3 (n^2+n)/2 // pivot row scaling line 9 \sum_{1<=i<=n-1}{i(i+1)} = \sum_{1<=i<=n-1}{i} + \sum_{1<=i<=n-1}{i^2} = (n^2-n)/2 + (2n^3-3n^2+n)/6 = (n^3-n)/3 // eliminations line 13 (n^2-n)/2 // back-substitutions ------------------- total (n^3+3n^2-n)/3 additions: line 9 \sum_{1<=i<=n-1}{i(i+1)} = (n^3-n)/3 // eliminations line 13 (n^2-n)/2 // back-substitutions ------------------- total (2n^3+3n^2-5n)/6 ==> asymptotically 2/3 of the work how might we parallelize this? problem: load balance inherent in problem induced by environment (# processors, competing jobs) banded v. cyclic task allocation is locality an issue? bag of tasks additional issue: pivoting for numerical stability, pivot row should be the one in which the leading non-zero element has the largest absolute value (and what do you do if the leading element *is* a zero?)