Difference in Matrix operations with very small values on Eigen (C++) vs Fortran 90

Question

Why am I getting slightly different results when I do the inverse of a matrix and multiply it by the original matrix? C++ (Eigen Library) vs Fortran (Using old Microsoft Developer Studio, Fortran PowerStation 4.0, and using MSIMSL module)

Code in C++

    #include <Eigen/Dense>
        #include <iostream>

        int main(){

        MatrixXd test_matrix
    {
        {1.3e-17, -3.2e-17, -4.1e-17, 0.2e-17},
        {-2.3e-17, 0.2e-17, 1.1e-17, -2.2e-17},
        {3.3e-17, -3.2e-17, -3.1e-17, -7.2e-17},
        {-4.3e-17, 0.2e-17, 6.1e-17, 0.2e-17},
    };
        
        cout << "Test Matrix: " << endl;
    cout << test_matrix << endl;
    cout << " " << endl;
    
    cout << "Inverse of Test Matrix: " << endl;
    cout << test_matrix.inverse() << endl;
    cout << " " << endl;

    cout << "Inverse * Original: " << endl;
    cout << test_matrix * test_matrix.inverse() << endl;
    
        cout << "Determinant of Test Matrix: " << endl;
    cout << test_matrix.determinant() << endl;
    
    return 0;
        }

Code in Fortran

    USE MSIMSL
    INTEGER M
    DOUBLE PRECISION A1(4,4),A2(4,4),B(4,4)
        OPEN(unit=1, file='output')

  1 FORMAT(16(E20.5,1X))


    M=4
    A1(1,1)=1.3e-17
    A1(1,2)=-3.2e-17
    A1(1,3)=-4.1e-17
    A1(1,4)=0.2e-17

    A1(2,1)=-2.3e-17
    A1(2,2)=0.2e-17
    A1(2,3)=1.1e-17
    A1(2,4)=-2.2e-17

    A1(3,1)=3.3e-17
    A1(3,2)=-3.2e-17
    A1(3,3)=-3.1e-17
    A1(3,4)=-7.2e-17

    A1(4,1)=-4.3e-17
    A1(4,2)=0.2e-17
    A1(4,3)=6.1e-17
    A1(4,4)=0.2e-17

C       returns inverse of matrix into A2
    CALL DLINRG(M,A1,M,A2,M)
C       B = A1 * A2
        CALL DMRRRR(M,M,A1,M,M,M,A2,M,M,M,B,M)

        DO 777 I=1,4
        WRITE(1,1)(A1(I,K),K=1,4)
 777    CONTINUE
        WRITE(1,*) '       '

        DO 888 I=1,4
        WRITE(1,1)(A2(I,K),K=1,4)
 888    CONTINUE
        WRITE(1,*) '       '

        DO 999 I=1,4
        WRITE(1,1)(B(I,K),K=1,4)
 999    CONTINUE
        WRITE(1,*) '       '

        END

Output in C++

Test Matrix:
 1.3e-17 -3.2e-17 -4.1e-17    2e-18
-2.3e-17    2e-18  1.1e-17 -2.2e-17
 3.3e-17 -3.2e-17 -3.1e-17 -7.2e-17
-4.3e-17    2e-18  6.1e-17    2e-18

Inverse of Test Matrix:
-1.46582e+16  -4.1296e+16  1.23181e+16  3.85461e+15
-2.41195e+16  1.97719e+16 -7.36473e+15 -2.35196e+16
-9.81172e+15 -2.92629e+16  9.21779e+15  1.97601e+16
 8.22593e+15 -1.51155e+16 -8.93865e+15  3.71206e+15

Inverse * Original:
           1  2.42861e-17            0  1.14492e-16
 1.38778e-16            1 -2.77556e-17  2.77556e-17
           0  4.44089e-16            1            0
 1.38778e-17  2.42861e-17            0            1

Determinant of Test Matrix:
-3.3674e-66

Output in Fortran

A1:
.13000E-16  -.32000E-16  -.41000E-16   .20000E-17
-.23000E-16   .20000E-17   .11000E-16  -.22000E-16
.33000E-16  -.32000E-16  -.31000E-16  -.72000E-16
-.43000E-16   .20000E-17   .61000E-16   .20000E-17
A2 = Inverse of A1:
-.14658E+17  -.41296E+17   .12318E+17   .38546E+16
-.24119E+17   .19772E+17  -.73647E+16  -.23520E+17
-.98117E+16  -.29263E+17   .92178E+16   .19760E+17
.82259E+16  -.15116E+17  -.89386E+16   .37121E+16
A2 * A1
.10000E+01   .21511E-15   .31225E-16  -.13878E-15
.55511E-16   .10000E+01  -.27756E-16   .13878E-16
.00000E+00  -.22204E-15   .10000E+01   .11102E-15
.33307E-15  -.69389E-17  -.19082E-15   .10000E+01

Calculating inverse if totally fine but: Why is there a difference when I multiply A1 to its inverse in C++ and Fortan. Does it relate to the fact that determinant of my test_matrix is too low (close to zero)?

I would expect that multiplication results will be relatively the same. Any suggestions? After some research i found this discussion that did not really help me, cause he is getting matrix with inf values.

The results are indeed very close up to the limit for double precision. The `A2 * A1` product both for C++ and for FORTRAN is very close to identity matrix. — Dima Chubarov, Jan 25 '23 at 18:05
Also the magnitude of the values is not of much influence on precision in your case since all values are more or less of the same magnitude. It would've resulted in rounding errors only if the elements of the original matrix were of different order of magnitude, i.e. some very small, while others close to 1. — Dima Chubarov, Jan 25 '23 at 18:11
Do you have a particular reason for `A1(1,1)=1.3e-17` instead of the (very different) `A1(1,1)=1.3d-17`? — francescalus, Jan 25 '23 at 18:48
No, I don't. Have almost no experience with Fortran and its syntax at all. — Yerlan Amir, Jan 25 '23 at 19:27
Indeed, got different results after using d-17 instead of e-17, but still mismatch between C++ and Fortran outputs. Thanks for the help, anyway! — Yerlan Amir, Jan 25 '23 at 19:39
Order of operations matters in floating point math. See [Is floating point math broken?](https://stackoverflow.com/questions/588004/is-floating-point-math-broken) The libraries may simply implement the order differently. Plus, depending on your compiler settings, it may change the order, too, or use fused multiply-add (FMA), see [Difference in gcc -ffp-contract options](https://stackoverflow.com/questions/43352510/difference-in-gcc-ffp-contract-options), meaning different compilers can cause different rounding errors, too. — Homer512, Jan 25 '23 at 21:03

Difference in Matrix operations with very small values on Eigen (C++) vs Fortran 90

0 Answers0