Template Gauss Elimination

First it will retrieve a possible solution for the augmented matrix. Then if valid() returns false, there is no solution. If it is required to distinguish whether there are infinite solutions or not, do an extra work checking whether there is a row where every element is 0.

const double EPS = 1e-3;
const double DINF = 1e9;

struct Gauss_Solve_Matrix{
    vvd mat;
    vd ans;
    Gauss_Solve_Matrix( const vvd &m ){
        mat = m;
        ans.resize( m[ 0 ].size() - 1 );
    }
    void Gauss_elimination(){
        for( int i = 0; i < min( mat.size(), mat[ 0 ].size() - 1 ); ++i ){
            if( mat[ i ][ i ] == 0 )
                for( int j = i + 1; j < mat.size(); ++j )
                    if( mat[ j ][ i ] ){
                        swap( mat[ i ], mat[ j ] );
                        break;
                    }
            if( mat[ i ][ i ] == 0 ) continue;

            for( int j = mat[ i ].size() - 1; j >= i; --j )
                mat[ i ][ j ] /= mat[ i ][ i ];

            for( int j = i + 1; j < mat.size(); ++j )
                for( int k = mat[ j ].size() - 1; k >= i; --k )
                    mat[ j ][ k ] -= mat[ j ][ i ] * mat[ i ][ k ];
        }
    }
    void back_substitution(){
        for( int i = mat.size() - 1; i >= 0; --i ){
            int t; // first non zero coefficient
            for( t = 0; t + 1 < mat[ i ].size(); ++t )
                if( mat[ i ][ t ] )
                    break;
            if( t >= ( int ) mat[ i ].size() - 1 ) continue;
            double s = 0;
            for( int j = t + 1; j + 1 < mat[ i ].size(); ++j )
                s += mat[ i ][ j ] * ans[ j ];
            ans[ t ] = ( mat[ i ].back() - s ) / mat[ i ][ t ];
        }
    }
};

bool valid( const vvd &mat, const vd &sol ){
    for( int i = 0; i < mat.size(); ++i ){
        double s = 0;
        for( int j = 0; j + 1 < mat[ i ].size(); ++j )
            s += mat[ i ][ j ] * sol[ j ];
        if( s > EPS + mat[ i ].back() ) return false;
    }
    return true;
}