/* Copyright (c) 1997-2010
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Darmstadt, Germany)
   http://www.polymake.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
--------------------------------------------------------------------------------
   $Project: polymake $$Id: totally_unimodular.h 9717 2010-07-02 13:13:53Z joswig $
*/

#ifndef _POLYMAKE_TOTALLY_UNIMODULAR_H
#define _POLYMAKE_TOTALLY_UNIMODULAR_H

#include <Matrix.h>
#include <PowerSet.h>
#include <linalg.h>

namespace pm {

// A matrix is totally unimodular if the determinant of each square submatrix equals 0, 1, or -1.

// This is the naive test (exponential in the size of the matrix).
// For polynomial time algorithms see
// Schrijver: Theory of Linear and Integer Programming, section 20.3.

bool totally_unimodular(const Matrix<int> M) {
   const int m=M.rows();
   const int n=M.cols();
   const int r=std::min(m,n);

   for (int k=1; k<=r; ++k)
      for (Entire< Subsets_of_k< sequence > >::const_iterator ri = entire(all_subsets_of_k(sequence(0,m),k)); !ri.at_end(); ++ri)
	 for (Entire< Subsets_of_k< sequence > >::const_iterator ci = entire(all_subsets_of_k(sequence(0,n),k)); !ci.at_end(); ++ci) {
	    const int d=det(M.minor(*ri,*ci));
	    if (d!=0 && d!=1 && d!=-1) return false;
	 }
   
   return true;
}

} // end namespace pm

namespace polymake {
   using pm::totally_unimodular;
}

#endif // _POLYMAKE_TOTALLY_UNIMODULAR_H

// Local Variables:
// mode:C++
// c-basic-offset:3
// End:
