0

So I'm trying to read a matrix A from a text file, which it does correctly. A vector B is entered by the user. Then I want to perform Gaussian Elimination (Ax = b) to get the solution vector x. The values I get for x are -1.#IND and I have no idea why...I'm guessing something is going wrong in SystemSolution?

 #include <iostream>
#include <vector>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>

using namespace std;

//this program does gaussian elimination for a matrix Ax=b

vector<double> SystemSolution(vector<vector<double>> A, vector<double> b)
{

    //Determine and test size of a matrix
    int n = A.size();
    for (int i = 0; i < n; i++)
        if (n != A[i].size())
            throw "Error! Number of rows and columns of matrix must be equal!";

    vector<double> x(b.size());
    //x is the vector of solutions


    for (int i = 0; i < n - 1; i++)
    {
        for (int j = i + 1; j < n; j++)
        {
            //Finding pivot
            double pivot = A[i][i];
            int index = i;
            for (int k = i + 1; k < n; k++)
            {
                if (pivot > abs(A[k][i]))
                {
                    index = k;
                    pivot = A[k][i];
                }
            }

            //Row exchange
            for (int k = 0; k < n; k++)
            {
                double tmp = A[i][k];
                A[i][k] = A[index][k];
                A[index][k] = tmp;
            }

            //Elimination
            double coefficient = -(A[j][i] / A[i][i]);
            for (int k = i; k < n; k++)
            {
                A[j][k] += coefficient*A[i][k];
            }

            b[j] += coefficient*b[i];
        }
    }

    //Back-substitution
    x[n - 1] = b[n - 1] / A[n - 1][n - 1];
    for (int i = n - 2; i >= 0; i--)
    {
        double sum = 0;
        for (int j = i; j < n; j++)
        {
            sum += x[j] * A[i][j];
        }
        x[i] = (b[i] - sum) / A[i][i];
    }

    return x;
}



void PrintVector(const vector<double> &b)
{
    for (int i = 0; i < b.size(); i++)
        cout << setiosflags(ios::showpoint | ios::fixed | ios::right)
        << setprecision(4)
        << setw(8) << b[i] << endl;
}

void PrintMatrix(const vector<vector<double> > &A)
{
    for (int i = 0; i < A.size(); i++)
    {
        for (int j = 0; j < A[i].size(); j++)
            cout << setiosflags(ios::showpoint | ios::fixed | ios::right)
            << setprecision(4)
            << setw(8) << A[i][j];
        cout << endl;
    }
}
int main()
{
    int n;
    cout << "Please enter the number of rows/columns:";

    cin >> n;

    ifstream matrixFile;

    matrixFile.open("matrix.txt");

    if (matrixFile.is_open()){

        //matrix A values
        vector<vector<double>> A(n, vector<double>(n));
        vector<double> b(n);
        string line;
        int col = 0;
        int row = 0;

        while (getline(matrixFile, line)){

            istringstream stream(line);

            int x;
            col = 0; //reset 

            while (stream >> x){
                A[row][col] = x;
                col++;
            }

            row++;

        }





        cout << "Please enter vector b:"<<endl;
        //vector b values
        for (int i = 0; i < row; i++)
        {
            cout << "b[" << i + 1 << "]= ";
            cin >> b[i];    
        }



        vector<double> x = SystemSolution(A, b);
        cout << "- SOLUTIONS -" << endl;
        cout << "Matrix:" << endl;
        PrintMatrix(A);
        cout << "\nVector x:" << endl;
        PrintVector(x);

    }
    else{
        cout << "File failed to open!";
    }
    matrixFile.close();

    return 0;
}
Xar
  • 41
  • 1
  • 1
  • 5
  • 2
    Hi. You say your I/O is working and don't document your input or exact output. I suggest creating an illustration of your unresolved problem that hardcodes a couple vectors minimally illustrating the problem, then listing the actual output of that program and what you'd expected, and how far you got diagnosing the difference. As is, "values I get for x are -1.#IND" is unintelligible. – Tony Delroy Nov 26 '14 at 08:41
  • 1
    Why have you tagged with `file-io` if that is not your issue? – Vinayak Garg Nov 26 '14 at 08:50
  • 1
    @Xar add debug print of your elimination (each step with comment what it does) and then visually check where the problem is. When you find specific problem update question... – Spektre Nov 26 '14 at 09:26
  • Plase show some sample input. – Codor Nov 26 '14 at 11:48
  • Here is a screenshot: http://i.imgur.com/p5be57R.jpg I enter the values for b. Matrix A is read from a txt file. And the values for x are -1.#IND – Xar Nov 26 '14 at 14:11
  • `-1.#IND` is a particular NaN encoding (INDEFINITE). This can be generated by a division 0/0, for example. Check the intermediate data in your code. I don't see a check for a pivot of 0, so a division of 0 by 0 seems a plausible working hypothesis. – njuffa Nov 26 '14 at 17:11

1 Answers1

0

There could be some divisions by zero in your code:

double coefficient = -(A[j][i] / A[i][i]);
/* .. */
x[n - 1] = b[n - 1] / A[n - 1][n - 1];
/* .. */
x[i] = (b[i] - sum) / A[i][i];

Check out Gauss-elimination here: Square Matrix Inversion in C

Review and debug yours.

Community
  • 1
  • 1
renonsz
  • 571
  • 1
  • 4
  • 17