带有 /arch:AVX2 /fp:fast breaks C++ 矩阵求逆算法的 MSVC 编译器

Question

我有矩阵求逆算法和一个用于测试的 19x19 矩阵。运行它在 Windows 的 Linux 子系统中使用 g++ 生成正确的倒矩阵。

在 Visual Studio 2017 年的 Windows 中，它曾经也能正常工作。但是，在我升级到 Visual Studio 2019 之后，我得到了一个包含所有零的矩阵。是 MSVC 坏了还是哪里出了问题？

到目前为止，我已经测试过用 greater/smaller 替换 ==0.0/!=0.0（在某些情况下由于四舍五入可能不安全）而不是公差值，但这确实也不行。

（我用 g++ 得到的）预期结果是 {5.26315784E-2,-1.25313282E-2, ... -1.25000000E-1}。我使用 Visual Studio 2019 v142，Windows SDK 10.0.19041.0，版本 x64 进行编译，得到 {0,0,...0}。使用 Debug x64 我得到一个不同的（也是错误的）结果：所有矩阵条目都是 -1.99839720E18。我使用 C++17 语言标准，/O2 /Oi /Ot /Qpar /arch:AVX2 /fp:fast /fp:except-。我有一个支持 AVX2 的 i7-8700K。

我在代码中Windows也标出了returns错误的地方。非常感谢任何帮助。

编辑：我刚刚发现错误行为的原因是 /arch:AVX2 与 /fp:fast 结合使用。但我不明白为什么。 /arch:SSE、/arch:SSE2 和 /arch:AVX with /fp:fast 可以正常工作，但 /arch:AVX2 不行。此处不同的舍入或操作顺序如何触发完全不同的行为？

#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());
}

Answer 1

似乎是与文字 0 的比较破坏了优化的算法。特别是第一个杀死它，因为正确的操作结果需要正好为 0。

一般来说，与其将浮点值与文字零进行比较，不如将绝对值与一个非常小的常量进行比较。

这似乎有效：

#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());
}

带有 /arch:AVX2 /fp:fast breaks C++ 矩阵求逆算法的 MSVC 编译器

MSVC compiler with /arch:AVX2 /fp:fast breaks C++ matrix inversion algorithm

c++

g++

matrix-inverse

visual-studio

visual-c++