LU is a direct method, not an iterative method (see Matlab code below). It forms matrices L, U such that PA = LU (exactly, in exact arithmetic). L has unit diagonal, U has non-zero diagonal (if it has a zero, then A is singular). There is no tolerance level epsilon in the algorithm.

Code: Select all

```
% LU factorization, but lacking pivoting.
for k = 1 : n-1
% multipliers
for i = k+1 : n
A( i, k ) = A( i, k ) / A( k, k );
end
% update trailing matrix
for j = k+1 : n
for i = k+1 : n
A( i, j ) = A( i, j ) - A( i, k ) * A( k, j );
end
end
end
% L is unit lower triangle of A,
% U is upper triangle of A.
```

Solving AX = B using LU is just forward and back-substitution, again a direct method that would be exact in exact arithmetic.

However, there is a backwards error check. For solving AX = B,

error = norm( AX - B ) / (n * norm(A) * norm(X)) is O(epsilon),

where epsilon is unit roundoff (1.11e-16 for double). By default in MAGMA and LAPACK, we check that the error is < tol*epsilon with tol = 30. This error check is done in the tester (magma/testing/testing_zgesv.cpp) after calling gesv to check its result. Frequently, the error is much smaller than this tolerance. For instance, these tests achieve errors around 1e-19, while tol*epsilon is 33e-16:

Code: Select all

```
magma/testing> ./testing_dgesv -n 100:1000:100 --check
% N NRHS CPU Gflop/s (sec) GPU Gflop/s (sec) ||B - AX|| / N*||A||*||X||
%===============================================================================
100 1 --- ( --- ) 0.13 ( 0.01) 1.40e-19 ok
200 1 --- ( --- ) 0.90 ( 0.01) 4.76e-19 ok
300 1 --- ( --- ) 1.93 ( 0.01) 6.07e-19 ok
400 1 --- ( --- ) 2.26 ( 0.02) 2.91e-19 ok
500 1 --- ( --- ) 3.15 ( 0.03) 3.20e-19 ok
600 1 --- ( --- ) 4.28 ( 0.03) 3.32e-19 ok
700 1 --- ( --- ) 5.08 ( 0.05) 3.79e-19 ok
800 1 --- ( --- ) 5.68 ( 0.06) 3.36e-19 ok
900 1 --- ( --- ) 5.88 ( 0.08) 2.24e-19 ok
1000 1 --- ( --- ) 6.55 ( 0.10) 3.80e-19 ok
```

-mark