Gauss-Jordan消元法求矩陣的逆

　　最近信息安全做業要寫hill密碼，涉及到求矩陣的逆。不想用一個函數就用那些亂七八糟的庫，大概找了找又沒現成的，只能本身擼了。

　　寫的時候才發現線代都忘光了QAQ

　　

 1 //Gauss-Jordan Elimination method to get inverse matrix
 2 template <typename T=double>
 3 std::vector<std::vector<double>> matrixInversion(const std::vector<std::vector<T>>& a)  4 {  5     const static double eps = 1e-6;  6     std::vector<std::vector<double>> augma; //augmented matrix
 7     for (auto& vec : a)  8  {  9         std::vector<double> tmp; 10         std::transform(vec.begin(), vec.end(), std::back_inserter(tmp), [](const T interger) 11  { 12                 return static_cast<double>(interger); 13  }); 14  augma.push_back(std::move(tmp)); 15  } 16     int rank = a.size(); 17     //init 
18     for (int i = 0; i < rank; ++i) 19  { 20         augma[i].resize(rank * 2); 21         augma[i][rank + i] = 1; 22  } 23     for (int i = 0; i < rank; ++i) 24  { 25         if (::abs(augma[i][i]) < eps) 26  { 27             int j; 28             for (j = i + 1; j < rank; ++j) 29  { 30                 if (::abs(augma[j][i]) > eps) break; 31  } 32             if (j == rank) return {}; 33             for (int r = i; r < 2 * rank; ++r) 34  { 35                 augma[i][r] += augma[j][r]; 36  } 37  } 38         double ep = augma[i][i]; 39         for (int r = i; r < 2 * rank; ++r) 40  { 41             augma[i][r] /= ep; 42  } 43 
44         for (int j = i + 1; j < rank; ++j) 45  { 46             double e = -1 * (augma[j][i] / augma[i][i]); 47             for (int r = i; r < 2 * rank; ++r) 48  { 49                 augma[j][r] += e * augma[i][r]; 50  } 51  } 52  } 53     for (int i = rank - 1; i >= 0; --i) 54  { 55         for (int j = i - 1; j >= 0; --j) 56  { 57             double e = -1 * (augma[j][i] / augma[i][i]); 58             for (int r = i; r < 2 * rank; ++r) 59  { 60                 augma[j][r] += e * augma[i][r]; 61  } 62  } 63  } 64 
65     for (auto& vec : augma) 66  { 67         vec.erase(vec.begin(), vec.begin() + rank); 68  } 69 
70     return augma; 71 }