C++ Blatt 5 A13.cc

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cstdlib>
#include <sstream>
#include <vector>
#include <valarray>
using namespace std;

void get_matrix(string filename, vector< valarray<double> > &A, vector<double> &x, int &dim) {
  string line;
  ifstream file(filename);

  for (file >> dim; A.size() < dim;) {
    getline(file, line);
    valarray<double> row(dim);
    istringstream stream(line);

    for (int i = 0; i < dim; i++)
      stream >> row[i];

    if (line.length() > 1)
      A.push_back(row);
  }

  for (double i; file >> i; x.push_back(i));
}

void pivot(vector< valarray<double> > &A, vector<double> &x, int n) {
  int pivot = n;
  for (int i = n + 1; i < A.size(); i++) {
    pivot = (abs(A[i][n]) > abs(A[pivot][n])) ? i : pivot;
  }
  A[n].swap(A[pivot]);
  swap(x[n], x[pivot]);
}

void print(vector< valarray<double> > &A, vector<double> &x) {
  for (int i = 0; i < A.size(); i++) {
    for (int j = 0; j < A.size(); j++) {
      printf("%4.0f", A[i][j]);
    }
    cout << "  |  " << x[i] << endl;
  }
}

int main(int argc, char* argv[]) {
  int dim = 0, usepivot = argc > 1;
  vector< valarray<double> > A;
  vector<double> x;

  cout << (usepivot ? "Using pivot" : "Run with --pivot to enable pivot") << endl << endl << "Input" << endl;

  get_matrix("a13-lgs1.dat", A, x, dim);
  print(A, x);

  cout << endl << "Gauss" << endl;
  for (int n = 0; n < dim; n++) {
    if (usepivot)
      pivot(A, x, n);

    for (int i = n + 1; i < dim; i++) {
      x[i] -= x[n] * (A[i][n]/A[n][n]);
      A[i] -= A[n] * (A[i][n]/A[n][n]);
    }
  }
  print(A, x);

  cout << endl << "Gauss-Jordan" << endl;
  for (int n = dim - 1; n >= 0; n--) {
    x[n] /= A[n][n];
    A[n] /= A[n][n];

    for (int i = 0; i < n; i++) {
      x[i] -= x[n] * A[i][n];
      A[i] -= A[n] * A[i][n];
    }
  }
  print(A, x);
}