一个完善的高斯消元法(C++)

基本上可以应付主元为0,不满秩等所有情况,同时可以获得零空间.
所用到的矩阵类长这样

 class matrix {
public:
 short r; //行
 short c; //列
 vector<vector<double> > x; //数据地址
 //构造函数 ,要求写出行数和列数
 matrix(const short ri = 1, const short ci = 1);
//重置大小
void matrix::reset(int ir, int ic) {
 if (ic > c) {
  auto temp = x.end();
  for (auto i = x.begin();i != temp;++i) {
   i->resize(ic);
  }
 }
 this->c = ic;
 if(ir>this->r){
  x.resize(ir);
  auto temp = x.end();
  for (auto i = x.begin();i != temp;++i) {
   i->resize(c);
  }
 }
 this->r = ir;
}
//其实他还有很多功能,这里只放我们要用的

这是高斯消元法的本体

//高斯消元法
//返回值为解空间的维数,0为唯一解,-1为无解
//保存一个解,如果有无穷解,保存零空间中的若干个线性无关向量
int GaussEliminate(matrix coeficient, vector<double>constant, vector<double>& a_solution, matrix& null_space) {
 if (constant.size() < (unsigned int)coeficient.r) { constant.resize(coeficient.r); }
 vector<int>pivot, freec;
 double tempmax, ratio, temp2;int temp;
 for (int ic = 0, curr = 0, swapr, icol;ic < coeficient.c;++ic) {
  //找出当前列主元绝对值最大的
  for (swapr = curr + 1, temp = curr, tempmax = abs(coeficient.x[curr][ic]);swapr < coeficient.r;++swapr) {
   temp2 = abs(coeficient.x[swapr][ic]);
   if (temp2 > tempmax) { temp = swapr;tempmax = temp2; }
  }
  if (temp != curr) {
   for (icol = ic;icol < coeficient.c;++icol) {
    std::swap(coeficient.x[curr][icol], coeficient.x[temp][icol]);
   }
   std::swap(constant[curr], constant[temp]);
  }
  //如果主元不为0就继续消
  if (coeficient.x[curr][ic] != 0) {
   pivot.push_back(ic);
   for (temp = curr + 1;temp < coeficient.r;++temp) {
    ratio = coeficient.x[temp][ic] / coeficient.x[curr][ic];
    if (ratio != 0) {
     for (icol = ic;icol < coeficient.c;++icol) {
      coeficient.x[temp][icol] -= coeficient.x[curr][icol] * ratio;
     }
     constant[temp] -= constant[curr] * ratio;
    }
   }
   ++curr;
  }
  else {
   freec.push_back(ic);
  }
 }
 //现在你得到了一个行阶梯
 //查看是否无解
 temp = pivot.size();
 if (temp < coeficient.c) {
  for (int ir = temp;ir < coeficient.r;++ir) {
   if (constant[ir] != 0) { return -1; }
  }
  //如果不无解,就给零空间留好位置
  temp = coeficient.c - temp;
  null_space.reset(coeficient.c, temp);
 }
 else {
  null_space.reset(0, 0);
 }
 //现在开始反向消元
 auto itc = pivot.end() - 1, init = pivot.begin();
 for (;itc - init >= 0;) {
  temp = itc - init;
  ratio = coeficient.x[temp][*itc];
  constant[temp] /= ratio;
  for (int ic = *itc;ic < coeficient.c;++ic) {
   coeficient.x[temp][ic] /= ratio;
  }
  for (int ir = temp - 1;ir >= 0;--ir) {
   ratio = coeficient.x[ir][*itc];
   if (ratio == 0) { continue; }
   for (int ic = *itc;ic < coeficient.c;++ic) {
    coeficient.x[ir][ic] -= coeficient.x[temp][ic] * ratio;
   }
   constant[ir] -= constant[temp] * ratio;
  }
  if (itc == init) { break; }
  else { --itc; }
 }
 //现在是行最简形
 //可以输出结果了
 a_solution = constant; a_solution.resize(coeficient.c, 0);
 if (freec.size()) {
  itc = freec.end() - 1, init = freec.begin();
  for (int cnt = 0;itc - init >= 0;++cnt) {
   for (int ir = 0;ir < coeficient.r;++ir) { null_space.x[ir][cnt] = coeficient.x[ir][*itc]; }
   null_space.x[*itc][cnt] = -1;
   if (itc == init) { break; }
   else { --itc; }
  }
 }
 return coeficient.c - pivot.size();
}

版权声明:本文为m0_52942601原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
原文链接:https://blog.csdn.net/m0_52942601/article/details/112253897