LU factorization numerical techniques

Open discussion for MAGMA library (Matrix Algebra on GPU and Multicore Architectures)
Post Reply
NASA_SimDeveloper
Posts: 8
Joined: Wed Apr 17, 2019 11:50 am

LU factorization numerical techniques

Post by NASA_SimDeveloper » Thu May 02, 2019 1:34 pm

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
Posts: 893
Joined: Fri Jan 06, 2012 2:13 pm

Re: LU factorization numerical techniques

Post by mgates3 » Sat May 04, 2019 1:28 am

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
Last edited by mgates3 on Sat May 04, 2019 1:33 am, edited 1 time in total.
Reason: MAGMA doesn't divide by epsilon as LAPACK does.

Post Reply