diff --git a/CMakeLists.txt b/CMakeLists.txt index f89d64e5..cb33f56e 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/doc/tu.md b/doc/tu.md index 6742f6f5..d6151e3e 100644 --- a/doc/tu.md +++ b/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 diff --git a/include/cmr/env.h b/include/cmr/env.h index 64105d4f..95e47213 100644 --- a/include/cmr/env.h +++ b/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. **/ diff --git a/include/cmr/linear_algebra.h b/include/cmr/linear_algebra.h new file mode 100644 index 00000000..4306140d --- /dev/null +++ b/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 +#include + +#include + +#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 */ diff --git a/include/cmr/matrix.h b/include/cmr/matrix.h index 4dedad25..6e6be78a 100644 --- a/include/cmr/matrix.h +++ b/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 diff --git a/include/cmr/tu.h b/include/cmr/tu.h index c54dddb7..e5107d28 100644 --- a/include/cmr/tu.h +++ b/include/cmr/tu.h @@ -18,8 +18,16 @@ extern "C" { #include #include +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 diff --git a/src/cmr/equimodular.c b/src/cmr/equimodular.c index c88a406b..b411a35f 100644 --- a/src/cmr/equimodular.c +++ b/src/cmr/equimodular.c @@ -9,7 +9,7 @@ #include #include "env_internal.h" -#include "linalg.h" +#include "linear_algebra_internal.h" CMR_ERROR CMRparamsEquimodularityInit(CMR_EQUIMODULAR_PARAMETERS* params) { diff --git a/src/cmr/linalg.c b/src/cmr/linear_algebra.c similarity index 93% rename from src/cmr/linalg.c rename to src/cmr/linear_algebra.c index 8b54c0e3..b451ac1a 100644 --- a/src/cmr/linalg.c +++ b/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,7 +632,7 @@ 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) { @@ -639,6 +640,8 @@ CMR_ERROR CMRintmatComputeUpperDiagonal(CMR* cmr, CMR_INTMAT* matrix, bool inver 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, ©) ); + CMR_CALL( CMRintmatDeterminant(cmr, copy, pdeterminant) ); + CMR_CALL( CMRintmatFree(cmr, ©) ); + + return CMR_OKAY; +} diff --git a/src/cmr/linalg.h b/src/cmr/linear_algebra_internal.h similarity index 91% rename from src/cmr/linalg.h rename to src/cmr/linear_algebra_internal.h index 266f9932..c7ddb0ed 100644 --- a/src/cmr/linalg.h +++ b/src/cmr/linear_algebra_internal.h @@ -1,8 +1,7 @@ #ifndef CMR_LINALG_INTERNAL_H #define CMR_LINALG_INTERNAL_H -#include -#include +#include #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). */ ); diff --git a/src/cmr/listmatrix.c b/src/cmr/listmatrix.c index 5a8ca807..518aba7e 100644 --- a/src/cmr/listmatrix.c +++ b/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; } diff --git a/src/cmr/matrix.c b/src/cmr/matrix.c index d9eb257e..95b244c8 100644 --- a/src/cmr/matrix.c +++ b/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); diff --git a/src/cmr/tu.c b/src/cmr/tu.c index fa26f82b..4f0db74c 100644 --- a/src/cmr/tu.c +++ b/src/cmr/tu.c @@ -1,6 +1,6 @@ -#include +#define CMR_DEBUG /* Uncomment to debug this file. */ -#include +#include #include "matrix_internal.h" #include "one_sum.h" @@ -16,6 +16,7 @@ CMR_ERROR CMRparamsTotalUnimodularityInit(CMR_TU_PARAMETERS* params) { assert(params); + params->algorithm = CMR_TU_ALGORITHM_DECOMPOSITION; CMR_CALL( CMRparamsRegularInit(¶ms->regular) ); return CMR_OKAY; @@ -97,6 +98,165 @@ CMR_ERROR tuTest( return CMR_OKAY; } + +/** + * \brief Recursively assigns +1 or -1 to each row that is part of the subset to test Ghouila-Houri. + * + * \return whether a feasible assignment was found. + */ + +static +bool testPartitionSearch( + CMR* cmr, /**< \ref CMR environment */ + CMR_CHRMAT* matrix, /**< Matrix \f$ M \f$. */ + int8_t* selection, /**< Array with selection. */ + size_t current, /**< Index to decide for selection. */ + int* columnSum /**< Array for computing column sums. */ +) +{ + assert(cmr); + assert(matrix); + assert(selection); + + while (current < matrix->numRows && selection[current] == 0) + ++current; + + if (current < matrix->numRows) + { + /* Recurse by keeping current row a +1. */ + bool found = testPartitionSearch(cmr, matrix, selection, current + 1, columnSum); + if (found) + return true; + + /* Recurse by making current row a -1. */ + size_t first = matrix->rowSlice[current]; + size_t beyond = matrix->rowSlice[current + 1]; + + selection[current] = -1; + for (size_t i = first; i < beyond; ++i) + columnSum[matrix->entryColumns[i]] -= 2 * matrix->entryValues[i]; + + found = testPartitionSearch(cmr, matrix, selection, current + 1, columnSum); + + selection[current] = +1; + for (size_t i = first; i < beyond; ++i) + columnSum[matrix->entryColumns[i]] += 2 * matrix->entryValues[i]; + + return found; + } + else + { + for (size_t column = 0; column < matrix->numColumns; ++column) + { + int sum = columnSum[column]; + if (sum < -1 || sum > +1) + return false; + } + return true; + } +} + +/** + * \brief Recursively selects a row subset and tests Ghouila-Houri for each. + * + * \return 1 if totally unimodular, 0 if not, and -1 if time limit was reached. + */ + +static +int testPartitionSubset( + CMR* cmr, /**< \ref CMR environment */ + CMR_CHRMAT* matrix, /**< Matrix \f$ M \f$. */ + int8_t* selection, /**< Array with selection. */ + size_t current, /**< Index to decide for selection. */ + int* columnSum, /**< Array for computing column sums. */ + clock_t startClock, /**< Clock for start for computation. */ + double timeLimit /**< Time limit for computation. */ +) +{ + assert(cmr); + assert(matrix); + assert(selection); + + if (current < matrix->numRows) + { + /* Recurse by not selecting a column. */ + selection[current] = 0; + int result = testPartitionSubset(cmr, matrix, selection, current + 1, columnSum, startClock, timeLimit); + if (result <= 0) + return result; + + /* Recurse by selecting a column unless we need to abort. */ + selection[current] = 1; + size_t first = matrix->rowSlice[current]; + size_t beyond = matrix->rowSlice[current + 1]; + for (size_t i = first; i < beyond; ++i) + columnSum[matrix->entryColumns[i]] += matrix->entryValues[i]; + + result = testPartitionSubset(cmr, matrix, selection, current + 1, columnSum, startClock, timeLimit); + + for (size_t i = first; i < beyond; ++i) + columnSum[matrix->entryColumns[i]] -= matrix->entryValues[i]; + + return result; + } + else + { + if (((clock() * 1.0 / CLOCKS_PER_SEC) - startClock) > timeLimit) + return -1; + + bool foundPartition = testPartitionSearch(cmr, matrix, selection, 0, columnSum); + return foundPartition ? 1 : 0; + } +} + +/** + * \brief Partition test based on Ghouila-Houri. + */ + +static +CMR_ERROR testPartition( + CMR* cmr, /**< \ref CMR environment */ + CMR_CHRMAT* matrix, /**< Matrix \f$ M \f$. */ + bool* pisTotallyUnimodular, /**< Pointer for storing whether \f$ M \f$ is totally unimodular. */ + double timeLimit /**< Time limit to impose. */ +) +{ + assert(cmr); + assert(matrix); + assert(pisTotallyUnimodular); + + CMR_ERROR error = CMR_OKAY; + + /* Consider transpose if this has fewer rows. */ + if (matrix->numRows > matrix->numColumns) + { + CMR_CHRMAT* transpose = NULL; + CMR_CALL( CMRchrmatTranspose(cmr, matrix, &transpose) ); + CMR_CALL( testPartition(cmr, transpose, pisTotallyUnimodular, timeLimit) ); + CMR_CALL( CMRchrmatFree(cmr, &transpose) ); + return CMR_OKAY; + } + + int8_t* selection = NULL; + CMR_CALL( CMRallocStackArray(cmr, &selection, matrix->numRows) ); + int* columnSum = NULL; + CMR_CALL( CMRallocStackArray(cmr, &columnSum, matrix->numColumns) ); + for (size_t column = 0; column < matrix->numColumns; ++column) + columnSum[column] = 0; + + clock_t startClock = clock(); + int result = testPartitionSubset(cmr, matrix, selection, 0, columnSum, startClock, timeLimit); + if (result < 0) + error = CMR_ERROR_TIMEOUT; + else + *pisTotallyUnimodular = (result > 0); + + CMR_CALL( CMRfreeStackArray(cmr, &columnSum) ); + CMR_CALL( CMRfreeStackArray(cmr, &selection) ); + + return error; +} + CMR_ERROR CMRtestTotalUnimodularity(CMR* cmr, CMR_CHRMAT* matrix, bool* pisTotallyUnimodular, CMR_DEC** pdec, CMR_SUBMAT** psubmatrix, CMR_TU_PARAMETERS* params, CMR_TU_STATISTICS* stats, double timeLimit) { @@ -128,17 +288,34 @@ CMR_ERROR CMRtestTotalUnimodularity(CMR* cmr, CMR_CHRMAT* matrix, bool* pisTotal return CMR_OKAY; } - // TODO: run regularity check with ternary = true. double remainingTime = timeLimit - ((clock() - totalClock) * 1.0 / CLOCKS_PER_SEC); - CMR_CALL( CMRtestRegular(cmr, matrix, false, pisTotallyUnimodular, pdec, NULL, ¶ms->regular, - stats ? &stats->regular : NULL, remainingTime) ); - if (!*pisTotallyUnimodular && psubmatrix) + if (params->algorithm == CMR_TU_ALGORITHM_DECOMPOSITION) + { + + // TODO: run regularity check with ternary = true. + CMR_CALL( CMRtestRegular(cmr, matrix, false, pisTotallyUnimodular, pdec, NULL, ¶ms->regular, + stats ? &stats->regular : NULL, remainingTime) ); + + if (!*pisTotallyUnimodular && psubmatrix) + { + assert(!*psubmatrix); + remainingTime = timeLimit - (clock() - totalClock) * 1.0 / CLOCKS_PER_SEC; + CMR_CALL( CMRtestHereditaryPropertySimple(cmr, matrix, tuTest, stats, psubmatrix, remainingTime) ); + } + } + else if (params->algorithm == CMR_TU_ALGORITHM_SUBMATRIX) + { + assert(!"Not implemented"); + } + else if (params->algorithm == CMR_TU_ALGORITHM_PARTITION) + { + CMR_CALL( testPartition(cmr, matrix, pisTotallyUnimodular, remainingTime) ); + } + else { - assert(!*psubmatrix); - remainingTime = timeLimit - (clock() - totalClock) * 1.0 / CLOCKS_PER_SEC; - CMR_CALL( CMRtestHereditaryPropertySimple(cmr, matrix, tuTest, stats, psubmatrix, remainingTime) ); + return CMR_ERROR_INVALID; } if (stats) diff --git a/src/main/tu_main.c b/src/main/tu_main.c index 10c76022..1590c5ea 100644 --- a/src/main/tu_main.c +++ b/src/main/tu_main.c @@ -6,6 +6,7 @@ #include #include +#include typedef enum { @@ -26,6 +27,7 @@ CMR_ERROR testTotalUnimodularity( bool printStats, /**< Whether to print statistics to stderr. */ bool directGraphicness, /**< Whether to use fast graphicness routines. */ bool seriesParallel, /**< Whether to allow series-parallel operations in the decomposition tree. */ + CMR_TU_ALGORITHM algorithm, /**< Algorithm to use for TU test. */ double timeLimit /**< Time limit to impose. */ ) { @@ -61,6 +63,7 @@ CMR_ERROR testTotalUnimodularity( CMR_SUBMAT* submatrix = NULL; CMR_TU_PARAMETERS params; CMR_CALL( CMRparamsTotalUnimodularityInit(¶ms) ); + params.algorithm = algorithm; params.regular.completeTree = outputTreeFileName; params.regular.matrices = outputTreeFileName ? CMR_DEC_CONSTRUCT_ALL : CMR_DEC_CONSTRUCT_NONE; params.regular.directGraphicness = directGraphicness; @@ -83,9 +86,17 @@ CMR_ERROR testTotalUnimodularity( if (submatrix && outputSubmatrixFileName) { + /* Extract submatrix to compute its determinant. */ + CMR_CHRMAT* violator = NULL; + int64_t determinant = 0; + CMR_CALL( CMRchrmatZoomSubmat(cmr, matrix, submatrix, &violator) ); + CMR_CALL( CMRchrmatDeterminant(cmr, violator, &determinant) ); + CMR_CALL( CMRchrmatFree(cmr, &violator) ); + bool outputSubmatrixToFile = strcmp(outputSubmatrixFileName, "-"); - fprintf(stderr, "Writing minimal non-totally-unimodular submatrix to %s%s%s.\n", outputSubmatrixToFile ? "file <" : "", - outputSubmatrixToFile ? outputSubmatrixFileName : "stdout", outputSubmatrixToFile ? ">" : ""); + fprintf(stderr, "Writing minimal non-totally-unimodular submatrix with absolute determinant %ld to %s%s%s.\n", determinant, + outputSubmatrixToFile ? "file <" : "", outputSubmatrixToFile ? outputSubmatrixFileName : "stdout", + outputSubmatrixToFile ? ">" : ""); CMR_CALL( CMRsubmatWriteToFile(cmr, submatrix, matrix->numRows, matrix->numColumns, outputSubmatrixFileName) ); } @@ -113,15 +124,20 @@ int printUsage(const char* program) fputs(" determines whether the matrix given in file IN-MAT is totally unimodular.\n\n", stderr); fputs("Options:\n", stderr); fputs(" -i FORMAT Format of file IN-MAT, among `dense' and `sparse'; default: dense.\n", stderr); - fputs(" -D OUT-DEC Write a decomposition tree of the underlying regular matroid to file OUT-DEC; default: skip computation.\n", stderr); - fputs(" -N NON-SUB Write a minimal non-totally-unimodular submatrix to file NON-SUB; default: skip computation.\n", stderr); + fputs(" -D OUT-DEC Write a decomposition tree of the underlying regular matroid to file OUT-DEC; " + "default: skip computation.\n", stderr); + fputs(" -N NON-SUB Write a minimal non-totally-unimodular submatrix to file NON-SUB; default: skip computation.\n", + stderr); fputs("Advanced options:\n", stderr); fputs(" --stats Print statistics about the computation to stderr.\n\n", stderr); fputs(" --time-limit LIMIT Allow at most LIMIT seconds for the computation.\n", stderr); fputs(" --no-direct-graphic Check only 3-connected matrices for regularity.\n", stderr); fputs(" --no-series-parallel Do not allow series-parallel operations in decomposition tree.\n\n", stderr); + fputs(" --algo ALGO Use algorithm from {decomposition, submatrix, partition}; default: decomposition.\n\n", + stderr); fputs("If IN-MAT is `-' then the matrix is read from stdin.\n", stderr); - fputs("If OUT-DEC or NON-SUB is `-' then the decomposition tree (resp. the submatrix) is written to stdout.\n", stderr); + fputs("If OUT-DEC or NON-SUB is `-' then the decomposition tree (resp. the submatrix) is written to stdout.\n", + stderr); return EXIT_FAILURE; } @@ -136,6 +152,7 @@ int main(int argc, char** argv) bool directGraphicness = true; bool seriesParallel = true; double timeLimit = DBL_MAX; + CMR_TU_ALGORITHM algorithm = CMR_TU_ALGORITHM_DECOMPOSITION; for (int a = 1; a < argc; ++a) { if (!strcmp(argv[a], "-h")) @@ -175,6 +192,21 @@ int main(int argc, char** argv) } ++a; } + else if (!strcmp(argv[a], "--algo") && (a+1 < argc)) + { + if (!strcmp(argv[a+1], "decomposition")) + algorithm = CMR_TU_ALGORITHM_DECOMPOSITION; + else if (!strcmp(argv[a+1], "submatrix")) + algorithm = CMR_TU_ALGORITHM_SUBMATRIX; + else if (!strcmp(argv[a+1], "partition")) + algorithm = CMR_TU_ALGORITHM_PARTITION; + else + { + fprintf(stderr, "Error: Invalid algorithm <%s> specified.\n\n", argv[a+1]); + return printUsage(argv[0]); + } + ++a; + } else if (!inputMatrixFileName) inputMatrixFileName = argv[a]; else @@ -192,7 +224,7 @@ int main(int argc, char** argv) CMR_ERROR error; error = testTotalUnimodularity(inputMatrixFileName, inputFormat, outputTree, outputSubmatrix, printStats, - directGraphicness, seriesParallel, timeLimit); + directGraphicness, seriesParallel, algorithm, timeLimit); switch (error) { diff --git a/test/test_tu.cpp b/test/test_tu.cpp index 55893d0c..7c896a67 100644 --- a/test/test_tu.cpp +++ b/test/test_tu.cpp @@ -833,3 +833,84 @@ TEST(TotallyUnimodular, ForbiddenSubmatrix) ASSERT_CMR_CALL( CMRfreeEnvironment(&cmr) ); } + +TEST(TotallyUnimodular, PartitionAlgorithm) +{ + CMR* cmr = NULL; + ASSERT_CMR_CALL( CMRcreateEnvironment(&cmr) ); + + { + CMR_CHRMAT* K_3_3 = NULL; + ASSERT_CMR_CALL( stringToCharMatrix(cmr, &K_3_3, "5 4 " + " 1 1 0 0 " + " 1 1 1 0 " + " 1 0 0 -1 " + " 0 1 1 1 " + " 0 0 1 1 " + ) ); + + CMR_CHRMAT* K_3_3_dual = NULL; + ASSERT_CMR_CALL( stringToCharMatrix(cmr, &K_3_3_dual, "4 5 " + " 1 1 1 0 0 " + " 1 1 0 1 0 " + " 0 1 0 1 1 " + " 0 0 -1 1 1 " + ) ); + + CMR_CHRMAT* twoSum = NULL; + ASSERT_CMR_CALL( CMRtwoSum(cmr, K_3_3, K_3_3_dual, CMRrowToElement(1), CMRcolumnToElement(1), &twoSum) ); + + size_t rowPermutations[] = { 4, 6, 5, 7, 0, 1, 2, 3 }; + CMR_CHRMAT* matrix = NULL; + ASSERT_CMR_CALL( CMRchrmatPermute(cmr, twoSum, rowPermutations, NULL, &matrix) ); + ASSERT_CMR_CALL( CMRchrmatFree(cmr, &twoSum) ); + + CMRchrmatPrintDense(cmr, matrix, stdout, '0', true); + + bool isTU; + CMR_DEC* dec = NULL; + CMR_TU_PARAMETERS params; + ASSERT_CMR_CALL( CMRparamsTotalUnimodularityInit(¶ms) ); + params.algorithm = CMR_TU_ALGORITHM_PARTITION; + ASSERT_CMR_CALL( CMRtestTotalUnimodularity(cmr, matrix, &isTU, &dec, NULL, ¶ms, NULL, DBL_MAX) ); + ASSERT_EQ( dec, (CMR_DEC*) NULL ); + + ASSERT_TRUE( isTU ); + + ASSERT_CMR_CALL( CMRchrmatFree(cmr, &matrix) ); + ASSERT_CMR_CALL( CMRchrmatFree(cmr, &K_3_3) ); + ASSERT_CMR_CALL( CMRchrmatFree(cmr, &K_3_3_dual) ); + } + + { + CMR_CHRMAT* matrix = NULL; + ASSERT_CMR_CALL( stringToCharMatrix(cmr, &matrix, "14 14 " + "1 1 1 0 1 0 1 0 1 1 1 1 1 1 " + "1 0 1 0 1 0 1 0 1 1 1 1 1 0 " + "0 1 1 0 0 0 0 0 0 0 0 0 0 0 " + "0 1 1 1 0 0 0 0 0 0 0 0 0 0 " + "0 1 1 1 1 0 0 0 0 0 0 0 0 0 " + "0 1 1 1 1 1 0 0 0 0 0 0 0 0 " + "0 1 1 1 1 1 1 0 0 0 0 0 0 0 " + "0 1 1 1 1 1 1 1 0 0 0 0 0 0 " + "0 1 1 1 1 1 1 1 1 0 0 0 0 0 " + "0 0 0 0 0 0 0 0 1 1 0 0 0 0 " + "0 0 0 0 0 0 0 0 0 1 1 0 0 0 " + "0 0 0 0 0 0 0 0 0 0 1 1 0 0 " + "0 0 0 0 0 0 0 0 0 0 0 1 1 0 " + "0 0 0 0 0 0 0 0 0 0 0 0 1 1 " + ) ); + + CMRchrmatPrintDense(cmr, matrix, stdout, '0', true); + + bool isTU; + CMR_TU_PARAMETERS params; + ASSERT_CMR_CALL( CMRparamsTotalUnimodularityInit(¶ms) ); + params.algorithm = CMR_TU_ALGORITHM_PARTITION; + ASSERT_CMR_CALL( CMRtestTotalUnimodularity(cmr, matrix, &isTU, NULL, NULL, ¶ms, NULL, DBL_MAX) ); + ASSERT_FALSE( isTU ); + ASSERT_CMR_CALL( CMRchrmatFree(cmr, &matrix) ); + } + + ASSERT_CMR_CALL( CMRfreeEnvironment(&cmr) ); +}