使用高斯消去法从文本文件中求逆矩阵.C++

Inverse matrix from text file using Gauss elimination. C++

本文关键字:C++ 文件 文本 高斯消      更新时间:2023-10-16

我需要编写一个程序,从逗号分隔的值文件中读取矩阵,然后使用高斯消去法计算逆,并将此逆写入一个新文件。

把它读进去没问题,把它写出来也没问题。我想我了解高斯消去法是如何工作的,并且能够通过使用数组来做到这一点。

#include <iostream>
#include <fstream>
#include <iomanip>
#include <vector>
#include <string>
using namespace std;
//  main program starts here
int main() {
//  create a vector of a vector to store original matrix
//  and matrices after each calculation
vector<vector<double>> data;
int n_com = 0;
//  try to read input file
ifstream readFile("test_data.txt");
if (readFile.is_open()){
    while (!readFile.eof()){
        int i;
        //  declare temporary line vector and line string
        vector<double> vline;
        string aLine;
        //  assign line of file to line string
        getline(readFile, aLine);
        //  count number of commas in line string
        n_com = count(aLine.begin(), aLine.end(), ',');
        //  define integer for start of each element
        int start = 0;
        //  loop over all but final element
        for (i=0; i < n_com; i++) {
            //  declare and find position of next comma
            int comma_pos;
            comma_pos = aLine.find(',', start);
            //  declare string for element and assign substring to it
            string elems;
            elems = aLine.substr(start, comma_pos - start);
            //  convert string to double
            double elemd = atof(elems.c_str());
            //  push back double to temporary vector
            vline.push_back(elemd);
            //  redefine start for next iteration
            start = comma_pos + 1;
        }
        //  assign final element to string
        string final_elems = aLine.substr(start, aLine.length() - start);
        //  convert final element to double and push back to vector
        double final_elemd = atof(final_elems.c_str());
        vline.push_back(final_elemd);
        //  push back line vector to data vector
        data.push_back(vline);
    }
}
else {
    //  print error if unable to open file
    printf("Error unable to open input file!n");
    //  exit program
    exit(1);
}
//  close input file
readFile.close();

这就是我对原始矩阵的解读。

这是我为阵列进行高斯消除的代码

//  calculate width and length of original data (no. of rows and columns)
int length = data.size();
int width = n_com + 1;
//  create new file to write to
ofstream writeFile ("tranpose.txt");
//  check outfile is open
if (writeFile.is_open()){
    //  declare indices
    //  width (columns)
    int i;
    // length (rows)0
    int j;
    //  k
    int k;
    //  declare a float?
    float data[10][10] = {0},d;
    //  identity matrix
    for (i=1; i <= length; i++){
        for (j=1; j <= 2 * length; j++){
            if (j == (i + length)){
                data[i][j] = 1;
            }
        }
    }
    //  partial pivoting
    for (i=length; i > 1; i--){
        if (data[i-1][1] < data[i][1]){
            for(j=1;j <= length * 2; j++){
                d = data[i][j];
                data[i][j] = data[i-1][j];
                data[i-1][j] = d;
            }
        }
    }
    cout<<"Augmented Matrix: "<<endl;
    for (i=1; i <= length; i++){
        for (j=1;j <= length * 2; j++){
            cout<<data[i][j]<<"  ";
        }
        cout<<endl;
    }
    //  reducing to diagonal matrix
    for (i=1; i <= length; i++){
        for (j=1; j <= length * 2; j++){
            if (j != i){
                d = data[j][i] / data [i][i];
                for (k=1; k<= length * 2; k++){
                    data[j][k] = data[j][k] - (data[i][k] * d);
                }
            }
        }
    }
    //  reducing to unit matrix
    for (i=1; i <= length; i++){
        d = data[i][i];
        for (j=1; j <= length * 2; j++){
            data[i][j] = data[i][j] / d;
        }
    }
    //  print inverse matrix in console
    cout<<"Inverse Matrix "<<endl;
    for (i=1; i <= length; i++){
        for (j = length + 1; j <= length * 2; j++){
            cout<<data[i][j]<<"  ";
        }
        cout<<endl;
    }
    //  loop over all rows
    for (i=1; i <= length; i++){
        //  loop over all columns
        for (j = length + 1; j <= length * 2; j++){
            //  print data in transposed positions excluding last value
            //  i.e. [j][i] instead of [i][j]
            writeFile << setw(4) << fixed << setprecision(2) << data[i][j] << ",";
        }
        // print onto new line
        int i = length - 1;
        writeFile << setw(4) << fixed << setprecision(2) << data[i][j] << "n";
    }
    //  close written file
    writeFile.close();

我如何编写它,使其使用从我的矢量或矩阵文件中存储的数据,而不是普通的数字数组?

我改进了整个过程,因为你的算法在我看来不合适。

using Matrix = std::vector<std::vector<double>>;
Matrix inverse(Matrix mat)
{
    // Use Gaussian elimination
    // Using two matrixs instead of one agumented
    // to improve peformance
    auto height = mat.size();
    auto width = mat[0].size();
    // Create an identity matrix
    Matrix result(height, Matrix::value_type(width));
    for (auto i = 0;i < width;++i) {
        result[i][i] = 1;
    }
    cout << "Augmented Matrix: " << endl;
    printTwo(mat, result);
    // reduce to Echelon form
    for (auto j = 0;j < width;++j) {
        // partial pivoting
        auto maxRow = j;
        for (auto i = j;i < height;++i) {
            maxRow = mat[i][j]>mat[maxRow][j] ? i : maxRow;
        }
        mat[j].swap(mat[maxRow]);
        result[j].swap(result[maxRow]);
        cout << "pivotted Matrix: " << endl;
        printTwo(mat, result);
        // Reduce row by row
        auto pivot = mat[j][j];
        auto& row1L = mat[j];
        auto& row1R = result[j];
        for (auto i = j + 1;i < height;++i) {
            auto& row2L = mat[i];
            auto& row2R = result[i];
            auto temp = row2L[j];
            for (auto k = 0;k < width;++k) {
                row2L[k] -= temp / pivot*row1L[k];
                row2R[k] -= temp / pivot*row1R[k];
            }
        }
        // Make diaganal elements to 1
        for (auto k = 0;k < width;++k) {
            row1L[k] /= pivot;
            row1R[k] /= pivot;
        }
        cout << "reduced Matrix: " << endl;
        printTwo(mat, result);
    }
    //back subsitution
    for (auto j = width - 1;;--j) {
        auto& row1L = mat[j];
        auto& row1R = result[j];
        for (auto i = 0;i < j;++i) {
            auto& row2L = mat[i];
            auto& row2R = result[i];
            auto temp = row2L[j];
            for (auto k = 0;k < width;++k) {
                row2L[k] -= temp*row1L[k];
                row2R[k] -= temp*row1R[k];
            }
        }
        cout << "subsituted Matrix: " << endl;
        printTwo(mat, result);
        if (j == 0) break;
    }
    return result;
}

这是代码中使用的助手:

void printTwo(const Matrix& lhs, const Matrix& rhs)
{
    for (auto i = 0;i < lhs.size();++i) {
        for (auto elm : lhs[i]) {
            cout << setw(4) << elm << ' ';
        }
        for (auto elm : rhs[i]) {
            cout << setw(4) << elm << ' ';
        }
        cout << endl;
    }
}