c++ - MSVC 编译器/arch :AVX2/fp:fast breaks C++ matrix inversion algorithm

标签 c++ visual-studio visual-c++ g++ matrix-inverse

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

在 Windows 上的 Visual Studio 2017 中，它也可以正常工作。但是，在升级到 Visual Studio 2019 后，我得到一个包含全零的矩阵。 MSVC 是否损坏或者还有什么问题？

到目前为止，我已经测试过将 ==0.0/!=0.0 替换为大于/小于公差值，但这也不起作用。

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

我还在代码中标记了Windows中错误返回的地方。非常感谢任何帮助。

编辑:我刚刚发现错误行为的原因是/arch:AVX2 与/fp:fast 的组合。但我不明白为什么。/arch:SSE、/arch:SSE2 和/arch:AVX 与/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());
}

最佳答案

似乎是与文字 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());
}

关于c++ - MSVC 编译器/arch :AVX2/fp:fast breaks C++ matrix inversion algorithm，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/69532196/

上一篇：python - 在 google ortools 中添加析取约束

下一篇：azure-devops - WhiteSource 无法扫描您的存储库。请导航到组织设置 > WhiteSource 选项卡连接到您的 WhiteSource 帐户

相关文章：

c++ - `==` 运算符的问题

.net - 调试 C# COMVisible 类方法时，Visual Studio 2010 不会在断点处停止

visual-studio - MSBUILD : error MSB1008: Only one project can be specified

c++ - 转换运算符 + 转换构造函数 = 不直观的行为？

c++ - 为什么结果输出为1

c++ - 如何在 C++ 中播放声音

c++ - 从被移出的对象 move

c++ - 使用 miniz 提取目录内的文件

c# - AXML 和 XAML 之间的区别？

C++:如何获取在同一类的不同函数中定义的元素？

©2024 IT工具网  联系我们