diff --git a/CMakeLists.txt b/CMakeLists.txt index 61cf305..4ac318f 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,20 +1,29 @@ cmake_minimum_required(VERSION 2.8) project(multirow) +set(CMAKE_CXX_STANDARD 11) + set(CMAKE_LIBRARY_PATH ${CMAKE_LIBRARY_PATH} /Users/axavier/Applications/IBM/ILOG/CPLEX_Studio1262/cplex/lib/x86-64_osx/static_pic) list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_LIST_DIR}/cmake") find_package(CPLEX REQUIRED) include_directories(${CPLEX_INCLUDE_DIR}) -include_directories(multirow/include) -include_directories(lifting/library/include) -include_directories(infinity/library/include) +find_package(GMP REQUIRED) +find_package(OpenMP REQUIRED) -add_subdirectory(googletest/googletest) include_directories(${gtest_SOURCE_DIR}/include) +include_directories(infinity/library/include) +include_directories(lifting/library/include) +include_directories(multirow/include) +include_directories(onerow/library/include) +include_directories(qxx/include) -add_subdirectory(multirow) -add_subdirectory(lifting/library) -add_subdirectory(lifting/benchmark) -add_subdirectory(infinity/library) +add_subdirectory(googletest/googletest) add_subdirectory(infinity/benchmark) +add_subdirectory(infinity/library) +add_subdirectory(lifting/benchmark) +add_subdirectory(lifting/library) +add_subdirectory(multirow) +add_subdirectory(onerow/benchmark) +add_subdirectory(onerow/library) +add_subdirectory(qxx) diff --git a/cmake/FindGMP.cmake b/cmake/FindGMP.cmake new file mode 100644 index 0000000..0956d38 --- /dev/null +++ b/cmake/FindGMP.cmake @@ -0,0 +1,23 @@ +# Try to find the GMP librairies +# GMP_FOUND - system has GMP lib +# GMP_INCLUDE_DIR - the GMP include directory +# GMP_LIBRARIES - Libraries needed to use GMP + +# Copyright (c) 2006, Laurent Montel, +# +# Redistribution and use is allowed according to the terms of the BSD license. +# For details see the accompanying COPYING-CMAKE-SCRIPTS file. + +if (GMP_INCLUDE_DIR AND GMP_LIBRARIES) + # Already in cache, be silent + set(GMP_FIND_QUIETLY TRUE) +endif (GMP_INCLUDE_DIR AND GMP_LIBRARIES) + +find_path(GMP_INCLUDE_DIR NAMES gmp.h ) +find_library(GMP_LIBRARIES NAMES gmp libgmp ) +find_library(GMPXX_LIBRARIES NAMES gmpxx libgmpxx ) + +include(FindPackageHandleStandardArgs) +FIND_PACKAGE_HANDLE_STANDARD_ARGS(GMP DEFAULT_MSG GMP_INCLUDE_DIR GMP_LIBRARIES) + +mark_as_advanced(GMP_INCLUDE_DIR GMP_LIBRARIES) diff --git a/infinity/benchmark/CMakeLists.txt b/infinity/benchmark/CMakeLists.txt index d936db5..241314f 100644 --- a/infinity/benchmark/CMakeLists.txt +++ b/infinity/benchmark/CMakeLists.txt @@ -1,2 +1,7 @@ add_executable(infinity-benchmark.run src/main.c) -target_link_libraries (infinity-benchmark.run LINK_PUBLIC lifting_static infinity_static multirow_static) +target_link_libraries (infinity-benchmark.run LINK_PUBLIC + lifting_static + infinity_static + multirow_static + m + pthread) diff --git a/infinity/library/src/greedy-2d.c b/infinity/library/src/greedy-2d.c index b155288..8a8fd12 100644 --- a/infinity/library/src/greedy-2d.c +++ b/infinity/library/src/greedy-2d.c @@ -14,8 +14,6 @@ * along with this program. If not, see . */ -#define _GNU_SOURCE - #include #include #include diff --git a/infinity/library/tests/greedy-2d-test.cpp b/infinity/library/tests/greedy-2d-test.cpp index 05ce119..5469eed 100644 --- a/infinity/library/tests/greedy-2d-test.cpp +++ b/infinity/library/tests/greedy-2d-test.cpp @@ -15,6 +15,7 @@ */ #include +using namespace std; #define TEST_SOURCE diff --git a/lifting/benchmark/CMakeLists.txt b/lifting/benchmark/CMakeLists.txt index 1219149..ddaa4b3 100644 --- a/lifting/benchmark/CMakeLists.txt +++ b/lifting/benchmark/CMakeLists.txt @@ -1,2 +1,2 @@ add_executable(lifting-benchmark.run src/main.c) -target_link_libraries (lifting-benchmark.run LINK_PUBLIC lifting_static multirow_static) +target_link_libraries (lifting-benchmark.run LINK_PUBLIC lifting_static multirow_static m pthread) diff --git a/onerow/benchmark/CMakeLists.txt b/onerow/benchmark/CMakeLists.txt new file mode 100644 index 0000000..2405d69 --- /dev/null +++ b/onerow/benchmark/CMakeLists.txt @@ -0,0 +1,2 @@ +add_executable(onerow-benchmark.run src/main.cpp) +target_link_libraries (onerow-benchmark.run LINK_PUBLIC qxx_static onerow_static m pthread ${GMP_LIBRARIES} ${CPLEX_LIBRARIES}) diff --git a/onerow/benchmark/Makefile b/onerow/benchmark/Makefile new file mode 100644 index 0000000..88d8d81 --- /dev/null +++ b/onerow/benchmark/Makefile @@ -0,0 +1,97 @@ +miplib2010 := \ + instances/30n20b8.pre.done \ + instances/aflow40b.pre.done \ + instances/air04.pre.done \ + instances/app1-2.pre.done \ + instances/bab5.pre.done \ + instances/beasleyC3.pre.done \ + instances/biella1.pre.done \ + instances/bienst2.pre.done \ + instances/binkar10_1.pre.done \ + instances/bley_xl1.pre.done \ + instances/core2536-691.pre.done \ + instances/cov1075.pre.done \ + instances/csched010.pre.done \ + instances/danoint.pre.done \ + instances/dfn-gwin-UUM.pre.done \ + instances/eil33-2.pre.done \ + instances/eilB101.pre.done \ + instances/enlight13.pre.done \ + instances/ex9.pre.done \ + instances/glass4.pre.done \ + instances/gmu-35-40.pre.done \ + instances/iis-100-0-cov.pre.done \ + instances/iis-bupa-cov.pre.done \ + instances/iis-pima-cov.pre.done \ + instances/lectsched-4-obj.pre.done \ + instances/macrophage.pre.done \ + instances/map18.pre.done \ + instances/map20.pre.done \ + instances/mcsched.pre.done \ + instances/mik-250-1-100-1.pre.done \ + instances/mine-166-5.pre.done \ + instances/mine-90-10.pre.done \ + instances/msc98-ip.pre.done \ + instances/mspp16.pre.done \ + instances/mzzv11.pre.done \ + instances/n3div36.pre.done \ + instances/n3seq24.pre.done \ + instances/n4-3.pre.done \ + instances/neos-1109824.pre.done \ + instances/neos-1337307.pre.done \ + instances/neos-1396125.pre.done \ + instances/neos13.pre.done \ + instances/neos-1601936.pre.done \ + instances/neos18.pre.done \ + instances/neos-476283.pre.done \ + instances/neos-686190.pre.done \ + instances/neos-916792.pre.done \ + instances/neos-934278.pre.done \ + instances/net12.pre.done \ + instances/netdiversion.pre.done \ + instances/newdano.pre.done \ + instances/noswot.pre.done \ + instances/ns1208400.pre.done \ + instances/ns1688347.pre.done \ + instances/ns1758913.pre.done \ + instances/ns1830653.pre.done \ + instances/opm2-z7-s2.pre.done \ + instances/pg5_34.pre.done \ + instances/pigeon-10.pre.done \ + instances/pw-myciel4.pre.done \ + instances/qiu.pre.done \ + instances/rail507.pre.done \ + instances/ran16x16.pre.done \ + instances/reblock67.pre.done \ + instances/rmatr100-p10.pre.done \ + instances/rmatr100-p5.pre.done \ + instances/rmine6.pre.done \ + instances/rocII-4-11.pre.done \ + instances/rococoC10-001000.pre.done \ + instances/roll3000.pre.done \ + instances/sp98ic.pre.done \ + instances/sp98ir.pre.done \ + instances/tanglegram1.pre.done \ + instances/tanglegram2.pre.done \ + instances/timtab1.pre.done \ + instances/triptim1.pre.done \ + instances/unitcal_7.pre.done \ + instances/vpphard.pre.done \ + instances/zib54-UUE.pre.done + +mir := $(patsubst instances/%, mir/%, $(patsubst %.done,%.yaml,$(miplib3-all-instances))) +bases := $(patsubst instances/%, bases/%, $(patsubst %.done,%.bas,$(miplib3-all-instances))) +solutions := $(patsubst instances/%, solutions/%, $(patsubst %.done,%.x,$(miplib3-all-instances))) + +benchmark := ../../build/onerow/benchmark/onerow-benchmark.run + +.SECONDARY: + +miplib2010: $(miplib2010) + +instances/%.done: $(benchmark) out/%.yaml + @touch $@ + +out/%.yaml: $(benchmark) instances/%.mps.gz + @$(benchmark) -mw -f instances/$*.mps.gz -s out/$*.yaml 2>&1 | tee out/$*.log + diff --git a/onerow/benchmark/generate_tables b/onerow/benchmark/generate_tables new file mode 100755 index 0000000..e39b901 --- /dev/null +++ b/onerow/benchmark/generate_tables @@ -0,0 +1,29 @@ +#!/bin/bash +for i in out/*yaml; do + IN=${i/.yaml/} + IN=${IN/out\//} + grep -q mip_value $i && continue + grep $IN instances/opt.tab | awk '{ print "mip_value:\n "$2 }' >> $i +done + +echo " TABLE 1 " +echo "--------------------------------------------------------------------------------" +printf "%s,%s,%s,%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw origgap mirperf wperf mircontrib wcontrib wimprov +for i in out/*.yaml; do + IN=${i/.yaml/} + IN=${IN/out\//} + printf "%s," $IN + scripts/gap.rb $i +done | sort -t',' -nrsk 8 | sed -e 's/,$//g' + +echo +echo +echo " TABLE 2 " +echo "--------------------------------------------------------------------------------" +printf "%s,%s,%s,%s,%s,%s\n" instance cutsmir cutsw mirt wedget avgm +for i in out/*.yaml; do + IN=${i/.yaml/} + IN=${IN/out\//} + printf "%s," $IN + scripts/speed.rb $i +done | sort | sed -e 's/,$//g' diff --git a/onerow/benchmark/instances/30n20b8.pre.mps.gz b/onerow/benchmark/instances/30n20b8.pre.mps.gz new file mode 100644 index 0000000..edad85d Binary files /dev/null and b/onerow/benchmark/instances/30n20b8.pre.mps.gz differ diff --git a/onerow/benchmark/instances/aflow40b.pre.mps.gz b/onerow/benchmark/instances/aflow40b.pre.mps.gz new file mode 100644 index 0000000..d9f205e Binary files /dev/null and b/onerow/benchmark/instances/aflow40b.pre.mps.gz differ diff --git a/onerow/benchmark/instances/air04.pre.mps.gz b/onerow/benchmark/instances/air04.pre.mps.gz new file mode 100644 index 0000000..476f4bb Binary files /dev/null and b/onerow/benchmark/instances/air04.pre.mps.gz differ diff --git a/onerow/benchmark/instances/app1-2.pre.mps.gz b/onerow/benchmark/instances/app1-2.pre.mps.gz new file mode 100644 index 0000000..479198b Binary files /dev/null and b/onerow/benchmark/instances/app1-2.pre.mps.gz differ diff --git a/onerow/benchmark/instances/bab5.pre.mps.gz b/onerow/benchmark/instances/bab5.pre.mps.gz new file mode 100644 index 0000000..5cfee60 Binary files /dev/null and b/onerow/benchmark/instances/bab5.pre.mps.gz differ diff --git a/onerow/benchmark/instances/beasleyC3.pre.mps.gz b/onerow/benchmark/instances/beasleyC3.pre.mps.gz new file mode 100644 index 0000000..0f43cca Binary files /dev/null and b/onerow/benchmark/instances/beasleyC3.pre.mps.gz differ diff --git a/onerow/benchmark/instances/biella1.pre.mps.gz b/onerow/benchmark/instances/biella1.pre.mps.gz new file mode 100644 index 0000000..bbcead9 Binary files /dev/null and b/onerow/benchmark/instances/biella1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/bienst2.pre.mps.gz b/onerow/benchmark/instances/bienst2.pre.mps.gz new file mode 100644 index 0000000..8d3b452 Binary files /dev/null and b/onerow/benchmark/instances/bienst2.pre.mps.gz differ diff --git a/onerow/benchmark/instances/binkar10_1.pre.mps.gz b/onerow/benchmark/instances/binkar10_1.pre.mps.gz new file mode 100644 index 0000000..3af092c Binary files /dev/null and b/onerow/benchmark/instances/binkar10_1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/bley_xl1.pre.mps.gz b/onerow/benchmark/instances/bley_xl1.pre.mps.gz new file mode 100644 index 0000000..f4409d5 Binary files /dev/null and b/onerow/benchmark/instances/bley_xl1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/core2536-691.pre.mps.gz b/onerow/benchmark/instances/core2536-691.pre.mps.gz new file mode 100644 index 0000000..bc88c3f Binary files /dev/null and b/onerow/benchmark/instances/core2536-691.pre.mps.gz differ diff --git a/onerow/benchmark/instances/cov1075.pre.mps.gz b/onerow/benchmark/instances/cov1075.pre.mps.gz new file mode 100644 index 0000000..4abeeeb Binary files /dev/null and b/onerow/benchmark/instances/cov1075.pre.mps.gz differ diff --git a/onerow/benchmark/instances/csched010.pre.mps.gz b/onerow/benchmark/instances/csched010.pre.mps.gz new file mode 100644 index 0000000..8719085 Binary files /dev/null and b/onerow/benchmark/instances/csched010.pre.mps.gz differ diff --git a/onerow/benchmark/instances/danoint.pre.mps.gz b/onerow/benchmark/instances/danoint.pre.mps.gz new file mode 100644 index 0000000..440ab0a Binary files /dev/null and b/onerow/benchmark/instances/danoint.pre.mps.gz differ diff --git a/onerow/benchmark/instances/dfn-gwin-UUM.pre.mps.gz b/onerow/benchmark/instances/dfn-gwin-UUM.pre.mps.gz new file mode 100644 index 0000000..7bf9bda Binary files /dev/null and b/onerow/benchmark/instances/dfn-gwin-UUM.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/infeasible/ash608gpia-3col.pre.mps.gz b/onerow/benchmark/instances/disabled/infeasible/ash608gpia-3col.pre.mps.gz new file mode 100644 index 0000000..f5f03bf Binary files /dev/null and b/onerow/benchmark/instances/disabled/infeasible/ash608gpia-3col.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/infeasible/enlight14.pre.mps.gz b/onerow/benchmark/instances/disabled/infeasible/enlight14.pre.mps.gz new file mode 100644 index 0000000..d11c500 Binary files /dev/null and b/onerow/benchmark/instances/disabled/infeasible/enlight14.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/infeasible/ns1766074.pre.mps.gz b/onerow/benchmark/instances/disabled/infeasible/ns1766074.pre.mps.gz new file mode 100644 index 0000000..a8e5f82 Binary files /dev/null and b/onerow/benchmark/instances/disabled/infeasible/ns1766074.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/no-gap/acc-tight5.pre.mps.gz b/onerow/benchmark/instances/disabled/no-gap/acc-tight5.pre.mps.gz new file mode 100644 index 0000000..a188ad7 Binary files /dev/null and b/onerow/benchmark/instances/disabled/no-gap/acc-tight5.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/no-gap/bnatt350.pre.mps.gz b/onerow/benchmark/instances/disabled/no-gap/bnatt350.pre.mps.gz new file mode 100644 index 0000000..f3da9cf Binary files /dev/null and b/onerow/benchmark/instances/disabled/no-gap/bnatt350.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/no-gap/m100n500k4r1.pre.mps.gz b/onerow/benchmark/instances/disabled/no-gap/m100n500k4r1.pre.mps.gz new file mode 100644 index 0000000..81a2a19 Binary files /dev/null and b/onerow/benchmark/instances/disabled/no-gap/m100n500k4r1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/disabled/no-gap/neos-849702.pre.mps.gz b/onerow/benchmark/instances/disabled/no-gap/neos-849702.pre.mps.gz new file mode 100644 index 0000000..7c14490 Binary files /dev/null and b/onerow/benchmark/instances/disabled/no-gap/neos-849702.pre.mps.gz differ diff --git a/onerow/benchmark/instances/eil33-2.pre.mps.gz b/onerow/benchmark/instances/eil33-2.pre.mps.gz new file mode 100644 index 0000000..c45cab8 Binary files /dev/null and b/onerow/benchmark/instances/eil33-2.pre.mps.gz differ diff --git a/onerow/benchmark/instances/eilB101.pre.mps.gz b/onerow/benchmark/instances/eilB101.pre.mps.gz new file mode 100644 index 0000000..8d5722d Binary files /dev/null and b/onerow/benchmark/instances/eilB101.pre.mps.gz differ diff --git a/onerow/benchmark/instances/enlight13.pre.mps.gz b/onerow/benchmark/instances/enlight13.pre.mps.gz new file mode 100644 index 0000000..e460b9a Binary files /dev/null and b/onerow/benchmark/instances/enlight13.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ex9.pre.mps.gz b/onerow/benchmark/instances/ex9.pre.mps.gz new file mode 100644 index 0000000..f1db4ea Binary files /dev/null and b/onerow/benchmark/instances/ex9.pre.mps.gz differ diff --git a/onerow/benchmark/instances/glass4.pre.mps.gz b/onerow/benchmark/instances/glass4.pre.mps.gz new file mode 100644 index 0000000..b66e06e Binary files /dev/null and b/onerow/benchmark/instances/glass4.pre.mps.gz differ diff --git a/onerow/benchmark/instances/gmu-35-40.pre.mps.gz b/onerow/benchmark/instances/gmu-35-40.pre.mps.gz new file mode 100644 index 0000000..5e2f817 Binary files /dev/null and b/onerow/benchmark/instances/gmu-35-40.pre.mps.gz differ diff --git a/onerow/benchmark/instances/iis-100-0-cov.pre.mps.gz b/onerow/benchmark/instances/iis-100-0-cov.pre.mps.gz new file mode 100644 index 0000000..d465b04 Binary files /dev/null and b/onerow/benchmark/instances/iis-100-0-cov.pre.mps.gz differ diff --git a/onerow/benchmark/instances/iis-bupa-cov.pre.mps.gz b/onerow/benchmark/instances/iis-bupa-cov.pre.mps.gz new file mode 100644 index 0000000..15b8610 Binary files /dev/null and b/onerow/benchmark/instances/iis-bupa-cov.pre.mps.gz differ diff --git a/onerow/benchmark/instances/iis-pima-cov.pre.mps.gz b/onerow/benchmark/instances/iis-pima-cov.pre.mps.gz new file mode 100644 index 0000000..5e1818a Binary files /dev/null and b/onerow/benchmark/instances/iis-pima-cov.pre.mps.gz differ diff --git a/onerow/benchmark/instances/lectsched-4-obj.pre.mps.gz b/onerow/benchmark/instances/lectsched-4-obj.pre.mps.gz new file mode 100644 index 0000000..352525e Binary files /dev/null and b/onerow/benchmark/instances/lectsched-4-obj.pre.mps.gz differ diff --git a/onerow/benchmark/instances/macrophage.pre.mps.gz b/onerow/benchmark/instances/macrophage.pre.mps.gz new file mode 100644 index 0000000..909c5b0 Binary files /dev/null and b/onerow/benchmark/instances/macrophage.pre.mps.gz differ diff --git a/onerow/benchmark/instances/map18.pre.mps.gz b/onerow/benchmark/instances/map18.pre.mps.gz new file mode 100644 index 0000000..774a795 Binary files /dev/null and b/onerow/benchmark/instances/map18.pre.mps.gz differ diff --git a/onerow/benchmark/instances/map20.pre.mps.gz b/onerow/benchmark/instances/map20.pre.mps.gz new file mode 100644 index 0000000..e6bd46a Binary files /dev/null and b/onerow/benchmark/instances/map20.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mcsched.pre.mps.gz b/onerow/benchmark/instances/mcsched.pre.mps.gz new file mode 100644 index 0000000..fe21fa8 Binary files /dev/null and b/onerow/benchmark/instances/mcsched.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mik-250-1-100-1.pre.mps.gz b/onerow/benchmark/instances/mik-250-1-100-1.pre.mps.gz new file mode 100644 index 0000000..44a2c6c Binary files /dev/null and b/onerow/benchmark/instances/mik-250-1-100-1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mine-166-5.pre.mps.gz b/onerow/benchmark/instances/mine-166-5.pre.mps.gz new file mode 100644 index 0000000..5672f2d Binary files /dev/null and b/onerow/benchmark/instances/mine-166-5.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mine-90-10.pre.mps.gz b/onerow/benchmark/instances/mine-90-10.pre.mps.gz new file mode 100644 index 0000000..221b985 Binary files /dev/null and b/onerow/benchmark/instances/mine-90-10.pre.mps.gz differ diff --git a/onerow/benchmark/instances/msc98-ip.pre.mps.gz b/onerow/benchmark/instances/msc98-ip.pre.mps.gz new file mode 100644 index 0000000..fa232bd Binary files /dev/null and b/onerow/benchmark/instances/msc98-ip.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mspp16.pre.mps.gz b/onerow/benchmark/instances/mspp16.pre.mps.gz new file mode 100644 index 0000000..b9b175b Binary files /dev/null and b/onerow/benchmark/instances/mspp16.pre.mps.gz differ diff --git a/onerow/benchmark/instances/mzzv11.pre.mps.gz b/onerow/benchmark/instances/mzzv11.pre.mps.gz new file mode 100644 index 0000000..f669584 Binary files /dev/null and b/onerow/benchmark/instances/mzzv11.pre.mps.gz differ diff --git a/onerow/benchmark/instances/n3div36.pre.mps.gz b/onerow/benchmark/instances/n3div36.pre.mps.gz new file mode 100644 index 0000000..b68dbab Binary files /dev/null and b/onerow/benchmark/instances/n3div36.pre.mps.gz differ diff --git a/onerow/benchmark/instances/n3seq24.pre.mps.gz b/onerow/benchmark/instances/n3seq24.pre.mps.gz new file mode 100644 index 0000000..c184461 Binary files /dev/null and b/onerow/benchmark/instances/n3seq24.pre.mps.gz differ diff --git a/onerow/benchmark/instances/n4-3.pre.mps.gz b/onerow/benchmark/instances/n4-3.pre.mps.gz new file mode 100644 index 0000000..a20efa4 Binary files /dev/null and b/onerow/benchmark/instances/n4-3.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-1109824.pre.mps.gz b/onerow/benchmark/instances/neos-1109824.pre.mps.gz new file mode 100644 index 0000000..b46d7d0 Binary files /dev/null and b/onerow/benchmark/instances/neos-1109824.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-1337307.pre.mps.gz b/onerow/benchmark/instances/neos-1337307.pre.mps.gz new file mode 100644 index 0000000..1f122dc Binary files /dev/null and b/onerow/benchmark/instances/neos-1337307.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-1396125.pre.mps.gz b/onerow/benchmark/instances/neos-1396125.pre.mps.gz new file mode 100644 index 0000000..a4db6b5 Binary files /dev/null and b/onerow/benchmark/instances/neos-1396125.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-1601936.pre.mps.gz b/onerow/benchmark/instances/neos-1601936.pre.mps.gz new file mode 100644 index 0000000..c1d10dc Binary files /dev/null and b/onerow/benchmark/instances/neos-1601936.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-476283.pre.mps.gz b/onerow/benchmark/instances/neos-476283.pre.mps.gz new file mode 100644 index 0000000..e468201 Binary files /dev/null and b/onerow/benchmark/instances/neos-476283.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-686190.pre.mps.gz b/onerow/benchmark/instances/neos-686190.pre.mps.gz new file mode 100644 index 0000000..d4ec957 Binary files /dev/null and b/onerow/benchmark/instances/neos-686190.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-916792.pre.mps.gz b/onerow/benchmark/instances/neos-916792.pre.mps.gz new file mode 100644 index 0000000..82bf309 Binary files /dev/null and b/onerow/benchmark/instances/neos-916792.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos-934278.pre.mps.gz b/onerow/benchmark/instances/neos-934278.pre.mps.gz new file mode 100644 index 0000000..bbbe3c7 Binary files /dev/null and b/onerow/benchmark/instances/neos-934278.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos13.pre.mps.gz b/onerow/benchmark/instances/neos13.pre.mps.gz new file mode 100644 index 0000000..d6ce32e Binary files /dev/null and b/onerow/benchmark/instances/neos13.pre.mps.gz differ diff --git a/onerow/benchmark/instances/neos18.pre.mps.gz b/onerow/benchmark/instances/neos18.pre.mps.gz new file mode 100644 index 0000000..b9161c6 Binary files /dev/null and b/onerow/benchmark/instances/neos18.pre.mps.gz differ diff --git a/onerow/benchmark/instances/net12.pre.mps.gz b/onerow/benchmark/instances/net12.pre.mps.gz new file mode 100644 index 0000000..08d610a Binary files /dev/null and b/onerow/benchmark/instances/net12.pre.mps.gz differ diff --git a/onerow/benchmark/instances/netdiversion.pre.mps.gz b/onerow/benchmark/instances/netdiversion.pre.mps.gz new file mode 100644 index 0000000..d06b611 Binary files /dev/null and b/onerow/benchmark/instances/netdiversion.pre.mps.gz differ diff --git a/onerow/benchmark/instances/newdano.pre.mps.gz b/onerow/benchmark/instances/newdano.pre.mps.gz new file mode 100644 index 0000000..a59860e Binary files /dev/null and b/onerow/benchmark/instances/newdano.pre.mps.gz differ diff --git a/onerow/benchmark/instances/noswot.pre.mps.gz b/onerow/benchmark/instances/noswot.pre.mps.gz new file mode 100644 index 0000000..298ff9e Binary files /dev/null and b/onerow/benchmark/instances/noswot.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ns1208400.pre.mps.gz b/onerow/benchmark/instances/ns1208400.pre.mps.gz new file mode 100644 index 0000000..4eb141d Binary files /dev/null and b/onerow/benchmark/instances/ns1208400.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ns1688347.pre.mps.gz b/onerow/benchmark/instances/ns1688347.pre.mps.gz new file mode 100644 index 0000000..3255012 Binary files /dev/null and b/onerow/benchmark/instances/ns1688347.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ns1758913.pre.mps.gz b/onerow/benchmark/instances/ns1758913.pre.mps.gz new file mode 100644 index 0000000..edf1094 Binary files /dev/null and b/onerow/benchmark/instances/ns1758913.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ns1830653.pre.mps.gz b/onerow/benchmark/instances/ns1830653.pre.mps.gz new file mode 100644 index 0000000..d65fe93 Binary files /dev/null and b/onerow/benchmark/instances/ns1830653.pre.mps.gz differ diff --git a/onerow/benchmark/instances/opm2-z7-s2.pre.mps.gz b/onerow/benchmark/instances/opm2-z7-s2.pre.mps.gz new file mode 100644 index 0000000..44d692f Binary files /dev/null and b/onerow/benchmark/instances/opm2-z7-s2.pre.mps.gz differ diff --git a/onerow/benchmark/instances/opt.tab b/onerow/benchmark/instances/opt.tab new file mode 100644 index 0000000..5aa0bce --- /dev/null +++ b/onerow/benchmark/instances/opt.tab @@ -0,0 +1,362 @@ +30_70_45_095_100 3.000000 +30n20b8 302.000000 +50v-10 3311.180000 +a1c1s1 11503.400000 +acc-tight4 0.000000 +acc-tight5 0.000000 +acc-tight6 0.000000 +aflow40b 1168.000000 +air04 56137.000000 +app1-2 -41.000000 +ash608gpia-3col 100000000000000000000.000000 +atlanta-ip 90.009900 +atm20-100 2463620.000000 +b2c1s1 25687.900000 +bab1 100000000000000000000.000000 +bab3 100000000000000000000.000000 +bab5 -106412.000000 +beasleyC3 754.000000 +berlin_5_8_0 62.000000 +bg512142 184203.000000 +biella1 3065010.000000 +bienst2 54.600000 +binkar10_1 6742.200000 +bley_xl1 190.000000 +blp-ar98 6205.210000 +blp-ic97 4025.020000 +bnatt350 0.000000 +bnatt400 1.000000 +buildingenergy 33283.900000 +cdma 100000000000000000000.000000 +circ10-3 100000000000000000000.000000 +co-100 2639940.000000 +core2536-691 689.000000 +core4872-1529 100000000000000000000.000000 +cov1075 20.000000 +csched007 351.000000 +csched008 173.000000 +csched010 408.000000 +d10200 12430.000000 +d20200 100000000000000000000.000000 +dano3mip 100000000000000000000.000000 +danoint 65.666700 +datt256 100000000000000000000.000000 +dc1c 1767900.000000 +dc1l 100000000000000000000.000000 +dfn-gwin-UUM 38752.000000 +dg012142 2300870.000000 +dolom1 6609250.000000 +ds-big 100000000000000000000.000000 +eil33-2 934.008000 +eilA101-2 880.920000 +eilB101 1216.920000 +enlight13 71.000000 +enlight14 100000000000000000000.000000 +enlight15 69.000000 +enlight16 100000000000000000000.000000 +enlight9 100000000000000000000.000000 +ex1010-pi 100000000000000000000.000000 +ex10 100.000000 +ex9 81.000000 +f2000 100000000000000000000.000000 +g200x740i 30086.000000 +ger50_17_trans 100000000000000000000.000000 +germanrr 47095900.000000 +germany50-DBM 473840.000000 +glass4 1200010000.000000 +gmu-35-40 -2406730.000000 +gmu-35-50 -2607960.000000 +gmut-75-50 100000000000000000000.000000 +gmut-77-40 100000000000000000000.000000 +go19 84.000000 +hanoi5 1931.000000 +harp2 -73899800.000000 +hawaiiv10-130 100000000000000000000.000000 +ic97_potential 3942.000000 +iis-100-0-cov 29.000000 +iis-bupa-cov 36.000000 +iis-pima-cov 33.000000 +in 58.000000 +ivu06-big 100000000000000000000.000000 +ivu52 100000000000000000000.000000 +janos-us-DDM 1492710.000000 +k16x240 10674.000000 +lectsched-1 0.000000 +lectsched-1-obj 100000000000000000000.000000 +lectsched-2 0.000000 +lectsched-3 0.000000 +lectsched-4-obj 4.000000 +leo1 404228000.000000 +leo2 404077000.000000 +liu 100000000000000000000.000000 +lotsize 1480200.000000 +lrsa120 100000000000000000000.000000 +m100n500k4r1 -25.000000 +macrophage 374.000000 +map06 -289.000000 +map10 -495.000000 +map14 -674.000000 +map18 -847.000000 +map20 -922.000000 +markshare_5_0 1.000000 +maxgasflow -44565800.000000 +mc11 11689.000000 +mcsched 211913.000000 +methanosarcina 100000000000000000000.000000 +mik-250-1-100-1 -66729.000000 +mine-166-5 -566396000.000000 +mine-90-10 -784302000.000000 +mining 100000000000000000000.000000 +mkc -563.846000 +momentum1 109143.000000 +momentum2 12314.100000 +momentum3 100000000000000000000.000000 +msc98-ip 19839500.000000 +mspp16 363.000000 +mzzv11 -21718.000000 +n15-3 100000000000000000000.000000 +n3-3 15915.000000 +n3700 100000000000000000000.000000 +n3705 100000000000000000000.000000 +n370a 100000000000000000000.000000 +n3div36 130800.000000 +n3seq24 52200.000000 +n4-3 8993.000000 +n9-3 14409.000000 +nag 100000000000000000000.000000 +nb10tb 100000000000000000000.000000 +neos-1109824 378.000000 +neos-1112782 572103000000.000000 +neos-1112787 565032000000.000000 +neos-1140050 100000000000000000000.000000 +neos-1171692 -273.000000 +neos-1171737 -195.000000 +neos-1224597 -428.000000 +neos-1225589 1231070000.000000 +neos-1311124 100000000000000000000.000000 +neos-1337307 -202319.000000 +neos-1396125 3000.050000 +neos13 -95.474800 +neos-1426635 -176.000000 +neos-1426662 -44.000000 +neos-1429212 100000000000000000000.000000 +neos-1436709 -128.000000 +neos-1440225 36.000000 +neos-1440460 -179.250000 +neos-1442119 -181.000000 +neos-1442657 -154.500000 +neos15 80598.400000 +neos-1601936 3.000000 +neos-1605061 12.000000 +neos-1605075 9.000000 +neos-1616732 159.000000 +neos-1620770 9.000000 +neos16 446.000000 +neos18 16.000000 +neos-476283 406.363000 +neos-506422 0.000000 +neos-506428 583780.000000 +neos-520729 -1385000.000000 +neos-555424 1286800.000000 +neos-631710 203.000000 +neos-686190 6730.000000 +neos-693347 234.000000 +neos6 83.000000 +neos-738098 -1099.000000 +neos-777800 -80.000000 +neos-785912 100000000000000000000.000000 +neos788725 100000000000000000000.000000 +neos-799711 -11170200.000000 +neos-807456 100000000000000000000.000000 +neos808444 0.000000 +neos-820146 100000000000000000000.000000 +neos-820157 100000000000000000000.000000 +neos-824661 33.000000 +neos-824695 31.000000 +neos-826650 29.000000 +neos-826694 58.000000 +neos-826812 58.011000 +neos-826841 29.008200 +neos-847302 4.000000 +neos-849702 0.000000 +neos858960 100000000000000000000.000000 +neos-859770 100000000000000000000.000000 +neos-885086 -243.000000 +neos-885524 12320.100000 +neos-911880 54.760000 +neos-916792 31.870400 +neos-932816 15376.000000 +neos-933638 276.000000 +neos-933966 318.000000 +neos-934278 260.000000 +neos-935627 2598.000000 +neos-935769 3010.000000 +neos-937511 3510.000000 +neos-937815 100000000000000000000.000000 +neos-941262 2791.000000 +neos-941313 9361.000000 +neos-942830 16.000000 +neos-948126 2607.000000 +neos-952987 100000000000000000000.000000 +neos-957389 1.500000 +neos-984165 100000000000000000000.000000 +net12 214.000000 +netdiversion 242.000000 +newdano 65.666700 +nobel-eu-DBE 608910.000000 +noswot -41.000000 +npmv07 104810000000.000000 +ns1111636 162.000000 +ns1116954 0.000000 +ns1158817 100000000000000000000.000000 +ns1208400 2.000000 +ns1456591 100000000000000000000.000000 +ns1606230 21.000000 +ns1631475 100000000000000000000.000000 +ns1644855 -1524.330000 +ns1663818 86.000000 +ns1685374 -13.000000 +ns1686196 100000000000000000000.000000 +ns1688347 27.000000 +ns1696083 45.000000 +ns1702808 100000000000000000000.000000 +ns1745726 100000000000000000000.000000 +ns1758913 -1454.670000 +ns1766074 100000000000000000000.000000 +ns1769397 100000000000000000000.000000 +ns1778858 100000000000000000000.000000 +ns1830653 20622.000000 +ns1853823 100000000000000000000.000000 +ns1854840 100000000000000000000.000000 +ns1856153 100000000000000000000.000000 +ns1904248 100000000000000000000.000000 +ns1905797 100000000000000000000.000000 +ns1905800 100000000000000000000.000000 +ns1952667 0.000000 +ns2017839 77030500000000.000000 +ns2081729 9.000000 +ns2118727 100000000000000000000.000000 +ns2122603 100000000000000000000.000000 +ns2124243 100000000000000000000.000000 +ns2137859 100000000000000000000.000000 +ns4-pr3 100000000000000000000.000000 +ns4-pr9 100000000000000000000.000000 +ns894236 100000000000000000000.000000 +ns894244 15.000000 +ns894786 100000000000000000000.000000 +ns894788 7.000000 +ns903616 100000000000000000000.000000 +ns930473 100000000000000000000.000000 +nsr8k 100000000000000000000.000000 +nu120-pr3 28130.000000 +nu60-pr9 24940.000000 +ofi 6155380000.000000 +opm2-z10-s2 -33826.000000 +opm2-z11-s8 -43485.000000 +opm2-z12-s14 -64291.000000 +opm2-z12-s7 -65514.000000 +opm2-z7-s2 -10280.000000 +p100x588b 47878.000000 +p2m2p1m1p0n100 100000000000000000000.000000 +p6b -63.000000 +p80x400b 39667.000000 +pb-simp-nonunif 100000000000000000000.000000 +pg5_34 -14339.400000 +pg -8674.340000 +pigeon-10 -9000.000000 +pigeon-11 -10000.000000 +pigeon-12 -11000.000000 +pigeon-13 -12000.000000 +pigeon-19 100000000000000000000.000000 +probportfolio 16.734200 +protfold -31.000000 +pw-myciel4 10.000000 +qiu -132.873000 +queens-30 -40.000000 +r80x800 5332.000000 +rail01 -70.570000 +rail02 -200.450000 +rail03 -867.094000 +rail507 174.000000 +ramos3 100000000000000000000.000000 +ran14x18-disj-8 3712.000000 +ran14x18 3712.000000 +ran16x16 3823.000000 +reblock166 -600052000.000000 +reblock354 -39280500.000000 +reblock420 -517793000.000000 +reblock67 -34630600.000000 +rmatr100-p10 423.000000 +rmatr100-p5 976.000000 +rmatr200-p10 2017.000000 +rmatr200-p20 837.000000 +rmatr200-p5 4521.000000 +rmine10 100000000000000000000.000000 +rmine14 100000000000000000000.000000 +rmine21 100000000000000000000.000000 +rmine25 100000000000000000000.000000 +rmine6 -457.186000 +rocII-4-11 -6.655640 +rocII-7-11 100000000000000000000.000000 +rocII-9-11 100000000000000000000.000000 +rococoB10-011000 19449.000000 +rococoC10-001000 11460.000000 +rococoC11-011100 20889.000000 +rococoC12-111000 100000000000000000000.000000 +roll3000 12890.000000 +rvb-sub 100000000000000000000.000000 +satellites1-25 -5.000000 +satellites2-60-fs -19.000000 +satellites2-60 -19.000000 +satellites3-40-fs -25.000000 +satellites3-40 -25.000000 +sct1 100000000000000000000.000000 +sct32 100000000000000000000.000000 +sct5 100000000000000000000.000000 +set3-10 100000000000000000000.000000 +set3-15 124886.000000 +set3-20 100000000000000000000.000000 +seymour 423.000000 +seymour-disj-10 287.000000 +shipsched 100000000000000000000.000000 +shs1023 13136.600000 +siena1 100000000000000000000.000000 +sing161 100000000000000000000.000000 +sing245 100000000000000000000.000000 +sing2 100000000000000000000.000000 +sing359 100000000000000000000.000000 +sp97ar 660706000.000000 +sp98ic 449145000.000000 +sp98ir 219677000.000000 +splan1 100000000000000000000.000000 +stockholm 100000000000000000000.000000 +stp3d 493.720000 +sts405 100000000000000000000.000000 +sts729 100000000000000000000.000000 +swath 467.407000 +t1717 100000000000000000000.000000 +t1722 100000000000000000000.000000 +tanglegram1 5182.000000 +tanglegram2 443.000000 +timtab1 764772.000000 +toll-like 610.000000 +transportmoment -3063100000.000000 +triptim1 22.868100 +triptim2 100000000000000000000.000000 +triptim3 100000000000000000000.000000 +tw-myciel4 10.000000 +uc-case11 100000000000000000000.000000 +uc-case3 7204.920000 +uct-subprob 314.000000 +umts 30090300.000000 +unitcal_7 19635600.000000 +usAbbrv-8-25_70 100000000000000000000.000000 +van 100000000000000000000.000000 +vpphard2 81.000000 +vpphard 5.000000 +wachplan -8.000000 +wnq-n100-mw99-14 259.000000 +zib01 100000000000000000000.000000 +zib02 100000000000000000000.000000 +zib54-UUE 10334000.000000 + diff --git a/onerow/benchmark/instances/pg5_34.pre.mps.gz b/onerow/benchmark/instances/pg5_34.pre.mps.gz new file mode 100644 index 0000000..46b4a52 Binary files /dev/null and b/onerow/benchmark/instances/pg5_34.pre.mps.gz differ diff --git a/onerow/benchmark/instances/pigeon-10.pre.mps.gz b/onerow/benchmark/instances/pigeon-10.pre.mps.gz new file mode 100644 index 0000000..0f9c7dd Binary files /dev/null and b/onerow/benchmark/instances/pigeon-10.pre.mps.gz differ diff --git a/onerow/benchmark/instances/pw-myciel4.pre.mps.gz b/onerow/benchmark/instances/pw-myciel4.pre.mps.gz new file mode 100644 index 0000000..0048db9 Binary files /dev/null and b/onerow/benchmark/instances/pw-myciel4.pre.mps.gz differ diff --git a/onerow/benchmark/instances/qiu.pre.mps.gz b/onerow/benchmark/instances/qiu.pre.mps.gz new file mode 100644 index 0000000..6a7836d Binary files /dev/null and b/onerow/benchmark/instances/qiu.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rail507.pre.mps.gz b/onerow/benchmark/instances/rail507.pre.mps.gz new file mode 100644 index 0000000..989d3d0 Binary files /dev/null and b/onerow/benchmark/instances/rail507.pre.mps.gz differ diff --git a/onerow/benchmark/instances/ran16x16.pre.mps.gz b/onerow/benchmark/instances/ran16x16.pre.mps.gz new file mode 100644 index 0000000..a0252d7 Binary files /dev/null and b/onerow/benchmark/instances/ran16x16.pre.mps.gz differ diff --git a/onerow/benchmark/instances/reblock67.pre.mps.gz b/onerow/benchmark/instances/reblock67.pre.mps.gz new file mode 100644 index 0000000..6be6864 Binary files /dev/null and b/onerow/benchmark/instances/reblock67.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rmatr100-p10.pre.mps.gz b/onerow/benchmark/instances/rmatr100-p10.pre.mps.gz new file mode 100644 index 0000000..c381742 Binary files /dev/null and b/onerow/benchmark/instances/rmatr100-p10.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rmatr100-p5.pre.mps.gz b/onerow/benchmark/instances/rmatr100-p5.pre.mps.gz new file mode 100644 index 0000000..9d45cc3 Binary files /dev/null and b/onerow/benchmark/instances/rmatr100-p5.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rmine6.pre.mps.gz b/onerow/benchmark/instances/rmine6.pre.mps.gz new file mode 100644 index 0000000..6198e3c Binary files /dev/null and b/onerow/benchmark/instances/rmine6.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rocII-4-11.pre.mps.gz b/onerow/benchmark/instances/rocII-4-11.pre.mps.gz new file mode 100644 index 0000000..c8b0b5a Binary files /dev/null and b/onerow/benchmark/instances/rocII-4-11.pre.mps.gz differ diff --git a/onerow/benchmark/instances/rococoC10-001000.pre.mps.gz b/onerow/benchmark/instances/rococoC10-001000.pre.mps.gz new file mode 100644 index 0000000..07f9faa Binary files /dev/null and b/onerow/benchmark/instances/rococoC10-001000.pre.mps.gz differ diff --git a/onerow/benchmark/instances/roll3000.pre.mps.gz b/onerow/benchmark/instances/roll3000.pre.mps.gz new file mode 100644 index 0000000..9d86f90 Binary files /dev/null and b/onerow/benchmark/instances/roll3000.pre.mps.gz differ diff --git a/onerow/benchmark/instances/sp98ic.pre.mps.gz b/onerow/benchmark/instances/sp98ic.pre.mps.gz new file mode 100644 index 0000000..bb9a65a Binary files /dev/null and b/onerow/benchmark/instances/sp98ic.pre.mps.gz differ diff --git a/onerow/benchmark/instances/sp98ir.pre.mps.gz b/onerow/benchmark/instances/sp98ir.pre.mps.gz new file mode 100644 index 0000000..7d21f68 Binary files /dev/null and b/onerow/benchmark/instances/sp98ir.pre.mps.gz differ diff --git a/onerow/benchmark/instances/tanglegram1.pre.mps.gz b/onerow/benchmark/instances/tanglegram1.pre.mps.gz new file mode 100644 index 0000000..cbbed62 Binary files /dev/null and b/onerow/benchmark/instances/tanglegram1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/tanglegram2.pre.mps.gz b/onerow/benchmark/instances/tanglegram2.pre.mps.gz new file mode 100644 index 0000000..ac1498b Binary files /dev/null and b/onerow/benchmark/instances/tanglegram2.pre.mps.gz differ diff --git a/onerow/benchmark/instances/timtab1.pre.mps.gz b/onerow/benchmark/instances/timtab1.pre.mps.gz new file mode 100644 index 0000000..97eaad0 Binary files /dev/null and b/onerow/benchmark/instances/timtab1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/triptim1.pre.mps.gz b/onerow/benchmark/instances/triptim1.pre.mps.gz new file mode 100644 index 0000000..2faf737 Binary files /dev/null and b/onerow/benchmark/instances/triptim1.pre.mps.gz differ diff --git a/onerow/benchmark/instances/unitcal_7.pre.mps.gz b/onerow/benchmark/instances/unitcal_7.pre.mps.gz new file mode 100644 index 0000000..b41e973 Binary files /dev/null and b/onerow/benchmark/instances/unitcal_7.pre.mps.gz differ diff --git a/onerow/benchmark/instances/vpphard.pre.mps.gz b/onerow/benchmark/instances/vpphard.pre.mps.gz new file mode 100644 index 0000000..c520700 Binary files /dev/null and b/onerow/benchmark/instances/vpphard.pre.mps.gz differ diff --git a/onerow/benchmark/instances/zib54-UUE.pre.mps.gz b/onerow/benchmark/instances/zib54-UUE.pre.mps.gz new file mode 100644 index 0000000..f591869 Binary files /dev/null and b/onerow/benchmark/instances/zib54-UUE.pre.mps.gz differ diff --git a/onerow/benchmark/out/.keep b/onerow/benchmark/out/.keep new file mode 100644 index 0000000..e69de29 diff --git a/onerow/benchmark/out/30n20b8.pre.log b/onerow/benchmark/out/30n20b8.pre.log new file mode 100644 index 0000000..e2208a2 --- /dev/null +++ b/onerow/benchmark/out/30n20b8.pre.log @@ -0,0 +1,7 @@ +[ 0.00] Reading input file: instances/30n20b8.pre.mps.gz... +[ 0.00] Fetched 406 rows, 5077 cols. +[ 0.00] Solving first relaxation... +[ 0.36] 93.274056 [optimal] +[ 0.36] Reading basis... + [ 0.41] Generating MIR cuts... +[ 0.41] Finding interesting rows... diff --git a/onerow/benchmark/out/bienst2.log b/onerow/benchmark/out/bienst2.log new file mode 100644 index 0000000..95598f5 --- /dev/null +++ b/onerow/benchmark/out/bienst2.log @@ -0,0 +1,32 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/bienst2.pre.mps.gz... +[ 0.00] Fetched 520 rows, 898 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 11.724138 [optimal] +[ 0.00] Reading basis... + [ 0.01] Generating MIR cuts... +[ 0.01] Finding interesting rows... +[ 3.22] 33 rows found +[ 3.22] Starting timer 1... +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/bienst2.pre.mps.gz... +[ 0.00] Fetched 520 rows, 898 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 11.724138 [optimal] +[ 0.00] Reading basis... + [ 0.01] Generating MIR cuts... +[ 0.01] Finding interesting rows... +[ 3.00] 33 rows found +[ 3.00] Starting timer 1... +[ 3.20] Ending timer 1: 0.20s + [ 3.20] Added 33 violated cuts... +[ 3.22] 14.442185 [optimal] +[ 3.22] Generating wedge cuts... +[ 3.22] Starting timer 2... +[ 3.59] Ending timer 2: 0.37s + [ 3.59] Added 0 violated cuts... +[ 3.59] 14.442185 [optimal] +[ 3.59] Writting stats: out/bienst2.yaml... +[ 3.60] Done. diff --git a/onerow/benchmark/out/bienst2.yaml b/onerow/benchmark/out/bienst2.yaml new file mode 100644 index 0000000..4c82b78 --- /dev/null +++ b/onerow/benchmark/out/bienst2.yaml @@ -0,0 +1,36 @@ +input_file: + instances/bienst2.pre.mps.gz +sol_value: + 0: 11.724138 + 1: 14.442185 + 2: 14.442185 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 33 + depth: + 0: 33 +n_generated_cuts: + total: + 41 + round: + 1: 33 + 2: 8 + depth: + 0: 41 +trivial_lifting: + max_m: 0 + average_m: 0.000000 + integral_coefficients: 0.002685 + slowdown: 0.000000 +timers: + 1: 0.2040 + 2: 0.3680 +cut_speed: + round_1: 0.0062 + round_2: 0.0460 +mip_value: + 54.600000 diff --git a/onerow/benchmark/out/csched010.log b/onerow/benchmark/out/csched010.log new file mode 100644 index 0000000..1d7f024 --- /dev/null +++ b/onerow/benchmark/out/csched010.log @@ -0,0 +1,17 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/csched010.pre.mps.gz... +[ 0.00] Fetched 271 rows, 1679 cols. +[ 0.00] Solving first relaxation... +[ 0.07] 332.422727 [optimal] +[ 0.07] Reading basis... + [ 0.08] Generating MIR cuts... +[ 0.08] Finding interesting rows... +[ 3.42] 84 rows found +[ 3.42] Starting timer 1... +[ 5.77] Ending timer 1: 2.35s + [ 5.79] Added 84 violated cuts... +[ 5.81] 335.363758 [optimal] +[ 5.81] Generating wedge cuts... +[ 5.81] Starting timer 2... + 2% ETA: 2017-04-05 11:24:18 2 / 84 2% ETA: 2017-04-05 11:25:00 2 / 84 2% ETA: 2017-04-05 11:25:42 2 / 84 2% ETA: 2017-04-05 11:26:24 2 / 84 2% ETA: 2017-04-05 11:27:06 2 / 84 2% ETA: 2017-04-05 11:27:48 2 / 84 2% ETA: 2017-04-05 11:28:30 2 / 84 2% ETA: 2017-04-05 11:29:12 2 / 84 2% ETA: 2017-04-05 11:29:54 2 / 84 2% ETA: 2017-04-05 11:30:36 2 / 84 2% ETA: 2017-04-05 11:31:18 2 / 84 2% ETA: 2017-04-05 11:32:00 2 / 84 5% ETA: 2017-04-05 11:28:09 4 / 84 5% ETA: 2017-04-05 11:28:30 4 / 84 5% ETA: 2017-04-05 11:28:51 4 / 84 6% ETA: 2017-04-05 11:28:04 5 / 84 7% ETA: 2017-04-05 11:27:34 6 / 84 7% ETA: 2017-04-05 11:27:48 6 / 84 7% ETA: 2017-04-05 11:28:02 6 / 84 7% ETA: 2017-04-05 11:28:16 6 / 84 7% ETA: 2017-04-05 11:28:30 6 / 84 7% ETA: 2017-04-05 11:28:44 6 / 84 7% ETA: 2017-04-05 11:28:58 6 / 84 7% ETA: 2017-04-05 11:29:12 6 / 84 7% ETA: 2017-04-05 11:29:26 6 / 84 7% ETA: 2017-04-05 11:29:40 6 / 84 7% ETA: 2017-04-05 11:29:54 6 / 84 8% ETA: 2017-04-05 11:29:12 7 / 84 8% ETA: 2017-04-05 11:29:24 7 / 84 8% ETA: 2017-04-05 11:29:36 7 / 84 8% ETA: 2017-04-05 11:29:48 7 / 84 8% ETA: 2017-04-05 11:30:00 7 / 84 8% ETA: 2017-04-05 11:30:12 7 / 84 10% ETA: 2017-04-05 11:29:33 8 / 84 11% ETA: 2017-04-05 11:29:02 9 / 84 11% ETA: 2017-04-05 11:29:12 9 / 84 12% ETA: 2017-04-05 11:28:46 10 / 84 13% ETA: 2017-04-05 11:28:26 11 / 84 13% ETA: 2017-04-05 11:28:33 11 / 84 13% ETA: 2017-04-05 11:28:41 11 / 84 13% ETA: 2017-04-05 11:28:49 11 / 84 13% ETA: 2017-04-05 11:28:56 11 / 84 13% ETA: 2017-04-05 11:29:04 11 / 84 14% ETA: 2017-04-05 11:28:44 12 / 84 15% ETA: 2017-04-05 11:28:26 13 / 84 15% ETA: 2017-04-05 11:28:33 13 / 84 15% ETA: 2017-04-05 11:28:39 13 / 84 15% ETA: 2017-04-05 11:28:46 13 / 84 15% ETA: 2017-04-05 11:28:52 13 / 84 17% ETA: 2017-04-05 11:28:36 14 / 84 17% ETA: 2017-04-05 11:28:42 14 / 84 17% ETA: 2017-04-05 11:28:48 14 / 84 18% ETA: 2017-04-05 11:28:32 15 / 84 19% ETA: 2017-04-05 11:28:19 16 / 84 19% ETA: 2017-04-05 11:28:24 16 / 84 19% ETA: 2017-04-05 11:28:30 16 / 84 19% ETA: 2017-04-05 11:28:35 16 / 84 20% ETA: 2017-04-05 11:28:22 17 / 84 \ No newline at end of file diff --git a/onerow/benchmark/out/danoint.log b/onerow/benchmark/out/danoint.log new file mode 100644 index 0000000..43bf66a --- /dev/null +++ b/onerow/benchmark/out/danoint.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/danoint.pre.mps.gz... +[ 0.00] Fetched 600 rows, 978 cols. +[ 0.00] Solving first relaxation... +[ 0.03] 62.637280 [optimal] +[ 0.03] Reading basis... + [ 0.04] Generating MIR cuts... +[ 0.04] Finding interesting rows... +[ 10.11] 52 rows found +[ 10.11] Starting timer 1... +[ 10.85] Ending timer 1: 0.74s + [ 10.86] Added 52 violated cuts... +[ 10.90] 62.690000 [optimal] +[ 10.90] Generating wedge cuts... +[ 10.90] Starting timer 2... +[ 13.15] Ending timer 2: 2.26s + [ 13.16] Added 0 violated cuts... +[ 13.16] 62.690000 [optimal] +[ 13.16] Writting stats: out/danoint.yaml... +[ 13.16] Done. diff --git a/onerow/benchmark/out/danoint.yaml b/onerow/benchmark/out/danoint.yaml new file mode 100644 index 0000000..01b1c0d --- /dev/null +++ b/onerow/benchmark/out/danoint.yaml @@ -0,0 +1,44 @@ +input_file: + instances/danoint.pre.mps.gz +sol_value: + 0: 62.637280 + 1: 62.690000 + 2: 62.690000 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 52 + depth: + 0: 52 +n_generated_cuts: + total: + 593 + round: + 1: 52 + 2: 541 + depth: + 0: 145 + 1: 93 + 2: 93 + 3: 84 + 4: 75 + 5: 58 + 6: 33 + 7: 10 + 8: 2 +trivial_lifting: + max_m: 1035 + average_m: 61.002773 + integral_coefficients: 0.005403 + slowdown: 0.329606 +timers: + 1: 0.7440 + 2: 2.2560 +cut_speed: + round_1: 0.0143 + round_2: 0.0042 +mip_value: + 65.666700 diff --git a/onerow/benchmark/out/dfn-gwin-UUM.log b/onerow/benchmark/out/dfn-gwin-UUM.log new file mode 100644 index 0000000..751a2e3 --- /dev/null +++ b/onerow/benchmark/out/dfn-gwin-UUM.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/dfn-gwin-UUM.pre.mps.gz... +[ 0.01] Fetched 156 rows, 984 cols. +[ 0.01] Solving first relaxation... +[ 0.02] 27467.257235 [optimal] +[ 0.02] Reading basis... + [ 0.02] Generating MIR cuts... +[ 0.02] Finding interesting rows... +[ 0.17] 45 rows found +[ 0.17] Starting timer 1... +[ 0.27] Ending timer 1: 0.10s + [ 0.28] Added 45 violated cuts... +[ 0.29] 32186.983731 [optimal] +[ 0.29] Generating wedge cuts... +[ 0.29] Starting timer 2... +[ 0.61] Ending timer 2: 0.32s + [ 0.62] Added 1 violated cuts... +[ 0.62] 32195.320369 [optimal] +[ 0.62] Writting stats: out/dfn-gwin-UUM.yaml... +[ 0.62] Done. diff --git a/onerow/benchmark/out/dfn-gwin-UUM.yaml b/onerow/benchmark/out/dfn-gwin-UUM.yaml new file mode 100644 index 0000000..4c51686 --- /dev/null +++ b/onerow/benchmark/out/dfn-gwin-UUM.yaml @@ -0,0 +1,37 @@ +input_file: + instances/dfn-gwin-UUM.pre.mps.gz +sol_value: + 0: 27467.257235 + 1: 32186.983731 + 2: 32195.320369 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 46 + depth: + 0: 45 + 1: 1 +n_generated_cuts: + total: + 201 + round: + 1: 45 + 2: 156 + depth: + 0: 90 + 1: 44 + 2: 30 + 3: 27 + 4: 9 + 5: 1 +timers: + 1: 0.1000 + 2: 0.3240 +cut_speed: + round_1: 0.0022 + round_2: 0.0021 +mip_value: + 38752.000000 diff --git a/onerow/benchmark/out/enlight13.log b/onerow/benchmark/out/enlight13.log new file mode 100644 index 0000000..0ce9ad0 --- /dev/null +++ b/onerow/benchmark/out/enlight13.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/enlight13.pre.mps.gz... +[ 0.00] Fetched 169 rows, 339 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 1.000000 [optimal] +[ 0.00] Reading basis... + [ 0.00] Generating MIR cuts... +[ 0.00] Finding interesting rows... +[ 0.03] 2 rows found +[ 0.04] Starting timer 1... +[ 0.04] Ending timer 1: 0.00s + [ 0.05] Added 2 violated cuts... +[ 0.05] 1.500000 [optimal] +[ 0.05] Generating wedge cuts... +[ 0.05] Starting timer 2... +[ 0.06] Ending timer 2: 0.01s + [ 0.08] Added 0 violated cuts... +[ 0.08] 1.500000 [optimal] +[ 0.08] Writting stats: out/enlight13.yaml... +[ 0.08] Done. diff --git a/onerow/benchmark/out/enlight13.yaml b/onerow/benchmark/out/enlight13.yaml new file mode 100644 index 0000000..bc3100a --- /dev/null +++ b/onerow/benchmark/out/enlight13.yaml @@ -0,0 +1,36 @@ +input_file: + instances/enlight13.pre.mps.gz +sol_value: + 0: 1.000000 + 1: 1.500000 + 2: 1.500000 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 2 + depth: + 0: 2 +n_generated_cuts: + total: + 10 + round: + 1: 2 + 2: 8 + depth: + 0: 10 +trivial_lifting: + max_m: 1 + average_m: 0.750000 + integral_coefficients: 0.750000 + slowdown: 0.562500 +timers: + 1: 0.0000 + 2: 0.0080 +cut_speed: + round_1: 0.0000 + round_2: 0.0010 +mip_value: + 71.000000 diff --git a/onerow/benchmark/out/glass4.log b/onerow/benchmark/out/glass4.log new file mode 100644 index 0000000..ca2b186 --- /dev/null +++ b/onerow/benchmark/out/glass4.log @@ -0,0 +1,28 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/glass4.pre.mps.gz... +[ 0.00] Fetched 392 rows, 674 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 800002400.000000 [optimal] +[ 0.00] Reading basis... +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/glass4.pre.mps.gz... +[ 0.00] Fetched 392 rows, 674 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 800002400.000000 [optimal] +[ 0.00] Reading basis... + [ 0.01] Generating MIR cuts... +[ 0.01] Finding interesting rows... +[ 0.08] 72 rows found +[ 0.08] Starting timer 1... +[ 0.10] Ending timer 1: 0.02s + [ 0.12] Added 72 violated cuts... +[ 0.12] 1269062192.116870 [optimal] +[ 0.12] Generating wedge cuts... +[ 0.12] Starting timer 2... +[ 0.19] Ending timer 2: 0.07s + [ 0.19] Added 0 violated cuts... +[ 0.20] 1269062192.116883 [optimal] +[ 0.20] Writting stats: out/glass4.yaml... +[ 0.20] Done. diff --git a/onerow/benchmark/out/glass4.yaml b/onerow/benchmark/out/glass4.yaml new file mode 100644 index 0000000..f292d85 --- /dev/null +++ b/onerow/benchmark/out/glass4.yaml @@ -0,0 +1,36 @@ +input_file: + instances/glass4.pre.mps.gz +sol_value: + 0: 800002400.000000 + 1: 1269062192.116870 + 2: 1269062192.116883 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 72 + depth: + 0: 72 +n_generated_cuts: + total: + 298 + round: + 1: 72 + 2: 226 + depth: + 0: 298 +trivial_lifting: + max_m: 0 + average_m: 0.000000 + integral_coefficients: 0.518489 + slowdown: 0.000000 +timers: + 1: 0.0240 + 2: 0.0720 +cut_speed: + round_1: 0.0003 + round_2: 0.0003 +mip_value: + 1200010000.000000 diff --git a/onerow/benchmark/out/mik-250-1-100-1.log b/onerow/benchmark/out/mik-250-1-100-1.log new file mode 100644 index 0000000..a7449e8 --- /dev/null +++ b/onerow/benchmark/out/mik-250-1-100-1.log @@ -0,0 +1,31 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/mik-250-1-100-1.pre.mps.gz... +[ 0.00] Fetched 100 rows, 352 cols. +[ 0.00] Solving first relaxation... +[ 0.00] -79842.423635 [optimal] +[ 0.00] Reading basis... + [ 0.00] Generating MIR cuts... +[ 0.00] Finding interesting rows... +[ 0.21] 100 rows found +[ 0.22] Starting timer 1... +[ 0.42] -72824.000000 [optimal] +[ 0.42] Ending timer 1: 0.20s + [ 0.43] Added 100 violated cuts... +[ 0.43] -72824.000000 [optimal] +[ 0.43] Generating wedge cuts... +[ 0.43] Starting timer 2... +[ 0.97] -72659.578424 [optimal] +[ 1.34] -72398.775758 [optimal] +[ 2.65] -71934.177977 [optimal] +[ 3.47] -71625.914521 [optimal] +[ 4.02] -71357.200706 [optimal] + 26% ETA: 2017-04-05 11:23:32 26 / 100 [ 5.10] -71103.868040 [optimal] +[ 7.52] -70668.685770 [optimal] + 53% ETA: 2017-04-05 11:23:32 53 / 100 [ 9.99] -70461.015542 [optimal] + 83% ETA: 2017-04-05 11:23:32 83 / 100 [ 11.95] -70278.553377 [optimal] +[ 13.81] Ending timer 2: 13.38s + [ 13.81] Added 955 violated cuts... +[ 13.82] -70219.174256 [optimal] +[ 13.82] Writting stats: out/mik-250-1-100-1.yaml... +[ 13.82] Done. diff --git a/onerow/benchmark/out/mik-250-1-100-1.yaml b/onerow/benchmark/out/mik-250-1-100-1.yaml new file mode 100644 index 0000000..59486eb --- /dev/null +++ b/onerow/benchmark/out/mik-250-1-100-1.yaml @@ -0,0 +1,55 @@ +input_file: + instances/mik-250-1-100-1.pre.mps.gz +sol_value: + 0: -79842.423635 + 1: -72824.000000 + 2: -70219.174256 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 1055 + depth: + 0: 100 + 1: 110 + 2: 306 + 3: 294 + 4: 157 + 5: 65 + 6: 20 + 7: 2 + 8: 1 +n_generated_cuts: + total: + 30321 + round: + 1: 100 + 2: 30221 + depth: + 0: 5197 + 1: 5082 + 2: 5042 + 3: 4686 + 4: 4146 + 5: 3132 + 6: 1897 + 7: 826 + 8: 256 + 9: 48 + 10: 8 + 11: 1 +trivial_lifting: + max_m: 1405 + average_m: 45.731028 + integral_coefficients: 0.944252 + slowdown: 43.181598 +timers: + 1: 0.2040 + 2: 13.3800 +cut_speed: + round_1: 0.0020 + round_2: 0.0004 +mip_value: + -66729.000000 diff --git a/onerow/benchmark/out/newdano.log b/onerow/benchmark/out/newdano.log new file mode 100644 index 0000000..da82d39 --- /dev/null +++ b/onerow/benchmark/out/newdano.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/newdano.pre.mps.gz... +[ 0.00] Fetched 520 rows, 898 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 11.724138 [optimal] +[ 0.00] Reading basis... + [ 0.02] Generating MIR cuts... +[ 0.02] Finding interesting rows... +[ 3.05] 52 rows found +[ 3.05] Starting timer 1... +[ 3.33] Ending timer 1: 0.28s + [ 3.34] Added 52 violated cuts... +[ 3.36] 14.442185 [optimal] +[ 3.36] Generating wedge cuts... +[ 3.36] Starting timer 2... +[ 3.79] Ending timer 2: 0.44s + [ 3.80] Added 0 violated cuts... +[ 3.80] 14.442185 [optimal] +[ 3.80] Writting stats: out/newdano.yaml... +[ 3.80] Done. diff --git a/onerow/benchmark/out/newdano.yaml b/onerow/benchmark/out/newdano.yaml new file mode 100644 index 0000000..df192f5 --- /dev/null +++ b/onerow/benchmark/out/newdano.yaml @@ -0,0 +1,36 @@ +input_file: + instances/newdano.pre.mps.gz +sol_value: + 0: 11.724138 + 1: 14.442185 + 2: 14.442185 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 52 + depth: + 0: 52 +n_generated_cuts: + total: + 72 + round: + 1: 52 + 2: 20 + depth: + 0: 72 +trivial_lifting: + max_m: 0 + average_m: 0.000000 + integral_coefficients: 0.002313 + slowdown: 0.000000 +timers: + 1: 0.2800 + 2: 0.4360 +cut_speed: + round_1: 0.0054 + round_2: 0.0218 +mip_value: + 65.666700 diff --git a/onerow/benchmark/out/noswot.log b/onerow/benchmark/out/noswot.log new file mode 100644 index 0000000..2d132e5 --- /dev/null +++ b/onerow/benchmark/out/noswot.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/noswot.pre.mps.gz... +[ 0.00] Fetched 172 rows, 293 cols. +[ 0.00] Solving first relaxation... +[ 0.00] -43.000000 [optimal] +[ 0.00] Reading basis... + [ 0.00] Generating MIR cuts... +[ 0.00] Finding interesting rows... +[ 0.11] 49 rows found +[ 0.11] Starting timer 1... +[ 0.14] Ending timer 1: 0.04s + [ 0.15] Added 49 violated cuts... +[ 0.15] -43.000000 [optimal] +[ 0.15] Generating wedge cuts... +[ 0.15] Starting timer 2... +[ 0.70] Ending timer 2: 0.54s + [ 0.70] Added 5 violated cuts... +[ 0.70] -43.000000 [optimal] +[ 0.70] Writting stats: out/noswot.yaml... +[ 0.70] Done. diff --git a/onerow/benchmark/out/noswot.yaml b/onerow/benchmark/out/noswot.yaml new file mode 100644 index 0000000..fa63cfe --- /dev/null +++ b/onerow/benchmark/out/noswot.yaml @@ -0,0 +1,46 @@ +input_file: + instances/noswot.pre.mps.gz +sol_value: + 0: -43.000000 + 1: -43.000000 + 2: -43.000000 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 54 + depth: + 0: 49 + 1: 5 +n_generated_cuts: + total: + 1089 + round: + 1: 49 + 2: 1040 + depth: + 0: 294 + 1: 207 + 2: 200 + 3: 142 + 4: 103 + 5: 72 + 6: 42 + 7: 19 + 8: 8 + 9: 2 +trivial_lifting: + max_m: 755 + average_m: 31.524828 + integral_coefficients: 0.571790 + slowdown: 18.025587 +timers: + 1: 0.0360 + 2: 0.5440 +cut_speed: + round_1: 0.0007 + round_2: 0.0005 +mip_value: + -41.000000 diff --git a/onerow/benchmark/out/pigeon-10.log b/onerow/benchmark/out/pigeon-10.log new file mode 100644 index 0000000..0ef23b9 --- /dev/null +++ b/onerow/benchmark/out/pigeon-10.log @@ -0,0 +1,22 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/pigeon-10.pre.mps.gz... +[ 0.00] Fetched 525 rows, 866 cols. +[ 0.00] Solving first relaxation... +[ 0.01] -10000.000000 [optimal] +[ 0.01] Reading basis... + [ 0.02] Generating MIR cuts... +[ 0.02] Finding interesting rows... +[ 0.12] 126 rows found +[ 0.12] Starting timer 1... +[ 0.17] -10000.000000 [optimal] +[ 0.18] Ending timer 1: 0.06s + [ 0.18] Added 112 violated cuts... +[ 0.19] -10000.000000 [optimal] +[ 0.19] Generating wedge cuts... +[ 0.19] Starting timer 2... +[ 0.90] Ending timer 2: 0.71s + [ 0.90] Added 0 violated cuts... +[ 0.90] -10000.000000 [optimal] +[ 0.90] Writting stats: out/pigeon-10.yaml... +[ 0.90] Done. diff --git a/onerow/benchmark/out/pigeon-10.yaml b/onerow/benchmark/out/pigeon-10.yaml new file mode 100644 index 0000000..f6f4cff --- /dev/null +++ b/onerow/benchmark/out/pigeon-10.yaml @@ -0,0 +1,42 @@ +input_file: + instances/pigeon-10.pre.mps.gz +sol_value: + 0: -10000.000000 + 1: -10000.000000 + 2: -10000.000000 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 112 + depth: + 0: 112 +n_generated_cuts: + total: + 2237 + round: + 1: 126 + 2: 2111 + depth: + 0: 1189 + 1: 495 + 2: 324 + 3: 135 + 4: 79 + 5: 11 + 6: 4 +trivial_lifting: + max_m: 910 + average_m: 8.150823 + integral_coefficients: 0.534935 + slowdown: 4.360160 +timers: + 1: 0.0600 + 2: 0.7080 +cut_speed: + round_1: 0.0005 + round_2: 0.0003 +mip_value: + -9000.000000 diff --git a/onerow/benchmark/out/qiu.log b/onerow/benchmark/out/qiu.log new file mode 100644 index 0000000..e1351a6 --- /dev/null +++ b/onerow/benchmark/out/qiu.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/qiu.pre.mps.gz... +[ 0.00] Fetched 1192 rows, 1901 cols. +[ 0.00] Solving first relaxation... +[ 0.03] -931.638854 [optimal] +[ 0.03] Reading basis... + [ 0.07] Generating MIR cuts... +[ 0.07] Finding interesting rows... +[ 7.61] 36 rows found +[ 7.61] Starting timer 1... +[ 8.18] Ending timer 1: 0.58s + [ 8.19] Added 36 violated cuts... +[ 8.26] -919.644028 [optimal] +[ 8.26] Generating wedge cuts... +[ 8.26] Starting timer 2... +[ 9.27] Ending timer 2: 1.02s + [ 9.28] Added 0 violated cuts... +[ 9.28] -919.644028 [optimal] +[ 9.28] Writting stats: out/qiu.yaml... +[ 9.28] Done. diff --git a/onerow/benchmark/out/qiu.yaml b/onerow/benchmark/out/qiu.yaml new file mode 100644 index 0000000..84101eb --- /dev/null +++ b/onerow/benchmark/out/qiu.yaml @@ -0,0 +1,30 @@ +input_file: + instances/qiu.pre.mps.gz +sol_value: + 0: -931.638854 + 1: -919.644028 + 2: -919.644028 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 36 + depth: + 0: 36 +n_generated_cuts: + total: + 36 + round: + 1: 36 + depth: + 0: 36 +timers: + 1: 0.5760 + 2: 1.0160 +cut_speed: + round_1: 0.0160 + round_2: inf +mip_value: + -132.873000 diff --git a/onerow/benchmark/out/ran16x16.log b/onerow/benchmark/out/ran16x16.log new file mode 100644 index 0000000..0706114 --- /dev/null +++ b/onerow/benchmark/out/ran16x16.log @@ -0,0 +1,21 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/ran16x16.pre.mps.gz... +[ 0.00] Fetched 288 rows, 769 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 3116.429512 [optimal] +[ 0.00] Reading basis... + [ 0.01] Generating MIR cuts... +[ 0.01] Finding interesting rows... +[ 0.12] 20 rows found +[ 0.12] Starting timer 1... +[ 0.19] Ending timer 1: 0.07s + [ 0.20] Added 20 violated cuts... +[ 0.21] 3238.332602 [optimal] +[ 0.21] Generating wedge cuts... +[ 0.21] Starting timer 2... +[ 1.51] Ending timer 2: 1.30s + [ 1.51] Added 52 violated cuts... +[ 1.52] 3244.082675 [optimal] +[ 1.52] Writting stats: out/ran16x16.yaml... +[ 1.52] Done. diff --git a/onerow/benchmark/out/ran16x16.yaml b/onerow/benchmark/out/ran16x16.yaml new file mode 100644 index 0000000..e6a582a --- /dev/null +++ b/onerow/benchmark/out/ran16x16.yaml @@ -0,0 +1,43 @@ +input_file: + instances/ran16x16.pre.mps.gz +sol_value: + 0: 3116.429512 + 1: 3238.332602 + 2: 3244.082675 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 72 + depth: + 0: 20 + 1: 41 + 2: 10 + 3: 1 +n_generated_cuts: + total: + 1716 + round: + 1: 20 + 2: 1696 + depth: + 0: 903 + 1: 471 + 2: 285 + 3: 54 + 4: 3 +trivial_lifting: + max_m: 45 + average_m: 3.522548 + integral_coefficients: 0.440125 + slowdown: 1.550363 +timers: + 1: 0.0680 + 2: 1.3000 +cut_speed: + round_1: 0.0034 + round_2: 0.0008 +mip_value: + 3823.000000 diff --git a/onerow/benchmark/out/timtab1.log b/onerow/benchmark/out/timtab1.log new file mode 100644 index 0000000..79491c7 --- /dev/null +++ b/onerow/benchmark/out/timtab1.log @@ -0,0 +1,22 @@ +version: benchmark.run +machine: mensa +[ 0.00] Reading input file: instances/timtab1.pre.mps.gz... +[ 0.00] Fetched 166 rows, 379 cols. +[ 0.00] Solving first relaxation... +[ 0.00] 28694.000000 [optimal] +[ 0.00] Reading basis... + [ 0.00] Generating MIR cuts... +[ 0.00] Finding interesting rows... +[ 0.06] 135 rows found +[ 0.06] Starting timer 1... +[ 0.09] 112659.532725 [optimal] +[ 0.10] Ending timer 1: 0.05s + [ 0.11] Added 121 violated cuts... +[ 0.11] 205801.284255 [optimal] +[ 0.11] Generating wedge cuts... +[ 0.11] Starting timer 2... +[ 0.19] Ending timer 2: 0.08s + [ 0.20] Added 0 violated cuts... +[ 0.20] 205801.284255 [optimal] +[ 0.20] Writting stats: out/timtab1.yaml... +[ 0.20] Done. diff --git a/onerow/benchmark/out/timtab1.yaml b/onerow/benchmark/out/timtab1.yaml new file mode 100644 index 0000000..579391d --- /dev/null +++ b/onerow/benchmark/out/timtab1.yaml @@ -0,0 +1,36 @@ +input_file: + instances/timtab1.pre.mps.gz +sol_value: + 0: 28694.000000 + 1: 205801.284255 + 2: 205801.284255 +sol_status: + 0: optimal + 1: optimal + 2: optimal +n_added_cuts: + total: + 121 + depth: + 0: 121 +n_generated_cuts: + total: + 250 + round: + 1: 135 + 2: 115 + depth: + 0: 250 +trivial_lifting: + max_m: 0 + average_m: 0.000000 + integral_coefficients: 0.088368 + slowdown: 0.000000 +timers: + 1: 0.0480 + 2: 0.0800 +cut_speed: + round_1: 0.0004 + round_2: 0.0007 +mip_value: + 764772.000000 diff --git a/onerow/benchmark/scripts/gap.rb b/onerow/benchmark/scripts/gap.rb new file mode 100755 index 0000000..4d86882 --- /dev/null +++ b/onerow/benchmark/scripts/gap.rb @@ -0,0 +1,45 @@ +#!/usr/bin/env ruby +require 'yaml' +yaml = YAML::load(File.open(ARGV[0])) + +int_format = "%d," +number_format = "%.2f," +string_format = "%s," + +presolve_opt = yaml['sol_value'][0] +mir_opt = yaml['sol_value'][1] +wedges_opt = yaml['sol_value'][2] +mip_opt = yaml['mip_value'] + +if presolve_opt == mip_opt then + print (string_format * 6 + "\n") % (["---"] * 6) +else + + orig_gap = (mip_opt-presolve_opt) / mip_opt.abs + mip_perf = (mir_opt-presolve_opt) / (mip_opt-presolve_opt) + w_perf = (wedges_opt-presolve_opt) / (mip_opt-presolve_opt) + + n_cuts_mir = yaml['n_added_cuts']['depth'][0] + n_cuts_w = yaml['n_added_cuts']['depth'].values.reduce(:+) - n_cuts_mir + + print int_format % n_cuts_mir + print int_format % n_cuts_w + print number_format % (100 * orig_gap) + print number_format % (100 * mip_perf) + print number_format % (100 * w_perf) + + if presolve_opt == wedges_opt + print (string_format * 2) % (["---"] * 2) + else + mir_contrib = (mir_opt-presolve_opt) / (wedges_opt-presolve_opt) + w_contrib = (wedges_opt-mir_opt) / (wedges_opt-presolve_opt) + + print number_format % (100 * mir_contrib) + print number_format % (100 * w_contrib) + end + + w_improv = (wedges_opt-mir_opt) / (mip_opt-mir_opt) + + # improvement + print (number_format + "\n") % (100 * w_improv.abs) +end diff --git a/onerow/benchmark/scripts/speed.rb b/onerow/benchmark/scripts/speed.rb new file mode 100755 index 0000000..a6ab3e3 --- /dev/null +++ b/onerow/benchmark/scripts/speed.rb @@ -0,0 +1,15 @@ +#!/usr/bin/env ruby +require 'yaml' +yaml = YAML::load(File.open(ARGV[0])) + +int_format = "%d," +number_format = "%.2f," +string_format = "%s," +blank = "---," + +print int_format % yaml['n_generated_cuts']['round'][1] rescue print blank +print int_format % yaml['n_generated_cuts']['round'][2] rescue print blank +print number_format % (yaml['cut_speed']['round_1'] * 1000) rescue print blank +print number_format % (yaml['cut_speed']['round_2'] * 1000) rescue print blank +print number_format % yaml['trivial_lifting']['average_m'] rescue print blank +print "\n" diff --git a/onerow/benchmark/src/main.cpp b/onerow/benchmark/src/main.cpp new file mode 100644 index 0000000..20a5d69 --- /dev/null +++ b/onerow/benchmark/src/main.cpp @@ -0,0 +1,169 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include +#include + +using namespace std; + +char *input_filename = 0; +char *stats_filename = 0; +char *sol_filename = 0; +bool enable_gomory_cuts = false; +bool enable_wedge_cuts = false; +bool enable_mir_cuts = false; + +int c; +extern char *optarg; +char usage[] = "usage: %s [-gmw] -f model.mps [-s stats.yaml]\n"; + +void read_params(int argc, char **argv) +{ + while ((c = getopt(argc, argv, "gwmf:s:x:")) != -1) + { + switch (c) + { + case 'g': + enable_gomory_cuts = true; + break; + case 'w': + enable_wedge_cuts = true; + break; + case 'm': + enable_mir_cuts = true; + break; + case 'f': + input_filename = optarg; + break; + case 's': + stats_filename = optarg; + break; + case 'x': + sol_filename = optarg; + break; + } + } + + if (!input_filename) + { + fprintf(stderr, "%s: missing model filename\n", argv[0]); + fprintf(stderr, usage, argv[0]); + exit(1); + } +} + +int main(int argc, char **argv) +{ + read_params(argc, argv); + + Stats::init(); + + int status; + CPXENVptr env = CPXopenCPLEX(&status); + CPXLPptr lp = CPXcreateprob(env, &status, ""); + CplexHelper cplexHelper(env, lp); + + CPXsetlogfile(env, NULL); + CPXsetintparam(env, CPX_PARAM_PREIND, CPX_OFF); // disable presolve + CPXsetintparam(env, CPX_PARAM_DATACHECK, CPX_ON); // check consistency + CPXsetintparam(env, CPX_PARAM_NUMERICALEMPHASIS, CPX_ON); // numerical precision + + Stats::set_input_filename(string(input_filename)); + + // reads input file + time_printf("Reading input file: %s...\n", input_filename); + status = CPXreadcopyprob(env, lp, input_filename, NULL); + if (status) + { + fprintf(stderr, "could not read input file (%d)\n", status); + return 1; + } + + cplexHelper.read_columns(); + + // read solution + if(sol_filename) + { + time_printf("Reading solution file: %s...\n", sol_filename); + FILE *f = fopen(sol_filename, "r"); + if(!f) + { + fprintf(stderr, "Could not open solution file (%s).", sol_filename); + return 1; + } + + double *solution = new double[cplexHelper.n_cols]; + for(int i=0; i(); + cplexHelper.solve(true); + } + + if (enable_mir_cuts) + { + time_printf("Generating MIR cuts...\n"); + cplexHelper.add_single_row_cuts(); + cplexHelper.solve(true); + } + + if (enable_wedge_cuts) + { + time_printf("Generating wedge cuts...\n"); + cplexHelper.add_single_row_cuts(); + cplexHelper.solve(true); + } + + if(stats_filename != 0) + { + time_printf("Writting stats: %s...\n", stats_filename); + Stats::write_stats(string(stats_filename)); + } + + time_printf("Done.\n"); + CPXcloseCPLEX(&env); + + return 0; + +} diff --git a/onerow/library/CMakeLists.txt b/onerow/library/CMakeLists.txt new file mode 100644 index 0000000..8ae7974 --- /dev/null +++ b/onerow/library/CMakeLists.txt @@ -0,0 +1,35 @@ +include_directories(include) + +set(COMMON_SOURCES + src/cplex_helper.cpp + src/geometry.cpp + src/gomory_cut_generator.cpp + src/knapsack2.cpp + src/mir_cut_generator.cpp + src/single_row_cut_generator.cpp + src/stats.cpp + src/wedge_cut_generator.cpp + include/onerow/cplex_helper.hpp + include/onerow/cplex_helper.tpp + include/onerow/geometry.hpp + include/onerow/gomory_cut_generator.hpp + include/onerow/knapsack2.hpp + include/onerow/mir_cut_generator.hpp + include/onerow/params.hpp + include/onerow/single_row_cut_generator.hpp + include/onerow/stats.hpp + include/onerow/wedge_cut_generator.hpp) + +set(TEST_SOURCES + tests/knapsack2_test.cpp + tests/rational_test.cpp + tests/single_row_generator_test.cpp + tests/wedge_cut_generator_test.cpp) + +add_library(onerow_static ${COMMON_SOURCES}) +set_target_properties(onerow_static PROPERTIES OUTPUT_NAME onerow) +target_include_directories (onerow_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) +target_link_libraries(onerow_static qxx_static ${GMP_LIBRARIES} ${CPLEX_LIBRARIES}) + +add_executable(onerow-test.run ${COMMON_SOURCES} ${TEST_SOURCES}) +target_link_libraries(onerow-test.run gtest_main qxx_static ${GMP_LIBRARIES} ${CPLEX_LIBRARIES}) diff --git a/onerow/library/include/onerow/cplex_helper.hpp b/onerow/library/include/onerow/cplex_helper.hpp new file mode 100644 index 0000000..b9a8c15 --- /dev/null +++ b/onerow/library/include/onerow/cplex_helper.hpp @@ -0,0 +1,141 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef CPLEX_HELPER_HPP_ +#define CPLEX_HELPER_HPP_ + +#include +#include +#include "single_row_cut_generator.hpp" +using std::set; +using std::vector; + +struct CplexRow { + int nz; + double pi_zero; + double *pi; + int *indices; + int depth; + int head; + double dynamism; + + bool operator<(const CplexRow &other) const; + double get_violation(double *x); + void print(double *x); +}; + +/** + * This class provides useful methods for dealing with CPLEX. + */ +class CplexHelper { +private: + const CPXENVptr env; + const CPXLPptr lp; + bool *is_integer; + int n_cuts; + +public: + /** + * Constructs a helper associated with the provided environment and problem. + * + * @param env Pointer to the ILOG CPLEX environment. + * @param lp Pointer to a CPLEX LP problem object. + */ + CplexHelper(CPXENVptr env, CPXLPptr lp); + + ~CplexHelper(); + + /** + * Adds the specified constraints to the model. + * + * @param cuts Set of constraints to add. + */ + void add_cut(Constraint *cut); + + /** + * For each fractional row of the current tableau, adds as many single row + * cuts as possible. The cuts are generated by the provided generator class. + * + * @tparam Generator Class used to generate the cuts. + * @returns The number of cuts added. + */ + template + int add_single_row_cuts(); + + /** + * Gets a single row from the current tableau. + * + * @param index Index of the row to fetch. + * @returns The selected tableau row. + */ + Row* get_tableau_row(int index); + + void read_basis(); + void read_columns(); + + /** + * Solves the current problem, and prints some solution information. + * + * @returns The solution status, as returned by CPXgetstat. + */ + int solve(bool save_stats = false); + + void dump_constraint(const Constraint &c, const char *msg = ""); + + void print_basis(); + + void flush_cuts(); + + void eta_print(); + void eta_reset(); + + CplexRow constraint_to_cplex_row(const Constraint &cut); + + Constraint cplex_row_to_constraint(const CplexRow &cplex_row); + + void print_solution(double *x); + + void find_good_rows(); + + int n_rows; + int n_cols; + + double *ub; + double *lb; + int *cstat; + + CplexRow *cplex_rows; + + int eta_count; + int eta_total; + time_t eta_start; + + set cut_buffer; + + double *first_solution; + double *current_solution; + double *optimal_solution; + + long total_cuts; + + int current_round; + + int n_good_rows; + int *good_rows; +}; + +#include "cplex_helper.tpp" + +#endif /* CPLEX_HELPER_HPP_ */ diff --git a/onerow/library/include/onerow/cplex_helper.tpp b/onerow/library/include/onerow/cplex_helper.tpp new file mode 100644 index 0000000..8e06989 --- /dev/null +++ b/onerow/library/include/onerow/cplex_helper.tpp @@ -0,0 +1,104 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef CPLEX_HELPER_TPP_ +#define CPLEX_HELPER_TPP_ +#include +#include +#include +#include +#include + +#include +#include +#include + +using std::cout; +using std::endl; + +#include +#include + + +template +int CplexHelper::add_single_row_cuts() +{ + total_cuts = 0; + + if(n_good_rows < 0) + find_good_rows(); + + eta_reset(); + eta_count = 0; + eta_total = n_good_rows; + std::thread eta(&CplexHelper::eta_print, this); + + Stats::start_timer(); + + #pragma omp parallel for schedule(dynamic) + for (int i = 0; i < n_good_rows; i++) + { + Row *row = get_tableau_row(good_rows[i]); + //Row *row = good_rows[i]; + + Generator generator(*row); + + while (generator.has_next()) + { + Constraint *cut = generator.next(); + + if (cut->pi.nz() == 0) + { + delete cut; + continue; + } + + #ifdef PRETEND_TO_ADD_CUTS + + delete(cut); + + #else + + #pragma omp critical(cplex) + { + add_cut(cut); + } + + #endif + + #ifdef ENABLE_EXTENDED_STATISTICS + Stats::add_generated_cut(current_round, cut->depth); + #endif + } + + #pragma omp critical + { + eta_count++; + } + + delete row; + } + + Stats::end_timer(); + + eta.join(); + + flush_cuts(); + time_printf("Added %d violated cuts...\n", total_cuts); + + return 0; +} + +#endif /* CPLEX_HELPER_TPP_ */ diff --git a/onerow/library/include/onerow/geometry.hpp b/onerow/library/include/onerow/geometry.hpp new file mode 100644 index 0000000..16663ce --- /dev/null +++ b/onerow/library/include/onerow/geometry.hpp @@ -0,0 +1,56 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef GEOMETRY_HPP_ +#define GEOMETRY_HPP_ + +#include +#include +typedef q::mpq rational; +typedef q::svec svec; + +/** + * Models a two-dimensional rational point (x,y). + */ +struct Point { + + rational x; + rational y; + + Point(); + Point(rational x, rational y); + + Point operator+(const Point& p) const; + Point operator-(const Point& p) const; + Point operator*(rational scale) const; + rational operator*(const Point& p) const; +}; + +/** + * Models a two-dimensional rational line, defined by a pair of points p1 and p2. + */ +struct Line { + Point p1; + Point p2; + + Line(); + Line(Point p1, Point p2); + Line(rational x1, rational y1, rational x2, rational y2); +}; + +std::ostream& operator<<(std::ostream& os, const Line &l); +std::ostream& operator<<(std::ostream& os, const Point &p); + +#endif /* GEOMETRY_HPP_ */ diff --git a/onerow/library/include/onerow/gomory_cut_generator.hpp b/onerow/library/include/onerow/gomory_cut_generator.hpp new file mode 100644 index 0000000..07dd6a0 --- /dev/null +++ b/onerow/library/include/onerow/gomory_cut_generator.hpp @@ -0,0 +1,37 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef GOMORY_CUT_GENERATOR_HPP_ +#define GOMORY_CUT_GENERATOR_HPP_ + +#include "single_row_cut_generator.hpp" + +/** + * This class can be used to generate classic (not fractional) Gomory cuts. + * The cuts are only valid for models with integral, non-negative variables. + */ +class GomoryCutGenerator: public SingleRowCutGenerator { +private: + bool finished; + +public: + GomoryCutGenerator(Row &row); + ~GomoryCutGenerator(); + + bool has_next(); + Constraint* next(); +}; + +#endif /* GOMORY_CUT_GENERATOR_HPP_ */ diff --git a/onerow/library/include/onerow/knapsack2.hpp b/onerow/library/include/onerow/knapsack2.hpp new file mode 100644 index 0000000..a3e0043 --- /dev/null +++ b/onerow/library/include/onerow/knapsack2.hpp @@ -0,0 +1,54 @@ +#ifndef KNAPSACK_2_HPP_ +#define KNAPSACK_2_HPP_ +#include +#include +#include "qxx/rational.hpp" +#include "geometry.hpp" +using std::vector; +using std::pair; + +// ************************************************************************** +// +// ************************************************************************** +#define KNAPSACK2_LEFT 0 +#define KNAPSACK2_RIGHT 1 +#define KNAPSACK2_BOTH 2 +#define KNAPSACK2_RAY 3 + +class PairRational { +public: + size_t operator()(const pair &k) const; +}; + +struct Knapsack2Vertex { + + int side; + q::dvec lower, upper, opposed; + +}; + + +// ************************************************************************** +// +// ************************************************************************** +class Knapsack2 { + +public: + Knapsack2(); + Knapsack2(rational f, rational r1); + ~Knapsack2(); + + void clear(); + void eval(rational f, rational r1); + +private: + void push(int side, + const q::dvec &l, const q::dvec &u, const q::dvec &o); + +public: + + vector list; + +}; + +#endif diff --git a/onerow/library/include/onerow/mir_cut_generator.hpp b/onerow/library/include/onerow/mir_cut_generator.hpp new file mode 100644 index 0000000..9acb729 --- /dev/null +++ b/onerow/library/include/onerow/mir_cut_generator.hpp @@ -0,0 +1,35 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef MIR_CUT_GENERATOR_HPP_ +#define MIR_CUT_GENERATOR_HPP_ + +#include "single_row_cut_generator.hpp" + +class MIRCutGenerator: public SingleRowCutGenerator { +private: + bool finished; + rational f(rational q1, rational q2); + rational h(rational q); + +public: + MIRCutGenerator(Row &row); + ~MIRCutGenerator(); + + bool has_next(); + Constraint* next(); +}; + +#endif /* MIR_CUT_GENERATOR_HPP_ */ diff --git a/onerow/library/include/onerow/params.hpp b/onerow/library/include/onerow/params.hpp new file mode 100644 index 0000000..7746378 --- /dev/null +++ b/onerow/library/include/onerow/params.hpp @@ -0,0 +1,41 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef PARAMS_HPP_ +#define PARAMS_HPP_ + +const double ZERO_CUTOFF = 1e-8; + +const double MAX_CUT_DYNAMISM = 1e6; +const double MIN_CUT_VIOLATION = 1e-6; + +const long REDUCE_FACTOR_RHS = 1000000; +const long REDUCE_FACTOR_R1 = 1000; +const long REDUCE_FACTOR_COEFFICIENT = 1000000; + +const int MAX_CUT_DEPTH = 200; + +const int ETA_UPDATE_INTERVAL = 1; +const unsigned int MAX_CUT_BUFFER_SIZE = 100; + +const double COEFFICIENT_SAFETY_MARGIN = 0.00001; + +#define INTERSECTION_CUT_USE_DOUBLE + +#define ENABLE_TRIVIAL_LIFTING +#define ENABLE_EXTENDED_STATISTICS +//#define PRETEND_TO_ADD_CUTS + +#endif /* PARAMS_HPP_ */ diff --git a/onerow/library/include/onerow/single_row_cut_generator.hpp b/onerow/library/include/onerow/single_row_cut_generator.hpp new file mode 100644 index 0000000..5c7bec7 --- /dev/null +++ b/onerow/library/include/onerow/single_row_cut_generator.hpp @@ -0,0 +1,107 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef SINGLE_ROW_CUT_GENERATOR_HPP_ +#define SINGLE_ROW_CUT_GENERATOR_HPP_ + +#include +#include +#include "geometry.hpp" +using std::vector; + +/** + * Models a linear constraint. + */ +struct Constraint { + + /** + * Vector that holds the coefficients of the variables. + */ + svec pi; + + /** + * The right hand side of the constraint. + */ + rational pi_zero; + + int depth; + + /** + * Comparator used to sort the constraints. + * + * @returns bool True if this constraint should come before the given + * constraint when sorting. + */ + bool operator<(const Constraint &other) const; + + bool operator==(const Constraint &other) const; + + rational get_violation(const rational *x); +}; + +/** + * Models a single row from the simplex tableau. + */ +struct Row { + /** + * Constraint + */ + Constraint c; + + /** + * Index of the basic variable this row corresponds to. + */ + int basic_var_index; + + bool* is_integer; +}; + +/** + * A single row cut generator receives a row from the simplex tableau and generates one + * or more cuts that invalidate the current basic solution. + */ +class SingleRowCutGenerator { +protected: + const Row& row; + +public: + + /** + * Constructs a new generator that will generate cuts from the provided tableau row. + */ + SingleRowCutGenerator(Row &r); + + /** + * Destructor. + */ + virtual ~SingleRowCutGenerator() {}; + + /** + * Returns true if more cuts can be generated from the tableau row. + */ + virtual bool has_next() = 0; + + /** + * Retrieves a cut generated from the tableau row. + * + * @throws std::out_of_bounds if no more cuts can be generated from the current tableau row. + */ + virtual Constraint* next() = 0; +}; + +std::ostream& operator<<(std::ostream& os, const Constraint &c); +std::ostream& operator<<(std::ostream& os, const Row &r); + +#endif diff --git a/onerow/library/include/onerow/stats.hpp b/onerow/library/include/onerow/stats.hpp new file mode 100644 index 0000000..99ef0dd --- /dev/null +++ b/onerow/library/include/onerow/stats.hpp @@ -0,0 +1,41 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef STATS_HPP_ +#define STATS_HPP_ + +#include +using std::string; + +namespace Stats { + void init(); + void add_cut(int depth); + void add_generated_cut(int round, int depth); + void set_solution(int round, double sol, string status); + void set_input_filename(string n); + + void add_trivial_lifting_m(unsigned long m); + void add_coefficient(bool integral); + + void write_stats(string filename); + + void start_timer(); + void end_timer(); +}; + +double get_current_time(void); +void time_printf(const char *fmt, ...); + +#endif /* STATS_HPP_ */ diff --git a/onerow/library/include/onerow/wedge_cut_generator.hpp b/onerow/library/include/onerow/wedge_cut_generator.hpp new file mode 100644 index 0000000..ff9c92f --- /dev/null +++ b/onerow/library/include/onerow/wedge_cut_generator.hpp @@ -0,0 +1,132 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef WEDGE_CUT_GENERATOR_HPP_ +#define WEDGE_CUT_GENERATOR_HPP_ + +#include "params.hpp" +#include "qxx/dmat.hpp" +#include "geometry.hpp" +#include "knapsack2.hpp" +#include "single_row_cut_generator.hpp" + +/** + * Models a two-dimensional intersection cut, and its associated convex set. + * The convex set is given by + * \f[ + * C = \{ x \in R^n : (d^i)^T(x-f) \leq 1 \}, + * \f] + * where f is a point in the interior of C. + */ +class IntersectionCut { +public: + Point f; + Point *d; + const int n_faces; + + /** + * Constructs a new intersection cut and its associated convex set from + * the data provided. + * + * @param f A fractional point in the interior of the convex set. + * @param n_faces The number of faces of the convex set. + * @param lines (Optional) Array of lines that support each face of the + * convex. + */ + IntersectionCut(Point f, int n_faces, const Line *lines = 0); + ~IntersectionCut(); + + /** + * Sets the desired face of the convex set. + * + * @param index Index of the face. + * @param line Line that supports the face. + */ + void set_face(int index, Line line); + + rational get_continuous_coefficient(rational rx, rational ry); + + double get_trivial_lifting_coefficient_double(double rx, double ry); + rational get_trivial_lifting_coefficient(rational rx, rational ry); + void pre_lifting(); + +private: + double d_r0x, d_p, d_d0x, d_d0y, d_d1x, d_d1y; + rational p, r0x; + bool pre_lifting_ready; +}; + +/** + * Models an intersection cut associated with a wedge. + */ +class WedgeCut : public IntersectionCut { +public: + /** + * Constructs a wedge cut from the provided data. + * + * @param f A fractional point in the interior of the wedge. + * @param left A point on the left face of the wedge. + * @param apex The apex of the wedge. + * @param right A point on the right face of the wedge. + */ + WedgeCut(Point f, Point left, Point apex, Point right); +}; + +/** + * Models an intersection cut associated with a split. + */ +class SplitCut : public IntersectionCut { +public: + /** + * Constructs a split cut from the provided data. + * + * @param f A fractional point in the interior of the split + * @param left A point on the left face of the split + * @param right A poit on the right face of the split + * @param direction The direction in which the split in unbounded. + */ + SplitCut(Point f, Point left, Point right, Point direction); +}; + +/** + * This class can be used to generate wedge cuts. + */ +class WedgeCutGenerator: public SingleRowCutGenerator { +private: + bool finished; + q::dvec f, r1; + int r1_offset; + int cur_facet; + Knapsack2 knapsack; + static int max_depth; + int n_knapsacks; + +private: + Constraint *cut; + void eval_next(); + static q::dvec intersection( + q::dvec a, q::dvec b, q::dvec c, q::dvec d); + +public: + WedgeCutGenerator(Row &row); + ~WedgeCutGenerator(); + + bool has_next(); + Constraint* next(); +}; + + + +#endif /* WEDGE_CUT_GENERATOR_HPP_ */ diff --git a/onerow/library/src/cplex_helper.cpp b/onerow/library/src/cplex_helper.cpp new file mode 100644 index 0000000..307e6f8 --- /dev/null +++ b/onerow/library/src/cplex_helper.cpp @@ -0,0 +1,533 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::string; + +static bool debug = false; + +CplexHelper::CplexHelper(CPXENVptr _env, CPXLPptr _lp) : + env(_env), lp(_lp), is_integer(0), n_cuts(0), n_rows(0), ub(0), lb(0), + cstat(0), cplex_rows(0), first_solution(0), current_solution(0), + optimal_solution(0), current_round(0), n_good_rows(-1) +{ +} + + +CplexHelper::~CplexHelper() +{ + if (cplex_rows) + { + for (int i = 0; i < n_rows; i++) + { + delete[] cplex_rows[i].pi; + delete[] cplex_rows[i].indices; + } + + delete[] cplex_rows; + } + + if(ub) delete[] ub; + if(lb) delete[] lb; + if(cstat) delete[] cstat; + + if (is_integer) + delete[] is_integer; +} + + +#define return_if_neq(a,b) if((a)<(b)) return true; if((a)>(b)) return false; + +bool CplexRow::operator<(const CplexRow &other) const +{ + return_if_neq(fabs(pi_zero-0.5), fabs(other.pi_zero-0.5)); + + return_if_neq(nz, other.nz); + + for (int i = 0; i < nz; i++) + { + return_if_neq(indices[i], other.indices[i]); + double diff = pi[i] - other.pi[i]; + if(diff < -ZERO_CUTOFF) return true; + } + + return false; +} + +void CplexRow::print(double *solution) +{ + for (int i = 0; i < nz; i++) + { + printf("%lf x%d ", pi[i], indices[i]); + + if (solution) + printf("(%lf) ", solution[indices[i]]); + } + + printf("<= %lf ", pi_zero); + + if (solution) + printf("violation=%lf", get_violation(solution)); + + printf("\n"); +} + +double CplexRow::get_violation(double *solution) +{ + double v = 0; + for (int i = 0; i < nz; i++) + { + v += pi[i] * solution[indices[i]]; + } + v -= pi_zero; + + return v; +} + + +Constraint CplexHelper::cplex_row_to_constraint(const CplexRow &cplex_row) +{ + Constraint constraint; + constraint.pi.resize(n_cols); + + rational pi_zero = rational(cplex_row.pi_zero); + + for (int j = 0; j < cplex_row.nz; j++) + { + int index = cplex_row.indices[j]; + rational pij = rational(cplex_row.pi[j]); + + if (cstat[index] == CPX_AT_LOWER) + pi_zero -= rational(lb[index]) * pij; + + if (cstat[index] == CPX_AT_UPPER) + { + pi_zero -= rational(ub[index]) * pij; + pij = -pij; + } + + constraint.pi.push_nz(cplex_row.indices[j], pij.reduce(REDUCE_FACTOR_COEFFICIENT)); + } + + constraint.pi_zero = rational(pi_zero).reduce(REDUCE_FACTOR_RHS); + + return constraint; +} + + +CplexRow CplexHelper::constraint_to_cplex_row(const Constraint &cut) +{ + CplexRow cplex_row; + cplex_row.pi = new double[n_cols]; + cplex_row.indices = new int[cut.pi.nz()]; + cplex_row.pi_zero = cut.pi_zero.get_double(); + cplex_row.depth = cut.depth; + cplex_row.nz = cut.pi.nz(); + + double max_pi = -INFINITY; + double min_pi = INFINITY; + + int cut_nz = cut.pi.nz(); + for (int j = 0; j < cut_nz; j++) + { + int index = cut.pi.index(j); + double pij = cut.pi.value(j).get_double(); + + if (fabs(pij) < ZERO_CUTOFF) + pij = 0; + + max_pi = std::max(max_pi, fabs(pij)); + + if (fabs(pij) > 0) + min_pi = std::min(min_pi, fabs(pij)); + + if (cstat[index] == CPX_AT_LOWER) + cplex_row.pi_zero += lb[index] * pij; + + if (cstat[index] == CPX_AT_UPPER) + { + pij = -pij; + cplex_row.pi_zero += ub[index] * pij; + } + + cplex_row.indices[j] = index; + cplex_row.pi[j] = pij; + } + + cplex_row.dynamism = max_pi / min_pi; + + return cplex_row; +} + + +void CplexHelper::add_cut(Constraint *cut) +{ + + CplexRow cplex_row = constraint_to_cplex_row(*cut); + double violation = cplex_row.get_violation(current_solution); + + if (optimal_solution) + assert(cplex_row.get_violation(optimal_solution) <= MIN_CUT_VIOLATION); + + if (first_solution) + assert(cplex_row.get_violation(first_solution) >= MIN_CUT_VIOLATION); + + #pragma omp critical + if (cplex_row.dynamism < MAX_CUT_DYNAMISM && violation >= MIN_CUT_VIOLATION) + { + auto p = cut_buffer.insert(cplex_row); + + // duplicate cut + if (!p.second) + { + delete[] cplex_row.pi; + delete[] cplex_row.indices; + } + + if (cut_buffer.size() >= MAX_CUT_BUFFER_SIZE) + { + flush_cuts(); + solve(false); + } + } + // rejected cut + else + { + delete[] cplex_row.pi; + delete[] cplex_row.indices; + } + + delete cut; +} + + +void CplexHelper::flush_cuts() +{ + int begin = 0; + char sense = 'L'; + + for (CplexRow cplex_row : cut_buffer) + { + Stats::add_cut(cplex_row.depth); + + total_cuts++; + + CPXaddrows(env, lp, 0, 1, cplex_row.nz, &cplex_row.pi_zero, &sense, + &begin, cplex_row.indices, cplex_row.pi, NULL, NULL); + + delete[] cplex_row.pi; + delete[] cplex_row.indices; + } + + cut_buffer.clear(); +} + + +Row* CplexHelper::get_tableau_row(int index) +{ + Row *row = new Row; + row->basic_var_index = cplex_rows[index].head; + row->is_integer = is_integer; + row->c = cplex_row_to_constraint(cplex_rows[index]); + + if (optimal_solution) + assert(cplex_rows[index].get_violation(optimal_solution) <= 0.001); + + if (debug) + { + printf("ROW %d\n", index); + cplex_rows[index].print(optimal_solution); + } + + return row; +} + +int CplexHelper::solve(bool should_end_round) +{ + // Optimize + int status = CPXlpopt(env, lp); + if (status) + { + fprintf(stderr, "Could not optimize (%d)\n", status); + throw std::runtime_error("CplexHelper::solve_mip"); + } + + // Get status + char buffer[512]; + double objval; + + status = CPXgetstat(env, lp); + + CPXgetstatstring(env, status, buffer); + CPXgetobjval(env, lp, &objval); + time_printf(" %.6lf [%s] \n", objval, buffer); + + assert(status == 1); + + // Store current solution + if (!current_solution) + current_solution = new double[n_cols]; + + CPXgetx(env, lp, current_solution, 0, n_cols - 1); + + // During first round, store basis and fractional solution + if (current_round == 0) + { + read_basis(); + + first_solution = new double[n_cols]; + for (int i = 0; i < n_cols; i++) + first_solution[i] = current_solution[i]; + } + + if (should_end_round) + Stats::set_solution(current_round++, objval, string(buffer)); + + return objval; +} + + +void CplexHelper::read_columns() +{ + n_rows = CPXgetnumrows(env, lp); + n_cols = CPXgetnumcols(env, lp); + + is_integer = new bool[n_cols]; + + char ctype[n_cols]; + CPXgetctype(env, lp, ctype, 0, n_cols - 1); + + for (int i = 0; i < n_cols; i++) + { + switch (ctype[i]) + { + case CPX_BINARY: + case CPX_INTEGER: + case CPX_SEMIINT: + is_integer[i] = true; + break; + + default: + is_integer[i] = false; + break; + } + } + + time_printf("Fetched %d rows, %d cols.\n", n_rows, n_cols); +} + + +void CplexHelper::read_basis() +{ + time_printf("Reading basis...\n"); + + ub = new double[n_cols]; + lb = new double[n_cols]; + cstat = new int[n_cols]; + + int *head = new int[n_rows]; + int *rstat = new int[n_rows]; + double *rhs = new double[n_rows]; + + CPXgetbhead(env, lp, head, rhs); + CPXgetub(env, lp, ub, 0, n_cols - 1); + CPXgetlb(env, lp, lb, 0, n_cols - 1); + CPXgetbase(env, lp, cstat, rstat); + + cplex_rows = new CplexRow[n_rows]; + assert(cplex_rows != 0); + + eta_reset(); + eta_count = 0; + eta_total = n_rows; + std::thread eta(&CplexHelper::eta_print, this); + + for (int i = 0; i < n_rows; i++) + { + int nz = 0; + double pi[n_cols]; + + CPXbinvarow(env, lp, i, pi); + + for (int j = 0; j < n_cols; j++) + { + if (fabs(pi[j]) < ZERO_CUTOFF) + continue; + + if (cstat[j] == CPX_AT_LOWER) + rhs[i] += lb[j] * pi[j]; + + if (cstat[j] == CPX_AT_UPPER) + rhs[i] += ub[j] * pi[j]; + + nz++; + } + + cplex_rows[i].nz = nz; + cplex_rows[i].depth = 0; + cplex_rows[i].pi = new double[nz]; + cplex_rows[i].indices = new int[nz]; + cplex_rows[i].pi_zero = rhs[i]; + cplex_rows[i].head = head[i]; + + if(fabs(cplex_rows[i].pi_zero) < ZERO_CUTOFF) + cplex_rows[i].pi_zero = 0; + + int k = 0; + for (int j = 0; j < n_cols; j++) + { + if (fabs(pi[j]) < ZERO_CUTOFF) + continue; + cplex_rows[i].pi[k] = pi[j]; + cplex_rows[i].indices[k++] = j; + } + + eta_count++; + } + + eta.join(); + + delete[] head; + delete[] rstat; + delete[] rhs; +} + + +void CplexHelper::print_basis() +{ + int n_rows = CPXgetnumrows(env, lp); + int n_cols = CPXgetnumcols(env, lp); + + double y[n_rows]; + double ub[n_cols]; + double lb[n_cols]; + double dj[n_cols]; + + CPXgetub(env, lp, ub, 0, n_cols - 1); + CPXgetlb(env, lp, lb, 0, n_cols - 1); + CPXgetdj(env, lp, dj, 0, n_cols - 1); + + cout << "basis inverse:" << endl; + for (int i = 0; i < n_rows; i++) + { + CPXbinvrow(env, lp, i, y); + for (int k = 0; k < n_rows; k++) + cout << rational(y[k]) << " "; + cout << endl; + } + + int cstat[n_cols]; + int rstat[n_rows]; + CPXgetbase(env, lp, cstat, rstat); + + cout << "column status:" << endl; + for (int i = 0; i < n_cols; i++) + { + cout << i << ": " << cstat[i] << " " << lb[i] << "..." << ub[i] << " " + << dj[i] << " " << (is_integer[i] ? "int" : "cont") << endl; + } + +} + + +void CplexHelper::print_solution(double *x) +{ + for (int i = 0; i < n_cols; i++) + if (fabs(x[i]) > ZERO_CUTOFF) + time_printf(" x%d = %.6lf\n", i, x[i]); +} + +void CplexHelper::eta_reset() +{ + time(&eta_start); +} + +void CplexHelper::eta_print() +{ + while (true) + { + sleep(ETA_UPDATE_INTERVAL); + + if (eta_count == 0) + { + printf("\r\r%3.0f%% ETA: unknown\r", 0.0); + fflush(stdout); + continue; + } + + if (eta_count >= eta_total) + break; + + time_t eta_now; + time(&eta_now); + + double diff = difftime(eta_now, eta_start); + double eta = diff / eta_count * eta_total; + + time_t eta_date = eta_start + eta; + tm *ttm = localtime(&eta_date); + + printf("\r%3.0f%% ", 100.0 * eta_count / eta_total); + printf("ETA: %04d-%02d-%02d %02d:%02d:%02d %d / %d\r", + ttm->tm_year + 1900, ttm->tm_mon + 1, ttm->tm_mday, + ttm->tm_hour, ttm->tm_min, ttm->tm_sec, + eta_count, eta_total); + fflush(stdout); + } + + printf("\r \r"); +} + +void CplexHelper::find_good_rows() +{ + n_good_rows = 0; + good_rows = new int[n_rows]; + + time_printf("Finding interesting rows...\n"); + + #pragma omp parallel for schedule(dynamic) + for (int i = 0; i < n_rows; i++) + { + Row *row = get_tableau_row(i); + + if (row->c.pi_zero.frac() != 0 && + row->is_integer[row->basic_var_index]) + { + good_rows[n_good_rows++] = i; + } + + delete row; + } + + time_printf(" %d rows found\n", n_good_rows, n_rows); +} diff --git a/onerow/library/src/geometry.cpp b/onerow/library/src/geometry.cpp new file mode 100644 index 0000000..4d9f558 --- /dev/null +++ b/onerow/library/src/geometry.cpp @@ -0,0 +1,77 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include + +Point::Point() : + x(0), y(0) +{ + +} + +Point::Point(rational _x, rational _y) : + x(_x), y(_y) +{ + +} + +Point Point::operator+(const Point& p) const +{ + return Point(x + p.x, y + p.y); +} + +Point Point::operator-(const Point& p) const +{ + return Point(x - p.x, y - p.y); +} + +rational Point::operator*(const Point& p) const +{ + return x * p.x + y * p.y; +} + +Point Point::operator*(rational scale) const +{ + return Point(scale * x, scale * y); +} + +Line::Line() +{ + +} + +Line::Line(Point _p1, Point _p2) : + p1(_p1), p2(_p2) +{ + +} + +Line::Line(rational x1, rational y1, rational x2, rational y2) : + p1(Point(x1, y1)), p2(Point(x2, y2)) +{ + +} + +std::ostream& operator<<(std::ostream& os, const Point &p) +{ + os << "(" << p.x << "," << p.y << ")"; + return os; +} + +std::ostream& operator<<(std::ostream& os, const Line &l) +{ + os << l.p1 << "--" << l.p2; + return os; +} diff --git a/onerow/library/src/gomory_cut_generator.cpp b/onerow/library/src/gomory_cut_generator.cpp new file mode 100644 index 0000000..62d7eae --- /dev/null +++ b/onerow/library/src/gomory_cut_generator.cpp @@ -0,0 +1,50 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include + +GomoryCutGenerator::GomoryCutGenerator(Row &r) : + SingleRowCutGenerator(r), finished(false) +{ + +} + +GomoryCutGenerator::~GomoryCutGenerator() +{ + +} + +bool GomoryCutGenerator::has_next() +{ + return !finished; +} + +Constraint* GomoryCutGenerator::next() +{ + if (!has_next()) + throw std::out_of_range(""); + + Constraint *cut = new Constraint; + int nz = row.c.pi.nz(); + cut->pi_zero = row.c.pi_zero.floor(); + + cut->pi.resize(row.c.pi.size()); + for (int i = 0; i < nz; i++) + cut->pi.push(row.c.pi.index(i), row.c.pi.value(i).floor()); + + finished = true; + return cut; +} diff --git a/onerow/library/src/knapsack2.cpp b/onerow/library/src/knapsack2.cpp new file mode 100644 index 0000000..fda6295 --- /dev/null +++ b/onerow/library/src/knapsack2.cpp @@ -0,0 +1,248 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include + +#include +#include + +#include +#include + +using std::endl; +using std::cout; + + +// ************************************************************************** +// +// ************************************************************************** +//#define DEBUG(s...) gmp_printf(s) +#define DEBUG(s...) + +// ************************************************************************** +// +// ************************************************************************** +Knapsack2::Knapsack2() +{ +} + +Knapsack2::Knapsack2(rational f, rational r1) +{ + eval(f, r1); +} + +Knapsack2::~Knapsack2() +{ +} + +// ************************************************************************** +// +// ************************************************************************** +void Knapsack2::clear() +{ + list.clear(); +} + + +// ************************************************************************** +// +// ************************************************************************** +void Knapsack2::push(int side, + const q::dvec &l, const q::dvec &u, const q::dvec &o) +{ + Knapsack2Vertex v; + + v.side = side; + v.lower = l; + v.upper = u; + v.opposed = o; + v.lower.resize(2); + v.upper.resize(2); + v.opposed.resize(2); + + list.push_back(v); +} + +//int hit_count = 0; +//int miss_count = 0; + +// ************************************************************************** +// +// ************************************************************************** +void Knapsack2::eval(rational fx, rational r1x) +{ + clear(); + + rational fr1frac, fr1floor, hslope; + q::dvec a(3), b(3), p(3), f(3), g(3), h(3), r1(3), wy(3), wa(3), wb(3); + q::dmat w(3, 3), u(3, 3), ta(3, 3), tb(3, 3), tx(3, 3); + + a[0] = 0; a[1] = 0; a[2] = 1; + b[0] = 1; b[1] = 0; b[2] = 1; + p[2] = 1; + f[0] = fx.frac(); f[1] = 0; f[2] = 1; + g[2] = 1; + h[2] = 1; + r1[0] = r1x; r1[1] = 1; r1[2] = 0; + + w.set_identity(); + w(0, 2) = fx.floor(); + tb(0, 1) = 1; tb(0, 2) = 1; tb(1, 2) = 1; + tb(2, 0) = tb(2, 1) = tb(2, 2) = 1; + ta(0, 1) = 1; ta(1, 2) = 1; + ta(2, 0) = ta(2, 1) = ta(2, 2) = 1; + + int it = 0; + + wa = w * a; + DEBUG("-: A vertex (0, 0)\t\t--> v A (%Qd, %Qd)\n", + wa[0].v, wa[1].v); + wb = w * b; + DEBUG("-: B vertex (1, 0)\t\t--> v B (%Qd, %Qd)\n", + wb[0].v, wb[1].v); + + while (1) { + DEBUG("%d: -----------------------------\n", it); + DEBUG("%d: f = (%Qd, %Qd, %Qd)\n", + it, f[0].v, f[1].v, f[2].v); + DEBUG("%d: r1 = (%Qd, %Qd, %Qd)\n", + it, r1[0].v, r1[1].v, r1[2].v); + DEBUG("%d: a = (%Qd, %Qd, %Qd)\n", + it, a[0].v, a[1].v, a[2].v); + DEBUG("%d: b = (%Qd, %Qd, %Qd)\n", + it, b[0].v, b[1].v, b[2].v); + + // Step 1 + fr1frac = (f[0] + r1[0]).frac(); + fr1floor = (f[0] + r1[0]).floor(); + + if (f[0] == fr1frac) { + wy = w * r1; + DEBUG("%d: AB ray (%Qd, %Qd)" + "\t\t--> r AB (%Qd, %Qd)\n", + it, r1[0].v, r1[1].v, wy[0].v, wy[1].v); + push(KNAPSACK2_RAY, wa, wy, wb); + break; + } else if (f[0] < fr1frac) { + DEBUG("%d: hit right\n", it); + + g[1] = (q::mpq(1) - f[0]) / (fr1frac - f[0]); + g[0] = q::mpq(1) + fr1floor * g[1]; + + b[1] = g[1].floor(); + b[0] = q::mpq(1) + fr1floor * b[1]; + + p[1] = b[1] + 1; + p[0] = q::mpq(1) + fr1floor * p[1]; + + if (b[1] == g[1]) { + wy = w * b; + DEBUG("%d: AB vertex (%Qd, %Qd)" + "\t\t--> v AB (%Qd, %Qd)\n", + it, b[0].v, b[1].v, wy[0].v, wy[1].v); + push(KNAPSACK2_BOTH, wa, wy, wb); + break; + } + + tx(0, 0) = 0; + tx(1, 0) = 0; + tx(2, 0) = 1; + tx.set_col(1, b); + tx.set_col(2, p); + + u = tb * tx.inv(); + + //tb.dump("tb"); + //tx.dump("tx"); + //u.dump("u"); + + wy = w * b; + + push(KNAPSACK2_RIGHT, wb, wy, wa); + + wb = wy; + + DEBUG("%d: B vertex (%Qd, %Qd)" + "\t\t--> v B (%Qd, %Qd)" + " opposed (%Qd, %Qd)\n", + it, b[0].v, b[1].v, wb[0].v, wb[1].v, + wa[0].v, wa[1].v); + + hslope = b[0] / b[1]; + h[1] = f[0] / (hslope - r1[0]); + h[0] = hslope * h[1]; + } else { + DEBUG("%d: hit left\n", it); + + g[1] = f[0] / (f[0] - fr1frac); + g[0] = fr1floor * g[1]; + + a[1] = g[1].floor(); + a[0] = fr1floor * a[1]; + + p[1] = a[1] + 1; + p[0] = fr1floor * p[1]; + + if (a[1] == g[1]) { + wy = w * a; + DEBUG("%d: AB vertex (%Qd, %Qd)" + "\t\t--> v AB (%Qd, %Qd)\n", + it, a[0].v, a[1].v, wy[0].v, wy[1].v); + push(KNAPSACK2_BOTH, wa, wy, wb); + break; + } + + tx.set_col(0, a); + tx(0, 1) = 1; + tx(1, 1) = 0; + tx(2, 1) = 1; + tx.set_col(2, p); + + u = ta * tx.inv(); + + //tb.dump("ta"); + //u.dump("u"); + + wy = w * a; + + push(KNAPSACK2_LEFT, wa, wy, wb); + + wa = wy; + + DEBUG("%d: A vertex (%Qd, %Qd)" + "\t\t--> v A (%Qd, %Qd)" + " opposed (%Qd, %Qd)\n", + it, a[0].v, a[1].v, wa[0].v, wa[1].v, + wb[0].v, wb[1].v); + + hslope = (a[0] - 1) / a[1]; + h[1] = (f[0] - 1) / (hslope - r1[0]); + h[0] = q::mpq(1) + (hslope * h[1]); + } + + a = u * a; + b = u * b; + f = u * h; + r1 = u * r1; + r1[0] /= r1[1]; + r1[1] = 1; + + w = w * u.inv(); + + it++; + } +} + diff --git a/onerow/library/src/main.cpp b/onerow/library/src/main.cpp new file mode 100644 index 0000000..5fd4133 --- /dev/null +++ b/onerow/library/src/main.cpp @@ -0,0 +1,169 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include + +#include "cplex_helper.hpp" +#include "cplex_helper.tpp" +#include "stats.hpp" + +#include +#include + +#include "gomory_cut_generator.hpp" +#include "mir_cut_generator.hpp" +#include "wedge_cut_generator.hpp" + +using namespace std; + +char *input_filename = 0; +char *stats_filename = 0; +char *sol_filename = 0; +bool enable_gomory_cuts = false; +bool enable_wedge_cuts = false; +bool enable_mir_cuts = false; + +int c; +extern char *optarg; +char usage[] = "usage: %s [-gmw] -f model.mps [-s stats.yaml]\n"; + +void read_params(int argc, char **argv) +{ + while ((c = getopt(argc, argv, "gwmf:s:x:")) != -1) + { + switch (c) + { + case 'g': + enable_gomory_cuts = true; + break; + case 'w': + enable_wedge_cuts = true; + break; + case 'm': + enable_mir_cuts = true; + break; + case 'f': + input_filename = optarg; + break; + case 's': + stats_filename = optarg; + break; + case 'x': + sol_filename = optarg; + break; + } + } + + if (!input_filename) + { + fprintf(stderr, "%s: missing model filename\n", argv[0]); + fprintf(stderr, usage, argv[0]); + exit(1); + } +} + +int main(int argc, char **argv) +{ + read_params(argc, argv); + + Stats::init(); + + int status; + CPXENVptr env = CPXopenCPLEX(&status); + CPXLPptr lp = CPXcreateprob(env, &status, ""); + CplexHelper cplexHelper(env, lp); + + CPXsetlogfile(env, NULL); + CPXsetintparam(env, CPX_PARAM_PREIND, CPX_OFF); // disable presolve + CPXsetintparam(env, CPX_PARAM_DATACHECK, CPX_ON); // check consistency + CPXsetintparam(env, CPX_PARAM_NUMERICALEMPHASIS, CPX_ON); // numerical precision + + Stats::set_input_filename(string(input_filename)); + + // reads input file + time_printf("Reading input file: %s...\n", input_filename); + status = CPXreadcopyprob(env, lp, input_filename, NULL); + if (status) + { + fprintf(stderr, "could not read input file (%d)\n", status); + return 1; + } + + cplexHelper.read_columns(); + + // read solution + if(sol_filename) + { + time_printf("Reading solution file: %s...\n", sol_filename); + FILE *f = fopen(sol_filename, "r"); + if(!f) + { + fprintf(stderr, "Could not open solution file (%s).", sol_filename); + return 1; + } + + double *solution = new double[cplexHelper.n_cols]; + for(int i=0; i(); + cplexHelper.solve(true); + } + + if (enable_mir_cuts) + { + time_printf("Generating MIR cuts...\n"); + cplexHelper.add_single_row_cuts(); + cplexHelper.solve(true); + } + + if (enable_wedge_cuts) + { + time_printf("Generating wedge cuts...\n"); + cplexHelper.add_single_row_cuts(); + cplexHelper.solve(true); + } + + if(stats_filename != 0) + { + time_printf("Writting stats: %s...\n", stats_filename); + Stats::write_stats(string(stats_filename)); + } + + time_printf("Done.\n"); + CPXcloseCPLEX(&env); + + return 0; + +} diff --git a/onerow/library/src/mir_cut_generator.cpp b/onerow/library/src/mir_cut_generator.cpp new file mode 100644 index 0000000..d54e4d9 --- /dev/null +++ b/onerow/library/src/mir_cut_generator.cpp @@ -0,0 +1,74 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include + +MIRCutGenerator::MIRCutGenerator(Row &r) : + SingleRowCutGenerator(r), finished(false) +{ + +} + +MIRCutGenerator::~MIRCutGenerator() +{ + +} + +bool MIRCutGenerator::has_next() +{ + return !finished; +} + +rational MIRCutGenerator::h(rational a) +{ + if (a > 0) return a; + else return 0; +} + +rational MIRCutGenerator::f(rational a, rational b) +{ + if (a.frac() <= b.frac()) + return b.frac() * a.floor() + a.frac(); + else + return b.frac() * a.ceil(); +} + +Constraint* MIRCutGenerator::next() +{ + if (!has_next()) + throw std::out_of_range(""); + + Constraint *cut = new Constraint; + int nz = row.c.pi.nz(); + cut->pi.resize(row.c.pi.size()); + cut->pi_zero = -row.c.pi_zero.frac() * row.c.pi_zero.ceil(); + + for (int i = 0; i < nz; i++) + { + int idx = row.c.pi.index(i); + + if (row.is_integer[idx]) + cut->pi.push(idx, -f(row.c.pi.value(i), row.c.pi_zero)); + else + cut->pi.push(idx, -h(row.c.pi.value(i))); + } + + cut->depth = 0; + + finished = true; + return cut; +} diff --git a/onerow/library/src/single_row_cut_generator.cpp b/onerow/library/src/single_row_cut_generator.cpp new file mode 100644 index 0000000..416fc82 --- /dev/null +++ b/onerow/library/src/single_row_cut_generator.cpp @@ -0,0 +1,72 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include + +SingleRowCutGenerator::SingleRowCutGenerator(Row &r) : row(r) +{ + +} + +#define return_if_neq(a,b) if((a)<(b)) return true; if((a)>(b)) return false; + +bool Constraint::operator<(const Constraint &other) const +{ + return_if_neq(pi.nz(), other.pi.nz()); + return_if_neq(pi_zero, other.pi_zero); + + int nz = pi.nz(); + for (int i = 0; i < nz; i++) + { + return_if_neq(pi.index(i), other.pi.index(i)); + return_if_neq(pi.value(i), other.pi.value(i)); + } + + return false; +} + +bool Constraint::operator==(const Constraint &other) const +{ + return !(operator<(other) || other.operator<(*this)); +} + +rational Constraint::get_violation(const rational *x) +{ + rational v(0); + + int nz = pi.nz(); + for (int i = 0; i < nz; i++) + v += pi.value(i) * x[pi.index(i)]; + v -= pi_zero; + + return v; +} + +std::ostream& operator<<(std::ostream& os, const Constraint &c) +{ + int nz = c.pi.nz(); + for (int k = 0; k < nz; k++) + { + os << c.pi.value(k) << " x" << c.pi.index(k) << " "; + } + os << "<= " << c.pi_zero; + return os; +} + +std::ostream& operator<<(std::ostream& os, const Row &r) +{ + os << "[" << r.basic_var_index << "] " << r.c; + return os; +} diff --git a/onerow/library/src/stats.cpp b/onerow/library/src/stats.cpp new file mode 100644 index 0000000..8d42ff2 --- /dev/null +++ b/onerow/library/src/stats.cpp @@ -0,0 +1,226 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +#include +#include + +#include +#include + +namespace Stats +{ + const int MAX_ROUNDS = 100; + const int MAX_TIMERS = 100; + + string input_filename; + + double opt_value[MAX_ROUNDS]; + string opt_sol_status[MAX_ROUNDS]; + + unsigned long n_cuts_total = 0; + unsigned long n_cuts_depth[MAX_CUT_DEPTH] = { 0 }; + + unsigned long n_generated_cuts_total = 0; + unsigned long n_generated_cuts_round[MAX_CUT_DEPTH] = { 0 }; + unsigned long n_generated_cuts_depth[MAX_CUT_DEPTH] = { 0 }; + + unsigned long trivial_lifting_m_count = 0; + unsigned long trivial_lifting_m_sum = 0; + unsigned long trivial_lifting_m_max = 0; + + unsigned long n_coefficients = 0; + unsigned long n_integral_coefficients = 0; + + int n_timers = 0; + double current_timer_start; + double timers[MAX_TIMERS] = {0}; + + void init() + { + input_filename = ""; + + for (int i = 0; i < MAX_ROUNDS; i++) + { + opt_value[i] = INFINITY; + opt_sol_status[i] = ""; + } + } + + void add_cut(int depth) + { + n_cuts_total++; + n_cuts_depth[depth]++; + } + + void add_generated_cut(int round, int depth) + { + n_generated_cuts_total++; + n_generated_cuts_round[round]++; + n_generated_cuts_depth[depth]++; + } + + void add_trivial_lifting_m(unsigned long m) + { + trivial_lifting_m_count++; + trivial_lifting_m_sum += m; + trivial_lifting_m_max = std::max(trivial_lifting_m_max, m); + } + + void add_coefficient(bool integral) + { + n_coefficients++; + if(integral) n_integral_coefficients++; + } + + void set_solution(int round, double sol, string status) + { + opt_value[round] = sol; + opt_sol_status[round] = status; + } + + void set_input_filename(string n) + { + input_filename = n; + } + + void write_stats(string f) + { + FILE *out = fopen(f.c_str(), "w"); + + fprintf(out, "input_file:\n %s\n", input_filename.c_str()); + + // solution value and status + fprintf(out, "sol_value:\n"); + for (int i = 0; i < MAX_ROUNDS; i++) + if (opt_value[i] < INFINITY) + fprintf(out, " %d: %-.6f\n", i, opt_value[i]); + + fprintf(out, "sol_status:\n"); + for (int i = 0; i < MAX_ROUNDS; i++) + if (opt_value[i] < INFINITY) + fprintf(out, " %d: %s\n", i, opt_sol_status[i].c_str()); + + // added cuts + if(n_cuts_total > 0) + { + fprintf(out, "n_added_cuts:\n"); + + fprintf(out, " total:\n %ld\n", n_cuts_total); + + fprintf(out, " depth:\n"); + for (int i = 0; i < MAX_CUT_DEPTH; i++) + if (n_cuts_depth[i] > 0) + fprintf(out, " %d: %ld\n", i, n_cuts_depth[i]); + } + + // generated cuts + if(n_generated_cuts_total > 0) + { + fprintf(out, "n_generated_cuts:\n"); + + fprintf(out, " total:\n %ld\n", n_generated_cuts_total); + + fprintf(out, " round:\n"); + for (int i = 0; i < MAX_ROUNDS; i++) + if (n_generated_cuts_round[i] > 0) + fprintf(out, " %d: %ld\n", i, n_generated_cuts_round[i]); + + fprintf(out, " depth:\n"); + for (int i = 0; i < MAX_CUT_DEPTH; i++) + if (n_generated_cuts_depth[i] > 0) + fprintf(out, " %d: %ld\n", i, n_generated_cuts_depth[i]); + } + + // trivial lifting + if(trivial_lifting_m_count > 0) + { + double integral_coefficients = ((double) n_integral_coefficients) / n_coefficients; + double average_m = ((double) (trivial_lifting_m_sum)) + / trivial_lifting_m_count; + double slowdown = integral_coefficients * average_m; + + fprintf(out, "trivial_lifting:\n"); + fprintf(out, " max_m: %ld\n", trivial_lifting_m_max); + fprintf(out, " average_m: %.6lf\n", average_m); + fprintf(out, " integral_coefficients: %.6lf\n", integral_coefficients); + fprintf(out, " slowdown: %.6lf\n", slowdown); + } + + if(n_timers > 0) + { + fprintf(out, "timers:\n"); + for(int i=0; i. + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using std::cout; +using std::endl; +using std::min; + +static bool debug = false; + +IntersectionCut::IntersectionCut(Point _f, int _n, const Line *l) : + f(_f), n_faces(_n), pre_lifting_ready(false) +{ + d = new Point[n_faces]; + if (l != 0) + for (int i = 0; i < n_faces; i++) + set_face(i, l[i]); +} + +IntersectionCut::~IntersectionCut() +{ + delete[] d; +} + +void IntersectionCut::set_face(int index, Line line) +{ + if (index < 0 || index >= n_faces) + throw std::out_of_range(""); + + rational x1 = line.p1.x; + rational y1 = line.p1.y; + rational x2 = line.p2.x; + rational y2 = line.p2.y; + + rational h[3] = { -y2 + y1, x2 - x1, x1 * y2 - x2 * y1 }; + rational rhs = -(h[2] + h[0] * f.x + h[1] * f.y); + + if (rhs == rational(0)) + throw std::invalid_argument(""); + + d[index].x = h[0] / rhs; + d[index].y = h[1] / rhs; +} + +rational IntersectionCut::get_continuous_coefficient(rational rx, rational ry) +{ + #ifdef ENABLE_EXTENDED_STATISTICS + Stats::add_coefficient(false); + #endif + + rational max_coeff(-100000, 1); + for (int i = 0; i < n_faces; i++) + { + rational coeff = rx * d[i].x + ry * d[i].y; + if (coeff > max_coeff) + max_coeff = coeff; + } + return max_coeff; +} + +void IntersectionCut::pre_lifting() +{ + rational r0y; + rational fx = f.x, fy = f.y; + rational d0x = d[0].x, d0y = d[0].y; + rational d1x = d[1].x, d1y = d[1].y; + rational apex_x, apex_y; + rational w1 = rational(1) + (d0x * fx + d0y * fy); + rational w2 = rational(1) + (d1x * fx + d1y * fy); + + rational m = (d0x*d1y-d0y*d1x); + + if (m == 0) + { + r0x = -d0y/d0x; + r0y = 1; + } + else + { + + apex_x = (d1y * w1 - d0y * w2) / m; + apex_y = (-d1x * w1 + d0x * w2) / m; + + r0x = apex_x - fx; + r0y = apex_y - fy; + + r0x /= r0y; + r0y = 1; + } + + assert(d0x < 0); + assert(d1x > 0); + + p = r0x*d0x + r0y*d0y; + assert(p == r0x*d1x + r0y*d1y); + assert(p >= 0); + + this->d_p = p.get_double(); + this->d_d1x = d1x.get_double(); + this->d_d1y = d1y.get_double(); + this->d_d0x = d0x.get_double(); + this->d_d0y = d0y.get_double(); + this->d_r0x = r0x.get_double(); + + pre_lifting_ready = true; +} + +double IntersectionCut::get_trivial_lifting_coefficient_double(double rx, double ry) +{ + if(!pre_lifting_ready) pre_lifting(); + + #ifdef ENABLE_EXTENDED_STATISTICS + Stats::add_coefficient(true); + #endif + + double a = 100000; + + unsigned long k2 = 0; + unsigned long M = 10000; + + while (k2 < M) + { + double b = (k2*d_r0x-rx); + + double a1 = d_d1x * (rx + ceil(b)) + d_d1y * k2; + double a2 = d_d0x * (rx + floor(b)) + d_d0y * k2; + + assert(rx + ceil(b) + 0.0001 >= d_r0x*k2); + assert(rx + floor(b) - 0.0001 <= d_r0x*k2); + + assert(a1 >= 0); + assert(a2 >= 0); + + if(a1 < a) a = a1; + if(a2 < a) a = a2; + + if(fabs(d_p) < ZERO_CUTOFF) break; + + M = ceil(a/d_p); + + k2++; + } + + + #ifdef ENABLE_EXTENDED_STATISTICS + Stats::add_trivial_lifting_m(k2); + #endif + + assert(a >= 0); + return a; +} + + +rational IntersectionCut::get_trivial_lifting_coefficient(rational rx, rational ry) +{ + if(!pre_lifting_ready) pre_lifting(); + + #ifdef ENABLE_EXTENDED_STATISTICS + Stats::add_coefficient(true); + #endif + + rational a = rational(100000, 1); + + unsigned long k2 = 0; + unsigned long M = 10000; + rational b = -rx; + + while (k2 < M) + { + rational r_k2((signed) k2); + rational a1 = d[1].x*(rx + b.ceil()) + d[1].y * r_k2; + rational a2 = d[0].x*(rx + b.floor()) + d[0].y * r_k2; + + assert(rx + b.ceil() >= r0x*r_k2); + assert(rx + b.floor() <= r0x*r_k2); + + assert(a1 >= 0); + assert(a2 >= 0); + + if(a1 < a) a = a1; + if(a2 < a) a = a2; + + if(p == 0) break; + M = (a/p).ceil().get_double(); + + k2++; + b += r0x; + } + + return a; +} + + +WedgeCut::WedgeCut(Point _f, Point left, Point apex, Point right) : + IntersectionCut(_f, 2) +{ + set_face(0, Line(left, apex)); + set_face(1, Line(apex, right)); +} + +SplitCut::SplitCut(Point _f, Point left, Point right, Point direction) : + IntersectionCut(_f, 2) +{ + set_face(0, Line(left, left + direction)); + set_face(1, Line(right, right + direction)); +} + +WedgeCutGenerator::WedgeCutGenerator(Row &r) : + SingleRowCutGenerator(r), finished(false), + f(2), r1(2), + r1_offset(-1), + cur_facet(-1), + n_knapsacks(0) +{ + f[0] = r.c.pi_zero.frac(); + f[1] = 0; + cut = new Constraint; + + eval_next(); +} + +WedgeCutGenerator::~WedgeCutGenerator() +{ + delete cut; +} + +bool WedgeCutGenerator::has_next() +{ + return !finished; +} + +Constraint* WedgeCutGenerator::next() +{ + if (!has_next()) + throw std::out_of_range(""); + + Constraint *old_cut = cut; + + cut = new Constraint; + eval_next(); + + return(old_cut); +} + +q::dvec WedgeCutGenerator::intersection( + q::dvec a, q::dvec b, q::dvec c, q::dvec d) +{ + q::dmat g(2,2); + + g.set_col(0, b - a); + g.set_col(1, c - d); + + q::dvec mult = g.inv() * (c - a); + + q::dvec x = a + (b - a) * mult[0]; + + assert(x == c + (d - c) * mult[1]); + + return(x); +} + +void WedgeCutGenerator::eval_next() +{ + while (true) + { + if (0 <= cur_facet && cur_facet < (int) knapsack.list.size() && cur_facet <= MAX_CUT_DEPTH) + break; + + r1_offset++; + + if (r1_offset >= row.c.pi.nz()) + { + finished = true; + return; + } + + int r1_index = row.c.pi.index(r1_offset); + r1[0] = -row.c.pi.value(r1_offset).reduce(REDUCE_FACTOR_R1); + r1[1] = 1; + + if (r1_index == row.basic_var_index) + continue; + + if (!row.is_integer[r1_index]) + continue; + + if (r1[0] == 0) + continue; + + knapsack.clear(); + knapsack.eval(f[0], r1[0]); + + n_knapsacks++; + cur_facet = 0; + } + + Line lines[2]; + int side = knapsack.list[cur_facet].side; + q::dvec &a1 = knapsack.list[cur_facet].lower; + q::dvec &a2 = knapsack.list[cur_facet].upper; + q::dvec &o = knapsack.list[cur_facet].opposed; + + if (side == KNAPSACK2_RAY) + { + a2 = o; + + q::dvec b1 = a1 + r1; + q::dvec b2 = a2 + r1; + + lines[0] = Line(a1[0], a1[1], b1[0], b1[1]); + lines[1] = Line(a2[0], a2[1], b2[0], b2[1]); + } + else + { + q::dvec apex = intersection(a1, a2, f, f + r1); + + if(side == KNAPSACK2_RIGHT) + { + lines[1] = Line(a1[0], a1[1], a2[0], a2[1]); + lines[0] = Line(o[0], o[1], apex[0], apex[1]); + } + else + { + lines[0] = Line(a1[0], a1[1], a2[0], a2[1]); + lines[1] = Line(o[0], o[1], apex[0], apex[1]); + } + } + + if(debug) + cout << endl << lines[0] << " " << lines[1] << endl; + + cur_facet++; + + IntersectionCut ic(Point(f[0], f[1]), 2, lines); + cut->pi.clear(); + cut->pi_zero = -1; + cut->pi.resize(row.c.pi.size()); + cut->depth = cur_facet-1; + + Point ray; + rational alpha, rx, ry; + + for (int l = 0; l < row.c.pi.nz(); l++) + { + int j = row.c.pi.index(l); + if (j == row.basic_var_index) + continue; + + rx = -row.c.pi.value(l); + ry = (l == r1_offset ? 1 : 0); + + if (row.is_integer[j] && l != r1_offset) + { + #ifdef INTERSECTION_CUT_USE_DOUBLE + cut->pi.push(j, -ic.get_trivial_lifting_coefficient_double( + rx.get_double(), ry.get_double())); + #else + cut->pi.push(j, -ic.get_trivial_lifting_coefficient(rx, ry)); + #endif + } + else + { + cut->pi.push(j, -ic.get_continuous_coefficient(rx, ry)); + } + + if(debug) + cout << "r" << j << ": (" << rx << ", " << ry << ") " << cut->pi[j] << endl; + } +} diff --git a/onerow/library/tests/knapsack2_test.cpp b/onerow/library/tests/knapsack2_test.cpp new file mode 100644 index 0000000..e7dcb3f --- /dev/null +++ b/onerow/library/tests/knapsack2_test.cpp @@ -0,0 +1,49 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +#include +typedef q::mpq rational; + +using std::cout; +using std::endl; + +TEST(Knapsack2Test, test_1) +{ + Knapsack2 k(rational(5,7), rational(3,5)); + +// cout << "list:\n"; +// for(auto v : k.list) +// cout << v.lower << " " << v.upper << " " << v.side << " " << v.opposed << endl; + +// cout << "has ray? " << k.has_ray << endl; +} + +TEST(Knapsack2Test, test_2) +{ + Knapsack2 k(rational(1,2), rational(1,2)); + +// cout << "left:\n"; +// each(k.left, v) +// cout << v->active << v->opposed << endl; +// +// cout << "right:\n"; +// each(k.right, v) +// cout << v->active << v->opposed << endl; +// +// cout << "has ray? " << k.has_ray << endl; +} diff --git a/onerow/library/tests/rational_test.cpp b/onerow/library/tests/rational_test.cpp new file mode 100644 index 0000000..9e43098 --- /dev/null +++ b/onerow/library/tests/rational_test.cpp @@ -0,0 +1,66 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include +typedef q::mpq rational; + +TEST(RationalTest, reduce_test) +{ + EXPECT_EQ(rational(5), rational(5).reduce(100)); + EXPECT_EQ(rational(0), rational(0).reduce(100)); + EXPECT_EQ(rational(1,3), rational(1/3.0).reduce(100)); + EXPECT_EQ(rational(1,2), rational(1,2).reduce(100)); + EXPECT_EQ(rational(-1,2), rational(-1,2).reduce(100)); +} + +TEST(RationalTest, minus_test) +{ + EXPECT_EQ(rational(-1,3), -rational(1,3)); +} + +TEST(RationalTest, plus_minus_times_test) +{ + rational a(47,17); + rational b(113,98); + + EXPECT_EQ(rational(6527,1666), a+b); + EXPECT_EQ(rational(2685,1666), a-b); + EXPECT_EQ(rational(-2685,1666), -a+b); + EXPECT_EQ(rational(5311,1666), a*b); +} + +TEST(RationalTest, floor_ceil_frac_test) +{ + EXPECT_EQ(rational(45,17).floor(), rational(2)); + EXPECT_EQ(rational(45,17).ceil(), rational(3)); + EXPECT_EQ(rational(45,17).frac(), rational(11,17)); + + EXPECT_EQ(rational(-45,17).floor(), rational(-3)); + EXPECT_EQ(rational(-45,17).ceil(), rational(-2)); + EXPECT_EQ(rational(-45,17).frac(), rational(6,17)); + + EXPECT_EQ(rational(1,2).frac(), rational(1,2)); + EXPECT_EQ(rational(-1,2).frac(), rational(1,2)); +} + +TEST(RationalTest, double_to_rational_test) +{ + EXPECT_FLOAT_EQ(rational(1/3.0).get_double(), 1/3.0); + EXPECT_FLOAT_EQ(rational(127/15.0).get_double(), 127/15.0); + EXPECT_FLOAT_EQ(rational(127/183.0).get_double(), 127/183.0); + EXPECT_FLOAT_EQ(rational(513/3577.0).get_double(), 513/3577.0); +} diff --git a/onerow/library/tests/single_row_generator_test.cpp b/onerow/library/tests/single_row_generator_test.cpp new file mode 100644 index 0000000..1875e63 --- /dev/null +++ b/onerow/library/tests/single_row_generator_test.cpp @@ -0,0 +1,19 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include + diff --git a/onerow/library/tests/wedge_cut_generator_test.cpp b/onerow/library/tests/wedge_cut_generator_test.cpp new file mode 100644 index 0000000..05453ba --- /dev/null +++ b/onerow/library/tests/wedge_cut_generator_test.cpp @@ -0,0 +1,267 @@ +/* + * 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 3 of the License, or + * (at your option) any later version. + * + * 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. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include +#include + +using std::set; +using std::cout; +using std::endl; + +TEST(IntersectionCutTest, set_face_test) +{ + IntersectionCut ic(Point(rational(1,2),rational(1,4)), 1); + + // should calculate d correctly + ic.set_face(0, Line(rational(127,64), rational(0), rational(67,64), rational(1))); + EXPECT_EQ(ic.d[0].x, rational(4,5)); + EXPECT_EQ(ic.d[0].y, rational(3,4)); + + // shoud work with any sequence of points + ic.set_face(0, Line(rational(67,64), rational(1), rational(127,64), rational(0))); + EXPECT_EQ(ic.d[0].x, rational(4,5)); + EXPECT_EQ(ic.d[0].y, rational(3,4)); + + + // should throw exception on badly defined lines + EXPECT_THROW(ic.set_face(0, Line(rational(1), rational(1), rational(1), rational(1))), std::invalid_argument); + + // should throw exception on lines crossing f + EXPECT_THROW(ic.set_face(0, Line(rational(1,2), rational(1,4), rational(1), rational(1))), std::invalid_argument); + + // should not throw exception on good lines + EXPECT_NO_THROW(ic.set_face(0, Line(rational(0), rational(0), rational(1), rational(1)))); + + // should throw exception on index out of range + EXPECT_THROW(ic.set_face(-1, Line()), std::out_of_range); + EXPECT_THROW(ic.set_face( 1, Line()), std::out_of_range); + EXPECT_THROW(ic.set_face(50, Line()), std::out_of_range); +} + +TEST(IntersectionCutTest, set_face_test_2) +{ + IntersectionCut ic(Point(rational(1,2),0), 1); + + ic.set_face(0, Line(1, 0, 0, 1)); + EXPECT_EQ(ic.d[0].x, rational(2)); + EXPECT_EQ(ic.d[0].y, rational(2)); + + ic.set_face(0, Line(0, 1, 1, 0)); + EXPECT_EQ(ic.d[0].x, rational(2)); + EXPECT_EQ(ic.d[0].y, rational(2)); +} + +TEST(IntersectionCutTest, get_continuous_coefficient_test) +{ + IntersectionCut ic(Point(rational(5,7),rational(0)), 2); + ic.set_face(0, Line(3, 4, rational(29,7), rational(40,7))); + ic.set_face(1, Line(rational(29,7), rational(40,7), rational(2), rational(2))); + + // should calculate d correctly + EXPECT_EQ(rational(-21, 8), ic.d[0].x); + EXPECT_EQ(rational( 7, 4), ic.d[0].y); + EXPECT_EQ(rational( 91,12), ic.d[1].x); + EXPECT_EQ(rational(-35, 8), ic.d[1].y); + + // should calculate continuous coefficients correctly + EXPECT_EQ(rational(91,12), ic.get_continuous_coefficient(rational(1),rational(0))); + EXPECT_EQ(rational( 7, 4), ic.get_continuous_coefficient(rational(0),rational(1))); + EXPECT_EQ(rational( 7,40), ic.get_continuous_coefficient(rational(3,5),1)); + EXPECT_EQ(rational(77,72), ic.get_continuous_coefficient(rational(1,3),rational(1,3))); +} + +TEST(IntersectionCutTest, get_trivial_lifting_coefficient_test) +{ + IntersectionCut ic(Point(rational(5,7),rational(0)), 2); + ic.set_face(0, Line(3, 4, rational(29,7), rational(40,7))); + ic.set_face(1, Line(rational(29,7), rational(40,7), rational(2), rational(2))); + + // should calculate d correctly + EXPECT_EQ(rational(-21, 8), ic.d[0].x); + EXPECT_EQ(rational( 7, 4), ic.d[0].y); + EXPECT_EQ(rational( 91,12), ic.d[1].x); + EXPECT_EQ(rational(-35, 8), ic.d[1].y); + + // should calculate lifted coefficients correctly + EXPECT_EQ(rational(399,704), ic.get_trivial_lifting_coefficient(rational(119,264),rational(0))); +} + +TEST(IntersectionCutTest, get_trivial_lifting_coefficient_test_2) +{ + IntersectionCut ic(Point(rational(1,2),rational(0)), 2); + ic.set_face(0, Line(0,0,0,1)); + ic.set_face(1, Line(1,0,0,1)); + + // should calculate d correctly + EXPECT_EQ(rational(-2), ic.d[0].x); + EXPECT_EQ( rational(0), ic.d[0].y); + EXPECT_EQ( rational(2), ic.d[1].x); + EXPECT_EQ( rational(2), ic.d[1].y); + + // should calculate lifted coefficients correctly + EXPECT_EQ(rational(1), ic.get_trivial_lifting_coefficient(rational(-1,2),0)); +} + +TEST(WedgeCutGenerator, generate_test_1) +{ + Row r; + r.basic_var_index = 999; + r.is_integer = new bool[3]; + r.c.pi.resize(3); + r.c.pi.push(0, rational(-3,5)); + r.c.pi.push(1, rational(-1)); + r.c.pi.push(2, rational(1)); + r.c.pi_zero = rational(5,7); + + std::set expected_cuts; + rational cut_matrix[] = { + rational(-14,25), rational(-7,2), rational(-7,5), + rational(-7,20), rational(-7,2), rational(-91,44), + rational(-7,40), rational(-91,12), rational(-21,8), + rational(0), rational(-35,3), rational(-35,4), +// rational(-21,10), rational(0), rational(-7,5), +// rational(-21,10), rational(-7,2), rational(0), + }; + + for(int i=0; i<3; i++) + r.is_integer[i] = (i == 0); + + for(int i=0; i<4; i++) + { + Constraint c; + c.pi.resize(3); + c.pi_zero = -1; + for(int j=0; j<3; j++) + c.pi.push(j, cut_matrix[i*3+j]); + expected_cuts.insert(c); + } + + std::set actual_cuts; + + WedgeCutGenerator wcg(r); + while(wcg.has_next()) + actual_cuts.insert(*wcg.next()); + +// cout << "EXPECTED:" << endl; +// for(Constraint c : expected_cuts) +// cout << c << endl; +// +// cout << "ACTUAL:" << endl; +// for(Constraint c : actual_cuts) +// cout << c << endl; + + EXPECT_TRUE(expected_cuts == actual_cuts); +} + +TEST(WedgeCutGenerator, generate_test_2) +{ + Row r; + r.basic_var_index = 999; + r.is_integer = new bool[3]; + r.c.pi.resize(3); + r.c.pi.push(0, rational(1,2)); + r.c.pi.push(1, rational(-1)); + r.c.pi.push(2, rational(1)); + r.c.pi_zero = rational(1,2); + + set expected_cuts; + rational cut_matrix[] = { + rational(-1), rational(-2), rational(-2) + }; + + for(int i=0; i<5; i++) + r.is_integer[i] = (i == 0); + + for(int i=0; i<1; i++) + { + Constraint c; + c.pi.resize(3); + c.pi_zero = -1; + for(int j=0; j<3; j++) + { + c.pi.push(j, cut_matrix[i*3+j]); + } + expected_cuts.insert(c); + } + + set actual_cuts; + WedgeCutGenerator wcg(r); + while(wcg.has_next()) + actual_cuts.insert(*wcg.next()); + +// cout << "EXPECTED:" << endl; +// for(Constraint c : expected_cuts) +// cout << c << endl; +// +// cout << "ACTUAL:" << endl; +// for(Constraint c : actual_cuts) +// cout << c << endl; + + EXPECT_TRUE(expected_cuts == actual_cuts); +} + +TEST(WedgeCutGenerator, generate_test_3) +{ + bool is_integer[] = { true, false, false, true, true, true, true, true, false, + false, false, false, false, false }; + Row r; + r.basic_var_index = 0; + r.is_integer = is_integer; + r.c.pi.resize(14); + r.c.pi.push(0, rational(1)); + r.c.pi.push(1, rational(7338,415411)); + r.c.pi.push(2, rational(1271,71708)); + r.c.pi.push(3, rational(19,577624)); + r.c.pi.push(4, rational(-19,577624)); + r.c.pi.push(5, rational(7,189263)); + r.c.pi.push(6, rational(-7,189263)); + r.c.pi.push(7, rational(-5,83513)); + r.c.pi.push(8, rational(157727,47951)); + r.c.pi.push(9, rational(-157727,47951)); + r.c.pi.push(10, rational(536435,145039)); + r.c.pi.push(11, rational(-536435,145039)); + r.c.pi.push(12, rational(237478,39665)); + r.c.pi.push(13, rational(237478,39665)); + r.c.pi_zero = rational(631975,707651); + + rational solution[14]; + solution[0] = rational(0); + solution[1] = rational(0); + solution[2] = rational(0); + solution[3] = rational(0); + solution[4] = rational(0); + solution[5] = rational(0); + solution[6] = rational(0); + solution[7] = rational(1); + solution[8] = rational(5531144031239411,576460752303423488); + solution[9] = rational(8350617490688559,1152921504606846976); + solution[10] = rational(6818955564856927,576460752303423488); + solution[11] = rational(6377917030130699,4611686018427387904); + solution[12] = rational(2547758494699757,18014398509481984); + solution[13] = rational(0); + + cout << r.c.get_violation(solution).get_double() << endl; + EXPECT_TRUE(r.c.get_violation(solution) < 0.001); + + WedgeCutGenerator wcg(r); + while(wcg.has_next()) + { + Constraint *c = wcg.next(); + cout << *c << endl; + cout << c->get_violation(solution).get_double() << endl; + EXPECT_TRUE(c->get_violation(solution) < 0.001); + } +} diff --git a/qxx/CMakeLists.txt b/qxx/CMakeLists.txt new file mode 100644 index 0000000..c9aad2d --- /dev/null +++ b/qxx/CMakeLists.txt @@ -0,0 +1,21 @@ +set(COMMON_SOURCES + src/dlu.cpp + src/dmat.cpp + src/dvec.cpp + src/rational.cpp + src/slu.cpp + src/smat.cpp + src/strprintf.cpp + src/svec.cpp + include/qxx/dlu.hpp + include/qxx/dmat.hpp + include/qxx/dvec.hpp + include/qxx/rational.hpp + include/qxx/slu.hpp + include/qxx/smat.hpp + include/qxx/strprintf.hpp + include/qxx/svec.hpp) + +add_library(qxx_static ${COMMON_SOURCES}) +set_target_properties(qxx_static PROPERTIES OUTPUT_NAME qxx) +target_include_directories (qxx_static PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/qxx/include/qxx/debug.hpp b/qxx/include/qxx/debug.hpp new file mode 100644 index 0000000..28a9c47 --- /dev/null +++ b/qxx/include/qxx/debug.hpp @@ -0,0 +1,54 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_DEBUG_H +#define QXX_DEBUG_H +#include +#include + +// ************************************************************************** +// +// ************************************************************************** +#define Q_RANGE_CHECK(i, lb, ub) \ + do { \ + if (((i) < (lb)) || ((i) > (ub))) { \ + fprintf(stderr, "%s:%d: in %s(): range check failed " \ + "(%s not in %s..%s | %d not in %d..%d)\n", \ + __FILE__, __LINE__, __func__, \ + #i, #lb, #ub, (i), (lb), (ub)); \ + assert(!"range check"); \ + } \ + } while(0) + +#define Q_MATCH(a, b) \ + do { \ + if ((a) != (b)) { \ + fprintf(stderr, "%s:%d: in %s(): match failed " \ + "(%s != %s | %d != %d)\n", \ + __FILE__, __LINE__, __func__, \ + #a, #b, (a), (b)); \ + assert(!"match"); \ + } \ + } while(0) + + +// ************************************************************************** +// +// ************************************************************************** + +#endif + + diff --git a/qxx/include/qxx/dlu.hpp b/qxx/include/qxx/dlu.hpp new file mode 100644 index 0000000..9e215c6 --- /dev/null +++ b/qxx/include/qxx/dlu.hpp @@ -0,0 +1,60 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_DENSE_LU_HPP +#define QXX_DENSE_LU_HPP +#include "rational.hpp" +#include "dmat.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +class dlu { + +public: + dlu(const dmat &a); + ~dlu(); + + dmat L() const; + dmat U() const; + dvec solve_Ax(const dvec &b) const; + dvec solve_xA(const dvec &b) const; + mpq det() const; + +private: + void pivot(int k); + void factorize(const dmat &a); + + int n; + std::vector row_bwd; + dmat g; + +}; + +// ************************************************************************** +// +// ************************************************************************** +}; + +#endif + + diff --git a/qxx/include/qxx/dmat.hpp b/qxx/include/qxx/dmat.hpp new file mode 100644 index 0000000..63fff78 --- /dev/null +++ b/qxx/include/qxx/dmat.hpp @@ -0,0 +1,104 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_DMAT_HPP +#define QXX_DMAT_HPP +#include +#include "rational.hpp" +#include "dvec.hpp" + + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +// ************************************************************************** +// +// ************************************************************************** +class dmat { + +public: + dmat(); + dmat(const dmat &b); + dmat(int m, int n); + + int rows() const; + int cols() const; + void resize(int m, int n); + void clear(); + + const mpq &get(int i, int j) const; + void set(int i, int j, const mpq &v); + void set(const mpq &v); + void set_identity(); + + dvec get_col(int j) const; + void set_col(int j, const dvec &v); + const dvec &get_row(int i) const; + void set_row(int i, const dvec &v); + + dvec &operator[] (int i); + const dvec &operator[] (int i) const; + mpq &operator() (int i, int j); + const mpq &operator() (int i, int j) const; + + const dmat &operator+() const; + dmat operator-() const; + dmat operator+(const dmat &b) const; + dmat operator-(const dmat &b) const; + dmat operator*(const dmat &b) const; + dvec operator*(const dvec &v) const; + dmat operator*(const mpq &v) const; + dmat operator/(const mpq &v) const; + dmat inv() const; + dmat t() const; + mpq det() const; + + bool operator==(const dmat &b) const; + bool operator!=(const dmat &b) const; + + dmat &operator=(const dmat &b); + dmat &operator=(const dvec &v); + + dmat &operator+=(const dmat &b); + dmat &operator-=(const dmat &b); + dmat &operator*=(const dmat &b); + dmat &operator*=(const mpq &v); + dmat &operator/=(const mpq &v); + + void fdump(FILE *f, const std::string &name = "") const; + void dump(const std::string &name = "") const; + +private: + int m, n; + std::vector d; + +}; + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const dmat &a); + + +// ************************************************************************** +// +// ************************************************************************** +} + +#endif + diff --git a/qxx/include/qxx/dvec.hpp b/qxx/include/qxx/dvec.hpp new file mode 100644 index 0000000..517dad2 --- /dev/null +++ b/qxx/include/qxx/dvec.hpp @@ -0,0 +1,97 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_DENSE_VECTOR_HPP +#define QXX_DENSE_VECTOR_HPP +#include +#include +#include "rational.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +class svec; + +// ************************************************************************** +// +// ************************************************************************** +class dvec { + +public: + dvec(); + dvec(const dvec &b); + dvec(int n); + dvec(const svec &b); + + int size() const; + void resize(int n); + void clear(); + mpq &operator[] (int i); + const mpq &operator[] (int i) const; + + void set(const mpq &v); + void set(int i0, int i1, const mpq &v); + void assign(int n, const mpq &v); + + dvec operator+(const dvec &b) const; + dvec operator-(const dvec &b) const; + mpq operator*(const dvec &b) const; + mpq operator*(const svec &b) const; + + const dvec &operator+() const; + dvec operator-() const; + dvec operator+(const mpq &v) const; + dvec operator-(const mpq &v) const; + dvec operator*(const mpq &v) const; + dvec operator/(const mpq &v) const; + + bool operator==(const dvec &b) const; + bool operator!=(const dvec &b) const; + + dvec &operator=(const dvec &b); + dvec &operator=(const svec &b); + + dvec &operator+=(const dvec &b); + dvec &operator-=(const dvec &b); + dvec &operator+=(const mpq &v); + dvec &operator-=(const mpq &v); + dvec &operator*=(const mpq &v); + dvec &operator/=(const mpq &v); + + void fdump(FILE *f, const std::string &name = "") const; + void dump(const std::string &name = "") const; + +private: + std::vector d; + +}; + + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const dvec &v); + +// ************************************************************************** +// +// ************************************************************************** +} + +#endif + + diff --git a/qxx/include/qxx/rational.hpp b/qxx/include/qxx/rational.hpp new file mode 100644 index 0000000..7520697 --- /dev/null +++ b/qxx/include/qxx/rational.hpp @@ -0,0 +1,114 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_TYPES_H +#define QXX_TYPES_H +#include +#include +#include + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +// ************************************************************************** +// +// ************************************************************************** +class mpq { + +public: + mpq(); + mpq(const mpq &b); + mpq(long num, long den); + mpq(int num, int den); + mpq(long i); + mpq(int i); + mpq(double d); + mpq(const char *s); + ~mpq(); + + const mpq &operator+() const; + mpq operator-() const; + mpq operator+(const mpq &b) const; + mpq operator-(const mpq &b) const; + mpq operator*(const mpq &b) const; + mpq operator/(const mpq &b) const; + + mpq &operator=(const mpq &b); + mpq &operator=(int i); + mpq &operator+=(const mpq &b); + mpq &operator-=(const mpq &b); + mpq &operator*=(const mpq &b); + mpq &operator/=(const mpq &b); + + bool operator<(const mpq &b) const; + bool operator>(const mpq &b) const; + bool operator<=(const mpq &b) const; + bool operator>=(const mpq &b) const; + bool operator==(const mpq &b) const; + bool operator!=(const mpq &b) const; + + void set_neg(); + void set_neg(const mpq &b); + void set_inv(); + void set_inv(const mpq &b); + void set_abs(); + void set_abs(const mpq &b); + + mpq neg() const; + mpq inv() const; + mpq abs() const; + int sign() const; + + mpq num() const; + mpq den() const; + int bits() const; + + mpq frac() const; + mpq floor() const; + mpq trunc() const; + mpq ceil() const; + + long get_long_num() const; + long get_long_den() const; + double get_double() const; + mpq reduce(mpq max_den) const; + + void fdump(FILE *f, const std::string &name = "") const; + void dump(const std::string &name = "") const; + +public: + mpq_t v; + +}; + + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const mpq &v); + + + +// ************************************************************************** +// +// ************************************************************************** +} + +#endif + + diff --git a/qxx/include/qxx/slu.hpp b/qxx/include/qxx/slu.hpp new file mode 100644 index 0000000..28cc127 --- /dev/null +++ b/qxx/include/qxx/slu.hpp @@ -0,0 +1,86 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_SLU_H +#define QXX_SLU_H +#include +#include "qxx/smat.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +class perm { + +public: + perm(); + ~perm(); + + int size() const; + void resize(int n); + void clear(); + + void id(); + void id(int n); + + void pivot(int k, int i); + +public: + std::vector fwd; + std::vector bwd; + +}; + +// ************************************************************************** +// +// ************************************************************************** +class slu { + +public: + slu(const smat &a); + ~slu(); + + dvec solve_Ax(const dvec &b) const; + dvec solve_xA(const dvec &b) const; + mpq det() const; + +public: + void pivot(int k); + void factorize(const smat &a); + dvec solve_gen(const smat &l, const smat &u, + const perm &row, const perm &col, const dvec &b) const; + + int n; + + std::vector cnz; + perm row, col; + smat w, lt, u; + smat l, ut; +}; + +// ************************************************************************** +// +// ************************************************************************** +}; + +#endif + + diff --git a/qxx/include/qxx/smat.hpp b/qxx/include/qxx/smat.hpp new file mode 100644 index 0000000..440e100 --- /dev/null +++ b/qxx/include/qxx/smat.hpp @@ -0,0 +1,97 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_SMAT_HPP +#define QXX_SMAT_HPP +#include +#include "rational.hpp" +#include "dvec.hpp" +#include "dmat.hpp" +#include "svec.hpp" + + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +// ************************************************************************** +// +// ************************************************************************** +class smat { + +public: + smat(); + smat(const smat &b); + smat(int m, int n); + + int rows() const; + int cols() const; + void resize(int m, int n); + void clear(); + + const mpq &get(int i, int j) const; + void set(int i, int j, const mpq &v); + void set_zero(); + void set_identity(); + + void gather(const dmat &src); + void spread(dmat &r) const; + dmat dense() const; + + svec get_col(int j) const; + void set_col(int j, const svec &v); + const svec &get_row(int i) const; + void set_row(int i, const svec &v); + + svec &operator[] (int i); + const svec &operator[] (int i) const; + svec_ref operator() (int i, int j); + const mpq &operator() (int i, int j) const; + + smat &operator=(const smat &b); + + smat operator*(const smat &b) const; + dvec operator*(const dvec &v) const; + svec operator*(const svec &v) const; + smat inv() const; + void set_transpose(const smat &a); + smat t() const; + mpq det() const; + + void fdump(FILE *f, const std::string &name = "") const; + void dump(const std::string &name = "") const; + +private: + int m, n; + std::vector d; + +}; + + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const smat &a); + + +// ************************************************************************** +// +// ************************************************************************** +} + +#endif + diff --git a/qxx/include/qxx/strprintf.hpp b/qxx/include/qxx/strprintf.hpp new file mode 100644 index 0000000..1cecf1c --- /dev/null +++ b/qxx/include/qxx/strprintf.hpp @@ -0,0 +1,31 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_STRPRINTF_HPP +#define QXX_STRPRINTF_HPP +#include + +namespace q { + + +void format(std::string &s, const char *fmt, ...); +std::string strprintf(const char *fmt, ...); + + +} // namespace + +#endif + diff --git a/qxx/include/qxx/svec.hpp b/qxx/include/qxx/svec.hpp new file mode 100644 index 0000000..595be2b --- /dev/null +++ b/qxx/include/qxx/svec.hpp @@ -0,0 +1,178 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#ifndef QXX_SPARSE_VECTOR_HPP +#define QXX_SPARSE_VECTOR_HPP +#include +#include "rational.hpp" +#include "dvec.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +class sparse_el; +class svec_ref; +class svec; + +// ************************************************************************** +// +// ************************************************************************** +class svec_ref { + +public: + svec_ref(svec &vec, int i); + svec_ref(const svec_ref &src); + ~svec_ref(); + + const mpq &get() const; + + operator mpq() const; + + const mpq &operator=(const svec_ref &src) const; + const mpq &operator=(const mpq &b) const; + const mpq &operator+=(const mpq &b) const; + const mpq &operator-=(const mpq &b) const; + const mpq &operator*=(const mpq &b) const; + const mpq &operator/=(const mpq &b) const; + + const mpq &operator+() const; + mpq operator-() const; + mpq operator+(const mpq &b) const; + mpq operator-(const mpq &b) const; + mpq operator*(const mpq &b) const; + mpq operator/(const mpq &b) const; + + bool operator<(const mpq &b) const; + bool operator>(const mpq &b) const; + bool operator<=(const mpq &b) const; + bool operator>=(const mpq &b) const; + bool operator==(const mpq &b) const; + bool operator!=(const mpq &b) const; + +public: + svec &vector; + int index; +}; + + +// ************************************************************************** +// +// ************************************************************************** +class svec_el { + +public: + svec_el(); + svec_el(int idx, const mpq &v); + svec_el(const svec_el &e); + ~svec_el(); + +public: + int index; + mpq value; +}; + + +// ************************************************************************** +// +// ************************************************************************** +class svec { + +public: + svec(); + svec(const svec &src); + svec(int n); + svec(const dvec &src); + ~svec(); + + void clear(); + + int size() const; + void resize(int n); + int nz() const; + + bool offs_active() const; + void offs_build(); + void offs_clear(); + + void set_zero(); + + void gather(const dvec &src); + void spread(dvec &r) const; + dvec dense() const; + + int locate(int idx) const; + void remove(int offset); + + int &index(int offset); + mpq &value(int offset); + int index(int offset) const; + const mpq &value(int offset) const; + + void push_nz(int idx, const mpq &v); + void push(int idx, const mpq &v); + + const mpq &get(int idx) const; + const mpq &set(int idx, const mpq &v); + + const mpq &operator[](int idx) const; + svec_ref operator[](int idx); + + const svec &operator+() const; + svec operator-() const; + + svec operator+(const svec &b) const; + svec operator-(const svec &b) const; + mpq operator*(const svec &b) const; + mpq operator*(const dvec &b) const; + + svec operator*(const mpq &v); + svec operator/(const mpq &v); + + svec &operator=(const svec &b); + svec &operator+=(const svec &b); + svec &operator-=(const svec &b); + svec &operator*=(const mpq &v); + svec &operator/=(const mpq &v); + + void sort(); + + void fdump(FILE *f, const std::string &name = "") const; + void dump(const std::string &name = "") const; + +private: + int n; + std::vector d; + std::vector offs; + + static const mpq zero; +}; + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const svec &v); + + +// ************************************************************************** +// +// ************************************************************************** +} + +#endif + + diff --git a/qxx/src/dlu.cpp b/qxx/src/dlu.cpp new file mode 100644 index 0000000..31b35c6 --- /dev/null +++ b/qxx/src/dlu.cpp @@ -0,0 +1,240 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include "qxx/dlu.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +dlu::dlu(const dmat &a) + : n(0) +{ + factorize(a); +} + +dlu::~dlu() +{ +} + +// ************************************************************************** +// +// ************************************************************************** +void dlu::pivot(int p) +{ + int best_v = -1; + int best_k = -1; + + //printf("pivot(%d):\n", p); + + for (int k = p; k < n; k++) { + int i = row_bwd[k]; + + //gmp_printf("\tchoice [%d] row %d: %Qd\n", + // k, i, g(i, p).v); + + if (g(i, p).sign()) { + int v = g(i, p).bits(); + + if ((best_v < 0) || (v < best_v)) { + best_k = k; + best_v = v; + } + } + } + + assert((best_k >= 0) && "Singular matrix"); + + int i = row_bwd[best_k]; + row_bwd[best_k] = row_bwd[p]; + row_bwd[p] = i; +} + +// ************************************************************************** +// +// ************************************************************************** +void dlu::factorize(const dmat &a) +{ + // check + assert(a.rows() == a.cols()); + + // init + n = a.rows(); + g = a; + row_bwd.resize(n); + for (int p = 0; p < n; p++) + row_bwd[p] = p; + + // main loop + mpq pv, f; + + for (int p = 0; p < n; p++) { + pivot(p); + + int pi = row_bwd[p]; + pv = g(pi, p); + + for (int k = p + 1; k < n; k++) { + int i = row_bwd[k]; + f = g(i, p) / pv; + + for (int j = p + 1; j < n; j++) + g(i, j) -= f * g(pi, j); + } + + //printf("g%d = ", p); + //g.dump(stdout); + } + + //printf("row_bwd = ["); + //for (int p = 0; p < n; p++) + // printf(" %d", row_bwd[p]); + //printf(" ];\n"); +} + +// ************************************************************************** +// +// ************************************************************************** +dmat dlu::L() const +{ + dmat l(n, n); + + for (int k = 0; k < n; k++) { + for (int j = 0; j < k; j++) + l(k, j) = g(row_bwd[k], j) / g(row_bwd[j], j); + l(k, k) = 1; + } + + return(l); +} + +dmat dlu::U() const +{ + dmat u(n, n); + + for (int k = 0; k < n; k++) { + for (int j = k; j < n; j++) + u(k, j) = g(row_bwd[k], j); + } + + return(u); +} + +// ************************************************************************** +// +// ************************************************************************** +dvec dlu::solve_Ax(const dvec &b) const +{ + assert(b.size() == n); + + dvec x(n); + mpq v; + + int i0 = row_bwd[0]; + x[0] = b[i0] / g(i0, 0); + + for (int k = 1; k < n; k++) { + int i = row_bwd[k]; + v = 0; + + for (int j = 0; j < k; j++) + v += g(i, j) * x[j]; + + x[k] = (b[i] - v) / g(i, k); + } + + for (int k = n - 2; k != -1; k--) { + int i = row_bwd[k]; + v = 0; + + for (int j = k + 1; j < n; j++) + v += g(i, j) * x[j]; + + x[k] -= v / g(i, k); + } + + return(x); +} + +// ************************************************************************** +// +// ************************************************************************** +dvec dlu::solve_xA(const dvec &b) const +{ + assert(b.size() == n); + + dvec x(n); + mpq v; + + int i0 = row_bwd[0]; + x[i0] = b[0] / g(i0, 0); + + for (int k = 1; k < n; k++) { + int i = row_bwd[k]; + v = 0; + + for (int l = 0; l < k; l++) { + int il = row_bwd[l]; + v += g(il, k) * x[il]; + } + + x[i] = (b[k] - v) / g(i, k); + } + + for (int k = n - 2; k != -1; k--) { + int i = row_bwd[k]; + v = 0; + + for (int l = k + 1; l < n; l++) { + int il = row_bwd[l]; + v += g(il, k) * x[il]; + } + + x[i] -= v / g(i, k); + } + + return(x); +} + +// ************************************************************************** +// +// ************************************************************************** +mpq dlu::det() const +{ + mpq v; + + v = 1; + for (int k = 0; k < n; k++) { + int i = row_bwd[k]; + + v *= g(i, k); + } + + return(v); +} + +// ************************************************************************** +// +// ************************************************************************** +}; + + diff --git a/qxx/src/dmat.cpp b/qxx/src/dmat.cpp new file mode 100644 index 0000000..4bd7b8c --- /dev/null +++ b/qxx/src/dmat.cpp @@ -0,0 +1,486 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include "qxx/debug.hpp" +#include "qxx/rational.hpp" +#include "qxx/dvec.hpp" +#include "qxx/dlu.hpp" +#include "qxx/dmat.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +dmat::dmat() + : m(0), n(0) +{ +} + +dmat::dmat(const dmat &b) + : m(b.m), n(b.n), d(b.d) +{ +} + +dmat::dmat(int m0, int n0) + : m(m0), n(n0), d(m0) +{ + for (int i = 0; i < m; i++) + d[i].resize(n); +} + + +int dmat::rows() const +{ + return(m); +} + +int dmat::cols() const +{ + return(n); +} + +void dmat::resize(int m0, int n0) +{ + m = m0; + n = n0; + + d.resize(m); + for (int i = 0; i < m; i++) + d[i].resize(n); +} + +void dmat::clear() +{ + m = n = 0; + d.clear(); +} + +// ************************************************************************** +// +// ************************************************************************** +const mpq &dmat::get(int i, int j) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + +void dmat::set(int i, int j, const mpq &v) +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + d[i][j] = v; +} + +void dmat::set(const mpq &v) +{ + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] = v; +} + +void dmat::set_identity() +{ + mpq zero(0), one(1); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] = (i == j) ? one : zero; +} + + +// ************************************************************************** +// +// ************************************************************************** +dvec dmat::get_col(int j) const +{ + Q_RANGE_CHECK(j, 0, n - 1); + + dvec r(m); + + for (int i = 0; i < m; i++) + r[i] = d[i][j]; + + return(r); +} + +void dmat::set_col(int j, const dvec &v) +{ + Q_RANGE_CHECK(j, 0, n - 1); + Q_MATCH(v.size(), m); + + for (int i = 0; i < m; i++) + d[i][j] = v[i]; +} + +const dvec &dmat::get_row(int i) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + + return(d[i]); +} + +void dmat::set_row(int i, const dvec &v) +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_MATCH(v.size(), n); + + d[i] = v; +} + + +// ************************************************************************** +// +// ************************************************************************** +dvec &dmat::operator[] (int i) +{ + Q_RANGE_CHECK(i, 0, m - 1); + + return(d[i]); +} + +const dvec &dmat::operator[] (int i) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + + return(d[i]); +} + +mpq &dmat::operator() (int i, int j) +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + +const mpq &dmat::operator() (int i, int j) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + + +// ************************************************************************** +// +// ************************************************************************** +const dmat &dmat::operator+() const +{ + return(*this); +} + +dmat dmat::operator-() const +{ + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = -d[i][j]; + + return(r); +} + +dmat dmat::operator+(const dmat &b) const +{ + Q_MATCH(m, b.m); + Q_MATCH(n, b.n); + + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = d[i][j] + b.d[i][j]; + + return(r); +} + +dmat dmat::operator-(const dmat &b) const +{ + Q_MATCH(m, b.m); + Q_MATCH(n, b.n); + + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = d[i][j] - b.d[i][j]; + + return(r); +} + +dmat dmat::operator*(const dmat &b) const +{ + Q_MATCH(n, b.m); + + dmat r(m, b.n); + mpq v; + + for (int i = 0; i < m; i++) { + for (int j = 0; j < b.n; j++) { + v = 0; + + for (int k = 0; k < n; k++) + v += d[i][k] * b.d[k][j]; + + r.d[i][j] = v; + } + } + + return(r); +} + +dvec dmat::operator*(const dvec &b) const +{ + Q_MATCH(n, b.size()); + + dvec r(m); + mpq v; + + for (int i = 0; i < m; i++) { + v = 0; + for (int j = 0; j < n; j++) + v += d[i][j] * b[j]; + r[i] = v; + } + + return(r); +} + +dmat dmat::operator*(const mpq &v) const +{ + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = d[i][j] * v; + + return(r); +} + +dmat dmat::operator/(const mpq &v) const +{ + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = d[i][j] / v; + + return(r); +} + + +// ************************************************************************** +// +// ************************************************************************** +dmat dmat::inv() const +{ + Q_MATCH(m, n); + + dvec e(n); + dmat r(n, n); + dlu lu(*this); + + for (int i = 0; i < n; i++) { + e[i] = 1; + + r.set_col(i, lu.solve_Ax(e)); + + //e.dump(stdout, "e"); + //r.get_col(i).dump(stdout, "x"); + + e[i] = 0; + } + + return(r); +} + +dmat dmat::t() const +{ + dmat r(m, n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + r.d[i][j] = d[j][i]; + + return(r); +} + +mpq dmat::det() const +{ + Q_MATCH(m, n); + + if (n == 1) + return(d[0][0]); + + if (n == 2) + return(d[0][0] * d[1][1] - d[1][0] * d[0][1]); + + if (n == 3) + return( d[0][0] * d[1][1] * d[2][2] + + d[0][1] * d[1][2] * d[2][0] + + d[0][2] * d[1][0] * d[2][1] + - d[2][0] * d[1][1] * d[0][2] + - d[2][1] * d[1][2] * d[0][0] + - d[2][2] * d[1][0] * d[0][1] ); + + dlu lu(*this); + return(lu.det()); +} + +// ************************************************************************** +// +// ************************************************************************** +bool dmat::operator==(const dmat &b) const +{ + Q_MATCH(m, b.m); + Q_MATCH(n, b.n); + + for (int i = 0; i < m; i++) { + if (d[i] != b.d[i]) + return(false); + } + + return(true); +} + +bool dmat::operator!=(const dmat &b) const +{ + return(!((*this) == b)); +} + +// ************************************************************************** +// +// ************************************************************************** +dmat &dmat::operator=(const dmat &b) +{ + m = b.m; + n = b.n; + d = b.d; + return(*this); +} + +dmat &dmat::operator=(const dvec &v) +{ + resize(v.size(), 1); + + for (int i = 0; i < m; i++) + d[i][0] = v[i]; + + return(*this); +} + + +dmat &dmat::operator+=(const dmat &b) +{ + Q_MATCH(m, b.m); + Q_MATCH(n, b.n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] += b.d[i][j]; + + return(*this); +} + +dmat &dmat::operator-=(const dmat &b) +{ + Q_MATCH(m, b.m); + Q_MATCH(n, b.n); + + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] -= b.d[i][j]; + + return(*this); +} + +dmat &dmat::operator*=(const dmat &b) +{ + dmat r; + + r = (*this) * b; + (*this) = r; + return(*this); +} + +dmat &dmat::operator*=(const mpq &v) +{ + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] *= v; + + return(*this); +} + +dmat &dmat::operator/=(const mpq &v) +{ + for (int i = 0; i < m; i++) + for (int j = 0; j < n; j++) + d[i][j] /= v; + + return(*this); +} + + +// ************************************************************************** +// +// ************************************************************************** +void dmat::fdump(FILE *f, const std::string &name) const +{ + if (name.length()) + gmp_fprintf(f, "%s = [\n", name.c_str()); + else + gmp_fprintf(f, "[\n"); + + for (int i = 0; i < m; i++) + d[i].fdump(f); + fprintf(f, "];\n"); +} + +void dmat::dump(const std::string &name) const +{ + fdump(stdout, name); +} + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const dmat &a) +{ + for (int i = 0; i < a.rows(); i++) + os << a[i] << std::endl; + + return(os); +} + +// ************************************************************************** +// +// ************************************************************************** +} + + + + diff --git a/qxx/src/dvec.cpp b/qxx/src/dvec.cpp new file mode 100644 index 0000000..badff22 --- /dev/null +++ b/qxx/src/dvec.cpp @@ -0,0 +1,366 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include "qxx/debug.hpp" +#include "qxx/rational.hpp" +#include "qxx/dvec.hpp" +#include "qxx/svec.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +dvec::dvec() +{ +} + +dvec::dvec(const dvec &b) + : d(b.d) +{ +} + +dvec::dvec(int n) + : d(n) +{ +} + +dvec::dvec(const svec &b) +{ + b.spread(*this); +} + + +// ************************************************************************** +// +// ************************************************************************** +int dvec::size() const +{ + return(d.size()); +} + +void dvec::resize(int n) +{ + d.resize(n); +} + +void dvec::clear() +{ + d.clear(); +} + +mpq &dvec::operator[] (int i) +{ + Q_RANGE_CHECK(i, 0, size() - 1); + return(d[i]); +} + +const mpq &dvec::operator[] (int i) const +{ + Q_RANGE_CHECK(i, 0, size() - 1); + return(d[i]); +} + +// ************************************************************************** +// +// ************************************************************************** +void dvec::set(const mpq &v) +{ + set(0, size() - 1, v); +} + +void dvec::set(int i0, int i1, const mpq &v) +{ + for (int i = i0; i <= i1; i++) + d[i] = v; +} + +void dvec::assign(int n, const mpq &v) +{ + resize(n); + set(v); +} + +// ************************************************************************** +// +// ************************************************************************** +dvec dvec::operator+(const dvec &b) const +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] + b.d[i]; + + return(r); +} + +dvec dvec::operator-(const dvec &b) const +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] - b.d[i]; + + return(r); +} + +mpq dvec::operator*(const dvec &b) const +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + mpq r; + + for (int i = 0; i < n; i++) + r += d[i] * b.d[i]; + + return(r); +} + +mpq dvec::operator*(const svec &b) const +{ + return(b.operator*(*this)); +} + + +// ************************************************************************** +// +// ************************************************************************** +const dvec &dvec::operator+() const +{ + return(*this); +} + +dvec dvec::operator-() const +{ + int n = size(); + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = -d[i]; + + return(r); +} + +dvec dvec::operator+(const mpq &v) const +{ + int n = size(); + + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] + v; + + return(r); +} + +dvec dvec::operator-(const mpq &v) const +{ + int n = size(); + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] - v; + + return(r); +} + +dvec dvec::operator*(const mpq &v) const +{ + int n = size(); + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] * v; + + return(r); +} + +dvec dvec::operator/(const mpq &v) const +{ + int n = size(); + dvec r(n); + + for (int i = 0; i < n; i++) + r.d[i] = d[i] / v; + + return(r); +} + + +// ************************************************************************** +// +// ************************************************************************** +bool dvec::operator==(const dvec &b) const +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + for (int i = 0; i < n; i++) { + if (d[i] != b.d[i]) + return(false); + } + + return(true); +} + +bool dvec::operator!=(const dvec &b) const +{ + return(!((*this) == b)); +} + +// ************************************************************************** +// +// ************************************************************************** +dvec &dvec::operator=(const dvec &b) +{ + d = b.d; + return(*this); +} + +dvec &dvec::operator=(const svec &b) +{ + b.spread(*this); + return(*this); +} + + +// ************************************************************************** +// +// ************************************************************************** +dvec &dvec::operator+=(const dvec &b) +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + for (int i = 0; i < n; i++) + d[i] += b.d[i]; + + return(*this); +} + +dvec &dvec::operator-=(const dvec &b) +{ + Q_MATCH(size(), b.size()); + + int n = size(); + + for (int i = 0; i < n; i++) + d[i] -= b.d[i]; + + return(*this); +} + +dvec &dvec::operator+=(const mpq &v) +{ + int n = size(); + + for (int i = 0; i < n; i++) + d[i] += v; + + return(*this); +} + +dvec &dvec::operator-=(const mpq &v) +{ + int n = size(); + + for (int i = 0; i < n; i++) + d[i] -= v; + + return(*this); +} + +dvec &dvec::operator*=(const mpq &v) +{ + int n = size(); + + for (int i = 0; i < n; i++) + d[i] *= v; + + return(*this); +} + +dvec &dvec::operator/=(const mpq &v) +{ + int n = size(); + + for (int i = 0; i < n; i++) + d[i] /= v; + + return(*this); +} + + +// ************************************************************************** +// +// ************************************************************************** +void dvec::fdump(FILE *f, const std::string &name) const +{ + int n = size(); + + if (name.length()) + gmp_fprintf(f, "%s = [", name.c_str()); + + for (int i = 0; i < n; i++) { + d[i].fdump(f); + } + + if (name.length()) + gmp_fprintf(f, " ]';\n"); + else + gmp_fprintf(f, ";\n"); +} + +void dvec::dump(const std::string &name) const +{ + fdump(stdout, name); +} + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const dvec &v) +{ + int n = v.size(); + + for (int i = 0; i < n; i++) + os << " " << v[i]; + + return(os); +} + +// ************************************************************************** +// +// ************************************************************************** + +} // namespace + diff --git a/qxx/src/rational.cpp b/qxx/src/rational.cpp new file mode 100644 index 0000000..c634bc1 --- /dev/null +++ b/qxx/src/rational.cpp @@ -0,0 +1,437 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include "qxx/rational.hpp" +#include "qxx/strprintf.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + +// ************************************************************************** +// +// ************************************************************************** +mpq::mpq() +{ + mpq_init(v); +} + +mpq::mpq(const mpq &b) +{ + mpq_init(v); + mpq_set(v, b.v); +} + +mpq::mpq(long num, long den) +{ + mpq_init(v); + mpq_set_si(v, num, den); + mpq_canonicalize(v); +} + +mpq::mpq(int num, int den) +{ + mpq_init(v); + mpq_set_si(v, num, den); + mpq_canonicalize(v); +} + +mpq::mpq(long i) +{ + mpq_init(v); + mpq_set_si(v, i, 1); +} + +mpq::mpq(int i) +{ + mpq_init(v); + mpq_set_si(v, i, 1); +} + +mpq::mpq(double d) +{ + mpq_init(v); + mpq_set_d(v, d); +} + +mpq::mpq(const char *s) +{ + mpq_init(v); + mpq_set_str(v, s, 0); + mpq_canonicalize(v); +} + +mpq::~mpq() +{ + mpq_clear(v); +} + +// ************************************************************************** +// +// ************************************************************************** +const mpq &mpq::operator+() const +{ + return(*this); +} + +mpq mpq::operator-() const +{ + mpq r; + mpq_neg(r.v, v); + return(r); +} + +mpq mpq::operator+(const mpq &b) const +{ + mpq r; + mpq_add(r.v, v, b.v); + return(r); +} + +mpq mpq::operator-(const mpq &b) const +{ + mpq r; + mpq_sub(r.v, v, b.v); + return(r); +} + +mpq mpq::operator*(const mpq &b) const +{ + mpq r; + mpq_mul(r.v, v, b.v); + return(r); +} + +mpq mpq::operator/(const mpq &b) const +{ + mpq r; + mpq_div(r.v, v, b.v); + return(r); +} + + +mpq &mpq::operator=(const mpq &b) +{ + mpq_set(v, b.v); + return(*this); +} + +mpq &mpq::operator=(int i) +{ + mpq_set_si(v, i, 1); + return(*this); +} + +mpq &mpq::operator+=(const mpq &b) +{ + mpq_add(v, v, b.v); + return(*this); +} + +mpq &mpq::operator-=(const mpq &b) +{ + mpq_sub(v, v, b.v); + return(*this); +} + +mpq &mpq::operator*=(const mpq &b) +{ + mpq_mul(v, v, b.v); + return(*this); +} + +mpq &mpq::operator/=(const mpq &b) +{ + mpq_div(v, v, b.v); + return(*this); +} + + +bool mpq::operator<(const mpq &b) const +{ + return(mpq_cmp(v, b.v) < 0); +} + +bool mpq::operator>(const mpq &b) const +{ + return(mpq_cmp(v, b.v) > 0); +} + +bool mpq::operator<=(const mpq &b) const +{ + return(mpq_cmp(v, b.v) <= 0); +} + +bool mpq::operator>=(const mpq &b) const +{ + return(mpq_cmp(v, b.v) >= 0); +} + +bool mpq::operator==(const mpq &b) const +{ + return(mpq_equal(v, b.v)); +} + +bool mpq::operator!=(const mpq &b) const +{ + return(!mpq_equal(v, b.v)); +} + +// ************************************************************************** +// +// ************************************************************************** +void mpq::set_neg() +{ + mpq_neg(v, v); +} + +void mpq::set_neg(const mpq &b) +{ + mpq_neg(v, b.v); +} + +void mpq::set_inv() +{ + mpq_inv(v, v); +} + +void mpq::set_inv(const mpq &b) +{ + mpq_inv(v, b.v); +} + +void mpq::set_abs() +{ + mpq_abs(v, v); +} + +void mpq::set_abs(const mpq &b) +{ + mpq_abs(v, b.v); +} + +mpq mpq::neg() const +{ + mpq r; + mpq_neg(r.v, v); + return(r); +} + +mpq mpq::inv() const +{ + mpq r; + mpq_inv(r.v, v); + return(r); +} + +mpq mpq::abs() const +{ + mpq r; + mpq_abs(r.v, v); + return(r); +} + +int mpq::sign() const +{ + return(mpq_sgn(v)); +} + +// ************************************************************************** +// +// ************************************************************************** +mpq mpq::num() const +{ + mpq r; + + mpz_set(mpq_numref(r.v), mpq_numref(v)); + mpz_set_ui(mpq_denref(r.v), 1); + + return(r); +} + +mpq mpq::den() const +{ + mpq r; + + mpz_set(mpq_numref(r.v), mpq_denref(v)); + mpz_set_ui(mpq_denref(r.v), 1); + + return(r); +} + +int mpq::bits() const +{ + return( + mpz_sizeinbase(mpq_numref(v), 2) + + mpz_sizeinbase(mpq_denref(v), 2) + ); +} + +mpq mpq::frac() const +{ + mpq r; + + mpz_fdiv_r(mpq_numref(r.v), mpq_numref(v), mpq_denref(v)); + mpz_set(mpq_denref(r.v), mpq_denref(v)); + + return(r); +} + +mpq mpq::floor() const +{ + mpq r; + + mpz_fdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v)); + mpz_set_ui(mpq_denref(r.v), 1); + + return(r); +} + +mpq mpq::trunc() const +{ + mpq r; + + mpz_tdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v)); + mpz_set_ui(mpq_denref(r.v), 1); + + return(r); +} + +mpq mpq::ceil() const +{ + mpq r; + + mpz_cdiv_q(mpq_numref(r.v), mpq_numref(v), mpq_denref(v)); + mpz_set_ui(mpq_denref(r.v), 1); + + return(r); +} + + +// ************************************************************************** +// +// ************************************************************************** +long mpq::get_long_num() const +{ + assert(mpz_fits_slong_p(mpq_numref(v))); + + return(mpz_get_si(mpq_numref(v))); +} + +long mpq::get_long_den() const +{ + assert(mpz_fits_slong_p(mpq_denref(v))); + + return(mpz_get_si(mpq_denref(v))); +} + +double mpq::get_double() const +{ + return(mpq_get_d(v)); +} + +// ************************************************************************** +// +// ************************************************************************** +mpq mpq::reduce(mpq max_coeff) const +{ + mpq p[3], q[3]; + mpq f, ff, r; + + p[0] = 0; + q[0] = 1; + p[1] = 1; + q[1] = 0; + f = *this; + + while (1) { + ff = f.floor(); + r = f - ff; + + p[2] = ff * p[1] + p[0]; + q[2] = ff * q[1] + q[0]; + + if ((p[2].abs() > max_coeff) || (q[2] > max_coeff)) + break; + + if (r.sign() == 0) { + p[1] = p[2]; + q[1] = q[2]; + break; + } + + f = r.inv(); + + p[0] = p[1]; + q[0] = q[1]; + p[1] = p[2]; + q[1] = q[2]; + } + + + if (q[1].sign() == 0) + q[1] = 1; + + return(p[1] / q[1]); +} + +// ************************************************************************** +// +// ************************************************************************** +void mpq::fdump(FILE *f, const std::string &name) const +{ + mpq_t t; + + mpq_init(t); + mpq_set(t, v); + mpq_canonicalize(t); + + if (mpq_cmp(t, v) != 0) { + gmp_fprintf(f, "\nBUG: not canon.: %Qd\n", + mpq_numref(v), + mpq_denref(v)); + } + + mpq_clear(t); + + if (name.length()) + gmp_fprintf(f, "%s = %Qd;\n", name.c_str(), v); + else + gmp_fprintf(f, " %Qd", v); +} + +void mpq::dump(const std::string &name) const +{ + fdump(stdout, name); +} + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const mpq &v) +{ + os << strprintf("%Qd", v.v); + return(os); +} + +// ************************************************************************** +// +// ************************************************************************** +} + + diff --git a/qxx/src/slu.cpp b/qxx/src/slu.cpp new file mode 100644 index 0000000..ca1c233 --- /dev/null +++ b/qxx/src/slu.cpp @@ -0,0 +1,307 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include "qxx/strprintf.hpp" +#include "qxx/rational.hpp" +#include "qxx/slu.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +perm::perm() +{ +} + +perm::~perm() +{ +} + +int perm::size() const +{ + return(fwd.size()); +} + +void perm::resize(int n) +{ + fwd.resize(n); + bwd.resize(n); +} + +void perm::clear() +{ + fwd.clear(); + bwd.clear(); +} + +void perm::id() +{ + int n = fwd.size(); + + for (int i = 0; i < n; i++) + fwd[i] = i; + + for (int k = 0; k < n; k++) + bwd[k] = k; +} + +void perm::id(int n) +{ + resize(n); + id(); +} + +void perm::pivot(int k, int i) +{ + int t = fwd[i]; + int u = bwd[k]; + + fwd[i] = k; + fwd[u] = t; + bwd[k] = i; + bwd[t] = u; +} + + +// ************************************************************************** +// +// ************************************************************************** +slu::slu(const smat &a) +{ + factorize(a); +} + +slu::~slu() +{ +} + +// ************************************************************************** +// +// ************************************************************************** +void slu::pivot(int p) +{ + double best_v = n * n; + int best_i = -1; + int best_j = -1; + + for (int k = p; k < n; k++) { + int i = row.bwd[k]; + int inz = w[i].nz(); + + for (int l = 0; l < inz; l++) { + int j = w[i].index(l); + int jnz = cnz[j]; + int kj = col.fwd[j]; + double v = (double)(inz - 1) * (jnz - 1); + + if ((kj >= p) && (v < best_v)) { + best_v = v; + best_i = i; + best_j = j; + } + } + } + + assert((best_i >= 0) && "Singular matrix"); + + //best_i = best_j = p; + row.pivot(p, best_i); + col.pivot(p, best_j); +} + +// ************************************************************************** +// +// ************************************************************************** +void slu::factorize(const smat &a) +{ + // check + assert(a.rows() == a.cols()); + + // init + n = a.rows(); + w = a; + lt.resize(n, n); + u.resize(n, n); + row.id(n); + col.id(n); + + cnz.assign(n, 0); + + for (int i = 0; i < n; i++) { + for (int l = 0; l < w[i].nz(); l++) + cnz[w[i].index(l)]++; + } + + // main loop + svec pr; + mpq pv, f; + std::vector hit; + + for (int p = 0; p < n; p++) { + //w.dense().dump(strprintf("w%d", p).c_str()); + + pivot(p); + + int pi = row.bwd[p]; + int pj = col.bwd[p]; + + // get pivot row + pr = w[pi]; + pr.offs_build(); + pv = pr[pj]; + for (int pl = 0; pl < pr.nz(); pl++) + cnz[pr.index(pl)]--; + + // update u, lt, w + u[p] = w[pi]; + lt[p] = w.get_col(pj) / pv; + w[pi].set_zero(); + + // elimination + for (int k = p + 1; k < n; k++) { + int i = row.bwd[k]; + svec &ir = w[i]; + + f = -ir[pj] / pv; + + if (f.sign() == 0) + continue; + + hit.assign(pr.nz(), false); + + for (int l = 0; l < ir.nz(); l++) { + int j = ir.index(l); + int pl = pr.locate(j); + + if (pl == -1) + continue; + + hit[pl] = true; + + ir.value(l) += f * pr.value(pl); + + if (ir.value(l).sign() == 0) { + cnz[j]--; + ir.remove(l); + l--; + } + } + + for (int pl = 0; pl < pr.nz(); pl++) { + if (hit[pl]) + continue; + int j = pr.index(pl); + + ir.push(j, f * pr.value(pl)); + cnz[j]++; + } + } + } + + // postprocessing + for (int p = 0; p < n; p++) { + svec &pr = u[p]; + + for (int l = 0; l < pr.nz(); l++) + pr.index(l) = col.fwd[pr.index(l)]; + } + + + for (int p = 0; p < n; p++) { + svec &pc = lt[p]; + + for (int l = 0; l < pc.nz(); l++) + pc.index(l) = row.fwd[pc.index(l)]; + } + + l.set_transpose(lt); + ut.set_transpose(u); +} + + +// ************************************************************************** +// +// ************************************************************************** +dvec slu::solve_gen(const smat &l, const smat &u, + const perm &row, const perm &col, const dvec &b) const +{ + dvec x(n), y(n); + mpq v, diag; + + // Ly = b + for (int k = 0; k < n; k++) { + v = b[row.bwd[k]]; + + const svec &lr = l[k]; + + for (int l = 0; l < lr.nz(); l++) { + int j = lr.index(l); + + if (j == k) + diag = lr.value(l); + else + v -= lr.value(l) * y[j]; + } + + y[k] = v / diag; + } + + // Ux = y + for (int k = n - 1; k != -1; k--) { + v = y[k]; + + const svec &ur = u[k]; + + for (int l = 0; l < ur.nz(); l++) { + int j = ur.index(l); + + if (j == k) + diag = ur.value(l); + else + v -= ur.value(l) * y[j]; + } + + y[k] = v / diag; + } + + // postprocessing + for (int j = 0; j < n; j++) + x[j] = y[col.fwd[j]]; + + return(x); +} + +dvec slu::solve_Ax(const dvec &b) const +{ + return(solve_gen(l, u, row, col, b)); +} + +dvec slu::solve_xA(const dvec &b) const +{ + return(solve_gen(ut, lt, col, row, b)); +} + +// ************************************************************************** +// +// ************************************************************************** +}; + diff --git a/qxx/src/smat.cpp b/qxx/src/smat.cpp new file mode 100644 index 0000000..47c4347 --- /dev/null +++ b/qxx/src/smat.cpp @@ -0,0 +1,339 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include "qxx/debug.hpp" +#include "qxx/rational.hpp" +#include "qxx/dvec.hpp" +#include "qxx/slu.hpp" +#include "qxx/smat.hpp" + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// +// ************************************************************************** +smat::smat() + : m(0), n(0) +{ +} + +smat::smat(const smat &b) + : m(b.m), n(b.n), d(b.d) +{ +} + +smat::smat(int m0, int n0) + : m(m0), n(n0), d(m0) +{ + for (int i = 0; i < m; i++) + d[i].resize(n); +} + + +int smat::rows() const +{ + return(m); +} + +int smat::cols() const +{ + return(n); +} + +void smat::resize(int m0, int n0) +{ + m = m0; + n = n0; + + d.resize(m); + for (int i = 0; i < m; i++) + d[i].resize(n); +} + +void smat::clear() +{ + m = n = 0; + d.clear(); +} + +// ************************************************************************** +// +// ************************************************************************** +const mpq &smat::get(int i, int j) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + +void smat::set(int i, int j, const mpq &v) +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + d[i][j] = v; +} + +void smat::set_zero() +{ + for (int i = 0; i < m; i++) + d[i].set_zero(); +} + +void smat::set_identity() +{ + mpq one(1); + + for (int i = 0; i < m; i++) { + d[i].set_zero(); + d[i][i] = one; + } +} + + +// ************************************************************************** +// +// ************************************************************************** +void smat::gather(const dmat &src) +{ + resize(src.rows(), src.cols()); + + for (int i = 0; i < m; i++) + d[i].gather(src[i]); +} + +void smat::spread(dmat &r) const +{ + r.resize(m, n); + + for (int i = 0; i < m; i++) + d[i].spread(r[i]); +} + +dmat smat::dense() const +{ + dmat r; + + spread(r); + + return(r); +} + +// ************************************************************************** +// +// ************************************************************************** +svec smat::get_col(int j) const +{ + svec r(m); + + for (int i = 0; i < m; i++) + r.push(i, d[i][j]); + + return(r); +} + +void smat::set_col(int j, const svec &v) +{ + Q_MATCH(v.size(), m); + + for (int i = 0; i < m; i++) + d[i][j] = v[i]; +} + +const svec &smat::get_row(int i) const +{ + return(d[i]); +} + +void smat::set_row(int i, const svec &v) +{ + Q_MATCH(v.size(), n); + + d[i] = v; +} + + +// ************************************************************************** +// +// ************************************************************************** +svec &smat::operator[] (int i) +{ + return(d[i]); +} + +const svec &smat::operator[] (int i) const +{ + return(d[i]); +} + +// ************************************************************************** +// +// ************************************************************************** +svec_ref smat::operator() (int i, int j) +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + +const mpq &smat::operator() (int i, int j) const +{ + Q_RANGE_CHECK(i, 0, m - 1); + Q_RANGE_CHECK(j, 0, n - 1); + + return(d[i][j]); +} + + +// ************************************************************************** +// +// ************************************************************************** +smat &smat::operator=(const smat &a) +{ + m = a.m; + n = a.n; + d = a.d; + + return(*this); +} + +// ************************************************************************** +// +// ************************************************************************** +dvec smat::operator*(const dvec &b) const +{ + Q_MATCH(n, b.size()); + + dvec x(m); + mpq v; + + for (int i = 0; i < m; i++) { + v = 0; + + const svec &r = d[i]; + for (int l = 0; l < r.nz(); l++) + v += r.value(l) * b[r.index(l)]; + + x[i] = v; + } + + return(x); +} + +// ************************************************************************** +// +// ************************************************************************** +void smat::set_transpose(const smat &a) +{ + resize(a.n, a.m); + + for (int i = 0; i < a.n; i++) { + const svec &ar = a[i]; + + for (int l = 0; l < ar.nz(); l++) + d[ar.index(l)].push(i, ar.value(l)); + } +} + +// ************************************************************************** +// +// ************************************************************************** +smat smat::t() const +{ + smat a; + + a.set_transpose(*this); + + return(a); +} + + +// ************************************************************************** +// +// ************************************************************************** +smat smat::inv() const +{ + Q_MATCH(m, n); + + slu lu(*this); + smat g; + dvec e; + svec x; + + g.resize(n, n); + e.assign(n, 0); + + for (int i = 0; i < n; i++) { + e[i] = 1; + + x = lu.solve_xA(e); + + g.set_row(i, x); + + e[i] = 0; + } + + return(g); +} + + +// ************************************************************************** +// +// ************************************************************************** +void smat::fdump(FILE *f, const std::string &name) const +{ + if (name.length()) + gmp_fprintf(f, "%s = [\n", name.c_str()); + else + gmp_fprintf(f, "[\n"); + + for (int i = 0; i < m; i++) + d[i].fdump(f); + + gmp_fprintf(f, "];\n"); +} + +void smat::dump(const std::string &name) const +{ + fdump(stdout, name); +} + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const smat &a) +{ + for (int i = 0; i < a.rows(); i++) + os << a[i] << std::endl; + + return(os); +} + +// ************************************************************************** +// +// ************************************************************************** +} + + + + diff --git a/qxx/src/strprintf.cpp b/qxx/src/strprintf.cpp new file mode 100644 index 0000000..4e19859 --- /dev/null +++ b/qxx/src/strprintf.cpp @@ -0,0 +1,85 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include +#include +#include +#include "qxx/strprintf.hpp" + + +namespace q { + +// *************************************************************************** +// +// *************************************************************************** +void format(std::string &s, const char *fmt, ...) +{ + va_list ap; + int l; + + va_start(ap, fmt); + l = gmp_vsnprintf(NULL, 0, fmt, ap); + va_end(ap); + + if (l <= 0) { + s = ""; + return; + } + + char *store = new char[l + 1]; + + va_start(ap, fmt); + gmp_vsnprintf(store, l + 1, fmt, ap); + va_end(ap); + + s = store; + + delete[] store; +} + + +// *************************************************************************** +// +// *************************************************************************** +std::string strprintf(const char *fmt, ...) +{ + va_list ap; + int l; + + va_start(ap, fmt); + l = gmp_vsnprintf(NULL, 0, fmt, ap); + va_end(ap); + + if (l <= 0) + return(std::string()); + + char *store = new char[l + 1]; + + va_start(ap, fmt); + gmp_vsnprintf(store, l + 1, fmt, ap); + va_end(ap); + + std::string s = store; + + delete[] store; + + return(s); +} + + + +} // namespace + diff --git a/qxx/src/svec.cpp b/qxx/src/svec.cpp new file mode 100644 index 0000000..526c596 --- /dev/null +++ b/qxx/src/svec.cpp @@ -0,0 +1,695 @@ +/* + This file is part of qxx -- matrix algebra in exact arithmetic + Copyright (C) 2013-2014 Laurent Poirrier + + libp 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, version 3 of the License. + + 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. + + You should have received a copy of the GNU General Public License + along with pxx. If not, see . +*/ +#include "qxx/debug.hpp" +#include "qxx/svec.hpp" +#include + +// ************************************************************************** +// +// ************************************************************************** +namespace q { + + +// ************************************************************************** +// svec_ref: construction and destruction +// ************************************************************************** +svec_ref::svec_ref(svec &vec, int i) + : vector(vec), index(i) +{ +} + +svec_ref::svec_ref(const svec_ref &src) + : vector(src.vector), index(src.index) +{ +} + +svec_ref::~svec_ref() +{ +} + + +// ************************************************************************** +// svec_ref: methods +// ************************************************************************** +const mpq &svec_ref::get() const +{ + return(vector.get(index)); +} + +// ************************************************************************** +// svec_ref: assignment operators +// ************************************************************************** +svec_ref::operator mpq() const +{ + return(get()); +} + +const mpq &svec_ref::operator=(const svec_ref &b) const +{ + return(vector.set(index, b.get())); +} + +const mpq &svec_ref::operator=(const mpq &b) const +{ + return(vector.set(index, b)); +} + +const mpq &svec_ref::operator+=(const mpq &b) const +{ + return(vector.set(index, vector.get(index) + b)); +} + +const mpq &svec_ref::operator-=(const mpq &b) const +{ + return(vector.set(index, vector.get(index) - b)); +} + +const mpq &svec_ref::operator*=(const mpq &b) const +{ + return(vector.set(index, vector.get(index) * b)); +} + +const mpq &svec_ref::operator/=(const mpq &b) const +{ + return(vector.set(index, vector.get(index) / b)); +} + +// ************************************************************************** +// svec_ref: read-only operators +// ************************************************************************** +const mpq &svec_ref::operator+() const +{ + return(get()); +} + +mpq svec_ref::operator-() const +{ + return(-get()); +} + +mpq svec_ref::operator+(const mpq &b) const +{ + return(get() + b); +} + +mpq svec_ref::operator-(const mpq &b) const +{ + return(get() - b); +} + +mpq svec_ref::operator*(const mpq &b) const +{ + return(get() * b); +} + +mpq svec_ref::operator/(const mpq &b) const +{ + return(get() / b); +} + + +bool svec_ref::operator<(const mpq &b) const +{ + return(get() < b); +} + +bool svec_ref::operator>(const mpq &b) const +{ + return(get() > b); +} + +bool svec_ref::operator<=(const mpq &b) const +{ + return(get() <= b); +} + +bool svec_ref::operator>=(const mpq &b) const +{ + return(get() >= b); +} + +bool svec_ref::operator==(const mpq &b) const +{ + return(get() == b); +} + +bool svec_ref::operator!=(const mpq &b) const +{ + return(get() != b); +} + +// ************************************************************************** +// svec_el +// ************************************************************************** +svec_el::svec_el() +{ +} + +svec_el::svec_el(int idx, const mpq &v) + : index(idx), value(v) +{ +} + +svec_el::svec_el(const svec_el &e) + : index(e.index), value(e.value) +{ +} + +svec_el::~svec_el() +{ +} + +// ************************************************************************** +// svec +// ************************************************************************** +const mpq svec::zero; + +svec::svec() + : n(0) +{ +} + +svec::svec(const svec &src) + : n(src.n), d(src.d) +{ +} + +svec::svec(int n0) + : n(n0) +{ +} + +svec::svec(const dvec &src) + : n(src.size()) +{ + gather(src); +} + +svec::~svec() +{ +} + + +// ************************************************************************** +// +// ************************************************************************** +void svec::clear() +{ + n = 0; + d.clear(); +} + + +int svec::size() const +{ + return(n); +} + +void svec::resize(int n0) +{ + n = n0; +} + +int svec::nz() const +{ + return(d.size()); +} + +// ************************************************************************** +// +// ************************************************************************** +bool svec::offs_active() const +{ + return(offs.size() != 0); +} + +void svec::offs_build() +{ + if (offs_active()) + return; + + offs.resize(n); + + for (int i = 0; i < n; i++) + offs[i] = 0; + + for (int l = 0; l < nz(); l++) + offs[d[l].index] = l + 1; +} + +void svec::offs_clear() +{ + offs.resize(0); +} + + +// ************************************************************************** +// +// ************************************************************************** +void svec::set_zero() +{ + d.clear(); + + if (offs_active()) { + for (int i = 0; i < n; i++) + offs[i] = 0; + } +} + +// ************************************************************************** +// +// ************************************************************************** +void svec::gather(const dvec &src) +{ + n = src.size(); + d.clear(); + for (int i = 0; i < n; i++) + push(i, src[i]); +} + +void svec::spread(dvec &r) const +{ + r.resize(n); + r.set(zero); + for (int l = 0; l < nz(); l++) + r[d[l].index] = d[l].value; +} + +dvec svec::dense() const +{ + dvec r; + + spread(r); + return(r); +} + +// ************************************************************************** +// +// ************************************************************************** +int svec::locate(int idx) const +{ + Q_RANGE_CHECK(idx, 0, n - 1); + + if (offs_active()) + return(offs[idx] - 1); + + for (int l = 0; l < nz(); l++) { + if (d[l].index == idx) + return(l); + } + + return(-1); +} + +// ************************************************************************** +// +// ************************************************************************** +void svec::remove(int l) +{ + Q_RANGE_CHECK(l, 0, nz() - 1); + + int k = nz() - 1; + + if (offs_active()) + offs[d[l].index] = 0; + + if (l < k) + d[l] = d[k]; + + d.resize(k); +} + +// ************************************************************************** +// +// ************************************************************************** +int &svec::index(int offset) +{ + Q_RANGE_CHECK(offset, 0, nz() - 1); + + return(d[offset].index); +} + +mpq &svec::value(int offset) +{ + Q_RANGE_CHECK(offset, 0, nz() - 1); + + return(d[offset].value); +} + + +int svec::index(int offset) const +{ + Q_RANGE_CHECK(offset, 0, nz() - 1); + + return(d[offset].index); +} + +const mpq &svec::value(int offset) const +{ + Q_RANGE_CHECK(offset, 0, nz() - 1); + + return(d[offset].value); +} + +// ************************************************************************** +// +// ************************************************************************** +void svec::push_nz(int idx, const mpq &v) +{ + Q_RANGE_CHECK(idx, 0, n - 1); + + d.push_back(svec_el(idx, v)); + + if (offs_active()) + offs[idx] = nz(); +} + +void svec::push(int idx, const mpq &v) +{ + Q_RANGE_CHECK(idx, 0, n - 1); + + if (v.sign()) + push_nz(idx, v); +} + +const mpq &svec::get(int idx) const +{ + Q_RANGE_CHECK(idx, 0, n - 1); + + int l = locate(idx); + if (l == -1) + return(zero); + return(d[l].value); +} + +const mpq &svec::set(int idx, const mpq &v) +{ + Q_RANGE_CHECK(idx, 0, n - 1); + + int l = locate(idx); + + if (v.sign()) { + if (l == -1) { + push_nz(idx, v); + return(d[nz() - 1].value); + } + + d[l].value = v; + return(d[l].value); + } else { + if (l != -1) + remove(l); + return(zero); + } + +} + +// ************************************************************************** +// +// ************************************************************************** +const mpq &svec::operator[](int idx) const +{ + return(get(idx)); +} + +svec_ref svec::operator[](int idx) +{ + return(svec_ref(*this, idx)); +} + + +// ************************************************************************** +// +// ************************************************************************** +const svec &svec::operator+() const +{ + return(*this); +} + +svec svec::operator-() const +{ + svec r(*this); + + for (int l = 0; l < nz(); l++) + r.d[l].value.set_neg(); + + return(r); +} + + +svec svec::operator+(const svec &b) const +{ + Q_MATCH(n, b.n); + + svec r(*this); + + r += b; + + return(r); +} + +svec svec::operator-(const svec &b) const +{ + Q_MATCH(n, b.n); + + svec r(*this); + + r -= b; + + return(r); +} + + +// ************************************************************************** +// +// ************************************************************************** +mpq svec::operator*(const svec &b) const +{ + Q_MATCH(n, b.n); + + mpq v; + + if (offs_active()) { + for (int l = 0; l < b.nz(); l++) { + int i = b.d[l].index; + if ((i < n) && (offs[i])) + v += d[offs[i] - 1].value * b.d[l].value; + } + } else if (b.offs_active()) { + for (int l = 0; l < nz(); l++) { + int i = d[l].index; + if ((i < b.n) && (b.offs[i])) + v += d[l].value * b.d[b.offs[i] - 1].value; + } + } else { + dvec db; + + b.spread(db); + + v = (*this) * db; + } + + return(v); +} + +mpq svec::operator*(const dvec &b) const +{ + Q_MATCH(n, b.size()); + + mpq v; + + for (int l = 0; l < nz(); l++) { + int j = d[l].index; + + if (j < b.size()) + v += d[l].value * b[j]; + } + + return(v); +} + + +// ************************************************************************** +// +// ************************************************************************** +svec svec::operator*(const mpq &v) +{ + svec r(n); + + if (!v.sign()) + return(r); + + for (int l = 0; l < nz(); l++) + r.push(d[l].index, d[l].value * v); + + return(r); +} + +svec svec::operator/(const mpq &v) +{ + svec r(n); + + for (int l = 0; l < nz(); l++) + r.push(d[l].index, d[l].value / v); + + return(r); +} + +// ************************************************************************** +// +// ************************************************************************** +svec &svec::operator=(const svec &b) +{ + n = b.n; + d = b.d; + + if ((offs_active()) && (b.offs_active())) + offs = b.offs; + else + offs_clear(); + + return(*this); +} + +svec &svec::operator+=(const svec &b) +{ + Q_MATCH(n, b.n); + + if (offs_active()) { + for (int l = 0; l < b.nz(); l++) { + int i = b.d[l].index; + + if ((i < n) && (offs[i])) + d[offs[i] - 1].value += b.d[l].value; + else + push(i, b.d[l].value); + } + } else { + dvec da, db; + + spread(da); + b.spread(db); + + da += db; + + gather(da); + } + + return(*this); +} + +svec &svec::operator-=(const svec &b) +{ + Q_MATCH(n, b.n); + + if (offs_active()) { + for (int l = 0; l < b.nz(); l++) { + int i = b.d[l].index; + + if ((i < n) && (offs[i])) + d[offs[i] - 1].value -= b.d[l].value; + else + push(i, -b.d[l].value); + } + } else { + dvec da, db; + + spread(da); + b.spread(db); + + da += db; + + gather(da); + } + + return(*this); +} + +svec &svec::operator*=(const mpq &v) +{ + if (!v.sign()) { + set_zero(); + return(*this); + } + + for (int l = 0; l < nz(); l++) + d[l].value *= v; + + return(*this); +} + +svec &svec::operator/=(const mpq &v) +{ + for (int l = 0; l < nz(); l++) + d[l].value /= v; + + return(*this); +} + +// ************************************************************************** +// +// ************************************************************************** +static bool comp(const svec_el &a, const svec_el &b) +{ + return(a.index < b.index); +} + +void svec::sort() +{ + std::sort(d.begin(), d.end(), comp); +} + +// ************************************************************************** +// +// ************************************************************************** +void svec::fdump(FILE *f, const std::string &name) const +{ + if (name.length()) + gmp_fprintf(f, "%s = ", name.c_str()); + + for (int l = 0; l < nz(); l++) { + gmp_fprintf(f, " %d:%Zd/%Zd", + d[l].index, + mpq_numref(d[l].value.v), + mpq_denref(d[l].value.v) + ); + } + + gmp_fprintf(f, "\n"); +} + +void svec::dump(const std::string &name) const +{ + fdump(stdout, name); +} + +// ************************************************************************** +// +// ************************************************************************** +std::ostream &operator<<(std::ostream &os, const svec &v) +{ + for (int l = 0; l < v.nz(); l++) + os << " " << v.index(l) << ":" << v.value(l); + + return(os); +} + + +// ************************************************************************** +// +// ************************************************************************** +}