I have matrix inversion algorithm and a 19x19 matrix for testing. Running it in Windows Subsystem for Linux with g++ produces the correct inverted matrix.
In Visual Studio 2017 on Windows it used to work as well. However after I upgraded to Visual Studio 2019, I get a matrix containing all zeros as result. Is MSVC broken or what else is wrong here?
So far I have tested replacing ==0.0/!=0.0 (that in some cases can be unsafe due to round-off) with greater/smaller than a tolerance value, but this does not work either.
The expected result (that I get with g++) is {5.26315784E-2,-1.25313282E-2, ... -1.25000000E-1}.
I compile with Visual Studio 2019 v142, Windows SDK 10.0.19041.0, Release x64, and I get {0,0,...0}. With Debug x64 I get a different (also wrong) result: all matrix entries are -1.99839720E18.
I use C++17 language standard, /O2 /Oi /Ot /Qpar /arch:AVX2 /fp:fast /fp:except-.
I have an i7-8700K that supports AVX2.
I also marked the place where it wrongly returns in Windows in the code. Any help is greatly appreciated.
EDIT: I just found out that the cause of the wrong behavior is /arch:AVX2 in combination with /fp:fast. But I don't understand why. /arch:SSE, /arch:SSE2 and /arch:AVX with /fp:fast work correctly, but not /arch:AVX2. How can different round-off or operation order here trigger entirely different behavior?
#include <iostream>
#include <string>
#include <vector>
using namespace std;
typedef unsigned int uint;
void println(const string& s="") {
cout << s << endl;
}
struct floatNxN {
uint N{}; // matrix size is NxN
vector<float> M; // matrix data
floatNxN(const uint N, const float* M) : N {N}, M(N*N) {
for(uint i=0; i<N*N; i++) this->M[i] = M[i];
}
floatNxN(const uint N, const float x=0.0f) : N {N}, M(N*N, x) { // create matrix filled with zeros
}
floatNxN() = default;
~floatNxN() = default;
floatNxN invert() const { // returns inverse matrix
vector<double> A(2*N*N); // calculating intermediate values as double is strictly necessary
for(uint i=0; i<N; i++) {
for(uint j=0; j< N; j++) A[2*N*i+j] = (double)M[N*i+j];
for(uint j=N; j<2*N; j++) A[2*N*i+j] = (double)(i+N==j);
}
for(uint k=0; k<N-1; k++) { // at iteration k==2, the content of A is already different in MSVC and g++
if(A[2*N*k+k]==0.0) {
for(uint i=k+1; i<N; i++) {
if(A[2*N*i+k]!=0.0) {
for(uint j=0; j<2*N; j++) {
const double t = A[2*N*k+j];
A[2*N*k+j] = A[2*N*i+j];
A[2*N*i+j] = t;
}
break;
} else if(i+1==N) {
return floatNxN(N);
}
}
}
for(uint i=k+1; i<N; i++) {
const double t = A[2*N*i+k]/A[2*N*k+k];
for(uint j=k; j<2*N; j++) A[2*N*i+j] -= A[2*N*k+j]*t;
}
}
double det = 1.0;
for(uint k=0; k<N; k++) det *= A[2*N*k+k];
if(det==0.0) {
return floatNxN(N);
}
for(int k=N-1; k>0; k--) {
for(int i=k-1; i>=0; i--) {
const double t = A[2*N*i+k]/A[2*N*k+k];
for(uint j=k; j<2*N; j++) A[2*N*i+j] -= A[2*N*k+j]*t;
}
}
floatNxN r = floatNxN(N);
for(uint i=0; i<N; i++) {
const double t = A[2*N*i+i];
for(uint j=0; j<N; j++) r.M[N*i+j] = (float)(A[2*N*i+N+j]/t);
}
return r;
}
string stringify() const { // converts matrix into string without spaces or newlines
string s = "{"+to_string(M[0]);
for(uint i=1; i<N*N; i++) s += ","+to_string(M[i]);
return s+"}";
}
};
int main() {
const float Md[19*19] = {
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
-30,-11,-11,-11,-11,-11,-11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0, 1, -1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,
0, -4, 4, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,
0, 0, 0, 1, -1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,
0, 0, 0, -4, 4, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,
0, 0, 0, 0, 0, 1, -1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1,
0, 0, 0, 0, 0, -4, 4, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1,
0, 2, 2, -1, -1, -1, -1, 1, 1, 1, 1,-2,-2, 1, 1, 1, 1,-2,-2,
0, -4, -4, 2, 2, 2, 2, 1, 1, 1, 1,-2,-2, 1, 1, 1, 1,-2,-2,
0, 0, 0, 1, 1, -1, -1, 1, 1,-1,-1, 0, 0, 1, 1,-1,-1, 0, 0,
0, 0, 0, -2, -2, 2, 2, 1, 1,-1,-1, 0, 0, 1, 1,-1,-1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,-1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0,-1, 1, 1,-1
};
floatNxN M(19, Md);
floatNxN Mm1 = M.invert();
println(Mm1.stringify());
}
It seems that it is the comparison against literal 0 that breaks the algorithm with the optimizations on. Particularly the 1st one kills it, because the result would need to be exactly 0 for the correct action.
Generally instead of comparing a floating point value against literal zero, it is better to compare the absolute value against a very small constant.
This seems to work:
#include <iostream>
#include <string>
#include <vector>
using namespace std;
typedef unsigned int uint;
void println(const string& s = "") {
cout << s << endl;
}
struct floatNxN {
const double epsilon = 1E-12;
uint N{}; // matrix size is NxN
vector<float> M; // matrix data
floatNxN(const uint N, const float* M) : N{ N }, M(N* N) {
for (uint i = 0; i < N * N; i++) this->M[i] = M[i];
}
floatNxN(const uint N, const float x = 0.0f) : N{ N }, M(N* N, x) { // create matrix filled with zeros
}
floatNxN() = default;
~floatNxN() = default;
floatNxN invert() const { // returns inverse matrix
vector<double> A(2 * N * N); // calculating intermediate values as double is strictly necessary
for (uint i = 0; i < N; i++) {
for (uint j = 0; j < N; j++) A[2 * N * i + j] = (double)M[N * i + j];
for (uint j = N; j < 2 * N; j++) A[2 * N * i + j] = (double)(i + N == j);
}
for (uint k = 0; k < N - 1; k++) { // at iteration k==2, the content of A is already different in MSVC and g++
if (fabs(A[2 * N * k + k]) < epsilon) { // comparing with 0 was the killer here
for (uint i = k + 1; i < N; i++) {
if (fabs(A[2 * N * i + k]) > epsilon) {
for (uint j = 0; j < 2 * N; j++) {
const double t = A[2 * N * k + j];
A[2 * N * k + j] = A[2 * N * i + j];
A[2 * N * i + j] = t;
}
break;
} else if (i + 1 == N) {
return floatNxN(N);
}
}
}
for (uint i = k + 1; i < N; i++) {
const double t = A[2 * N * i + k] / A[2 * N * k + k];
for (uint j = k; j < 2 * N; j++) A[2 * N * i + j] -= A[2 * N * k + j] * t;
}
}
double det = 1.0;
for (uint k = 0; k < N; k++) det *= A[2 * N * k + k];
if (fabs(det) < epsilon) {
return floatNxN(N);
}
for (int k = N - 1; k > 0; k--) {
for (int i = k - 1; i >= 0; i--) {
const double t = A[2 * N * i + k] / A[2 * N * k + k];
for (uint j = k; j < 2 * N; j++) A[2 * N * i + j] -= A[2 * N * k + j] * t;
}
}
floatNxN r = floatNxN(N);
for (uint i = 0; i < N; i++) {
const double t = A[2 * N * i + i];
for (uint j = 0; j < N; j++) r.M[N * i + j] = (float)(A[2 * N * i + N + j] / t);
}
return r;
}
string stringify() const { // converts matrix into string without spaces or newlines
string s = "{" + to_string(M[0]);
for (uint i = 1; i < N * N; i++) s += "," + to_string(M[i]);
return s + "}";
}
};
int main() {
const float Md[19 * 19] = {
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
-30,-11,-11,-11,-11,-11,-11, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8,
12, -4, -4, -4, -4, -4, -4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
0, 1, -1, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,
0, -4, 4, 0, 0, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1, 1,-1, 0, 0,
0, 0, 0, 1, -1, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,
0, 0, 0, -4, 4, 0, 0, 1,-1, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,
0, 0, 0, 0, 0, 1, -1, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1,
0, 0, 0, 0, 0, -4, 4, 0, 0, 1,-1, 1,-1, 0, 0,-1, 1,-1, 1,
0, 2, 2, -1, -1, -1, -1, 1, 1, 1, 1,-2,-2, 1, 1, 1, 1,-2,-2,
0, -4, -4, 2, 2, 2, 2, 1, 1, 1, 1,-2,-2, 1, 1, 1, 1,-2,-2,
0, 0, 0, 1, 1, -1, -1, 1, 1,-1,-1, 0, 0, 1, 1,-1,-1, 0, 0,
0, 0, 0, -2, -2, 2, 2, 1, 1,-1,-1, 0, 0, 1, 1,-1,-1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0,-1,-1, 0, 0,
0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0, 1,-1,-1, 1, 0, 0,
0, 0, 0, 0, 0, 0, 0,-1, 1, 0, 0, 1,-1, 1,-1, 0, 0, 1,-1,
0, 0, 0, 0, 0, 0, 0, 0, 0, 1,-1,-1, 1, 0, 0,-1, 1, 1,-1
};
floatNxN M(19, Md);
floatNxN Mm1 = M.invert();
println(Mm1.stringify());
}