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?