this is the same as f03_exercise, but with C++.
like the following can be performed.
Matrix A(size), B(size), C(size);
A.set(-1.0, 1.0, 5555);
B.set(-2.0, 2.0, 7777);
A = A.transpose();
B = B.transpose();
C = 0.0;
std::cout << "A:\n";
A.show();
std::cout << "B:\n";
B.show();
C += A;
C -= A;
std::cout << "C: \n";
C.show();
C = A*B;
std::cout << "C = A*B\n";
C.show();
C /= B; // C = C * B^-1
std::cout << "C /= B(= A)\n";
C.show();
if (C == A)
std::cout << "C == A\n";
if (A != B)
std::cout << "A != B\n";
Matrix Ainv = A.inverse();
std::cout << "Ainv = A^-1\n";
Ainv.show();
$ make
$ ./a.out
A:
6.71751e-01 -5.11412e-01 5.03268e-01 -5.15220e-01
-1.42823e-01 -4.56335e-01 1.01004e-01 -7.82845e-01
-7.78005e-01 9.79638e-01 7.61025e-01 -6.15416e-01
-8.88567e-01 -5.99848e-01 3.96229e-01 -8.35100e-01
B:
0.00000e+00 -1.32831e-01 -3.12948e-01 1.46513e+00
-1.44351e+00 0.00000e+00 1.67454e-01 -1.87368e+00
1.24615e-01 -1.73444e+00 0.00000e+00 1.59954e-03
-7.00869e-01 -3.26451e-01 1.82880e+00 0.00000e+00
C:
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
C = A*B
1.16204e+00 -7.93924e-01 -1.23809e+00 1.94323e+00
1.21998e+00 9.93457e-02 -1.46338e+00 6.45931e-01
-8.87953e-01 -1.01571e+00 -7.17952e-01 -2.97419e+00
1.50056e+00 -2.96589e-01 -1.34960e+00 -1.77311e-01
C /= B(= A)
6.71751e-01 -5.11412e-01 5.03268e-01 -5.15220e-01
-1.42823e-01 -4.56335e-01 1.01004e-01 -7.82845e-01
-7.78005e-01 9.79638e-01 7.61025e-01 -6.15416e-01
-8.88567e-01 -5.99848e-01 3.96229e-01 -8.35100e-01
C == A
maximum error is: 2.71408e+00
A != B
Ainv = A^-1
5.24346e-01 5.24421e-01 1.30336e-02 -8.24710e-01
-3.50690e-01 7.29534e-01 6.14144e-01 -9.20109e-01
1.20070e+00 -2.07358e+00 2.73936e-01 1.00117e+00
2.63679e-01 -2.06586e+00 -3.25030e-01 8.15982e-01
A*Ainv (should equal to unit matrix)
1.00000e+00 -2.22045e-16 5.55112e-17 2.77556e-16
2.77556e-17 1.00000e+00 0.00000e+00 0.00000e+00
8.32667e-17 2.22045e-16 1.00000e+00 0.00000e+00
5.55112e-17 0.00000e+00 0.00000e+00 1.00000e+00
tr(A) = 1.41341e-01
det(A) = 6.91491e-01
det(Ainv) = 1.44615e+00
det(A)*det(Ainv) (should equal to one) = 1.00000e+00
trivial test: C = C - C (=0)
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
0.00000e+00 0.00000e+00 0.00000e+00 0.00000e+00
C == 0.0
maximum error is: 2.00000e+00
D != 2.0
performance comparison of C = A/B between Fortran 2003 vs C++:
matrix size: 1024x1024
compiler: intel compiler 2018u0
time measurement function: C++: std::chrono::system_clock::now(), Fortran 2003: system_clock()
-- Fortran 2003 --
$ OMP_NUM_THREADS=16 KMP_AFFINITY=compact ./a.out 1024
size: 1024
time[s]: 9.99890000000000
-- C++ --
$ OMP_NUM_THREADS=16 KMP_AFFINITY=compact ./a.out 1024
size: 1024
time[s]: 9.08929