2
$\begingroup$

I have transferred my MATLAB Lanczos solver for symmetric eigenvalue solvers to C++ with the help of Intel MKL and MTL4 libraries. I have some wrapper templates for MKL routines. However during the iterations, the results of my C++ implementation starts to deviate from the values that MATLAB finds. Up to some point in the iterations, the differences seem to be small but after about 9-10 iterations the results start to deviate.

I suspect that this is due to the full re-orthogonalization that I used in the Lanczos solver. The strange thing is that I am using a simple Gram-Schmidt orthogonalization in MATLAB, however use of the same operations in C++ results in the difference in the \alpha and \beta coefficients of my problem. MATLAB uses double precision as far as I understand and in C++ I also use double precision by default. Intel MKL also uses double precision. However, it seems that the round-off errors introduced in the intermediate computations severely affect the performance of the C++ code, at least this is my guess apart from orthogonalizations. Any ideas to recover this strange problem?

$\endgroup$
  • 1
    $\begingroup$ Are the MKL results the same after every iteration? Are you using a parallel version? Rounding errors depend on the order of operations and threaded versions of certain algorithms will give you different results each time you run them. $\endgroup$ – Lagerbaer Jan 20 '12 at 17:36
  • 1
    $\begingroup$ Sorry if this is obvious: you are using modified Gram-Schmidt, not ordinary Gram-Schmidt, right? The more straightforward version is numerically unstable; see en.wikipedia.org/wiki/Gram-Schmidt_process#Numerical_stability $\endgroup$ – David Ketcheson Jan 21 '12 at 12:03
3
$\begingroup$

I believe that this is the same question that you recently posed to the MUMPS list here, so I am assuming that the C++ code you refer to is a wrapper for MUMPS.

To be honest, what language you're computing in is irrelevant as long as you are using the same floating point standard and the same algorithms. It is also not clear to me that your final answer is changing significantly; if you are worried about variations in the tridiagonal matrix, then make sure that every algorithmic detail is the same between your two implementations.

| cite | improve this answer | |
$\endgroup$
1
$\begingroup$

I ran into a similar sort of issue when comparing results between SciPy/NumPy and MATLAB. (I was using the results from SciPy/NumPy as reference results in unit tests.) MATLAB uses its own version of LAPACK and BLAS. In order to get bit-for-bit agreement, you'd have to use the same libraries as MATLAB, especially since algorithms can change slightly (bug fixes, changes in ordering conventions) from version to version of LAPACK.

I doubt this is the exact problem you're facing (my use case was a QR factorization), but it should explain the minor discrepancies you're seeing, and probably not the magnification of those discrepancies unless you're otherwise using the exact same algorithm in both implementations, and that algorithm is numerically stable.

| cite | improve this answer | |
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.