c++ - NxM 矩阵的高斯消除

Question

/* Program to demonstrate gaussian <strong class="highlight">elimination</strong>
   on a set of linear simultaneous equations
 */

#include <iostream>
#include <cmath>
#include <vector>

using namespace std;

const double eps = 1.e-15;

/*Preliminary pivoting strategy
  Pivoting function
 */
double pivot(vector<vector<double> > &a, vector<double> &b, int i)
 {
     int n = a.size();
     int j=i;
     double t=0;

     for(int k=i; k<n; k+=1)
     {
         double aki = fabs(a[k][i]);
         if(aki>t)
         {
             t=aki;
             j=k;
         }
     }
     if(j>i)
     {
         double dummy;
         for(int L=0; L<n; L+=1)
         {
             dummy = a[i][L];
             a[i][L]= a[j][L];
             a[j][L]= dummy;
         }
         double temp = b[j];
         b[i]=b[j];
         b[j]=temp;
     }
     return a[i][i];
 }        



/* Forward <strong class="highlight">elimination</strong> */
void triang(vector<vector<double> > &a, vector<double> &b) 
{
    int n = a.size();
    for(int i=0; i<n-1; i+=1)
    {
        double diag = pivot(a,b,i);
        if(fabs(diag)<eps)
        {
            cout<<"zero det"<<endl;
            return;
        }
        for(int j=i+1; j<n; j+=1)
        {
            double mult = a[j][i]/diag;
            for(int k = i+1; k<n; k+=1)
            {
                a[j][k]-=mult*a[i][k];
            }
            b[j]-=mult*b[i];
        }
    }
}

/*
DOT PRODUCT OF TWO VECTORS
*/
double dotProd(vector<double> &u, vector<double> &v, int k1,int k2)
{
  double sum = 0;
  for(int i = k1; i <= k2; i += 1)
  {
     sum += u[i] * v[i];
  }
  return sum;
}

/*
BACK SUBSTITUTION STEP
*/
void backSubst(vector<vector<double> > &a, vector<double> &b, vector<double> &x)
{
    int n = a.size();
  for(int i =  n-1; i >= 0; i -= 1)
  {
    x[i] = (b[i] - dotProd(a[i], x, i + 1,  n-1))/ a[i][i];
  }
}
/*
REFINED GAUSSIAN <strong class="highlight">ELIMINATION</strong> PROCEDURE
*/
void gauss(vector<vector<double> > &a, vector<double> &b, vector<double> &x)
{
   triang(a, b);
   backSubst(a, b, x);
}                               

// EXAMPLE MAIN PROGRAM
int main()
{
    int n;
    cin >> n;
    vector<vector<double> > a;
    vector<double> x;
    vector<double> b;
    for (int i = 0; i < n; i++) {
        vector<double> temp;
        for (int j = 0; j < n; j++) {
            int no;
            cin >> no;
            temp.push_back(no);
        }
        a.push_back(temp);
        b.push_back(0);
        x.push_back(0);
    }
    /* 
    for (int i = 0; i < n; i++) {
        int no;
        cin >> no;
        b.push_back(no);
        x.push_back(0);
    }
    */

  gauss(a, b, x);
  for (size_t i = 0; i < x.size(); i++) {
      cout << x[i] << endl;
  }
   return 0;
}

上述高斯消除算法在 NxN 矩阵上运行良好。但我需要它在 NxM 矩阵上工作。任何人都可以帮我做吗？我不太擅长数学。我在某个网站上获得了这段代码，但我一直坚持下去。

Answer 1

1. (可选)了解this .在纸上做一些例子。
2. 不要自己编写高斯消元代码。如果不小心，朴素的高斯旋转是不稳定的。您必须缩放线条并注意使用最大元素进行旋转，起点是 there .请注意，此建议适用于大多数线性代数算法。
3. 如果您想求解方程组，LU decomposition , QR decomposition (比 LU 稳定，但速度较慢)，Cholesky decomposition (在系统对称的情况下)或 SVD (在系统不是正方形的情况下)几乎总是更好的选择。然而，高斯消元法最适合计算行列式。
4. 使用 LAPACK 中的算法对于需要高斯消除的问题(例如求解系统或计算行列式)。真的。不要自己动手。由于您正在使用 C++，因此您可能对 Armadillo 感兴趣它会为您处理很多事情。
5. 如果您出于教学原因必须自己动手，请先查看Numerical Recipes，第 3 版。如果您的预算不足/无法访问图书馆，可以免费在线找到第 2 版。
6. 作为一般性建议，不要编写您不理解的算法。