一个完善的高斯消元法(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 版权协议,转载请附上原文出处链接和本声明。