leq_solver.h

/****************************************************************************
  FileName     [ leq_solver.h ]
  Synopsis     [ Linear Equation Solver ]
  Author       [ Jui-Yang (Tony) Hsu]
  Copyright    [ Copyleft(c) 2017 NTUEE, Taiwan ]
****************************************************************************/
#ifndef LEQ_SOLVER_H
#define LEQ_SOLVER_H
#include <vector>
#include <algorithm>
#define print(x) std::cout << x << std::endl
template<typename T>
std::ostream& operator <<(std::ostream &s, const std::valarray<T> &c) {
//s<<"[ ";
for (auto it : c) s << it << " ";
//s<<"]";
return s;
}
namespace leq{
  //typedef double double;
  typedef std::valarray<std::valarray<double> > Matrix;
  const double eps = 1e-10;
  bool is_zero(std::valarray<double>& v){
    double s=std::abs(v).sum();
    if(s<eps){
      return true;
    }
    else
      return false;
  }

  //Gaussian Elimination
  void gauss_eliminate(Matrix& mat){
    int n = mat.size();
    for (int i = 0; i<n; i++) {

      //swap the upper row if 1st coeff. is 0
      if(fabs(mat[i][i]) < eps){
        for(int j=i+1;j<n;++j){
          if(fabs(mat[j][i]) > eps){
            std::swap(mat[i],mat[j]);
            break;
          }
        }
      }

      if(fabs(mat[i][i]) < eps) continue; // still 0

      mat[i] /= mat[i][i];
      
      //eliminate the following rows
      for(int j=i+1;j<n;++j){
        if(fabs(mat[j][i]) > eps){
          std::valarray<double> tmp  = mat[j] - mat[i] * mat[j][i];
          mat[j] = tmp;
        }
      }
    }
    //swap zero rows down
    for(int i=0,j=n-1;i<=j;){
      if(is_zero(mat[i])){
        std::swap(mat[i],mat[j]);
        --j;
      }
      else
        ++i;
    }
  }

  std::valarray<double> solve(Matrix& A,std::valarray<double>& b,int& status){
    int n = A.size();
    Matrix augment_A(std::valarray<double>(n+1),n);
    for(std::size_t i=0;i<n;++i){
      for(std::size_t j=0;j<n;++j){
        augment_A[i][j] = A[i][j];
      }
      augment_A[i][n] = b[i];
    }
    gauss_eliminate(augment_A);

    std::valarray<double> x(n);

    int next_var = n-1; //[solved_var ,n-1]-th var solved
    for(int i=n-1;i>=0 && status!=-1;--i){
      //print(i);
      if(!is_zero(augment_A[i])){
        double t=0;
        for(std::size_t k=next_var+1;k<=n-1;++k)
          t += augment_A[i][k] * x[k];
        while(fabs(augment_A[i][next_var]) < eps && next_var >i){
          --next_var;
        }
        if(next_var == i){
          if(fabs(augment_A[i][next_var]) < eps){
            status = (fabs(augment_A[i][n] - t) < eps)?1:-1;
          }
          else{
            x[next_var] = (augment_A[i][n] - t) / augment_A[i][next_var];
            next_var --;
            //print();
          }
        }
        else{
          x[next_var] = (augment_A[i][n] - t) / augment_A[i][next_var];
          next_var --;
        }
      }
      else
       status = 1; 
    }
    return x;
  }
}
#endif

main.cpp

#pragma GCC optimize ("O2")
#include<bits/stdc++.h>
#include<unistd.h>
#include "leq_solver.h"
using namespace std;
#define ALL(x) begin(x),end(x)
#define IOS ios_base::sync_with_stdio(0); cin.tie(0)
template<typename A, typename B>
ostream& operator <<(ostream &s, const pair<A,B> &p) {
return s<<"("<<p.first<<","<<p.second<<")";
}
//template<typename T>
//ostream& operator <<(ostream &s, const valarray<T> &c) {
//s<<"[ ";
//for (auto it : c) s << it << " ";
//s<<"]";
//return s;
//}

int main(){
  int n;
  while(cin >> n){
    int status = 0;
    if(n == 0) break;
    leq::Matrix A(valarray<double>(n),n);
    for(size_t i=0;i<n;++i){
      for(size_t j=0;j<n;++j)
        cin >> A[i][j];
    }
    valarray<double> b(n);
    for(size_t i=0;i<n;++i)
      cin >> b[i];
    valarray<double>ans = leq::solve(A,b,status);
    if(status == 1) cout << "multiple" << endl;
    else if(status == -1) cout << "inconsistent" << endl;
    else cout << ans << endl;
  }
}