c++ 如何计算浮点型的大型矩阵的逆？

我使用高斯消元法来求逆矩阵，矩阵是12 x 12。
当找到一个小矩阵时，结果是正确的，但当矩阵变得更大时，结果是不正确的。
我想这是一个浮点小数精度问题。
有没有一种方法可以使用单精度浮点小数来计算精确的逆矩阵？

template<typename T = float>
inline bool inverse(const MATRIX<Mtx_Rt, Mtx_Rt, T>& matrix, MATRIX<Mtx_Rt, Mtx_Rt, T>& out, void* pTmp) {
       if (matrix.GetNumOfRow() != matrix.GetNumOfCol())
        return false;
    MATRIX<Mtx_Rt, Mtx_Rt, T> augmentedMatrix = MATRIX<Mtx_Rt, Mtx_Rt, T>(matrix.GetNumOfRow(), matrix.GetNumOfCol() * 2, (T*)pTmp);
    memset(pTmp, 0, augmentedMatrix._sizeof());
    //[matrix | identity]
    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        for (int j = 0; j < matrix.GetNumOfCol(); j++)
            augmentedMatrix(i, j) = matrix(i, j);
        augmentedMatrix(i, i + matrix.GetNumOfRow()) = 1.0;
    }
    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        int pivotRow = i;
        for (int j = i + 1; j < matrix.GetNumOfRow(); j++)
            if (augmentedMatrix(j, i) > augmentedMatrix(pivotRow, i))
                pivotRow = j;
        for (int k = 0; k < 2 * matrix.GetNumOfRow(); k++) {
            T tmp = augmentedMatrix(i, k);
            augmentedMatrix(i, k) = augmentedMatrix(pivotRow, k);
            augmentedMatrix(pivotRow, k) = tmp;
        }
        for (int j = 0; j < matrix.GetNumOfRow(); j++)
            if (j != i) {
                T ratio;
                if (augmentedMatrix(i, i) != 0.0)
                    ratio = augmentedMatrix(j, i) / augmentedMatrix(i, i);
                else
                    return false;
                for (int k = 0; k < 2 * matrix.GetNumOfRow(); k++)
                    augmentedMatrix(j, k) -= ratio * augmentedMatrix(i, k);
            }
    }
    for (int i = 0; i < matrix.GetNumOfRow(); i++) {
        T diagonalElement = augmentedMatrix(i, i);
        for (int j = 0; j < 2 * matrix.GetNumOfRow(); j++)
            augmentedMatrix(i, j) /= diagonalElement;
    }
    for (int i = 0; i < matrix.GetNumOfRow(); i++)
        for (int j = 0; j < matrix.GetNumOfCol(); j++)
            out(i, j) = augmentedMatrix(i, j + matrix.GetNumOfRow());
    return true;
}
//The following is the test data
0x000001DEBFA44B70     0.00901737064    0.00851101335    0.00903018098    0.00852382369   1.71711508e-06   1.07878850e-05   1.07854212e-05  -5.71605915e-06       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44BA0     0.00851101428    0.00921893492    0.00835073180    0.00905865245   0.000142186589   6.97691357e-05   6.96772186e-05   0.000119451797       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44BD0     0.00903018005    0.00835073087    0.00916963257    0.00849018339  -0.000121859746  -5.83902292e-05  -5.83116416e-05  -0.000103386286       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44C00     0.00852382462    0.00905865245    0.00849018339    0.00902501121   1.86097168e-05   5.91026037e-07   5.80158485e-07   2.17815614e-05       0.00000000       0.00000000       0.00000000       0.00000000  
0x000001DEBFA44C30    1.71710417e-06   0.000142186575  -0.000121859761   1.86097132e-05      0.141640529       0.00000000       0.00000000       0.00000000    0.00899084471    0.00772444298   -0.00682374742   -0.00809014961  
0x000001DEBFA44C60    1.07878877e-05   6.97691430e-05  -5.83902220e-05   5.91027856e-07       0.00000000      0.122075945      0.122108623     0.0703755319     0.0109803220    0.00996958558   -0.00884070992   -0.00985144638  
0x000001DEBFA44C90    1.07854303e-05   6.96772186e-05  -5.83116343e-05   5.80155756e-07       0.00000000      0.122108623      0.122141339     0.0704018399     0.0109902890    0.00998063106   -0.00885062385   -0.00986028183  
0x000001DEBFA44CC0   -5.71606142e-06   0.000119451797  -0.000103386286   2.17815686e-05       0.00000000     0.0703755319     0.0704018399      0.110758379    0.00697547942    0.00992908329   -0.00901818927   -0.00606458541  
0x000001DEBFA44CF0        0.00000000       0.00000000       0.00000000       0.00000000    0.00899084471     0.0109803220     0.0109902890    0.00697547942      0.157142207      0.148548171      0.145359069      0.136765048  
0x000001DEBFA44D20        0.00000000       0.00000000       0.00000000       0.00000000    0.00772444159    0.00996958464    0.00998063106    0.00992908329      0.148548156      0.157674491      0.136261418      0.145387754  
0x000001DEBFA44D50        0.00000000       0.00000000       0.00000000       0.00000000   -0.00682374695   -0.00884070992   -0.00885062385   -0.00901818834      0.145359069      0.136261433      0.156348825      0.147251189  
0x000001DEBFA44D80        0.00000000       0.00000000       0.00000000       0.00000000   -0.00809014961   -0.00985144544   -0.00986028090   -0.00606458541      0.136765033      0.145387754      0.147251174      0.155873895

字符串

使用模板数值工具包[TNT] https://math.nist.gov/tnt/只需要很少的包含（没有库来构建）和下面的行代码来使用LU分解计算逆。它是一个模板类，所以你可以使用它与浮点。我用它来从10多年的两倍矩阵大小到2000。TNT还包括QR，特征值求解器和SVD。

#include <tnt/jama_lu.h>
template <class T>
TNT::Array2D<T> inverse(const TNT::Array2D<T>& M)
{
    TNT::Array2D<T> LU(M);
    int nline = dim1();
    int ncol = dim2();
    TNT::Array2D<T> input(nline, ncol);
    for (int i = 0; i < nline; i++)
    {
        for (int j = 0; j < ncol; j++)
        {
            if (i == j) input[i][j] = 1.;
            else input[i][j] = 0.;
        }
    }
    return LU.solve(input);
}

字符串

展开查看全部

c++ 如何计算浮点型的大型矩阵的逆？

1条答案

相关问题

热门标签

最新问答