LU factorization numerical techniques

[NASA_SimDeveloper]

It is my understanding that a LU factorization numerically runs until the diagonal reaches values less than epsilon. I am not a math expert so correct me if I am wrong. Is there a way for me to get this epsilon? I am running my simulation on a cpu using a different library and I want the results to be calculated the same and close. I am profiling our simulation on the cpu and gpu code to compare timing. I plan on using magma_dgetfr_gpu(). I looked through the src code for epsilon but I could not find it.

[mgates3]

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