LU factorization numerical techniques

 Posts: 8
 Joined: Wed Apr 17, 2019 11:50 am
LU factorization numerical techniques
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.
Re: LU factorization numerical techniques
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 nonzero diagonal (if it has a zero, then A is singular). There is no tolerance level epsilon in the algorithm.
Solving AX = B using LU is just forward and backsubstitution, 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.11e16 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 1e19, while tol*epsilon is 33e16:
mark
Code: Select all
% LU factorization, but lacking pivoting.
for k = 1 : n1
% 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.
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.11e16 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 1e19, while tol*epsilon is 33e16:
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.40e19 ok
200 1  (  ) 0.90 ( 0.01) 4.76e19 ok
300 1  (  ) 1.93 ( 0.01) 6.07e19 ok
400 1  (  ) 2.26 ( 0.02) 2.91e19 ok
500 1  (  ) 3.15 ( 0.03) 3.20e19 ok
600 1  (  ) 4.28 ( 0.03) 3.32e19 ok
700 1  (  ) 5.08 ( 0.05) 3.79e19 ok
800 1  (  ) 5.68 ( 0.06) 3.36e19 ok
900 1  (  ) 5.88 ( 0.08) 2.24e19 ok
1000 1  (  ) 6.55 ( 0.10) 3.80e19 ok
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.
Reason: MAGMA doesn't divide by epsilon as LAPACK does.