Merge branch '41-remove-boost-dependency'

CMakeLists.txt

@@ -55,7 +55,7 @@ add_library(cmr
   src/cmr/hashtable.c
   src/cmr/heap.c
   src/cmr/equimodular.c
-  src/cmr/linalg.c
+  src/cmr/linear_algebra.c
   src/cmr/listmatrix.c
   #src/cmr/logger.cpp
   #src/cmr/matrix.cpp

doc/tu.md

@@ -26,9 +26,9 @@ determines whether the matrix given in file `IN-MAT` is totally unimodular.
 If `IN-MAT` is `-` then the matrix is read from stdin.
 If `OUT-DEC` or `NON-SUB` is `-` then the decomposition tree (resp. the submatrix) is written to stdout.
 
-## Algorithm ##
+## Algorithms ##
 
-The implemented recognition algorithm is based on [Implementation of a unimodularity test](https://doi.org/10.1007/s12532-012-0048-x) by Matthias Walter and Klaus Truemper (Mathematical Programming Computation, 2013).
+The implemented default recognition algorithm is based on [Implementation of a unimodularity test](https://doi.org/10.1007/s12532-012-0048-x) by Matthias Walter and Klaus Truemper (Mathematical Programming Computation, 2013).
 It first runs \ref camion to reduce the question to that of [recognizing regular matroids](\ref regular).
 Please cite the paper in case the implementation contributed to your research:
 
@@ -45,6 +45,11 @@ Please cite the paper in case the implementation contributed to your research:
       publisher = {Springer-Verlag},
     }
 
+In order to repeat experiments described in the paper above, the function can be parameterized as to use exponential-time algorithms.
+
+  - The first is based on the criterion of Ghouila-Houri and runs in time \f$ \mathcal{O}( (m + n) \cdot 3^{\min(m, n)}) \f$.
+  - The second enumerates square [Eulerian submatrices](https://www.ams.org/journals/proc/1965-016-05/S0002-9939-1965-0180568-2/) and runs in time \f$ \mathcal{O}( (m+n) \cdot 2^{ m + n } ) \f$.
+
 ## C Interface ##
 
 The corresponding function in the library is

include/cmr/env.h

@@ -23,7 +23,6 @@ extern "C" {
 /* Macro for intended non-use of variables. */
 #define CMR_UNUSED(x) (void)(x)
 
-
 /**
  * \brief Type for return codes of library functions.
  **/

include/cmr/linear_algebra.h

@@ -0,0 +1,48 @@
+#ifndef CMR_LINEAR_ALGEBRA_H
+#define CMR_LINEAR_ALGEBRA_H
+
+/**
+ * \file linear_algebra.h
+ *
+ * \author Matthias Walter
+ */
+
+#include <cmr/env.h>
+#include <cmr/matrix.h>
+
+#include <stdio.h>
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/**
+ * \brief Computes the determinant of an 8-bit integer matrix.
+ */
+
+CMR_EXPORT
+CMR_ERROR CMRchrmatDeterminant(
+  CMR* cmr,               /**< \ref CMR environment. */
+  CMR_CHRMAT* matrix,     /**< Matrix. */
+  int64_t* pdeterminant   /**< Pointer for storing the determinant. */
+);
+
+/**
+ * \brief Computes the determinant of an int matrix.
+ */
+
+CMR_EXPORT
+CMR_ERROR CMRintmatDeterminant(
+  CMR* cmr,               /**< \ref CMR environment. */
+  CMR_INTMAT* matrix,     /**< Matrix. */
+  int64_t* pdeterminant   /**< Pointer for storing the determinant. */
+);
+
+
+/**@}*/
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif /* CMR_LINEAR_ALGEBRA_H */

include/cmr/matrix.h

@@ -1005,10 +1005,21 @@ CMR_ERROR CMRchrmatSignedSupport(
   CMR_CHRMAT** presult  /**< Pointer for storing the signed support matrix of \p matrix. */
 );
 
+/**
+ * \brief Converts a char matrix to an int matrix.
+ */
+
+CMR_EXPORT
+CMR_ERROR CMRchrmatToInt(
+  CMR* cmr,             /**< \ref CMR environment. */
+  CMR_CHRMAT* matrix,   /**< Input matrix. */
+  CMR_INTMAT** presult  /**< Pointer for storing the output matrix. */
+);
+
 /**
  * \brief Converts an int matrix to a char matrix.
  *
- * \returns \ref CMR_ERROR_INPUT in case of overflow.
+ * \returns \ref CMR_ERROR_OVERFLOW in case of overflow.
  */
 
 CMR_EXPORT

include/cmr/tu.h

@@ -18,8 +18,16 @@ extern "C" {
 #include <cmr/matrix.h>
 #include <cmr/camion.h>
 
+typedef enum
+{
+  CMR_TU_ALGORITHM_DECOMPOSITION = 0, /**< \brief Algorithm based on Seymour's decomposition of regular matroids. */
+  CMR_TU_ALGORITHM_SUBMATRIX = 1,     /**< \brief Enumeration algorithm based on submatrices. */
+  CMR_TU_ALGORITHM_PARTITION = 2      /**< \brief Enumeration algorithm based on criterion of Ghouila-Houri. */
+} CMR_TU_ALGORITHM;
+
 typedef struct
 {
+  CMR_TU_ALGORITHM algorithm;     /**< \brief Algorithm to use. */
   CMR_REGULAR_PARAMETERS regular; /**< \brief Parameters for regularity test. */
 } CMR_TU_PARAMETERS;
 
@@ -90,6 +98,7 @@ CMR_ERROR CMRtestTotalUnimodularity(
   double timeLimit            /**< Time limit to impose. */
 );
 
+
 #ifdef __cplusplus
 }
 #endif

src/cmr/equimodular.c

@@ -9,7 +9,7 @@
 #include <inttypes.h>
 
 #include "env_internal.h"
-#include "linalg.h"
+#include "linear_algebra_internal.h"
 
 CMR_ERROR CMRparamsEquimodularityInit(CMR_EQUIMODULAR_PARAMETERS* params)
 {

src/cmr/linalg.c → src/cmr/linear_algebra.c

@@ -1,6 +1,6 @@
 // #define CMR_DEBUG /* Uncomment to debug this file. */
 
-#include "linalg.h"
+#include "linear_algebra_internal.h"
 #include "listmatrix.h"
 #include "sort.h"
 
@@ -571,7 +571,6 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(CMR* cmr, CMR_INTMAT* matrix, bool inver
 {
   assert(cmr);
   assert(matrix);
-  assert(ppermutations);
   assert(prank);
 
   bool isIntTooSmall = false;
@@ -581,8 +580,10 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(CMR* cmr, CMR_INTMAT* matrix, bool inver
   CMRintmatPrintDense(cmr, matrix, stdout, '0', true);
 #endif /* CMR_DEBUG */
 
-  CMR_CALL( CMRsubmatCreate(cmr, matrix->numRows, matrix->numColumns, ppermutations) );
-  CMR_SUBMAT* permutations = *ppermutations;
+  CMR_SUBMAT* permutations = NULL;
+  CMR_CALL( CMRsubmatCreate(cmr, matrix->numRows, matrix->numColumns, &permutations) );
+  if (ppermutations)
+    *ppermutations = permutations;
   assert(permutations);
 
   for (size_t row = 0; row < matrix->numRows; ++row)
@@ -631,14 +632,16 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(CMR* cmr, CMR_INTMAT* matrix, bool inver
     for (size_t permutedRow = *prank; permutedRow < matrix->numRows; ++permutedRow)
     {
       size_t row = permutations->rows[permutedRow];
-      CMRdbgMsg(4, "Permuted row %ld is row %ld\n", permutedRow, row);
+      CMRdbgMsg(4, "Permuted row %ld is original row %ld\n", permutedRow, row);
       for (ListMat64Nonzero* nz = listmatrix->rowElements[row].head.right; nz != &listmatrix->rowElements[row].head;
         nz = nz->right)
       {
         int64_t x = llabs(nz->value);
         if (x > llabs(minEntryValue))
           continue;
 
+        CMRdbgMsg(6, "Area for potential pivot row %ld with column %ld (abs value %ld) is %ldx%ld.\n", row, nz->column, x,
+          (listmatrix->rowElements[row].numNonzeros - 1), (listmatrix->columnElements[nz->column].numNonzeros - 1));
         size_t area = (listmatrix->rowElements[row].numNonzeros - 1)
           * (listmatrix->columnElements[nz->column].numNonzeros - 1);
         if (x < llabs(minEntryValue) || (x == llabs(minEntryValue) && area < minEntryArea))
@@ -940,11 +943,84 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(CMR* cmr, CMR_INTMAT* matrix, bool inver
   if (isIntTooSmall)
   {
     CMR_CALL( CMRintmatComputeUpperDiagonalGMP(cmr, matrix, invert, prank, *ppermutations, presult, ptranspose) );
-
-    return CMR_OKAY;
+    isIntTooSmall = false;
   }
 
 #endif /* CMR_WITH_GMP */
 
+  if (!ppermutations)
+    CMR_CALL( CMRsubmatFree(cmr, &permutations) );
+
+#if defined(CMR_DEBUG)
+
+  if (presult)
+  {
+    CMRdbgMsg(0, "CMRintmatComputeUpperDiagonal computed the following transformed matrix:\n");
+    CMRintmatPrintDense(cmr, *presult, stdout, '0', true);
+  }
+  else if (ptranspose)
+  {
+    CMRdbgMsg(0, "CMRintmatComputeUpperDiagonal computed the following transposed transformed matrix:\n");
+    CMRintmatPrintDense(cmr, *ptranspose, stdout, '0', true);
+  }
+
+#endif /* CMR_DEBUG */
+
   return isIntTooSmall ? CMR_ERROR_OVERFLOW : CMR_OKAY;
 }
+
+CMR_ERROR CMRintmatDeterminant(CMR* cmr, CMR_INTMAT* matrix, int64_t* pdeterminant)
+{
+  assert(cmr);
+  assert(matrix);
+  assert(pdeterminant);
+  CMRconsistencyAssert(CMRintmatConsistency(matrix));
+
+  if (matrix->numRows != matrix->numColumns)
+    return CMR_ERROR_INPUT;
+
+  CMR_INTMAT* transformed = NULL;
+  size_t rank = SIZE_MAX;
+  CMR_CALL( CMRintmatComputeUpperDiagonal(cmr, matrix, false, &rank, NULL, &transformed, NULL) );
+  if (rank < matrix->numRows)
+    *pdeterminant = 0;
+  else
+  {
+    int64_t old;
+    int64_t det = 1;
+    for (size_t row = 0; row < transformed->numRows; ++row)
+    {
+      int64_t x = transformed->entryValues[transformed->rowSlice[row]];
+      old = det;
+      det *= x;
+      if (det / x != old)
+      {
+        CMR_CALL( CMRintmatFree(cmr, &transformed) );
+        return CMR_ERROR_OVERFLOW;
+      }
+    }
+    *pdeterminant = det;
+  }
+
+  CMR_CALL( CMRintmatFree(cmr, &transformed) );
+
+  return CMR_OKAY;
+}
+
+
+CMR_ERROR CMRchrmatDeterminant(CMR* cmr, CMR_CHRMAT* matrix, int64_t* pdeterminant)
+{
+  assert(cmr);
+  assert(matrix);
+  CMRconsistencyAssert(CMRchrmatConsistency(matrix));
+
+  if (matrix->numRows != matrix->numColumns)
+    return CMR_ERROR_INPUT;
+
+  CMR_INTMAT* copy = NULL;
+  CMR_CALL( CMRchrmatToInt(cmr, matrix, &copy) );
+  CMR_CALL( CMRintmatDeterminant(cmr, copy, pdeterminant) );
+  CMR_CALL( CMRintmatFree(cmr, &copy) );
+
+  return CMR_OKAY;
+}

src/cmr/linalg.h → src/cmr/linear_algebra_internal.h

@@ -1,8 +1,7 @@
 #ifndef CMR_LINALG_INTERNAL_H
 #define CMR_LINALG_INTERNAL_H
 
-#include <cmr/env.h>
-#include <cmr/matrix.h>
+#include <cmr/linear_algebra.h>
 
 #ifdef __cplusplus
 extern "C" {
@@ -23,7 +22,8 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(
   CMR_INTMAT* matrix,         /**< A matrix */
   bool invert,                /**< Whether the transformed basis columns shall be strictly diagonally dominant. */
   size_t* prank,              /**< Pointer for storing the rank of the basis matrix. */
-  CMR_SUBMAT** ppermutations, /**< Pointer for storing the row- and column permutations applied to \p matrix. */
+  CMR_SUBMAT** ppermutations, /**< Pointer for storing the row- and column permutations applied to \p matrix
+                               **  (may be \c NULL). */
   CMR_INTMAT** presult,       /**< Pointer for storing the resulting int matrix (may be \c NULL). */
   CMR_INTMAT** ptranspose     /**< Pointer for storing the transpose of the result int matrix (may be \c NULL). */
 );

src/cmr/listmatrix.c

@@ -1460,7 +1460,7 @@ CMR_ERROR CMRlistmat64PrintDense(CMR* cmr, ListMat64* listmatrix, FILE* stream)
       fprintf(stream, " %" PRId64 "", dense[column]);
       dense[column] = 0;
     }
-    fprintf(stream, "\n");
+    fprintf(stream, " (%ld nonzeros)\n", listmatrix->rowElements[row].numNonzeros);
   }
   fflush(stream);
 
@@ -1834,6 +1834,10 @@ CMR_ERROR CMRlistmat8Delete(CMR* cmr, ListMat8* listmatrix, ListMat8Nonzero* nz)
   assert(nz != &listmatrix->rowElements[nz->row].head);
   assert(nz != &listmatrix->columnElements[nz->column].head);
 
+  listmatrix->numNonzeros--;
+  listmatrix->rowElements[nz->row].numNonzeros--;
+  listmatrix->columnElements[nz->column].numNonzeros--;
+
   nz->left->right = nz->right;
   nz->right->left = nz->left;
   nz->above->below = nz->below;
@@ -1857,13 +1861,16 @@ CMR_ERROR CMRlistmat64Delete(CMR* cmr, ListMat64* listmatrix, ListMat64Nonzero*
   assert(nz != &listmatrix->rowElements[nz->row].head);
   assert(nz != &listmatrix->columnElements[nz->column].head);
 
+  listmatrix->numNonzeros--;
+  listmatrix->rowElements[nz->row].numNonzeros--;
+  listmatrix->columnElements[nz->column].numNonzeros--;
+
   nz->left->right = nz->right;
   nz->right->left = nz->left;
   nz->above->below = nz->below;
   nz->below->above = nz->above;
   nz->right = listmatrix->firstFreeNonzero;
   listmatrix->firstFreeNonzero = nz;
-  listmatrix->numNonzeros--;
 
   return CMR_OKAY;
 }
@@ -1882,13 +1889,16 @@ CMR_ERROR CMRlistmatGMPDelete(CMR* cmr, ListMatGMP* listmatrix, ListMatGMPNonzer
   assert(nz != &listmatrix->rowElements[nz->row].head);
   assert(nz != &listmatrix->columnElements[nz->column].head);
 
+  listmatrix->numNonzeros--;
+  listmatrix->rowElements[nz->row].numNonzeros--;
+  listmatrix->columnElements[nz->column].numNonzeros--;
+
   nz->left->right = nz->right;
   nz->right->left = nz->left;
   nz->above->below = nz->below;
   nz->below->above = nz->above;
   nz->right = listmatrix->firstFreeNonzero;
   listmatrix->firstFreeNonzero = nz;
-  listmatrix->numNonzeros--;
 
   return CMR_OKAY;
 }

src/cmr/matrix.c

@@ -2645,6 +2645,27 @@ CMR_ERROR CMRchrmatSignedSupport(CMR* cmr, CMR_CHRMAT* matrix, CMR_CHRMAT** pres
   return CMR_OKAY;
 }
 
+CMR_ERROR CMRchrmatToInt(CMR* cmr, CMR_CHRMAT* matrix, CMR_INTMAT** presult)
+{
+  assert(cmr);
+  CMRconsistencyAssert( CMRchrmatConsistency(matrix) );
+  assert(presult);
+
+  CMR_CALL( CMRintmatCreate(cmr, presult, matrix->numRows, matrix->numColumns, matrix->numNonzeros) );
+  CMR_INTMAT* result = *presult;
+
+  for (size_t row = 0; row <= matrix->numRows; ++row)
+    result->rowSlice[row] = matrix->rowSlice[row];
+
+  for (size_t e = 0; e < matrix->numNonzeros; ++e)
+  {
+    result->entryColumns[e] = matrix->entryColumns[e];
+    result->entryValues[e] = matrix->entryValues[e];
+  }
+
+  return CMR_OKAY;
+}
+
 CMR_ERROR CMRintmatToChr(CMR* cmr, CMR_INTMAT* matrix, CMR_CHRMAT** presult)
 {
   assert(cmr);