diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e72ecddc..d7cefc16 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -65,27 +65,58 @@ jobs: steps: - uses: actions/checkout@v4 + # set up OpenMp for ColPack + - name: Install OpenMP (Linux) + if: startsWith(matrix.os, 'ubuntu') + run: sudo apt-get install libomp-dev + + - name: Install libomp (macOS) + if: startsWith(matrix.os, 'macos') + run: | + brew install libomp + - name: Set reusable strings - # Turn repeated input strings (such as the build output directory) into step outputs. These step outputs can be used throughout the workflow file. + # Turn repeated input strings (such as the build output directory) into step outputs. These step outputs can be used throughout the workflow file. id: strings shell: bash run: | echo "build-output-dir=${{ github.workspace }}/build" >> "$GITHUB_OUTPUT" - - name: Configure and build ADOL-C (Linux, MacOs) - if: matrix.c_compiler != 'msvc' + - name: Configure and build ADOL-C (Linux) + if: startsWith(matrix.os, 'ubuntu') + run: > + cmake -B ${{ steps.strings.outputs.build-output-dir }} + -DCMAKE_CXX_COMPILER=${{ matrix.cpp_compiler }} + -DCMAKE_C_COMPILER=${{ matrix.c_compiler }} + -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} + -DENABLE_DOCEXA=1 + -DENABLE_ADDEXA=1 + -DENABLE_SPARSE=1 + -DBUILD_SHARED_LIBS=${{ matrix.build_shared_libs }} + -S ${{ github.workspace }} + + + - name: Configure and build ADOL-C (MacOs) + if: startsWith(matrix.os, 'macos') run: > cmake -B ${{ steps.strings.outputs.build-output-dir }} -DCMAKE_CXX_COMPILER=${{ matrix.cpp_compiler }} -DCMAKE_C_COMPILER=${{ matrix.c_compiler }} -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} + -DOpenMP_C_FLAGS="-Xclang -fopenmp -I/opt/homebrew/opt/libomp/include" + -DOpenMP_C_LIB_NAMES="omp" + -DOpenMP_CXX_FLAGS="-Xclang -fopenmp -I/opt/homebrew/opt/libomp/include" + -DOpenMP_CXX_LIB_NAMES="omp" + -DOpenMP_omp_LIBRARY="/opt/homebrew/opt/libomp/lib/libomp.dylib" -DENABLE_DOCEXA=1 -DENABLE_ADDEXA=1 + -DENABLE_SPARSE=1 -DBUILD_SHARED_LIBS=${{ matrix.build_shared_libs }} -S ${{ github.workspace }} + # we don't use sparse drivers on windows! - name: Configure and build ADOL-C (Windows) - if: matrix.c_compiler == 'msvc' + if: startsWith(matrix.os, 'windows') shell: powershell run: # Setup Visual Studio build environment diff --git a/.github/workflows/coverage.yml b/.github/workflows/coverage.yml index bc57ad38..5918e87f 100644 --- a/.github/workflows/coverage.yml +++ b/.github/workflows/coverage.yml @@ -37,6 +37,7 @@ jobs: -DCMAKE_C_COMPILER=gcc-13 \ -DCMAKE_BUILD_TYPE=Debug \ -DBUILD_TESTS_WITH_COV=ON \ + -DENABLE_SPARSE=1 \ -S .. \ -B . cmake --build . @@ -62,5 +63,4 @@ jobs: files: build/coverage.info token: ${{ secrets.CODECOV_TOKEN }} flags: unittests - name: code-coverage-report - + name: code-coverage-report \ No newline at end of file diff --git a/.github/workflows/min_compiler_ubuntu.yml b/.github/workflows/min_compiler_ubuntu.yml index 5caee600..16fde669 100644 --- a/.github/workflows/min_compiler_ubuntu.yml +++ b/.github/workflows/min_compiler_ubuntu.yml @@ -1,4 +1,4 @@ -name: Build ubuntu-latest with g++11 and clang++-13 +name: minimal compiler version ubuntu-latest on: push: branches: ["master"] diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index a42c0a5b..840a0e8e 100644 --- a/.github/workflows/tests.yml +++ b/.github/workflows/tests.yml @@ -48,7 +48,8 @@ jobs: -DCMAKE_CXX_COMPILER=${{ matrix.cpp_compiler }} -DCMAKE_C_COMPILER=${{ matrix.c_compiler }} -DCMAKE_BUILD_TYPE=${{ matrix.build_type }} - -DBUILD_TESTS=ON \ + -DBUILD_TESTS=ON + -DENABLE_SPARSE=1 -S ${{ github.workspace }} - name: Build ADOL-C diff --git a/.gitignore b/.gitignore index 32e03b1f..dddad4c9 100644 --- a/.gitignore +++ b/.gitignore @@ -58,7 +58,6 @@ install_manifest.txt *.dll build/ lib/ -include/ # Nix build results result* diff --git a/ADOL-C/boost-test/CMakeLists.txt b/ADOL-C/boost-test/CMakeLists.txt index 8130221e..115b8f2d 100644 --- a/ADOL-C/boost-test/CMakeLists.txt +++ b/ADOL-C/boost-test/CMakeLists.txt @@ -48,6 +48,12 @@ list(APPEND SOURCE_FILES ${VALUETAPE_FILES}) file(GLOB VALUETAPE_FILES "ho_rev/*.cpp") list(APPEND SOURCE_FILES ${VALUETAPE_FILES}) +if (ENABLE_SPARSE) +# Add all source files from sparse +file(GLOB VALUETAPE_FILES "sparse/*.cpp") +list(APPEND SOURCE_FILES ${VALUETAPE_FILES}) +endif() + add_executable(boost-test-adolc ${SOURCE_FILES}) target_include_directories(boost-test-adolc PRIVATE "${ADOLC_INCLUDE_DIR}") target_link_libraries(boost-test-adolc PRIVATE diff --git a/ADOL-C/boost-test/sparse/sparse_abs_normal.cpp b/ADOL-C/boost-test/sparse/sparse_abs_normal.cpp new file mode 100644 index 00000000..48d6b760 --- /dev/null +++ b/ADOL-C/boost-test/sparse/sparse_abs_normal.cpp @@ -0,0 +1,459 @@ +/* +This file defines eight abs-smooth test functions f_V, V = 1,...,8, +each represented by a templated taping routine `taping()`. +All functions use ADOL-C's ABS-normal form for studying piecewise +differentiability, switching structures, and Jacobian sparsity. + +-------------------------------------------------------------------- +VERSION 1 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 2 +Input x = (x0, x1) + +f1(x) = ( f1_0(x), f1_1(x) ) where + f1_0(x) = x0 + |x0 - x1| + |x0 - |x1||, + f1_1(x) = x1. + +sparsity pattern (crs): +[3, 0, 2, 4] +[1, 1] +[2, 0, 1] +[1, 1] +[3, 0, 3] + +-------------------------------------------------------------------- +VERSION 2 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 2 +Input x = (x0, x1) + +f2_0(x) = x0 + |x0 - x1| + |x0 - |x1||, +f2_1(x) = |x1 - 5|. + +sparsity pattern (crs): +[3, 0, 2, 4] +[1, 5] +[2, 0, 1] +[1, 1] +[2, 0, 3] +[1, 1] + +-------------------------------------------------------------------- +VERSION 3 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 1 +Input x = (x0, x1) + +f3(x) = max( 0, x1² - max(0, x0) ). + +sparsity pattern (crs): +[4, 0, 1, 2, 3] +[1, 0] +[3, 0, 1, 2] + +-------------------------------------------------------------------- +VERSION 4 +-------------------------------------------------------------------- +dimIn = 5, dimOut = 1 +Input x = (x0, …, x4) + +Let + z1 = | sum_{j=0..4} x_j / (j+1) |. + +For each i = 1..4 define + z2(i) = | sum_{j=0..4} x_j / (i + j + 1) |. + +Then + f4(x) = max( z1, z2(1), z2(2), z2(3), z2(4) ). + +sparsity pattern (crs): +[9, 5, 6, 7, 8, 9, 10, 11, 12, 13] +[5, 0, 1, 2, 3, 4] +[5, 0, 1, 2, 3, 4] +[2, 5, 6] +[5, 0, 1, 2, 3, 4] +[4, 5, 6, 7, 8] +[5, 0, 1, 2, 3, 4] +[7, 5, 6, 7, 8, 9, 10] +[5, 0, 1, 2, 3, 4] +[8, 5, 6, 7, 8, 9, 10, 11, 12] +-------------------------------------------------------------------- +VERSION 5 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 1 +Input x = (x0, x1) + +f5(x) = max + 0, 2x0 - 5x1, 3x0 + 2x1, 2x0 + 5x1, 3x0 - 2x1, -100). + + +sparsity patter (crs): +[7, 0, 1, 2, 3, 4, 5, 6] +[2, 0, 1] +[3, 0, 1, 2] +[4, 0, 1, 2, 3] +[5, 0, 1, 2, 3, 4] +[6, 0, 1, 2, 3, 4, 5] +-------------------------------------------------------------------- +VERSION 6 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 1 +Input x = (x0, x1) + + f6(x) = ( z1 + |x1² - x0 - 1 - (x0 + |-x0|)/2| ) / 2 . + +sparsity pattern (crs): +[4, 0, 1, 2, 3] +[1, 0] +[3, 0, 1, 2] +-------------------------------------------------------------------- +VERSION 7 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 2 +Input x = (x0, x1) + +Let z = |x1|. Then + + f7_0(x) = max( -100, 3x0 + 2z, 2x0 + 5z ), + f7_1(x) = max( -100, 3x0 + 2z, 2x0 + 5z ). + +sparsity patter (crs): +[4, 0, 2, 3, 4] +[4, 0, 2, 5, 6] +[1, 1] +[2, 0, 2] +[3, 0, 2, 3] +[2, 0, 2] +[3, 0, 2, 5] + + +-------------------------------------------------------------------- +VERSION 8 +-------------------------------------------------------------------- +dimIn = 2, dimOut = 1 +Input x = (x0, x1) + +Let z = |x1|. Then + + f8(x) = max( -100, 3x0 + 2z, 2x0 + 5z ). + +sparsity pattern (crs): +[3, 0, 2, 3, 4] +[1, 1] +[2, 0, 2] +[3, 0, 2, 3] +-------------------------------------------------------------------- +*/ +#include +#define BOOST_TEST_DYN_LINK +#include +#include +#include + +BOOST_AUTO_TEST_SUITE(test_sparse_abs_normal) + +template struct version_trait {}; + +template <> struct version_trait<1> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {3.0, 2.0}; + static constexpr size_t dimOut = 2; + static constexpr std::array, 5> crs = { + {{3, 0, 2, 4}, {1, 1, 0, 0}, {2, 0, 1, 0}, {1, 1, 0, 0}, {2, 0, 3, 0}}}; +}; + +template <> struct version_trait<2> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {3.0, 2.0}; + static constexpr size_t dimOut = 2; + static constexpr std::array, 6> crs = {{{3, 0, 2, 4}, + {1, 5, 0, 0}, + {2, 0, 1, 0}, + {1, 1, 0, 0}, + {2, 0, 3, 0}, + {1, 1, 0, 0}}}; +}; + +template <> struct version_trait<3> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {1.0, 1.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 3> crs = { + {{4, 0, 1, 2, 3}, {1, 0, 0, 0, 0}, {3, 0, 1, 2, 0}}}; +}; + +template <> struct version_trait<4> { + static constexpr size_t dimIn = 5; + static constexpr std::array init = {1.0, 1.0, 1.0, 1.0, 1.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 10> crs = { + {{9, 5, 6, 7, 8, 9, 10, 11, 12, 13}, + {5, 0, 1, 2, 3, 4, 0, 0, 0, 0}, + {5, 0, 1, 2, 3, 4, 0, 0, 0, 0}, + {2, 5, 6, 0, 0, 0, 0, 0, 0, 0}, + {5, 0, 1, 2, 3, 4, 0, 0, 0, 0}, + {4, 5, 6, 7, 8, 0, 0, 0, 0, 0}, + {5, 0, 1, 2, 3, 4, 0, 0, 0, 0}, + {6, 5, 6, 7, 8, 9, 10, 0, 0, 0}, + {5, 0, 1, 2, 3, 4, 0, 0, 0, 0}, + {8, 5, 6, 7, 8, 9, 10, 11, 12, 0}}}; +}; + +template <> struct version_trait<5> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {0.0, 0.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 6> crs = { + {{7, 0, 1, 2, 3, 4, 5, 6}, + {2, 0, 1, 0, 0, 0, 0, 0}, + {3, 0, 1, 2, 0, 0, 0, 0}, + {4, 0, 1, 2, 3, 0, 0, 0}, + {5, 0, 1, 2, 3, 4, 0, 0}, + {6, 0, 1, 2, 3, 4, 5, 0}}}; +}; + +template <> struct version_trait<6> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {0.0, -1.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 3> crs = { + {{4, 0, 1, 2, 3}, {1, 0, 0, 0, 0}, {3, 0, 1, 2, 0}}}; +}; + +template <> struct version_trait<7> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {2.0, 2.0}; + static constexpr size_t dimOut = 2; + static constexpr std::array, 7> crs = { + {{4, 0, 2, 3, 4}, + {4, 0, 2, 5, 6}, + {1, 1, 0, 0, 0}, + {2, 0, 2, 0, 0}, + {3, 0, 2, 3, 0}, + {2, 0, 2, 0, 0}, + {3, 0, 2, 5, 0}}}; +}; + +template <> struct version_trait<8> { + static constexpr size_t dimIn = 2; + static constexpr std::array init = {2.0, 2.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 4> crs = { + {{4, 0, 2, 3, 4}, {1, 1, 0, 0, 0}, {2, 0, 2, 0, 0}, {3, 0, 2, 3, 0}}}; +}; + +struct UnAllocated {}; +struct Allocated {}; + +template struct ANFProblem; + +template struct ANFProblem { + static constexpr auto init = version_trait::init; + static constexpr size_t dimIn = version_trait::dimIn; + static constexpr size_t dimOut = version_trait::dimOut; + + std::array in; + std::array out; + + static constexpr short tapeId = 716 + Version; + + constexpr ANFProblem() : in(init) {} + ANFProblem allocateBuffers(short numSwitchingVars) { + return ANFProblem(numSwitchingVars, *this); + } +}; + +template struct ANFProblem { + static constexpr auto init = version_trait::init; + static constexpr size_t dimIn = version_trait::dimIn; + static constexpr size_t dimOut = version_trait::dimOut; + + static constexpr auto crs = version_trait::crs; + + std::array in; + std::array out; + + static constexpr short tapeId = 716 + Version; + + int numSwitchingVars; + std::vector z; + std::vector cz; + std::vector cy; + + double **J, **Y; + double **Z, **L; + + std::vector d; + std::vector g; + double **gradz; + + std::vector sigma_x; + std::vector sigma_g; + + ANFProblem(ANFProblem &anfProblem) = delete; + ANFProblem(ANFProblem &&anfProblem) = delete; + + ANFProblem operator=(ANFProblem &anfProblem) = delete; + ANFProblem operator=(ANFProblem &&anfProblem) = delete; + + ~ANFProblem() { + myfree2(Z); + myfree2(L); + myfree2(J); + myfree2(Y); + myfree2(gradz); + } + constexpr ANFProblem(short numSVars, + const ANFProblem &base) + : in(base.init), out(base.out) { + + numSwitchingVars = numSVars; + sigma_x.resize(numSwitchingVars); + sigma_g.resize(numSwitchingVars); + z.resize(numSwitchingVars); + + cz.resize(numSwitchingVars); + cy.resize(dimOut); + Z = myalloc2(numSwitchingVars, dimIn); + L = myalloc2(numSwitchingVars, numSwitchingVars); + J = myalloc2(dimOut, numSwitchingVars); + Y = myalloc2(dimOut, dimIn); + + d.resize(dimIn); + g.resize(dimIn); + gradz = myalloc2(numSwitchingVars, dimIn); + } +}; + +template +static void computeANF(ANFProblem &anfProblemAlloc) { + + zos_pl_forward(anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, + anfProblemAlloc.dimIn, 1, anfProblemAlloc.in.data(), + anfProblemAlloc.out.data(), anfProblemAlloc.z.data()); + + abs_normal(anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, + anfProblemAlloc.dimIn, anfProblemAlloc.numSwitchingVars, + anfProblemAlloc.in.data(), anfProblemAlloc.out.data(), + anfProblemAlloc.z.data(), anfProblemAlloc.cz.data(), + anfProblemAlloc.cy.data(), anfProblemAlloc.Y, anfProblemAlloc.J, + anfProblemAlloc.Z, anfProblemAlloc.L); +} + +template +static void taping(ANFProblem &anfProblem) { + std::vector x(anfProblem.dimIn); + std::vector y(anfProblem.dimOut); + if constexpr (Version == 1) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + y[0] = x[0] + fabs(x[0] - x[1]) + fabs(x[0] - fabs(x[1])); + y[1] = x[1]; + y[0] >>= anfProblem.out[0]; + y[1] >>= anfProblem.out[1]; + } else if (Version == 2) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + y[0] = x[0] + fabs(x[0] - x[1]) + fabs(x[0] - fabs(x[1])); + y[1] = fabs(x[1] - 5); + y[0] >>= anfProblem.out[0]; + y[1] >>= anfProblem.out[1]; + } else if (Version == 3) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + y[0] = fmax(0, x[1] * x[1] - fmax(0, x[0])); + y[0] >>= anfProblem.out[0]; + } else if (Version == 4) { + { + for (int i = 0; i < anfProblem.dimIn; i++) { + x[i] <<= anfProblem.in[i]; + } + int j = 0; + adouble z1 = + fabs(std::accumulate(x.begin(), x.end(), adouble{0.0}, + [&](adouble &&res, adouble val) -> adouble { + res += val / (j + 1); + j++; + return res; + })); + for (int i = 1; i < anfProblem.dimIn; i++) { + adouble z2 = 0; + for (int j = 0; j < anfProblem.dimIn; j++) { + z2 += x[j] / (i + j + 1); + } + z2 = fabs(z2); + z1 = fmax(z1, z2); + } + y[0] = z1; + y[0] >>= anfProblem.out[0]; + } + } else if (Version == 5) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + y[0] = 0; + y[0] = fmax(y[0], 2 * x[0] - 5 * x[1]); + y[0] = fmax(y[0], 3 * x[0] + 2 * x[1]); + y[0] = fmax(y[0], 2 * x[0] + 5 * x[1]); + y[0] = fmax(y[0], 3 * x[0] - 2 * x[1]); + y[0] = fmax(y[0], static_cast(-100)); + y[0] >>= anfProblem.out[0]; + } else if (Version == 6) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + adouble z1 = x[1] * x[1] - x[0] - 1.0 - (x[0] + fabs(-x[0])) / 2.0; + y[0] = (z1 + fabs(-z1)) / 2.0; + y[0] >>= anfProblem.out[0]; + } else if (Version == 7) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + adouble z1 = fabs(x[1]); + y[0] = fmax(fmax(-100, 3 * x[0] + 2 * z1), 2 * x[0] + 5 * z1); + y[1] = fmax(fmax(-100, 3 * x[0] + 2 * z1), 2 * x[0] + 5 * z1); + y[0] >>= anfProblem.out[0]; + y[1] >>= anfProblem.out[1]; + } else if (Version == 8) { + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + adouble z1 = fabs(x[1]); + y[0] = fmax(fmax(-100, 3 * x[0] + 2 * z1), 2 * x[0] + 5 * z1); + y[0] >>= anfProblem.out[0]; + } +} +template static void problem() { + ANFProblem anfProblem{}; + createNewTape(anfProblem.tapeId); + setCurrentTape(anfProblem.tapeId); + currentTape().enableMinMaxUsingAbs(); + trace_on(anfProblem.tapeId); + taping(anfProblem); + trace_off(); + + int numSwitchingVars = get_num_switches(anfProblem.tapeId); + auto anfProblemAlloc = anfProblem.allocateBuffers(numSwitchingVars); + computeANF(anfProblemAlloc); + + std::vector crs(anfProblemAlloc.dimOut + + anfProblemAlloc.numSwitchingVars); + std::span crsSpan(crs); + ADOLC::Sparse::absnormal_jac_pat( + anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, anfProblemAlloc.dimIn, + anfProblemAlloc.numSwitchingVars, anfProblemAlloc.in.data(), crsSpan); + + for (int row{0}; row < anfProblemAlloc.crs.size(); row++) { + BOOST_TEST(crs[row][0] == anfProblemAlloc.crs[row][0]); + for (int col{1}; col < anfProblemAlloc.crs[row][0]; col++) { + BOOST_TEST(crs[row][col] == anfProblemAlloc.crs[row][col]); + } + }; +} + +BOOST_AUTO_TEST_CASE(SparseANFPatTest) { + constexpr std::array problems = { + problem<1>, problem<2>, problem<3>, problem<4>, + problem<5>, problem<6>, problem<7>, problem<8>}; + for (auto problem : problems) { + problem(); + } +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/ADOL-C/boost-test/sparse/sparse_hessian.cpp b/ADOL-C/boost-test/sparse/sparse_hessian.cpp new file mode 100644 index 00000000..351a50b5 --- /dev/null +++ b/ADOL-C/boost-test/sparse/sparse_hessian.cpp @@ -0,0 +1,544 @@ +#include +#include +#include +#define BOOST_TEST_DYN_LINK +#include +namespace tt = boost::test_tools; + +#include "../const.h" +#include +#include + +BOOST_AUTO_TEST_SUITE(test_sparse_hessian) + +using ADOLC::Sparse::ControlFlowMode; +using ADOLC::Sparse::RecoveryMethod; + +/* + This test constructs a simple function whose Hessian is strictly diagonal. + It is designed to verify that ADOL-C’s sparse Hessian driver correctly + identifies and recovers only the diagonal nonzero entries. + + The function is: + f(x) = Σ_{i=0}^{4} [ sin(x_i) + (x_{i+10})^2 ] + + Analytical Hessian (expected): + For i = 0,... , 4: + (i,i) = -sin(x_i) + For i = 10,... , 14: + (i,i) = 2 + All off-diagonal entries are zero. + + Therefore, the Hessian is diagonal with 10 nonzero entries in total, + all located on the diagonal. The sparse driver should only report those + diagonal positions (rowIndices[i] == columnIndices[i]). + + Total unique non-zero positions: 10 (purely diagonal). +*/ + +template +static void testSparseHessWithDiagonal(short tapeId) { + createNewTape(tapeId); + constexpr int dimOut = 5; + constexpr int dimIn = 20; + + std::array in; + int i = 0; + std::for_each(in.begin(), in.end(), [&i](double &v) { + i++; + v = i * 1.0 / (i + 2.0); + }); + + trace_on(tapeId); + { + std::array indeps; + for (int i{0}; i < indeps.size(); ++i) + indeps[i] <<= in[i]; + + adouble deps; + for (int i{0}; i < dimOut; ++i) + deps += sin(indeps[i]) + indeps[i + 10] * indeps[i + 10]; + + double out; + deps >>= out; + } + trace_off(); + + std::array hess; + for (auto &h : hess) { + h = new double[dimIn]; + for (int j{0}; j < dimIn; j++) { + h[j] = 0.0; + } + } + hessian(tapeId, dimIn, in.data(), hess.data()); + unsigned int *rowIndices = nullptr; /* row indices */ + unsigned int *columnIndices = nullptr; /* column indices */ + int numberOfNonzeros = 0; + double *sparseValues = nullptr; + + ADOLC::Sparse::sparse_hess(tapeId, dimIn, 0, in.data(), + &numberOfNonzeros, &rowIndices, + &columnIndices, &sparseValues); + std::array trueSolution = {-std::sin(in[0]), + -std::sin(in[1]), + -std::sin(in[2]), + -std::sin(in[3]), + -std::sin(in[4]), + 2, + 2, + 2, + 2, + 2}; + + BOOST_TEST(numberOfNonzeros == trueSolution.size()); + + for (int i{0}; i < numberOfNonzeros; i++) { + // test how it should be + BOOST_TEST(sparseValues[i] == trueSolution[i]); + + // sanity check with non-sparse hess + if (i < 5) + BOOST_TEST(hess[i][i] == sparseValues[i]); + else if (i >= 5) + BOOST_TEST(hess[i + 5][i + 5] == sparseValues[i]); + + // we only have diagonal elements + BOOST_TEST(rowIndices[i] == columnIndices[i]); + } + + // test with rowIndices, columnIndices and sparseValues != nullptr + ADOLC::Sparse::sparse_hess(tapeId, dimIn, 0, in.data(), + &numberOfNonzeros, &rowIndices, + &columnIndices, &sparseValues); + + BOOST_TEST(numberOfNonzeros == trueSolution.size()); + + for (int i{0}; i < numberOfNonzeros; i++) { + // test how it should be + BOOST_TEST(sparseValues[i] == trueSolution[i]); + + // sanity check with non-sparse hess + if (i < 5) + BOOST_TEST(hess[i][i] == sparseValues[i]); + else if (i >= 5) + BOOST_TEST(hess[i + 5][i + 5] == sparseValues[i]); + + // we only have diagonal elements + BOOST_TEST(rowIndices[i] == columnIndices[i]); + } + + // test with repeat == 1 + ADOLC::Sparse::sparse_hess(tapeId, dimIn, 1, in.data(), + &numberOfNonzeros, &rowIndices, + &columnIndices, &sparseValues); + + BOOST_TEST(numberOfNonzeros == trueSolution.size()); + + for (int i{0}; i < numberOfNonzeros; i++) { + // test how it should be + BOOST_TEST(sparseValues[i] == trueSolution[i]); + + // sanity check with non-sparse hess + if (i < 5) + BOOST_TEST(hess[i][i] == sparseValues[i]); + else if (i >= 5) + BOOST_TEST(hess[i + 5][i + 5] == sparseValues[i]); + + // we only have diagonal elements + BOOST_TEST(rowIndices[i] == columnIndices[i]); + } + delete sparseValues; + for (auto &h : hess) + delete[] h; +} + +BOOST_AUTO_TEST_CASE(SparseHessSafeDirect) { + const short tapeId = 632; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessSafeIndirect) { + const short tapeId = 633; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessTightDirect) { + const short tapeId = 634; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessTightIndirect) { + const short tapeId = 635; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessOldSafeDirect) { + const short tapeId = 636; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessOldSafeIndirect) { + const short tapeId = 637; + testSparseHessWithDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessOldTightDirect) { + const short tapeId = 638; + testSparseHessWithDiagonal( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessOldTightIndirect) { + const short tapeId = 639; + testSparseHessWithDiagonal(tapeId); +} + +/* + This test constructs a function with mixed-product terms so that the Hessian + contains off-diagonal nonzeros. + + f(x) = sin(x0) + sin(x1) + + x2*x3 + + 0.5*x0*x2 + + x4*x5 + + x4*x4 + + Analytical Hessian nonzeros (expected): + (0,0) = -sin(x0) + (1,1) = -sin(x1) + (2,3) = 1 + (3,2) = 1 // symmetric of (2,3) + (0,2) = 0.5 + (2,0) = 0.5 // symmetric of (0,2) + (4,5) = 1 + (5,4) = 1 // symmetric of (4,5) + (4,4) = 2 + Total unique non-zero positions: 9 (but sparse_hess return only lower triangle + like hessian driver). +*/ + +template +static void testSparseHessWithOffDiagonals(short tapeId) { + createNewTape(tapeId); + constexpr int dimIn = 6; + + std::array in; + int i = 0; + std::for_each(in.begin(), in.end(), [&i](double &v) { + i++; + v = 0.1 * (i + 1); + }); + + trace_on(tapeId); + { + std::array indeps; + for (int i{0}; i < indeps.size(); ++i) + indeps[i] <<= in[i]; + + adouble deps = 0.0; + deps += sin(indeps[0]); // contributes (0,0) -> -sin(x0) + deps += sin(indeps[1]); // contributes (1,1) -> -sin(x1) + deps += indeps[2] * indeps[3]; // contributes (2,3) and (3,2) -> 1 + deps += 0.5 * indeps[0] * indeps[2]; // contributes (0,2) and (2,0) -> 0.5 + deps += indeps[4] * indeps[5]; // contributes (4,5) and (5,4) -> 1 + deps += indeps[4] * indeps[4]; // contributes (4,4) -> 2 + + double out; + deps >>= out; + } + trace_off(); + + // compute dense hessian with ADOL-C for sanity comparison + std::array hess; + for (auto &h : hess) { + h = new double[dimIn]; + for (int j{0}; j < dimIn; j++) { + h[j] = 0.0; + } + } + hessian(tapeId, dimIn, in.data(), hess.data()); + + // call sparse_hess + unsigned int *rowIndices = nullptr; /* row indices */ + unsigned int *columnIndices = nullptr; /* column indices */ + double *sparseValues = nullptr; /* values */ + int numberOfNonzeros = 0; + + ADOLC::Sparse::sparse_hess(tapeId, dimIn, 0, in.data(), + &numberOfNonzeros, &rowIndices, + &columnIndices, &sparseValues); + + // build the expected set of nonzeros (list of triples) + struct SparseResult { + unsigned r, c; + double v; + }; + std::vector expected{{0, 0, -std::sin(in[0])}, + {1, 1, -std::sin(in[1])}, + {3, 2, 1.0}, + {2, 0, 0.5}, + {5, 4, 1.0}, + {4, 4, 2.0}}; + + BOOST_TEST(numberOfNonzeros == expected.size()); + for (int i{0}; i < numberOfNonzeros; i++) { + bool found = false; + + auto row = std::max(rowIndices[i], columnIndices[i]); + auto col = std::min(rowIndices[i], columnIndices[i]); + + // search if the computed row and column indices are present + for (auto &expect : expected) { + if (expect.r == row && expect.c == col) { + found = true; + BOOST_TEST(expect.v == sparseValues[i]); + BOOST_TEST(hess[row][col] == sparseValues[i], tt::tolerance(tol)); + } + } + if (!found) { + auto s = "Non-zero at (" + std::to_string(rowIndices[i]) + "," + + std::to_string(columnIndices[i]) + ") not expected!"; + BOOST_CHECK(found && s.c_str()); + } + } + + // cleanup + delete[] rowIndices; + delete[] columnIndices; + delete[] sparseValues; + for (auto &ptr : hess) + delete[] ptr; +} + +BOOST_AUTO_TEST_CASE(SparseHessNonDiagSafeDirect) { + const short tapeId = 700; + testSparseHessWithOffDiagonals( + tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagSafeIndirect) { + const short tapeId = 701; + testSparseHessWithOffDiagonals(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagTightDirect) { + const short tapeId = 702; + testSparseHessWithOffDiagonals(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagTightIndirect) { + const short tapeId = 703; + testSparseHessWithOffDiagonals(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessNonDiagOldSafeDirect) { + const short tapeId = 704; + testSparseHessWithOffDiagonals(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagOldSafeIndirect) { + const short tapeId = 705; + testSparseHessWithOffDiagonals(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagOldTightDirect) { + const short tapeId = 706; + testSparseHessWithOffDiagonals(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseHessNonDiagOldTightIndirect) { + const short tapeId = 707; + testSparseHessWithOffDiagonals(tapeId); +} + +template +static void testSparseHessPatWithDiagonal(short tapeId) { + createNewTape(tapeId); + constexpr int dimOut = 5; + constexpr int dimIn = 20; + + std::array in; + int i = 0; + std::for_each(in.begin(), in.end(), [&i](double &v) { + i++; + v = i * 1.0 / (i + 2.0); + }); + + trace_on(tapeId); + { + std::array indeps; + for (int i{0}; i < indeps.size(); ++i) + indeps[i] <<= in[i]; + + adouble deps; + for (int i{0}; i < dimOut; ++i) + deps += sin(indeps[i]) + indeps[i + 10] * indeps[i + 10]; + + double out; + deps >>= out; + } + trace_off(); + std::array compressedRowStorage; + std::span compressedRowStorageSpan(compressedRowStorage); + ADOLC::Sparse::hess_pat(tapeId, dimIn, in.data(), + compressedRowStorageSpan); + + // first hessian row only has first entry + BOOST_TEST(compressedRowStorage[0][0] == 1); + BOOST_TEST(compressedRowStorage[0][1] == 0); + + // etc... + BOOST_TEST(compressedRowStorage[1][0] == 1); + BOOST_TEST(compressedRowStorage[1][1] == 1); + + BOOST_TEST(compressedRowStorage[2][0] == 1); + BOOST_TEST(compressedRowStorage[2][1] == 2); + + BOOST_TEST(compressedRowStorage[3][0] == 1); + BOOST_TEST(compressedRowStorage[3][1] == 3); + + BOOST_TEST(compressedRowStorage[4][0] == 1); + BOOST_TEST(compressedRowStorage[4][1] == 4); + + BOOST_TEST(compressedRowStorage[5][0] == 0); + BOOST_TEST(compressedRowStorage[6][0] == 0); + BOOST_TEST(compressedRowStorage[7][0] == 0); + BOOST_TEST(compressedRowStorage[8][0] == 0); + BOOST_TEST(compressedRowStorage[9][0] == 0); + + BOOST_TEST(compressedRowStorage[10][0] == 1); + BOOST_TEST(compressedRowStorage[10][1] == 10); + + BOOST_TEST(compressedRowStorage[11][0] == 1); + BOOST_TEST(compressedRowStorage[11][1] == 11); + + BOOST_TEST(compressedRowStorage[12][0] == 1); + BOOST_TEST(compressedRowStorage[12][1] == 12); + + BOOST_TEST(compressedRowStorage[13][0] == 1); + BOOST_TEST(compressedRowStorage[13][1] == 13); + + BOOST_TEST(compressedRowStorage[14][0] == 1); + BOOST_TEST(compressedRowStorage[14][1] == 14); + + BOOST_TEST(compressedRowStorage[15][0] == 0); + BOOST_TEST(compressedRowStorage[16][0] == 0); + BOOST_TEST(compressedRowStorage[17][0] == 0); + BOOST_TEST(compressedRowStorage[18][0] == 0); + BOOST_TEST(compressedRowStorage[19][0] == 0); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatDiagSafe) { + const short tapeId = 708; + testSparseHessPatWithDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatDiagTight) { + const short tapeId = 709; + testSparseHessPatWithDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatDiagOldSafe) { + const short tapeId = 710; + testSparseHessPatWithDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatDiagOldTight) { + const short tapeId = 711; + testSparseHessPatWithDiagonal(tapeId); +} + +template +static void testSparseHessPatWithOffDiagonal(short tapeId) { + createNewTape(tapeId); + constexpr int dimIn = 6; + + std::array in; + int i = 0; + std::for_each(in.begin(), in.end(), [&i](double &v) { + i++; + v = 0.1 * (i + 1); + }); + + trace_on(tapeId); + { + std::array indeps; + for (int i{0}; i < indeps.size(); ++i) + indeps[i] <<= in[i]; + + adouble deps = 0.0; + deps += sin(indeps[0]); // contributes (0,0) -> -sin(x0) + deps += sin(indeps[1]); // contributes (1,1) -> -sin(x1) + deps += indeps[2] * indeps[3]; // contributes (2,3) and (3,2) -> 1 + deps += 0.5 * indeps[0] * indeps[2]; // contributes (0,2) and (2,0) -> 0.5 + deps += indeps[4] * indeps[5]; // contributes (4,5) and (5,4) -> 1 + deps += indeps[4] * indeps[4]; // contributes (4,4) -> 2 + + double out; + deps >>= out; + } + trace_off(); + + std::array compressedRowStorage; + std::span compressedRowStorageSpan(compressedRowStorage); + ADOLC::Sparse::hess_pat(tapeId, dimIn, in.data(), + compressedRowStorageSpan); + + // first hessian row has first entry and third + BOOST_TEST(compressedRowStorage[0][0] == 2); + BOOST_TEST(compressedRowStorage[0][1] == 0); + BOOST_TEST(compressedRowStorage[0][2] == 2); + + // second hessian row has second entry + BOOST_TEST(compressedRowStorage[1][0] == 1); + BOOST_TEST(compressedRowStorage[1][1] == 1); + + // third hessian row has first and fourth + BOOST_TEST(compressedRowStorage[2][0] == 2); + BOOST_TEST(compressedRowStorage[2][1] == 0); + BOOST_TEST(compressedRowStorage[2][2] == 3); + + // fourth hessian row has third + BOOST_TEST(compressedRowStorage[3][0] == 1); + BOOST_TEST(compressedRowStorage[3][1] == 2); + + // fifth hessian row has fifth and sixth + BOOST_TEST(compressedRowStorage[4][0] == 2); + BOOST_TEST(compressedRowStorage[4][1] == 4); + BOOST_TEST(compressedRowStorage[4][2] == 5); + + // sixth hessian row has fifth + BOOST_TEST(compressedRowStorage[5][0] == 1); + BOOST_TEST(compressedRowStorage[5][1] == 4); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatOffDiagSafe) { + const short tapeId = 712; + testSparseHessPatWithOffDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatOffDiagTight) { + const short tapeId = 713; + testSparseHessPatWithOffDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatOffDiagOldSafe) { + const short tapeId = 714; + testSparseHessPatWithOffDiagonal(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseHessPatOffDiagOldTight) { + const short tapeId = 715; + testSparseHessPatWithOffDiagonal(tapeId); +} + +BOOST_AUTO_TEST_SUITE_END() \ No newline at end of file diff --git a/ADOL-C/boost-test/sparse/sparse_jacobian.cpp b/ADOL-C/boost-test/sparse/sparse_jacobian.cpp new file mode 100644 index 00000000..cae17f8c --- /dev/null +++ b/ADOL-C/boost-test/sparse/sparse_jacobian.cpp @@ -0,0 +1,230 @@ + +#include +#include +#define BOOST_TEST_DYN_LINK +#include +namespace tt = boost::test_tools; + +#include +#include + +#include "../const.h" + +BOOST_AUTO_TEST_SUITE(test_sparse_jacobian) + +using ADOLC::Sparse::BitPatternPropagationDirection; +using ADOLC::Sparse::CompressionMode; +using ADOLC::Sparse::ControlFlowMode; +using ADOLC::Sparse::SparseMethod; + +template < + SparseMethod SM, CompressionMode CM, ControlFlowMode CFM, + BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto> +static void testSparseJac(short tapeId) { + createNewTape(tapeId); + constexpr int dimOut = 10; + constexpr int dimIn = 20; + + std::array in; + in.fill(std::rand()); + trace_on(tapeId); + { + std::array indeps; + for (int i{0}; i < indeps.size(); ++i) + indeps[i] <<= in[i]; + + std::array deps; + for (int i{0}; i < deps.size(); ++i) + deps[i] = sin(indeps[i]) + indeps[i + 10] * indeps[i + 10]; + + std::array out; + for (int i{0}; i < deps.size(); ++i) + deps[i] >>= out[i]; + } + trace_off(); + std::array jac; + for (auto &j : jac) + j = new double[dimIn]; + jacobian(tapeId, dimOut, dimIn, in.data(), jac.data()); + + unsigned int *rowIndices = nullptr; /* row indices */ + unsigned int *columnIndices = nullptr; /* column indices */ + double *nonzeroValues = nullptr; /* values */ + int numberOfNonzeros = 0; + + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 0, in.data(), &numberOfNonzeros, &rowIndices, + &columnIndices, &nonzeroValues); + BOOST_TEST(numberOfNonzeros == 20); + // go through all entries: either 0 or == sparse_jac + int nonzeroCounter = 0; + for (int row{0}; row < dimOut; ++row) { + for (int col{0}; col < dimIn; ++col) { + if (jac[row][col] != 0) { + BOOST_TEST(jac[row][col] == nonzeroValues[nonzeroCounter++], + tt::tolerance(tol)); + } + } + } + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 0, in.data(), &numberOfNonzeros, &rowIndices, + &columnIndices, &nonzeroValues); + BOOST_TEST(numberOfNonzeros == 20); + nonzeroCounter = 0; + // go through all entries: either 0 or == sparse_jac + for (int row{0}; row < dimOut; ++row) { + for (int col{0}; col < dimIn; ++col) { + if (jac[row][col] != 0) { + BOOST_TEST(jac[row][col] == nonzeroValues[nonzeroCounter++], + tt::tolerance(tol)); + } + } + } + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 1, in.data(), &numberOfNonzeros, &rowIndices, + &columnIndices, &nonzeroValues); + BOOST_TEST(numberOfNonzeros == 20); + nonzeroCounter = 0; + // go through all entries: either 0 or == sparse_jac + for (int row{0}; row < dimOut; ++row) { + for (int col{0}; col < dimIn; ++col) { + if (jac[row][col] != 0) { + BOOST_TEST(jac[row][col] == nonzeroValues[nonzeroCounter++], + tt::tolerance(tol)); + } + } + } + for (auto &j : jac) + delete[] j; +} + +BOOST_AUTO_TEST_CASE(SparseJacIndexColSafe) { + const short tapeId = 565; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacIndexColTight) { + const short tapeId = 566; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacIndexRowSafe) { + const short tapeId = 567; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacIndexRowTight) { + const short tapeId = 568; + testSparseJac(tapeId); +} +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropTight) { + const short tapeId = 569; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropSafe) { + const short tapeId = 570; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropForwardTight) { + const short tapeId = 571; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropForwardSafe) { + const short tapeId = 572; + testSparseJac( + tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropReverseTight) { + const short tapeId = 573; + testSparseJac(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacBitPatterPropReverseSafe) { + const short tapeId = 574; + testSparseJac( + tapeId); +} + +template +static void testSparseJacPat(short tapeId) { + createNewTape(tapeId); + constexpr int dimIn = 7; + constexpr int dimOut = 2; + + std::array in; + int i = 0; + std::for_each(in.begin(), in.end(), [&i](double &v) { v = (i++ * 2.3); }); + + trace_on(tapeId); + { + std::array indeps; + int i = 0; + std::for_each(indeps.begin(), indeps.end(), + [&i, &in](auto &&v) { v <<= in[i++]; }); + adouble out1 = pow(indeps[0], 3); + out1 += indeps[1] * indeps[4]; + + double dummyOut = 0.0; + out1 >>= dummyOut; + + adouble out2 = exp(indeps[6] + indeps[5]); + out2 -= indeps[2]; + out2 >>= dummyOut; + } + trace_off(); + + std::array compressedRowStorage; + std::span compressedRowStorageSpan(compressedRowStorage); + ADOLC::Sparse::jac_pat(tapeId, dimOut, dimIn, in.data(), + compressedRowStorageSpan); + + // first output depends on 0, 1 and 4 + BOOST_TEST(compressedRowStorage[0][0] == 3); + BOOST_TEST(compressedRowStorage[0][1] == 0); + BOOST_TEST(compressedRowStorage[0][2] == 1); + BOOST_TEST(compressedRowStorage[0][3] == 4); + + // second output depends on 2, 5 and 6 + BOOST_TEST(compressedRowStorage[1][0] == 3); + BOOST_TEST(compressedRowStorage[1][1] == 2); + BOOST_TEST(compressedRowStorage[1][2] == 5); + BOOST_TEST(compressedRowStorage[1][3] == 6); +} + +BOOST_AUTO_TEST_CASE(SparseJacPatIndexSafe) { + short tapeId = 575; + testSparseJacPat(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacPatBitPatternSafe) { + short tapeId = 576; + testSparseJacPat(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacPatIndexTight) { + short tapeId = 577; + testSparseJacPat(tapeId); +} + +BOOST_AUTO_TEST_CASE(SparseJacPatBitPatternTight) { + short tapeId = 578; + testSparseJacPat(tapeId); +} + +BOOST_AUTO_TEST_SUITE_END() diff --git a/ADOL-C/examples/additional_examples/CMakeLists.txt b/ADOL-C/examples/additional_examples/CMakeLists.txt index 5f499c1c..9688e63b 100644 --- a/ADOL-C/examples/additional_examples/CMakeLists.txt +++ b/ADOL-C/examples/additional_examples/CMakeLists.txt @@ -1,6 +1,16 @@ add_executable(checkpointing checkpointing/checkpointing.cpp) target_link_libraries(checkpointing PRIVATE adolc) - add_executable(helm_auto_example helm/helm-auto-exam.cpp) -target_link_libraries(helm_auto_example PRIVATE adolc) \ No newline at end of file +target_link_libraries(helm_auto_example PRIVATE adolc) + +if (ENABLE_SPARSE) + add_executable(abs_normal sparse/abs_normal.cpp) + target_link_libraries(abs_normal PRIVATE adolc) + + add_executable(sparse_jacobian sparse/sparse_jacobian.cpp) + target_link_libraries(sparse_jacobian PRIVATE adolc) + + add_executable(sparse_hessian sparse/sparse_hessian.cpp) + target_link_libraries(sparse_hessian PRIVATE adolc) +endif() diff --git a/ADOL-C/examples/additional_examples/sparse/abs_normal.cpp b/ADOL-C/examples/additional_examples/sparse/abs_normal.cpp new file mode 100644 index 00000000..e869d4f5 --- /dev/null +++ b/ADOL-C/examples/additional_examples/sparse/abs_normal.cpp @@ -0,0 +1,170 @@ +/* +This file defines an abs-smooth example function f and computes its jacobian +pattern + +dimIn = 2, dimOut = 1 +Input x = (x0, x1) + +Let z = |x1|. Then + + f8(x) = max( -100, 3x0 + 2z, 2x0 + 5z ). + +sparsity pattern (crs): +[3, 0, 2, 3, 4] +[1, 1] +[2, 0, 2] +[3, 0, 2, 3] +*/ + +#include +#include +#include + +struct UnAllocated {}; +struct Allocated {}; + +template struct ANFProblem; +template <> struct ANFProblem; +template <> struct ANFProblem { + static constexpr size_t dimIn = 2; + static constexpr std::array in = {2.0, 2.0}; + static constexpr size_t dimOut = 1; + static constexpr std::array, 4> crs = { + {{4, 0, 2, 3, 4}, {1, 1, 0, 0, 0}, {2, 0, 2, 0, 0}, {3, 0, 2, 3, 0}}}; + + std::array out; + + static constexpr short tapeId = 1; + + constexpr ANFProblem() = default; + ANFProblem allocateBuffers(short numSwitchingVars); +}; + +template <> struct ANFProblem { + static constexpr size_t dimIn = 2; + static constexpr size_t dimOut = 1; + static constexpr std::array, 4> crs = { + {{4, 0, 2, 3, 4}, {1, 1, 0, 0, 0}, {2, 0, 2, 0, 0}, {3, 0, 2, 3, 0}}}; + + std::array in; + std::array out; + + static constexpr short tapeId = 1; + + int numSwitchingVars; + std::vector z; + std::vector cz; + std::vector cy; + + double **J, **Y; + double **Z, **L; + + std::vector d; + std::vector g; + double **gradz; + + std::vector sigma_x; + std::vector sigma_g; + + ANFProblem(ANFProblem &anfProblem) = delete; + ANFProblem(ANFProblem &&anfProblem) = delete; + + ANFProblem operator=(ANFProblem &anfProblem) = delete; + ANFProblem operator=(ANFProblem &&anfProblem) = delete; + + ~ANFProblem() { + myfree2(Z); + myfree2(L); + myfree2(J); + myfree2(Y); + myfree2(gradz); + } + constexpr ANFProblem(short numSVars, const ANFProblem &base) + : in(base.in), out(base.out) { + + numSwitchingVars = numSVars; + sigma_x.resize(numSwitchingVars); + sigma_g.resize(numSwitchingVars); + z.resize(numSwitchingVars); + + cz.resize(numSwitchingVars); + cy.resize(dimOut); + Z = myalloc2(numSwitchingVars, dimIn); + L = myalloc2(numSwitchingVars, numSwitchingVars); + J = myalloc2(dimOut, numSwitchingVars); + Y = myalloc2(dimOut, dimIn); + + d.resize(dimIn); + g.resize(dimIn); + gradz = myalloc2(numSwitchingVars, dimIn); + } +}; + +ANFProblem +ANFProblem::allocateBuffers(short numSwitchingVars) { + return ANFProblem(numSwitchingVars, *this); +} +void computeANF(ANFProblem &anfProblemAlloc) { + zos_pl_forward(anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, + anfProblemAlloc.dimIn, 1, anfProblemAlloc.in.data(), + anfProblemAlloc.out.data(), anfProblemAlloc.z.data()); + + abs_normal(anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, + anfProblemAlloc.dimIn, anfProblemAlloc.numSwitchingVars, + anfProblemAlloc.in.data(), anfProblemAlloc.out.data(), + anfProblemAlloc.z.data(), anfProblemAlloc.cz.data(), + anfProblemAlloc.cy.data(), anfProblemAlloc.Y, anfProblemAlloc.J, + anfProblemAlloc.Z, anfProblemAlloc.L); +} + +void taping(ANFProblem &anfProblem) { + std::vector x(anfProblem.dimIn); + std::vector y(anfProblem.dimOut); + x[0] <<= anfProblem.in[0]; + x[1] <<= anfProblem.in[1]; + adouble z1 = fabs(x[1]); + y[0] = fmax(fmax(-100, 3 * x[0] + 2 * z1), 2 * x[0] + 5 * z1); + y[0] >>= anfProblem.out[0]; +} + +void problem() { + ANFProblem anfProblem{}; + createNewTape(anfProblem.tapeId); + currentTape().enableMinMaxUsingAbs(); + trace_on(anfProblem.tapeId); + taping(anfProblem); + trace_off(); + + int numSwitchingVars = get_num_switches(anfProblem.tapeId); + auto anfProblemAlloc = anfProblem.allocateBuffers(numSwitchingVars); + computeANF(anfProblemAlloc); + + std::vector crs(anfProblemAlloc.dimOut + + anfProblemAlloc.numSwitchingVars); + std::span crsSpan(crs); + ADOLC::Sparse::absnormal_jac_pat( + anfProblemAlloc.tapeId, anfProblemAlloc.dimOut, anfProblemAlloc.dimIn, + anfProblemAlloc.numSwitchingVars, anfProblemAlloc.in.data(), crsSpan); + + std::cout << "Jacobian Pattern (ANF): " << std::endl; + for (int i = 0; i < anfProblemAlloc.dimOut + anfProblemAlloc.numSwitchingVars; + i++) { + std::cout << "dependent variable: " << i << ": " + << "number of non-zeros: " << crs[i][0] << " (should be " + << anfProblemAlloc.crs[i][0] << ")\n"; + std::cout << "non-zero indices: "; + for (int j = 1; j <= crs[i][0]; j++) + std::cout << crs[i][j] << ", "; + + std::cout << " (should be "; + for (int j = 1; j <= anfProblemAlloc.crs[i][0]; j++) + std::cout << anfProblemAlloc.crs[i][j] << ", "; + std::cout << ")\n"; + } +} +//****************************************************************** + +int main() { + problem(); + return 0; +} diff --git a/ADOL-C/examples/additional_examples/sparse/jacpatexam.cpp b/ADOL-C/examples/additional_examples/sparse/jacpatexam.cpp index 9dd562a8..f83ec434 100644 --- a/ADOL-C/examples/additional_examples/sparse/jacpatexam.cpp +++ b/ADOL-C/examples/additional_examples/sparse/jacpatexam.cpp @@ -15,31 +15,24 @@ /****************************************************************************/ /* INCLUDES */ -#include "../clock/myclock.h" #include -#include #include #include #include -using namespace std; - /****************************************************************************/ /* DEFINES */ -#define TAG 11 - /****************************************************************************/ /* EVALUATION FUNCTIONS */ /*--------------------------------------------------------------------------*/ const unsigned int N = 5, M = 6; -void eval_small(short tag, double *xp, double *yp) { +void eval_small(short tapeId, double *xp, double *yp) { unsigned int i, j; - - trace_on(tag); + trace_on(tapeId); adouble *x, *y; x = new adouble[N]; @@ -66,17 +59,17 @@ void eval_small(short tag, double *xp, double *yp) { trace_off(1); // force a numbered tape file to be written - cout << "\nproblem definition in " << __FILE__ << ", lines " << PD1B - << " - " << PD1E << "\n"; + std::cout << "\nproblem definition in " << __FILE__ << ", lines " << PD1B + << " - " << PD1E << "\n"; } /*--------------------------------------------------------------------------*/ const unsigned int NM = 961; // PQ_STRIPMINE_MAX * 8*sizeof(size_t) + 1 -void eval_arrow_like_matrix(short tag, double *xp, double *yp) { +void eval_arrow_like_matrix(short tapeId, double *xp, double *yp) { unsigned int i, j; - trace_on(tag); + trace_on(tapeId); adouble *x, *y; x = new adouble[NM]; @@ -104,14 +97,15 @@ void eval_arrow_like_matrix(short tag, double *xp, double *yp) { trace_off(1); // force a numbered tape file to be written - cout << "\nproblem definition in " << __FILE__ << ", lines " << PD2B - << " - " << PD2E << "\n"; + std::cout << "\nproblem definition in " << __FILE__ << ", lines " << PD2B + << " - " << PD2E << "\n"; } /****************************************************************************/ /* MAIN PROGRAM */ int main(void) { - short tag; + const short tapeId = 1; + createNewTape(tapeId); int ret_c = -1, choice; int oper, op_buf_size, loc_buf_size, con_buf_size, maxlive, deaths; unsigned int depen, indep; @@ -129,22 +123,25 @@ int main(void) { double basepoint; int ctrl[3]; - cout << "----------------------------------------------------------------\n"; - cout << "\n Jacobian Pattern Example\n\n"; - cout << "Tape identification tag ( [-4..-1] for standart examples ) : ?\b"; - cin >> choice; + std::cout + << "----------------------------------------------------------------\n"; + std::cout << "\n Jacobian Pattern Example\n\n"; + std::cout + << "Tape identification tag ( [-4..-1] for standart examples ) : ?\b"; + std::cin >> choice; - cout << "\n\nOutput Jacobian pattern? (y/n) ?\b"; - cin >> outp; + std::cout << "\n\nOutput Jacobian pattern? (y/n) ?\b"; + std::cin >> outp; - cout << "\n\nCompare with the full Jacobian calculation? (y/n) ?\b"; - cin >> full_jac; + std::cout << "\n\nCompare with the full Jacobian calculation? (y/n) ?\b"; + std::cin >> full_jac; if ((full_jac == 'y') || (full_jac == 'Y')) full_jac = 1; else full_jac = 0; - cout << "----------------------------------------------------------------\n"; + std::cout + << "----------------------------------------------------------------\n"; /*--------------------------------------------------------------------------*/ @@ -156,13 +153,10 @@ int main(void) { base[i] = i + 1; value = new double[M]; + eval_small(tapeIdSmall, base, value); - tag = TAG; - eval_small(tag, base, value); - - cout << "\n\nCreated ADOL-C tape with identification tag " << tag - << ".\n\n"; - cout.flush(); + std::cout << "\n\nCreated ADOL-C tape with identification tag " + << tapeIdSmall << ".\n\n"; } else // ( choice == -4 ) ----------------------------------------------- { @@ -172,64 +166,22 @@ int main(void) { base[i] = i; value = new double[NM]; + eval_arrow_like_matrix(tapeId, base, value); - tag = TAG; - eval_arrow_like_matrix(tag, base, value); - - cout << "\n\nCreated ADOL-C tape with identification tag " << tag - << ".\n\n"; - cout.flush(); + std::cout << "\n\nCreated ADOL-C tape with identification tag " << tapeId + << ".\n\n"; } - tapestats(tag, tape_stats); /* Reading of tape statistics */ - indep = tape_stats[TapeInfos::NUM_INDEPENDENTS]; - depen = tape_stats[TapeInfos::NUM_DEPENDENTS]; - op_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - loc_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - con_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - oper = tape_stats[TapeInfos::NUM_OPERATIONS]; - deaths = tape_stats[TapeInfos::TAY_STACK_SIZE]; - maxlive = tape_stats[TapeInfos::NUM_MAX_LIVES]; - - cout << "\nTape " << tag << ": \n"; - cout << " independents " << indep << "\n"; - cout << " dependents " << depen << "\n"; - cout << " operations " << oper << "\n"; - cout << " operations buffer size " << op_buf_size << "\n"; - cout << " locations buffer size " << loc_buf_size << "\n"; - cout << " constants buffer size " << con_buf_size << "\n"; - cout << " maxlive " << maxlive << "\n"; - cout << " valstack size " << deaths << "\n"; - cout << "\n"; + printTapeStats(tapeId); } else // ( choice >= 0 ) : Take a written tape ------------------------------ { - tag = choice; - - cout << "\nproblem definition in tape " << tag << "\n"; - - tapestats(tag, tape_stats); - indep = tape_stats[TapeInfos::NUM_INDEPENDENTS]; - depen = tape_stats[TapeInfos::NUM_DEPENDENTS]; - op_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - loc_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - con_buf_size = tape_stats[TapeInfos::OP_BUFFER_SIZE]; - oper = tape_stats[TapeInfos::NUM_OPERATIONS]; - deaths = tape_stats[TapeInfos::TAY_STACK_SIZE]; - maxlive = tape_stats[TapeInfos::NUM_MAX_LIVES]; - - cout << "\nTape " << tag << ": \n"; - cout << " independents " << indep << "\n"; - cout << " dependents " << depen << "\n"; - cout << " operations " << oper << "\n"; - cout << " operations buffer size " << op_buf_size << "\n"; - cout << " locations buffer size " << loc_buf_size << "\n"; - cout << " constants buffer size " << con_buf_size << "\n"; - cout << " maxlive " << maxlive << "\n"; - cout << " valstack size " << deaths << "\n"; - - cout << "\n\nbasepoint[0.." << indep << "] = ?\b"; - cin >> basepoint; + tapeId = choice; + + std::cout << "\nproblem definition in tape " << tapeId << "\n"; + + printTapeStats(tapeId); + std::cin >> basepoint; base = new double[indep]; value = new double[depen]; @@ -238,12 +190,12 @@ int main(void) { base[i] = basepoint; } - tape_doc(tag, depen, indep, base, value); // write a tape into a tex-file + tape_doc(tapeId, depen, indep, base, value); // write a tape into a tex-file /*--------------------------------------------------------------------------*/ /* Jacobian pattern by index domains, safe */ - cout << "\nJacobian Pattern by index domain propagation, safe ...\n"; + std::cout << "\nJacobian Pattern by index domain propagation, safe ...\n"; jacpat = new unsigned int *[depen]; @@ -252,17 +204,17 @@ int main(void) { ctrl[2] = 0; // safe z1 = myclock(); - ret_c = jac_pat(tag, depen, indep, base, jacpat, ctrl); + ret_c = ADOLC::Sparse::jac_pat(tapeId, depen, indep, base, jacpat, ctrl); z2 = myclock(); if ((outp == 'y') || (outp == 'Y')) { for (i = 0; i < depen; i++) { - cout << "depen[" << i << "], " << jacpat[i][0] << " non-zero entries :\n"; + std::cout << "depen[" << i << "], " << jacpat[i][0] + << " non-zero entries :\n"; for (j = 1; j <= jacpat[i][0]; j++) - cout << jacpat[i][j] << " "; - cout << "\n"; + std::cout << jacpat[i][j] << " "; + std::cout << "\n"; } - cout.flush(); } t0 = z2 - z1; @@ -278,10 +230,10 @@ int main(void) { maxnz = jacpat[i][0]; } nz_rel = (int)ceil(100 * nz / ((double)depen * indep)); - cout << nz << " non-zero Jacobian elements of total " << depen * indep - << " elements <= " << nz_rel << "%\n"; - cout << "min " << minnz << " non-zeros per row; max " << maxnz - << " non-zeros per row;\n"; + std::cout << nz << " non-zero Jacobian elements of total " << depen * indep + << " elements <= " << nz_rel << "%\n"; + std::cout << "min " << minnz << " non-zeros per row; max " << maxnz + << " non-zeros per row;\n"; nzref = nz; for (i = 0; i < depen; i++) @@ -292,7 +244,7 @@ int main(void) { /*--------------------------------------------------------------------------*/ /* Jacobian pattern by index domains, tight */ - cout << "\nJacobian Pattern by index domain propagation, tight ...\n"; + std::cout << "\nJacobian Pattern by index domain propagation, tight ...\n"; jacpat = new unsigned int *[depen]; diff --git a/ADOL-C/examples/additional_examples/sparse/sparse_hessian.cpp b/ADOL-C/examples/additional_examples/sparse/sparse_hessian.cpp index 9721226b..0e61f4ce 100644 --- a/ADOL-C/examples/additional_examples/sparse/sparse_hessian.cpp +++ b/ADOL-C/examples/additional_examples/sparse/sparse_hessian.cpp @@ -13,24 +13,64 @@ ---------------------------------------------------------------------------*/ +#include +#include #include #include -#include -#include -#include +/***************************************************************************/ + +double feval(double *x) { + double res; + + res = 0.5 * (x[0] - 1) * (x[0] - 1) + 0.8 * (x[1] - 2) * (x[1] - 2) + + 0.9 * (x[2] - 3) * (x[2] - 3); + res += 5 * x[0] * x[1]; + res += cos(x[3]); + res += sin(x[4]) * pow(x[1], 2); + res += exp(x[5]) * x[2]; + res += sin(x[4] * x[5]); + + return res; +} + +/***************************************************************************/ + +adouble feval_ad(adouble *x) { + adouble res; + + res = 0.5 * (x[0] - 1) * (x[0] - 1) + 0.8 * (x[1] - 2) * (x[1] - 2) + + 0.9 * (x[2] - 2) * (x[2] - 2); + res += 5 * x[0] * x[1]; + res += cos(x[3]); + res += sin(x[4]) * x[1] * x[1]; + res += exp(x[5]) * x[2]; + res += sin(x[4] * x[5]); -#define tag 1 + return res; +} -double feval(double *x); -adouble feval_ad(adouble *x); +/***************************************************************************/ -void printmat(const char *kette, int n, int m, double **M); +void printmat(const char *name, int m, int n, double **M) { + int i, j; + + printf("%s \n", name); + for (i = 0; i < m; i++) { + printf("\n %d: ", i); + for (j = 0; j < n; j++) + printf(" %10.4f ", M[i][j]); + } + printf("\n"); +} int main() { - int n = 6; - double f, x[6]; - adouble fad, xad[6]; + + const short tapeId = 1; + createNewTape(tapeId); + constexpr size_t dim = 6; + double f, x[dim]; + adouble fad, xad[dim]; int i, j; @@ -38,13 +78,13 @@ int main() { /******* function evaluation ***************/ /****************************************************************************/ - for (i = 0; i < n; i++) + for (i = 0; i < dim; i++) x[i] = log(1.0 + i); /* Tracing of function f(x) */ - trace_on(tag); - for (i = 0; i < n; i++) + trace_on(tapeId); + for (i = 0; i < dim; i++) xad[i] <<= x[i]; fad = feval_ad(xad); @@ -59,11 +99,11 @@ int main() { /****************************************************************************/ double **H; - H = myalloc2(n, n); + H = myalloc2(dim, dim); - hessian(tag, n, x, H); + hessian(tapeId, dim, x, H); - printmat(" H", n, n, H); + printmat(" H", dim, dim, H); printf("\n"); /****************************************************************************/ @@ -75,12 +115,9 @@ int main() { unsigned int *cind = NULL; double *values = NULL; int nnz; - int options[2]; - - options[0] = 0; /* safe mode (default) */ - options[1] = 0; /* indirect recovery (default) */ - - sparse_hess(tag, n, 0, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_hess( + tapeId, dim, 0, x, &nnz, &rind, &cind, &values); printf("In sparse format:\n"); for (i = 0; i < nnz; i++) @@ -93,10 +130,9 @@ int main() { free(values); values = NULL; - options[0] = 0; /* safe mode (default) */ - options[1] = 1; /* direct recovery */ - - sparse_hess(tag, n, 0, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_hess( + tapeId, dim, 0, x, &nnz, &rind, &cind, &values); printf("In sparse format:\n"); for (i = 0; i < nnz; i++) @@ -113,19 +149,21 @@ int main() { /* change value of x, but not the sparsity pattern */ /*--------------------------------------------------------------------------*/ - for (i = 0; i < n; i++) + for (i = 0; i < dim; i++) x[i] = 2.0 * i; /* For comparisons: Full Hessian: */ - hessian(tag, n, x, H); + hessian(tapeId, dim, x, H); - printmat(" H", n, n, H); + printmat(" H", dim, dim, H); printf("\n"); /* repeated call of sparse_hess with same sparsity pattern => repeat = 1 */ - sparse_hess(tag, n, 0, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_hess( + tapeId, dim, 0, x, &nnz, &rind, &cind, &values); printf("In sparse format:\n"); for (i = 0; i < nnz; i++) @@ -146,17 +184,14 @@ int main() { /* sparsity pattern Hessian */ /*--------------------------------------------------------------------------*/ - unsigned int **HP = NULL; /* compressed block row storage */ - int ctrl; - - HP = (unsigned int **)malloc(n * sizeof(unsigned int *)); - ctrl = 0; - - hess_pat(tag, n, x, HP, ctrl); + std::vector HP(dim); /* compressed block row storage */ + std::span HP_(HP); + ADOLC::Sparse::hess_pat(tapeId, dim, x, + HP_); printf("\n"); printf("Sparsity pattern of Hessian: \n"); - for (i = 0; i < n; i++) { + for (i = 0; i < dim; i++) { printf(" %d: ", i); for (j = 1; j <= (int)HP[i][0]; j++) printf(" %d ", HP[i][j]); @@ -170,14 +205,11 @@ int main() { double **Seed; int p; - int option = 1; - /* option = 0 indirect recovery (default), - option = 1 direct recovery */ + ADOLC::Sparse::generate_seed_hess( + dim, HP_, &Seed, &p); - generate_seed_hess(n, HP, &Seed, &p, option); - - printmat(" Seed matrix", n, p, Seed); + printmat(" Seed matrix", dim, p, Seed); printf("\n"); /*--------------------------------------------------------------------------*/ @@ -185,86 +217,40 @@ int main() { /*--------------------------------------------------------------------------*/ double **Hcomp; - Hcomp = myalloc2(n, p); + Hcomp = myalloc2(dim, p); - hess_mat(tag, n, p, x, Seed, Hcomp); + hess_mat(tapeId, dim, p, x, Seed, Hcomp); - printmat("compressed H:", n, p, Hcomp); + printmat("compressed H:", dim, p, Hcomp); printf("\n"); /*--------------------------------------------------------------------------*/ /* change value of x, but not the sparsity pattern */ /*--------------------------------------------------------------------------*/ - for (i = 0; i < n; i++) + for (i = 0; i < dim; i++) x[i] = 2.0 * i; /* For comparisons: Full Hessian */ - hessian(tag, n, x, H); + hessian(tapeId, dim, x, H); - printmat(" H", n, n, H); + printmat(" H", dim, dim, H); printf("\n"); - hess_mat(tag, n, p, x, Seed, Hcomp); + hess_mat(tapeId, dim, p, x, Seed, Hcomp); - printmat("compressed H:", n, p, Hcomp); + printmat("compressed H:", dim, p, Hcomp); printf("\n"); - for (i = 0; i < n; i++) - free(HP[i]); - free(HP); - + for (i = 0; i < dim; i++) { + delete[] HP[i]; + HP[i] = nullptr; + } myfree2(H); myfree2(Hcomp); - for (i = 0; i < n; i++) + for (i = 0; i < dim; i++) delete[] Seed[i]; delete[] Seed; } - -/***************************************************************************/ - -double feval(double *x) { - double res; - - res = 0.5 * (x[0] - 1) * (x[0] - 1) + 0.8 * (x[1] - 2) * (x[1] - 2) + - 0.9 * (x[2] - 3) * (x[2] - 3); - res += 5 * x[0] * x[1]; - res += cos(x[3]); - res += sin(x[4]) * pow(x[1], 2); - res += exp(x[5]) * x[2]; - res += sin(x[4] * x[5]); - - return res; -} - -/***************************************************************************/ - -adouble feval_ad(adouble *x) { - adouble res; - - res = 0.5 * (x[0] - 1) * (x[0] - 1) + 0.8 * (x[1] - 2) * (x[1] - 2) + - 0.9 * (x[2] - 2) * (x[2] - 2); - res += 5 * x[0] * x[1]; - res += cos(x[3]); - res += sin(x[4]) * x[1] * x[1]; - res += exp(x[5]) * x[2]; - res += sin(x[4] * x[5]); - - return res; -} - -/***************************************************************************/ - -void printmat(const char *name, int m, int n, double **M) { - int i, j; - - printf("%s \n", name); - for (i = 0; i < m; i++) { - printf("\n %d: ", i); - for (j = 0; j < n; j++) - printf(" %10.4f ", M[i][j]); - } - printf("\n"); -} diff --git a/ADOL-C/examples/additional_examples/sparse/sparse_jacobian.cpp b/ADOL-C/examples/additional_examples/sparse/sparse_jacobian.cpp index fdcfda8b..7e09a692 100644 --- a/ADOL-C/examples/additional_examples/sparse/sparse_jacobian.cpp +++ b/ADOL-C/examples/additional_examples/sparse/sparse_jacobian.cpp @@ -13,24 +13,51 @@ ---------------------------------------------------------------------------*/ +#include #include #include #include #include -#include -#define tag 1 +/***************************************************************************/ + +void ceval(double *x, double *c) { + c[0] = 2 * x[0] + x[1] - 2.0; + c[1] = x[2] * x[2] + x[3] * x[3] - 2.0; + c[2] = 3 * x[4] * x[5] - 3.0; +} -void ceval_ad(adouble *x, adouble *c); -void ceval(double *x, double *c); +/***************************************************************************/ + +void ceval_ad(adouble *x, adouble *c) { + c[0] = 2 * x[0] + x[1] - 2.0; + c[0] += cos(x[3]) * sin(x[4]); + c[1] = x[2] * x[2] + x[3] * x[3] - 2.0; + c[2] = 3 * x[4] * x[5] - 3.0 + sin(x[4] * x[5]); +} + +/***************************************************************************/ + +void printmat(const char *name, int m, int n, double **M) { + int i, j; -void printmat(const char *name, int n, int m, double **M); + printf("%s \n", name); + for (i = 0; i < m; i++) { + printf("\n %d: ", i); + for (j = 0; j < n; j++) + printf(" %10.4f ", M[i][j]); + } + printf("\n"); +} int main() { - int n = 6, m = 3; - double x[6], c[3]; - adouble xad[6], cad[3]; + const short tapeId = 1; + createNewTape(tapeId); + constexpr size_t dimIn = 6; + constexpr size_t dimOut = 3; + double x[dimIn], c[dimOut]; + adouble xad[dimIn], cad[dimOut]; int i, j; @@ -38,23 +65,23 @@ int main() { /******* function evaluation ***************/ /****************************************************************************/ - for (i = 0; i < n; i++) + for (i = 0; i < dimIn; i++) x[i] = log(1.0 + i); /* Tracing of function c(x) */ - trace_on(tag); - for (i = 0; i < n; i++) + trace_on(tapeId); + for (i = 0; i < dimIn; i++) xad[i] <<= x[i]; ceval_ad(xad, cad); - for (i = 0; i < m; i++) + for (i = 0; i < dimOut; i++) cad[i] >>= c[i]; trace_off(); printf("\n c = "); - for (j = 0; j < m; j++) + for (j = 0; j < dimOut; j++) printf(" %e ", c[j]); printf("\n"); @@ -63,11 +90,11 @@ int main() { /****************************************************************************/ double **J; - J = myalloc2(m, n); + J = myalloc2(dimOut, dimIn); - jacobian(tag, m, n, x, J); + jacobian(tapeId, dimOut, dimIn, x, J); - printmat(" J", m, n, J); + printmat(" J", dimOut, dimIn, J); printf("\n"); /****************************************************************************/ @@ -79,14 +106,9 @@ int main() { unsigned int *cind = NULL; /* column indices */ double *values = NULL; /* values */ int nnz; - int options[4]; - - options[0] = 0; /* sparsity pattern by index domains (default) */ - options[1] = 0; /* safe mode (default) */ - options[2] = 0; /* not required if options[0] = 0 */ - options[3] = 0; /* column compression (default) */ - sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_jac(tapeId, dimOut, dimIn, 0, x, &nnz, &rind, &cind, + &values); printf("In sparse format:\n"); for (i = 0; i < nnz; i++) @@ -102,11 +124,10 @@ int main() { /* same approach but using row compression */ /*--------------------------------------------------------------------------*/ - options[3] = 1; /* row compression => reverse mode, */ - /* sometimes better than forward mode */ - /* due to sparsity structure */ - - sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 0, x, &nnz, &rind, &cind, &values); printf("In sparse format (using row compression): \n"); for (i = 0; i < nnz; i++) @@ -122,39 +143,46 @@ int main() { /* change value of x, but not the sparsity pattern */ /*--------------------------------------------------------------------------*/ - for (i = 0; i < n; i++) + for (i = 0; i < dimIn; i++) x[i] = 2.0 * i; /* For comparisons: Full Jacobian */ - jacobian(tag, m, n, x, J); + jacobian(tapeId, dimOut, dimIn, x, J); - printmat(" J", m, n, J); + printmat(" J", dimOut, dimIn, J); printf("\n"); /* repeated call of sparse_jac with same sparsity pattern => repeat = 1 */ - sparse_jac(tag, m, n, 1, x, &nnz, &rind, &cind, &values, options); + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 1, x, &nnz, &rind, &cind, &values); printf("In sparse format:\n"); for (i = 0; i < nnz; i++) printf("%2d %2d %10.6f\n\n", rind[i], cind[i], values[i]); - free(rind); - rind = NULL; - free(cind); - cind = NULL; - free(values); - values = NULL; /*--------------------------------------------------------------------------*/ /* same approach but using row compression */ /*--------------------------------------------------------------------------*/ - options[3] = 1; /* row compression => reverse mode, */ - /* sometimes better than forward mode */ - /* due to sparsity structure */ + ADOLC::Sparse::sparse_jac( + tapeId, dimOut, dimIn, 0, x, &nnz, &rind, &cind, &values); + + printf("In sparse format (using row compression): \n"); + for (i = 0; i < nnz; i++) + printf("%2d %2d %10.6f\n\n", rind[i], cind[i], values[i]); - sparse_jac(tag, m, n, 0, x, &nnz, &rind, &cind, &values, options); + constexpr int dimOutCT = 3; + constexpr int dimInCT = 6; + ADOLC::Sparse::sparse_jac( + tapeId, dimOutCT, dimInCT, 0, x, &nnz, &rind, &cind, &values); printf("In sparse format (using row compression): \n"); for (i = 0; i < nnz; i++) @@ -174,19 +202,15 @@ int main() { /* sparsity pattern Jacobian */ /*--------------------------------------------------------------------------*/ - unsigned int **JP = NULL; /* compressed block row storage */ - int ctrl[3]; - - JP = (unsigned int **)malloc(m * sizeof(unsigned int *)); - ctrl[0] = 0; - ctrl[1] = 0; - ctrl[2] = 0; - - jac_pat(tag, m, n, x, JP, ctrl); + std::vector JP(dimOut); /* compressed block row storage */ + std::span JP_(JP); + ADOLC::Sparse::jac_pat(tapeId, dimOut, + dimIn, x, JP_); printf("\n"); printf("Sparsity pattern of Jacobian: \n"); - for (i = 0; i < m; i++) { + for (i = 0; i < dimOut; i++) { printf(" %d: ", i); for (j = 1; j <= (int)JP[i][0]; j++) printf(" %d ", JP[i][j]); @@ -200,15 +224,11 @@ int main() { double **Seed; int p; - int option = 0; - - /* option = 0 column compression (default), - option = 1 rom compression */ - - generate_seed_jac(m, n, JP, &Seed, &p, option); + ADOLC::Sparse::generate_seed_jac( + dimOut, dimIn, JP_, &Seed, &p); printf(" p_J = %d \n", p); - printmat(" Seed matrix", n, p, Seed); + printmat(" Seed matrix", dimIn, p, Seed); printf("\n"); /*--------------------------------------------------------------------------*/ @@ -216,69 +236,37 @@ int main() { /*--------------------------------------------------------------------------*/ double **Jcomp; - Jcomp = myalloc2(m, p); + Jcomp = myalloc2(dimOut, p); - fov_forward(tag, m, n, p, x, Seed, c, Jcomp); - printmat("compressed J:", m, p, Jcomp); + fov_forward(tapeId, dimOut, dimIn, p, x, Seed, c, Jcomp); + printmat("compressed J:", dimOut, p, Jcomp); printf("\n"); /*--------------------------------------------------------------------------*/ /* change value of x, but not the sparsity pattern */ /*--------------------------------------------------------------------------*/ - for (i = 0; i < n; i++) + for (i = 0; i < dimIn; i++) x[i] = 2.0 * i; /* For comparisons: Full Jacobian */ - jacobian(tag, m, n, x, J); + jacobian(tapeId, dimOut, dimIn, x, J); - printmat(" J", m, n, J); + printmat(" J", dimOut, dimIn, J); printf("\n"); - fov_forward(tag, m, n, p, x, Seed, c, Jcomp); - printmat("compressed J:", m, p, Jcomp); + fov_forward(tapeId, dimOut, dimIn, p, x, Seed, c, Jcomp); + printmat("compressed J:", dimOut, p, Jcomp); printf("\n"); - for (i = 0; i < m; i++) - free(JP[i]); - free(JP); + for (i = 0; i < dimOut; i++) + delete[] JP[i]; myfree2(J); - for (i = 0; i < n; i++) + for (i = 0; i < dimIn; i++) delete[] Seed[i]; delete[] Seed; myfree2(Jcomp); } - -/***************************************************************************/ - -void ceval(double *x, double *c) { - c[0] = 2 * x[0] + x[1] - 2.0; - c[1] = x[2] * x[2] + x[3] * x[3] - 2.0; - c[2] = 3 * x[4] * x[5] - 3.0; -} - -/***************************************************************************/ - -void ceval_ad(adouble *x, adouble *c) { - c[0] = 2 * x[0] + x[1] - 2.0; - c[0] += cos(x[3]) * sin(x[4]); - c[1] = x[2] * x[2] + x[3] * x[3] - 2.0; - c[2] = 3 * x[4] * x[5] - 3.0 + sin(x[4] * x[5]); -} - -/***************************************************************************/ - -void printmat(const char *name, int m, int n, double **M) { - int i, j; - - printf("%s \n", name); - for (i = 0; i < m; i++) { - printf("\n %d: ", i); - for (j = 0; j < n; j++) - printf(" %10.4f ", M[i][j]); - } - printf("\n"); -} diff --git a/ADOL-C/examples/traceless_vector_indo.cpp b/ADOL-C/examples/traceless_vector_indo.cpp index f4e7129b..f992aeb0 100644 --- a/ADOL-C/examples/traceless_vector_indo.cpp +++ b/ADOL-C/examples/traceless_vector_indo.cpp @@ -38,7 +38,7 @@ int main(int argc, char *argv[]) { adtl::setNumDir(n); my_function fun; - my_function fun_indo; + my_function fun_indo; adtl::adouble x[n], y[m]; for (int i = 0; i < n; ++i) // Initialize x_i @@ -82,8 +82,8 @@ int main(int argc, char *argv[]) { int nnz; unsigned int *rind = NULL, *cind = NULL; double *values = NULL; - retVal = ::ADOLC_get_sparse_jacobian(&fun, &fun_indo, n, m, 0, basepoints, - &nnz, &rind, &cind, &values); + retVal = ADOLC::Sparse::adtl_indo::ADOLC_get_sparse_jacobian( + &fun, &fun_indo, n, m, 0, basepoints, &nnz, &rind, &cind, &values); cout << endl; cout << "Checking results with ADOLC_get_sparse_jacobian functionality..." diff --git a/ADOL-C/include/adolc/CMakeLists.txt b/ADOL-C/include/adolc/CMakeLists.txt index ebd57ece..cbd58276 100644 --- a/ADOL-C/include/adolc/CMakeLists.txt +++ b/ADOL-C/include/adolc/CMakeLists.txt @@ -5,7 +5,6 @@ install(FILES adolcexport.h adolc.h adolc_openmp.h - adolc_sparse.h adoublecuda.h adtl.h adtl_hov.h @@ -35,5 +34,6 @@ install(FILES add_subdirectory(drivers) add_subdirectory(internal) add_subdirectory(lie) +add_subdirectory(sparse) add_subdirectory(tapedoc) add_subdirectory(valuetape) diff --git a/ADOL-C/include/adolc/adolc.h b/ADOL-C/include/adolc/adolc.h index 95f059db..b85ce0ec 100644 --- a/ADOL-C/include/adolc/adolc.h +++ b/ADOL-C/include/adolc/adolc.h @@ -53,8 +53,8 @@ #include /*--------------------------------------------------------------------------*/ -/* interfaces to SPARSE package */ -#if defined(SPARSE_DRIVERS) +/* interfaces to sparse drivers */ +#ifdef ADOLC_SPARSE #include #include #endif diff --git a/ADOL-C/include/adolc/adolc_sparse.h b/ADOL-C/include/adolc/adolc_sparse.h deleted file mode 100644 index 39cb3f3c..00000000 --- a/ADOL-C/include/adolc/adolc_sparse.h +++ /dev/null @@ -1,16 +0,0 @@ -/*---------------------------------------------------------------------------- - ADOL-C -- Automatic Differentiation by Overloading in C++ - File: adolc_sparse.h - Revision: $Id$ - Contents: Provides C/C++ interfaces of ADOL-C sprase drivers. - - Copyright (c) Andrea Walther - - This file is part of ADOL-C. This software is provided as open source. - Any use, reproduction, or distribution of the software constitutes - recipient's acceptance of the terms of the accompanying license file. - -----------------------------------------------------------------------------*/ - -#include -#include diff --git a/ADOL-C/include/adolc/adolcerror.h b/ADOL-C/include/adolc/adolcerror.h index 992853aa..ae0903ec 100644 --- a/ADOL-C/include/adolc/adolcerror.h +++ b/ADOL-C/include/adolc/adolcerror.h @@ -233,8 +233,9 @@ enum class ADOLC_API ErrorType : size_t { SPARSE_HESS_IND, SPARSE_CRS, SPARSE_JAC_MALLOC, - SPARSE_JAC_NO_BP + SPARSE_JAC_NO_BP, + NOT_IMPLEMENTED }; // wrapper for information of errors diff --git a/ADOL-C/include/adolc/adtl_hov.h b/ADOL-C/include/adolc/adtl_hov.h index fd565535..8587d00a 100644 --- a/ADOL-C/include/adolc/adtl_hov.h +++ b/ADOL-C/include/adolc/adtl_hov.h @@ -2598,4 +2598,4 @@ inline size_t adouble::get_pattern_size() const { } } // namespace adtl_hov -#endif +#endif \ No newline at end of file diff --git a/ADOL-C/include/adolc/adtl_indo.h b/ADOL-C/include/adolc/adtl_indo.h index 7b3e7d05..c651741a 100644 --- a/ADOL-C/include/adolc/adtl_indo.h +++ b/ADOL-C/include/adolc/adtl_indo.h @@ -18,11 +18,14 @@ #define ADOLC_ADTL_INDO_H #include +#include +#include +#include +#include #include #include #include - -#include +#include using std::istream; using std::list; @@ -34,13 +37,12 @@ template class func_ad { virtual int operator()(int n, T *x, int m, T *y) = 0; }; -namespace adtl_indo { +namespace ADOLC::Sparse::adtl_indo { class adouble; ADOLC_API int ADOLC_Init_sparse_pattern(adouble *a, int n, unsigned int start_cnt); ADOLC_API int ADOLC_get_sparse_pattern(const adouble *const b, int m, - unsigned int **&pat); -} // namespace adtl_indo + std::vector &pat); ADOLC_API int ADOLC_get_sparse_jacobian(func_ad *const func, @@ -49,8 +51,6 @@ ADOLC_get_sparse_jacobian(func_ad *const func, unsigned int **rind, unsigned int **cind, double **values); -namespace adtl_indo { - ADOLC_API double makeNaN(); ADOLC_API double makeInf(); @@ -210,7 +210,7 @@ class ADOLC_API adouble { ADOLC_API friend int ADOLC_Init_sparse_pattern(adouble *a, int n, unsigned int start_cnt); ADOLC_API friend int ADOLC_get_sparse_pattern(const adouble *const b, int m, - unsigned int **&pat); + std::vector &pat); /******************* i/o operations *********************************/ ADOLC_API friend ostream &operator<<(ostream &, const adouble &); ADOLC_API friend istream &operator>>(istream &, adouble &); @@ -220,14 +220,6 @@ class ADOLC_API adouble { list pattern; }; -} // namespace adtl_indo - -#include -#include -#include - -namespace adtl_indo { - #if defined(HAVE_BUILTIN_EXPECT) && HAVE_BUILTIN_EXPECT #define likely(x) __builtin_expect(!!(x), 1) #define unlikely(x) __builtin_expect(!!(x), 0) @@ -881,6 +873,6 @@ inline size_t adouble::get_pattern_size() const { return s; } -} // namespace adtl_indo +} // namespace ADOLC::Sparse::adtl_indo #endif diff --git a/ADOL-C/include/adolc/interfaces.h b/ADOL-C/include/adolc/interfaces.h index 0e61c9a3..ac0e3bac 100644 --- a/ADOL-C/include/adolc/interfaces.h +++ b/ADOL-C/include/adolc/interfaces.h @@ -49,11 +49,6 @@ #include #include -#if defined(SPARSE) -#include -#include -#endif - /****************************************************************************/ /****************************************************************************/ /* Now the C++ THINGS */ @@ -295,15 +290,16 @@ ADOLC_API fint hov_wk_forward_(const fint *, const fint *, const fint *, /* INT_FOR, SAFE */ /* int_forward_safe(tag, m, n, p, X[n][p], Y[m][p]) */ -ADOLC_API int int_forward_safe(short, int, int, int, const size_t *const *, - size_t **); +ADOLC_API int int_forward_safe(short, int, int, int, const bitword_t *const *, + bitword_t **); /*--------------------------------------------------------------------------*/ /* INT_FOR, TIGHT */ /* int_forward_tight(tag, m, n, p, x[n], X[n][p], y[m], Y[m][p]) */ ADOLC_API int int_forward_tight(short, int, int, int, const double *, - const size_t *const *, double *, size_t **); + const bitword_t *const *, double *, + bitword_t **); /****************************************************************************/ /* INDEX DOMAIN UTILITIES */ @@ -311,15 +307,13 @@ ADOLC_API int int_forward_tight(short, int, int, int, const double *, /* INDOPRO, SAFE */ /* indopro_forward_safe(tag, m, n, p, x[n], *Y[m]) */ -ADOLC_API int indopro_forward_safe(short, int, int, const double *, - unsigned int **); +ADOLC_API int indopro_forward_safe(short, int, int, const double *, uint **); /*--------------------------------------------------------------------------*/ /* INDOPRO, TIGHT */ /* indopro_forward_tight(tag, m, n, x[n], *Y[m]) */ -ADOLC_API int indopro_forward_tight(short, int, int, const double *, - unsigned int **); +ADOLC_API int indopro_forward_tight(short, int, int, const double *, uint **); /****************************************************************************/ /* NONLINEAR INDEX DOMAIN UTILITIES */ @@ -327,28 +321,26 @@ ADOLC_API int indopro_forward_tight(short, int, int, const double *, /* INDOPRO, SAFE */ /* nonl_ind_forward_safe(tag, m, n, p, x[n], *Y[m]) */ -ADOLC_API int nonl_ind_forward_safe(short, int, int, const double *, - unsigned int **); +ADOLC_API int nonl_ind_forward_safe(short, int, int, const double *, uint **); /*--------------------------------------------------------------------------*/ /* INDOPRO, TIGHT */ /* nonl_ind_forward_tight(tag, m, n, x[n], *Y[m]) */ -ADOLC_API int nonl_ind_forward_tight(short, int, int, const double *, - unsigned int **); +ADOLC_API int nonl_ind_forward_tight(short, int, int, const double *, uint **); /*--------------------------------------------------------------------------*/ /* INDOPRO, SAFE */ /* nonl_ind_old_forward_safe(tag, m, n, p, x[n], *Y[m]) */ ADOLC_API int nonl_ind_old_forward_safe(short, int, int, const double *, - unsigned int **); + uint **); /*--------------------------------------------------------------------------*/ /* INDOPRO, TIGHT */ /* nonl_ind_old_forward_tight(tag, m, n, x[n], *Y[m]) */ ADOLC_API int nonl_ind_old_forward_tight(short, int, int, const double *, - unsigned int **); + uint **); /****************************************************************************/ /* REVERSE MODE */ @@ -438,15 +430,15 @@ ADOLC_API fint hov_ti_reverse_(const fint *, const fint *, const fint *, /* INT_REV, TIGHT */ /* int_reverse_tight(tag, m, n, q, U[q][m], Z[q][n]) */ -ADOLC_API int int_reverse_tight(short, int, int, int, const size_t *const *, - size_t **); +ADOLC_API int int_reverse_tight(short, int, int, int, const bitword_t *const *, + bitword_t **); /*--------------------------------------------------------------------------*/ /* INT_REV, SAFE */ /* int_reverse_safe(tag, m, n, q, U[q][m], Z[q][n]) */ -ADOLC_API int int_reverse_safe(short, int, int, int, const size_t *const *, - size_t **); +ADOLC_API int int_reverse_safe(short, int, int, int, const bitword_t *const *, + bitword_t **); /*--------------------------------------------------------------------------*/ ADOLC_API int get_num_switches(short); @@ -469,7 +461,7 @@ ADOLC_API int fov_pl_sig_forward(short, int, int, int, const double *, const short *, double *, double **, double *, double **, short *); ADOLC_API int indopro_forward_absnormal(short, int, int, int, const double *, - unsigned int **); + uint **); /*--------------------------------------------------------------------------*/ ADOLC_API int fos_pl_reverse(short, int, int, int, int, double *); ADOLC_API int fos_pl_sig_reverse(short, int, int, int, const short *, diff --git a/ADOL-C/include/adolc/internal/adolc_settings.h.in b/ADOL-C/include/adolc/internal/adolc_settings.h.in index 4cf2d571..797b581b 100644 --- a/ADOL-C/include/adolc/internal/adolc_settings.h.in +++ b/ADOL-C/include/adolc/internal/adolc_settings.h.in @@ -37,7 +37,7 @@ typedef double revreal; /*--------------------------------------------------------------------------*/ /* Sparse drivers have been compiled */ -@SPARSE_DRIVERS@ +@ADOLC_SPARSE@ /*--------------------------------------------------------------------------*/ /* Use Boost Library Pool allocator */ diff --git a/ADOL-C/include/adolc/internal/common.h b/ADOL-C/include/adolc/internal/common.h index ed069f4f..84ead07e 100644 --- a/ADOL-C/include/adolc/internal/common.h +++ b/ADOL-C/include/adolc/internal/common.h @@ -15,20 +15,27 @@ #if !defined(ADOLC_COMMON_H) #define ADOLC_COMMON_H 1 - -#include /*--------------------------------------------------------------------------*/ /* standard includes */ #if !defined(__cplusplus) #include #include #else -#include +#include #include #endif /*--------------------------------------------------------------------------*/ /* type definitions */ + +// define bit vector propagation type used in sparse drivers +#if ADOLC_BITWORD_BITS == 32 +typedef uint32_t bitword_t; +#elif ADOLC_BITWORD_BITS == 64 +typedef uint64_t bitword_t; +#else +typedef uint32_t bitword_t; +#endif typedef unsigned int uint; /*--------------------------------------------------------------------------*/ @@ -78,6 +85,7 @@ typedef unsigned int uint; #include #if defined(__cplusplus) + #define BEGIN_C_DECLS extern "C" { #define END_C_DECLS } #else diff --git a/ADOL-C/include/adolc/sparse/CMakeLists.txt b/ADOL-C/include/adolc/sparse/CMakeLists.txt new file mode 100644 index 00000000..e94bbd19 --- /dev/null +++ b/ADOL-C/include/adolc/sparse/CMakeLists.txt @@ -0,0 +1,4 @@ +install(FILES + sparse_fo_rev.h + sparsedrivers.h + DESTINATION "include/adolc/sparse") diff --git a/ADOL-C/include/adolc/sparse/sparse_fo_rev.h b/ADOL-C/include/adolc/sparse/sparse_fo_rev.h index 5de4e5dc..c323d3ef 100644 --- a/ADOL-C/include/adolc/sparse/sparse_fo_rev.h +++ b/ADOL-C/include/adolc/sparse/sparse_fo_rev.h @@ -13,13 +13,13 @@ package. recipient's acceptance of the terms of the accompanying license file. ----------------------------------------------------------------------------*/ -#if !defined(ADOLC_SPARSE_SPARSE_H) -#define ADOLC_SPARSE_SPARSE_H 1 - +#ifndef ADOLC_SPARSE_FO_REV_H +#define ADOLC_SPARSE_FO_REV_H #include #include +#include -#if defined(__cplusplus) +namespace ADOLC::Sparse { /****************************************************************************/ /* FORWARD MODE, overloaded calls */ /* */ @@ -37,38 +37,40 @@ package. /* */ /* forward(tag, m, n, p, x[n], X[n][p], y[m], Y[m][p], mode) : intfov */ -ADOLC_API int forward(short, int, int, int, double *, size_t **, double *, - size_t **, char = 0); +ADOLC_API int forward(short tag, int m, int n, int p, double *x, + std::vector &X, double *y, + std::vector &Y, char mode); /*--------------------------------------------------------------------------*/ /* Bit pattern propagation call, d = 1, safe version (no x[] and y[]) */ /* */ /* forward(tag, m, n, p, X[n][p], Y[m][p], mode) : intfov */ -ADOLC_API int forward(short, int, int, int, size_t **, size_t **, char = 0); +ADOLC_API int forward(short tag, int m, int n, int p, + std::vector &X, std::vector &Y, + char mode); /****************************************************************************/ -/* REVERSE MODE, overloaded calls */ +/* REVERSE MODE, overloaded calls + */ /* */ -/* nBV = number of Boolean Vectors to be packed */ -/* (see Chapter Dependence Analysis, ADOL-C Documentation) */ +/* nBV = number of Boolean Vectors to be packed */ +/* (see Chapter Dependence Analysis, ADOL-C Documentation) */ /* bits_per_long = 8*sizeof(size_t) */ -/* q = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 ) */ +/* q = nBV / bits_per_long + ( (nBV % bits_per_long) != 0 ) */ /* */ -/* For the full Jacobian matrix set */ -/* q = depen / bits_per_long + ((depen % bits_per_long) != 0) */ -/* and pass a bit pattern version of the identity matrix as an argument */ +/* For the full Jacobian matrix set */ +/* q = depen / bits_per_long + ((depen % bits_per_long) != 0) */ +/* and pass a bit pattern version of the identity matrix as an argument */ /* */ /*--------------------------------------------------------------------------*/ /* */ -/* Bit pattern propagation call, d = 0, tight & safe version */ +/* Bit pattern propagation call, d = 0, tight & safe version */ /* */ -/* reverse(tag, m, n, q, U[q][m], Z[q][n], mode) : intfov */ - -ADOLC_API int reverse(short, int, int, int, size_t **, size_t **, char = 0); - -#endif - -/****************************************************************************/ +/* reverse(tag, m, n, q, U[q][m], Z[q][n], mode) : intfov */ -#endif +ADOLC_API int reverse(short tag, int m, int n, int q, + std::vector &U, std::vector &Z, + char mode); +} // namespace ADOLC::Sparse +#endif // ADOLC_SPARSE_FO_REV_H \ No newline at end of file diff --git a/ADOL-C/include/adolc/sparse/sparsedrivers.h b/ADOL-C/include/adolc/sparse/sparsedrivers.h index 903b7ce0..31adb87a 100644 --- a/ADOL-C/include/adolc/sparse/sparsedrivers.h +++ b/ADOL-C/include/adolc/sparse/sparsedrivers.h @@ -12,93 +12,1222 @@ package. recipient's acceptance of the terms of the accompanying license file. ----------------------------------------------------------------------------*/ -#if !defined(ADOLC_SPARSE_SPARSE_H) -#define ADOLC_SPARSE_SPARSE_H 1 - -#include +#ifndef ADOLC_SPARSE_DRIVERS_H +#define ADOLC_SPARSE_DRIVERS_H +#include +#include +#include #include +#include +#include +#include +#include +#include +#include -BEGIN_C_DECLS +// Max. number of unsigned ints to store the seed / jacobian matrix strips. +// Reduce this value to x if your system happens to run out of memory. x < 10 +// makes no sense. x = 50 or 100 is better. x stays for (x * sizeof(size_t) * 8) +// (block) variables at once +#define PQ_STRIPMINE_MAX 30 -/****************************************************************************/ +namespace ADOLC::Sparse { -/*--------------------------------------------------------------------------*/ -/* jacobian pattern */ -/* jac_pat(tag, m, n, argument, */ -/* crs[] [ crs[][0] = non-zero independent blocks per row ], */ -/* options[3]) */ -/* */ - -ADOLC_API int jac_pat(short, int, int, const double *, unsigned int **, int *); - -/*--------------------------------------------------------------------------*/ -/* abs-normal jacobian pattern */ -/* absnormal_jac_pat(tag, m, n, s, argument, */ -/* crs[] [ crs[][0] = non-zero independent blocks per row ]) */ -/* */ -ADOLC_API int absnormal_jac_pat(short, int, int, int, const double *, - unsigned int **); -/*--------------------------------------------------------------------------*/ -/* seed matrix for sparse jacobian */ -/* generate_seed_jac(m, n, crs, &seed, &p, option); */ - -ADOLC_API void generate_seed_jac(int, int, unsigned int **, double ***, int *, - int); - -/*--------------------------------------------------------------------------*/ -/* sparse jacobian */ -/* int sparse_jac(tag, m, n, repeat, x, &nnz, &row_ind, &col_ind, &values, */ -/* options[3]); */ - -ADOLC_API int sparse_jac(short, int, int, int, const double *, int *, - unsigned int **, unsigned int **, double **, int *); - -/*--------------------------------------------------------------------------*/ -/* hessian pattern */ -/* hess_pat(tag, n, x[n], crs[n][*], option) */ -/* */ -/* crs[i][ crs[i][0] = non-zero entries per row ] */ -/* */ - -ADOLC_API int hess_pat(short, int, const double *, unsigned int **, int); - -/*--------------------------------------------------------------------------*/ -/* seed matrix for sparse hessian */ -/* generate_seed_hess(n, crs, &seed, &p, option); */ - -ADOLC_API void generate_seed_hess(int, unsigned int **, double ***, int *, int); - -/*--------------------------------------------------------------------------*/ -/* sparse hessian */ -/* int sparse_hess(tag, n, repeat, x, &nnz, &row_ind, &col_ind, &values, */ -/* options[2]); */ - -ADOLC_API int sparse_hess(short, int, int, const double *, int *, - unsigned int **, unsigned int **, double **, int *); - -ADOLC_API void set_HP(short tag, /* tape identification */ - int indep, /* number of independent variables */ - unsigned int **HP); - -ADOLC_API void get_HP(short tag, /* tape identification */ - int indep, /* number of independent variables */ - unsigned int ***HP); - -/*--------------------------------------------------------------------------*/ -/* JACOBIAN BLOCK PATTERN */ - -/* Max. number of unsigned ints to store the seed / jacobian matrix strips. - Reduce this value to x if your system happens to run out of memory. - x < 10 makes no sense. x = 50 or 100 is better - x stays for ( x * sizeof(size_t) * 8 ) - (block) variables at once */ +enum class SparseMethod { + IndexDomains, + BitPattern, +}; -#define PQ_STRIPMINE_MAX 30 +enum class ControlFlowMode { Safe, Tight, OldSafe, OldTight }; + +enum class BitPatternPropagationDirection { + Auto, + Forward, + Reverse, +}; + +enum class CompressionMode { + Column, + Row, +}; + +enum class RecoveryMethod { + Indirect, + Direct, +}; + +namespace detail { +// a word represents an instance of bitword_t every letter of the word gives a +// bit. This bit represents a depence of the input i (= location of the bit) to +// the output that corresponds to the word. +static constexpr size_t BITS_PER_WORD = 8 * sizeof(bitword_t); +static constexpr bitword_t MOST_SIGNIFICANT_BIT = static_cast(1) + << (BITS_PER_WORD - 1); + +/** + * @brief Internal data structure managing buffers and parameters for bit-vector + * propagation. + * + * This struct encapsulates the state and storage used during the computation + * of Jacobian sparsity patterns using bit-vector propagation (BVP). + * + * It manages seed matrices, temporary buffers, and strip-mining batches, + * depending on the chosen bit pattern propagation direction. + * + * @tparam BPPD + * Direction of bit pattern propagation — Forward or Reverse. + */ +template struct BvpData { + int rc_{0}; ///< Return code from ADOL-C drivers. + int indep_{-1}; ///< Number of independent (input) variables. + int depen_{-1}; ///< Number of dependent (output) variables. + size_t wordsPerBatch_{0}; ///< Number of bit words per batch. + size_t bitsPerStrip_{0}; ///< Number of bits represented by one batch. + size_t numStripmineBatches_{ + 0}; ///< Number of batches needed to cover all variables. + std::vector + indepWordHasNonzero_; ///< Flags marking independent variables that appear + ///< in nonzero entries. + std::vector valuepoint_; ///< Evaluation point values (used in tight + ///< control-flow mode). + std::vector seed_; ///< Seed matrix partitions. + std::vector + jacobianBitPattern_; ///< Jacobian bit pattern partitions. + + BvpData(const BvpData &) = delete; + BvpData &operator=(const BvpData &) = delete; + BvpData(BvpData &&) = delete; + BvpData &operator=(BvpData &&) = delete; + + /** + * @brief Destructor — frees dynamically allocated memory for bit patterns. + */ + ~BvpData() { + for (auto &s : seed_) + delete[] s; + for (auto &j : jacobianBitPattern_) + delete[] j; + } + + /** + * @brief Constructs a bit-vector propagation data container. + * + * Initializes memory and parameters for strip-mining based on the + * propagation direction and number of variables. + * + * @param depen Number of dependent variables. + * @param indep Number of independent variables. + */ + BvpData(int depen, int indep) : depen_(depen), indep_(indep) { + size_t numWordsFullSeed = 0; + if constexpr (BPPD == BitPatternPropagationDirection::Forward) + numWordsFullSeed = + indep_ / BITS_PER_WORD + ((indep_ % BITS_PER_WORD) != 0); + else if (BPPD == BitPatternPropagationDirection::Reverse) + numWordsFullSeed = depen / BITS_PER_WORD + ((depen % BITS_PER_WORD) != 0); + + wordsPerBatch_ = (numWordsFullSeed <= PQ_STRIPMINE_MAX) ? numWordsFullSeed + : PQ_STRIPMINE_MAX; + bitsPerStrip_ = wordsPerBatch_ * BITS_PER_WORD; + numStripmineBatches_ = + (numWordsFullSeed <= PQ_STRIPMINE_MAX) + ? 1 + : numWordsFullSeed / PQ_STRIPMINE_MAX + + ((numWordsFullSeed % PQ_STRIPMINE_MAX) != 0); + + if constexpr (BPPD == BitPatternPropagationDirection::Forward) { + seed_.resize(indep_); + for (auto &s : seed_) + s = new bitword_t[wordsPerBatch_]; + + jacobianBitPattern_.resize(depen_); + for (auto &j : jacobianBitPattern_) + j = new bitword_t[wordsPerBatch_]; + + indepWordHasNonzero_.resize(bitsPerStrip_); + + } else if (BPPD == BitPatternPropagationDirection::Reverse) { + seed_.resize(wordsPerBatch_); + for (auto &s : seed_) + s = new bitword_t[depen_]; + + jacobianBitPattern_.resize(wordsPerBatch_); + for (auto &j : jacobianBitPattern_) + j = new bitword_t[indep_]; + + indepWordHasNonzero_.resize(indep_); + } + } +}; + +/** + * @brief Verifies input data consistency for bit-vector propagation. + * + * @tparam CFM Control flow mode (Safe, Tight, etc.). + * @param basepoint Pointer to the independent variable values (may be null in + * Safe mode). + */ +template void checkBVPInput(const double *basepoint) {}; + +/** + * @brief Resets all entries in the compressed row storage pointer array to + * nullptr. + * + * @param compressedRowStorage Span over the row storage pointer array to be + * reset. + */ +inline void resetInput(std::span &compressedRowStorage) { + for (auto &row : compressedRowStorage) { + row = nullptr; + } +} + +/** + * @brief Sets all seed matrix entries to zero for reuse in the next strip-mined + * batch. + * + * @tparam BPPD Bit pattern propagation direction. + * @param data Bit-vector propagation state data. + * @param dim2 Inner dimension to clear (depends on direction). + */ +template +void resetOldSeed(BvpData &data, size_t dim2) { + for (auto &s : data.seed_) + for (int j = 0; j < dim2; ++j) + s[j] = 0; +} + +/** + * @brief Builds a strip-mined partition of the seed matrix for a given batch. + * + * @tparam BPPD Bit pattern propagation direction. + * @param indep Number of independent variables. + * @param stripIdx Current strip-mined batch index. + * @param data Bit-vector propagation state data. + */ +template +void prepareSeedMatrix(int indep, size_t stripIdx, BvpData &data) { + if constexpr (BPPD == BitPatternPropagationDirection::Forward) { + resetOldSeed(data, data.wordsPerBatch_); + createSeed(data, stripIdx, data.indep_); + } else if (BPPD == BitPatternPropagationDirection::Reverse) { + resetOldSeed(data, data.depen_); + createSeed(data, stripIdx, data.depen_); + } +} + +/** + * @brief Sets a single bit in the seed matrix corresponding to a specific + * variable index. + * + * @tparam BPPD Bit pattern propagation direction. + * @param data Bit-vector propagation state data. + * @param idx Global variable index. + * @param currentStripWordIdx Starting index of the current strip. + */ +template +void setSeed(BvpData &data, size_t idx, size_t currentStripWordIdx) { + if constexpr (BPPD == BitPatternPropagationDirection::Forward) + data.seed_[idx][(idx - currentStripWordIdx) / BITS_PER_WORD] = + MOST_SIGNIFICANT_BIT >> ((idx - currentStripWordIdx) % BITS_PER_WORD); + + else if (BPPD == BitPatternPropagationDirection::Reverse) + data.seed_[(idx - currentStripWordIdx) / BITS_PER_WORD][idx] = + MOST_SIGNIFICANT_BIT >> ((idx - currentStripWordIdx) % BITS_PER_WORD); +} + +/** + * @brief Creates a seed matrix segment by invoking setSeed() for each variable + * in the strip. + * + * @tparam BPPD Bit pattern propagation direction. + * @param data Bit-vector propagation data container. + * @param stripIdx Index of the strip-mined batch. + * @param end End index (exclusive) of the variable range for this batch. + */ +template +void createSeed(BvpData &data, size_t stripIdx, size_t end) { + size_t currentStripWordIdx = stripIdx * data.bitsPerStrip_; + size_t nextStripWordIdx = (stripIdx + 1) * data.bitsPerStrip_; + if (nextStripWordIdx > end) + nextStripWordIdx = end; + + for (size_t idx = currentStripWordIdx; idx < nextStripWordIdx; idx++) + setSeed(data, idx, currentStripWordIdx); +} + +/** + * @brief Extracts the sparsity pattern for one dependent variable (row) in + * forward bit-vector propagation mode. + * + * This routine constructs or extends a row in compressed row storage format + * based on the current bit pattern flags in @p indepWordHasNonzero_. + * + * Each bit set in the propagated Jacobian pattern indicates a dependency of + * the current dependent variable on an independent variable. + * + * @param wordIdx Index of the dependent variable (row in the Jacobian). + * @param stripIdx Index of the current strip-mined batch. + * @param currentStripWordIdx Starting global index of the strip. + * @param data Bit-vector propagation data container (Forward mode). + * @param[in,out] compressedRowStorage + * Span of pointers to CRS rows (updated in-place). + * + * @note + * - The function merges new nonzeros with previously existing entries + * when processing subsequent strips. + * - Memory for each row is reallocated as needed and must be freed by the + * caller. + */ +void extract(size_t wordIdx, size_t stripIdx, size_t currentStripWordIdx, + BvpData &data, + std::span &compressedRowStorage); + +/** + * @brief Extracts the sparsity pattern for one dependent variable (row) + * in reverse bit-vector propagation mode. + * + * This function builds the corresponding CRS row for + * each output, using the flags stored in @p indepWordHasNonzero_. + * + * Each bit set in the propagated Jacobian pattern indicates a dependency of + * the current dependent variable on an independent variable. + + * @param wordIdx Index of the dependent variable (row in the Jacobian). + * @param stripIdx Index of the current strip-mined batch. + * @param data Bit-vector propagation data container (Reverse mode). + * @param[in,out] compressedRowStorage + * Span of pointers to CRS rows (updated in-place). + * + * @note + * - Memory for each row is reallocated every time the strip is processed. + * - The caller is responsible for freeing the allocated CRS memory. + */ +void extract(size_t wordIdx, size_t stripIdx, + BvpData &data, + std::span &compressedRowStorage); + +/** + * @brief Converts Jacobian bit patterns to compressed row storage (CRS) format + * for one strip. + * + * @tparam BPPD Bit pattern propagation direction. + * @param stripIdx Index of the strip-mined batch. + * @param compressedRowStorage Output CRS array (span over row pointers). + * @param data Bit-vector propagation data container. + */ +template +void extractCompressedRowStorage(size_t stripIdx, + std::span &compressedRowStorage, + BvpData &data) { + size_t currentStripWordIdx = stripIdx * data.bitsPerStrip_; + size_t nextStripWordIdx = (stripIdx + 1) * data.bitsPerStrip_; + if constexpr (BPPD == BitPatternPropagationDirection::Forward) + for (size_t wordIdx = 0; wordIdx < data.depen_; wordIdx++) { + size_t idx = 0; + bitword_t currentBit = MOST_SIGNIFICANT_BIT; + for (size_t i = 0; i < data.bitsPerStrip_; ++i) { + if (currentBit == 0) { + currentBit = MOST_SIGNIFICANT_BIT; + idx++; + } + if (currentBit & data.jacobianBitPattern_[wordIdx][idx]) + data.indepWordHasNonzero_[i] = 1; + + currentBit = currentBit >> 1; + } + extract(wordIdx, stripIdx, currentStripWordIdx, data, + compressedRowStorage); + } + else if (BPPD == BitPatternPropagationDirection::Reverse) { + if (nextStripWordIdx > data.depen_) + nextStripWordIdx = data.depen_; + int idx = 0; + bitword_t currentBit = MOST_SIGNIFICANT_BIT; + for (size_t wordIdx = currentStripWordIdx; wordIdx < nextStripWordIdx; + wordIdx++) { + if (currentBit == 0) { + currentBit = MOST_SIGNIFICANT_BIT; + idx++; + } + for (int i = 0; i < data.indep_; i++) { + if (currentBit & data.jacobianBitPattern_[idx][i]) { + data.indepWordHasNonzero_[i] = 1; + } + } + currentBit = currentBit >> 1; + extract(wordIdx, stripIdx, data, compressedRowStorage); + } + } +} + +/** + * @brief Prepares reverse call using zos_forward. + * + * @tparam BPPD Bit pattern propagation direction. + * @tparam CFM Control flow mode. + * @param tapeId ADOL-C tape identifier. + * @param basepoint Pointer to basepoint values (used in tight control-flow + * mode). + * @param data Bit-vector propagation data container. + */ +template +void prepareADCalls(short tapeId, const double *basepoint, + BvpData &data) { + if constexpr (BPPD == BitPatternPropagationDirection::Reverse && + CFM == ControlFlowMode::Tight) + data.rc_ = zos_forward(tapeId, data.depen_, data.indep_, 1, basepoint, + data.valuepoint_.data()); +} + +/** + * @brief Executes the ADOL-C internal bit-vector propagation calls (forward or + * reverse). + * + * @tparam BPPD Bit pattern propagation direction. + * @tparam CFM Control flow mode. + * @param tapeId ADOL-C tape identifier. + * @param basepoint Pointer to independent variable values (required for forward + * + tight). + * @param data Bit-vector propagation state container. + */ +template +void ADCalls(short tapeId, const double *basepoint, BvpData &data) { + if constexpr (BPPD == BitPatternPropagationDirection::Forward) { + if constexpr (CFM == ControlFlowMode::Tight) { + data.valuepoint_.resize(data.depen_); + data.rc_ = int_forward_tight(tapeId, data.depen_, data.indep_, + data.wordsPerBatch_, basepoint, + data.seed_.data(), data.valuepoint_.data(), + data.jacobianBitPattern_.data()); + } else if (CFM == ControlFlowMode::Safe) + data.rc_ = int_forward_safe(tapeId, data.depen_, data.indep_, + data.wordsPerBatch_, data.seed_.data(), + data.jacobianBitPattern_.data()); + + } else if (BPPD == BitPatternPropagationDirection::Reverse) { + if constexpr (CFM == ControlFlowMode::Tight) + data.rc_ = int_reverse_tight(tapeId, data.depen_, data.indep_, + data.wordsPerBatch_, data.seed_.data(), + data.jacobianBitPattern_.data()); + else if (CFM == ControlFlowMode::Safe) + data.rc_ = int_reverse_safe(tapeId, data.depen_, data.indep_, + data.wordsPerBatch_, data.seed_.data(), + data.jacobianBitPattern_.data()); + } +} + +/** + * @brief Main internal implementation for bit-vector propagation. + * + * Iteratively builds and propagates seed matrices, performs ADOL-C bit-vector + * propagation calls, and extracts the sparsity pattern into CRS format. + * + * @tparam CFM Control flow mode. + * @tparam BPPD Bit pattern propagation direction. + * @param tapeId ADOL-C tape identifier. + * @param depen Number of dependent variables. + * @param indep Number of independent variables. + * @param basepoint Pointer to basepoint values. + * @param compressedRowStorage Output CRS sparsity pattern. + * @return Return code from ADOL-C driver. + */ +template +int bitVectorPropagation(short tapeId, int depen, int indep, + const double *basepoint, + std::span &compressedRowStorage) { + using namespace detail; + checkBVPInput(basepoint); + BvpData data{depen, indep}; + prepareADCalls(tapeId, basepoint, data); + for (size_t stripIdx = 0; stripIdx < data.numStripmineBatches_; stripIdx++) { + prepareSeedMatrix(indep, stripIdx, data); + ADCalls(tapeId, basepoint, data); + extractCompressedRowStorage(stripIdx, compressedRowStorage, data); + } + return data.rc_; +} + +} // namespace detail + +/** + * @brief Computes the sparsity pattern of a Jacobian using bit-vector + * propagation. + * + * This routine determines the sparsity pattern of the Jacobian of a taped + * function by propagating bit patterns through the recorded computational + * graph. It supports both forward and reverse propagation strategies, with + * optional control-flow sensitivity. + * + * @param tag + * Tape identification number for the recorded function. + * @param depen + * Number of dependent (output) variables. + * @param indep + * Number of independent (input) variables. + * @param basepoint + * Pointer to an array of length `indep` specifying the values of the + * independent variables at which the sparsity pattern is evaluated. + * Required if `options.cfmode_ == ControlFlowMode::Tight`, otherwise + * may be `nullptr`. + * @param compressedRowStorage + * Output: Compressed Row Storage (compressedRowStorage) representation + * of the Jacobian sparsity pattern. + * - `compressedRowStorage[i][0]` → number of nonzero entries in row + * `i` + * - `compressedRowStorage[i][1..]` → column indices (0-based) of + * nonzero entries in row `i`. Memory for `compressedRowStorage[i]` is + * (re)allocated within this routine. + * @param options + * DriverOptions controlling propagation and control-flow handling: + * - `options.bpdir_` (BitPatternPropagationDirection): + * - `Automatic` → heuristic selection (forward if `depen >= + * indep/2`, else reverse) + * - `Forward` → forward mode propagation + * - `Reverse` → reverse mode propagation + * - `options.cfmode_` (ControlFlowMode): + * - `Safe` → control flow branches may be ignored + * - `Tight` → control flow branches are tested at `basepoint` + * + * @return + * Return code from the underlying forward/reverse driver: + * - `0` → success, no warnings + * - `1` → warnings occurred + * - `2` → error occurred during propagation + * - `3` → default initialization value (no propagation done) + * + * @note + * - Uses "strip-mining" to split the bit pattern matrix into manageable + * blocks if the number of bits per `size_t` is insufficient to store the full + * pattern. + * - In forward mode, seeds are built per independent variable block; + * in reverse mode, seeds are built per dependent variable block. + * - Memory allocated for `compressedRowStorage[i]` must be freed by the + * caller when no longer needed. + * + * @see jac_pat + * @see hess_pat + */ +template +ADOLC_API int bit_vector_propagation(short tapeId, int depen, int indep, + const double *basepoint, + std::span &compressedRowStorage) { + using namespace detail; + if constexpr (BPPD == BitPatternPropagationDirection::Auto) { + // we allow here run-time selection + if (depen >= indep / 2) + return bitVectorPropagation( + tapeId, depen, indep, basepoint, compressedRowStorage); + + else + return bitVectorPropagation( + tapeId, depen, indep, basepoint, compressedRowStorage); + } else + return bitVectorPropagation(tapeId, depen, indep, basepoint, + compressedRowStorage); +} + +/** + * @brief Compute Jacobian sparsity pattern using the selected sparse method. + * + * Dispatches to either index-domain propagation or bit-vector propagation + * depending on the template parameter @p SM and control-flow mode @p CFM. + * + * @tparam SM SparseMethod used to compute pattern: + * - SparseMethod::IndexDomains : index-domain propagation + * - SparseMethod::BitPattern : bit-vector propagation + * @tparam CFM ControlFlowMode controlling control-flow sensitivity: + * - ControlFlowMode::Safe : conservative (may ignore branches) + * - ControlFlowMode::Tight : evaluate control-flow at @p basepoint + * @tparam BPPD BitPatternPropagationDirection used when SM==BitPattern. + * If BitPatternPropagationDirection::Auto a runtime heuristic + * selects forward/reverse direction. + * + * @param tag Tape identifier (as used by ADOL-C tracing). + * @param depen Number of dependent (output) variables (rows). + * @param indep Number of independent (input) variables (cols). + * @param basepoint Pointer to an array of length @p indep. Required + * when @p CFM == ControlFlowMode::Tight; may be + * nullptr in Safe mode. + * @param[out] compressedRowStorage + * Span of length @p depen where each element is a + * pointer to a Compressed Row Storage (CRS) row: + * - entry [i][0] holds the count of nonzeros in + * row i + * - entry [i][1..] holds the 0-based column + * indices that represent the influencing + * independs + * + * @return Non-negative: number of nonzeros found; Negative: error code from + * the underlying propagation routine. + * + * @note Ownership: allocated row buffers in compressedRowStorage are + * transferred to the caller / tape infrastructure and must be freed + * appropriately. + */ +template < + SparseMethod SM, ControlFlowMode CFM, + BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto> +ADOLC_API int jac_pat(short tag, int depen, int indep, const double *basepoint, + std::span &compressedRowStorage) { + detail::resetInput(compressedRowStorage); + if constexpr (SM == SparseMethod::IndexDomains && + CFM == ControlFlowMode::Tight) + return indopro_forward_tight(tag, depen, indep, basepoint, + compressedRowStorage.data()); + else if (SM == SparseMethod::IndexDomains && CFM == ControlFlowMode::Safe) + return indopro_forward_safe(tag, depen, indep, basepoint, + compressedRowStorage.data()); + else if (SM == SparseMethod::BitPattern) + return bit_vector_propagation(tag, depen, indep, basepoint, + compressedRowStorage); +}; + +/** + * @brief Generate a seed matrix (coloring) for compressed Jacobian recovery. + * + * Uses ColPack's Bipartite graph coloring to produce the seed matrix suitable + * for recovering the full Jacobian from compressed directional evaluations. + * This function uses ColPack's unmanaged API and returns pointers via @p Seed + * and the compressed dimension via @p p. + * + * @tparam CM CompressionMode selecting orientation of compression: + * - CompressionMode::Row : seeds correspond to row-compression + * - CompressionMode::Column : seeds correspond to column-compression + * + * @param m Number of dependent variables (rows of JP). + * @param n Number of independent variables (columns of JP). + * @param JP Compressed Row Storage, storing the dependencies on + * independents. + * @param[out] Seed Output pointer to a 2D array representing the seed matrix. + * Memory is provided by ColPack (unmanaged) and returned to + * caller. + * @param[out] p On return: compressed dimension (number of seed columns + * when Row compression, or seed rows when Column + * compression). + * + * @note The caller is responsible for managing the returned @p Seed memory + * according to ColPack's unmanaged API semantics. + */ +template +ADOLC_API void generate_seed_jac(int m, int n, const std::span JP, + double ***Seed, int *p) { + if constexpr (CM == CompressionMode::Row) + generateSeedJac(m, n, JP, Seed, p, "ROW_PARTIAL_DISTANCE_TWO"); + else if (CM == CompressionMode::Column) + generateSeedJac(m, n, JP, Seed, p, "COLUMN_PARTIAL_DISTANCE_TWO"); +} + +namespace detail { + +/** + * @brief Allocates and computes pattern of jacobian and computes a seed to + * compute a sparse Jacobian. + * + * This routine computes the Jacobian sparsity pattern (via `jac_pat`) and then + * uses ColPack to generate a seed matrix (graph coloring) suitable for + * compressed Jacobian recovery. It stores results inside the tape's + * `sJInfos()` structure so subsequent `compute_sparse_jac`/repeat calls can + * reuse the seed and recovery objects. + * + * @tparam SM + * Sparsity-detection method used by `jac_pat`: + * - `SparseMethod::IndexDomains` : index domain propagation + * - `SparseMethod::BitPattern` : bit-vector propagation + * @tparam CM + * Compression orientation for ColPack seed generation: + * - `CompressionMode::Row` : row-compression (recover by rows) + * - `CompressionMode::Column` : column-compression (recover by columns) + * @tparam CFM + * Control-flow mode used when computing the sparsity pattern. + * @tparam BPPD + * Bit pattern propagation direction when `SM == BitPattern`. When + * `BitPatternPropagationDirection::Auto` a runtime heuristic is used. + * + * @param tag ADOL-C tape identifier. + * @param depen Number of dependent (output) variables (rows of Jacobian). + * @param indep Number of independent (input) variables (columns of + * Jacobian). + * @param basepoint Pointer to an array of length `indep` with the basepoint + * values. May be `nullptr` if control-flow mode does not + * require a basepoint. + * @param[out] nnz Pointer to an integer where the computed number of + * nonzeros in the Jacobian pattern will be stored. + * + * @return + * - >= 0 : Number of seed columns (when CM==Row) or seed rows (when + * CM==Column) produced by ColPack (i.e., the compressed dimension). + * - < 0 : Error code forwarded from `jac_pat` or other internal failures. + * + * @note + * - The function allocates and stores the Compressed Row Storage (CRS) + * pattern inside `tape.sJInfos().JP_`. Ownership of those row buffers is + * transferred to the tape (they will be freed/managed by the tape + * lifecycle). + * - After success, `tape.sJInfos().g_` and `tape.sJInfos().jr1d_` are created + * and kept in the tape for later use by `compute_sparse_jac`. + * - `*nnz` is set to the number of nonzero entries found in the Jacobian. + */ +template < + SparseMethod SM, CompressionMode CM, ControlFlowMode CFM, + BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto> +int buildJacPatternAndSeed(short tag, int depen, int indep, + const double *basepoint, int *nnz) { + + ValueTape &tape = findTape(tag); + tape.sJInfos().setJP(std::vector(depen)); + std::span JPSpan = tape.sJInfos().getJP(); + int ret_val = jac_pat(tag, depen, indep, basepoint, JPSpan); + + if (ret_val < 0) { + printf(" ADOL-C error in sparse_jac() \n"); + return ret_val; + } + + tape.sJInfos().depen_ = depen; + tape.sJInfos().nnzIn_ = 0; + for (int i = 0; i < depen; i++) { + for (int j = 1; j <= tape.sJInfos().JP_[i][0]; j++) + tape.sJInfos().nnzIn_++; + } + + *nnz = tape.sJInfos().nnzIn_; + tape.sJInfos().initColoring(depen, indep); + + if constexpr (CM == CompressionMode::Row) { + tape.sJInfos().generateSeedJac("ROW_PARTIAL_DISTANCE_TWO"); + tape.sJInfos().seedClms_ = indep; + ret_val = tape.sJInfos().seedRows_; + } else if (CM == CompressionMode::Column) { + tape.sJInfos().generateSeedJac("COLUMN_PARTIAL_DISTANCE_TWO"); + tape.sJInfos().seedRows_ = depen; + ret_val = tape.sJInfos().seedClms_; + } + return ret_val; +} + +/** + * @brief Compute sparse Jacobian using precomputed seed matrix and recovery + * information. + * + * This routine performs the actual compressed Jacobian computation using the + * seed matrix and graph/recovery structures stored in the tape's `sJInfos()`. + * It supports both row- and column-oriented compression strategies and both + * user-allocated and internally allocated (unmanaged) output buffers. + * + * @tparam CM + * Compression orientation (Row or Column) chosen when the seed was + * generated. + * + * @param tag ADOL-C tape identifier. + * @param depen Number of dependent variables (rows). + * @param indep Number of independent variables (columns). + * @param basepoint Pointer to basepoint values (length `indep`) used for + * evaluation in tight control-flow or other drivers. + * @param nnz Expected number of nonzeros (must match + * `tape.sJInfos().nnzIn_`). On exit: unchanged. + * @param[out] rind Pointer to pointer which will receive the row-index array + * of the coordinate-format Jacobian (allocated by this + * function if `*rind == nullptr`, otherwise user-provided). + * @param[out] cind Pointer to pointer which will receive the column-index + * array of the coordinate-format Jacobian (allocated by + * this function if `*cind == nullptr`, otherwise + * user-provided). + * @param[out] values Pointer to pointer which will receive the nonzero values + * (allocated by this function if `*values == nullptr`, + * otherwise user-provided). + * + * @return + * - `0` on success (or driver-specific non-negative codes), + * - negative error code on failure (propagated from internal ADOL-C calls). + * + * @pre + * - `buildJacPatternAndSeed` must have been called previously for this tape to + * populate `tape.sJInfos()` with a valid seed (`Seed_`), graph (`g_`) and + * recovery object (`jr1d_`). + * - `*nnz` must equal the number of nonzeros discovered earlier and stored + * in `tape.sJInfos().nnzIn_`. If not, the function returns error `-3`. + * + * @post + * - On success, `rind`, `cind`, and `values` point to coordinate-format + * arrays describing the sparse Jacobian. If the user supplied non-null + * pointers, user memory is used; otherwise memory is allocated by the + * recovery routine and must be freed by the caller. + * + * @note + * - The function allocates intermediate storage `B_` and `y_` inside the + * tape; these are freed with the tape lifecycle or overwritten on subsequent + * calls. + * - For `CompressionMode::Row`, a forward evaluation `zos_forward` is used to + * compute `y_` followed by `fov_reverse` to obtain `B_`. For + * `CompressionMode::Column`, `fov_forward` is used. + * - The ColPack `RecoverD2*` functions are invoked to transform the compressed + * representation to coordinate format; they accept both user-managed and + * unmanaged memory usage patterns (usermem vs unmanaged). + */ +template +int computeSparseJac(short tag, int depen, int indep, const double *basepoint, + int *nnz, unsigned int **rind, unsigned int **cind, + double **values) { + ValueTape &tape = findTape(tag); + int ret_val = 0; + tape.sJInfos().B_ = + myalloc2(tape.sJInfos().seedRows_, tape.sJInfos().seedClms_); + tape.sJInfos().y_ = myalloc1(depen); + + if (tape.sJInfos().nnzIn_ != *nnz) { + printf(" ADOL-C error in sparse_jac():" + " Number of nonzeros not consistent," + " repeat call with repeat = 0 \n"); + return -3; + } + + if constexpr (CM == CompressionMode::Row) { + ret_val = zos_forward(tag, depen, indep, 1, basepoint, tape.sJInfos().y_); + if (ret_val < 0) + return ret_val; + MINDEC(ret_val, fov_reverse(tag, depen, indep, tape.sJInfos().seedRows_, + tape.sJInfos().Seed_, tape.sJInfos().B_)); + } else if (CM == CompressionMode::Column) + ret_val = + fov_forward(tag, depen, indep, tape.sJInfos().seedClms_, basepoint, + tape.sJInfos().Seed_, tape.sJInfos().y_, tape.sJInfos().B_); + + if (values != nullptr && *values != nullptr && rind != nullptr && + *rind != nullptr && cind != nullptr && *cind != nullptr) { + // everything is preallocated, we assume correctly + // call usermem versions + if (CM == CompressionMode::Row) + tape.sJInfos().recoverRowFormat(rind, cind, values); + else if (CM == CompressionMode::Column) + tape.sJInfos().recoverColFormat(rind, cind, values); + } else { + // at least one of rind cind values is not allocated, deallocate others + // and call unmanaged versions + if (values != nullptr && *values != nullptr) + free(*values); + if (rind != nullptr && *rind != nullptr) + free(*rind); + if (cind != nullptr && *cind != nullptr) + free(*cind); + if (CM == CompressionMode::Row) { + tape.sJInfos().recoverRowFormat(rind, cind, values); + } else if (CM == CompressionMode::Column) { + tape.sJInfos().recoverColFormat(rind, cind, values); + } + } + return ret_val; +} + +} // namespace detail + +/** + * @brief High-level API to compute a sparse Jacobian (single call or repeat). + * + * This function coordinates sparsity detection, seed generation (via ColPack), + * compressed AD evaluations, and recovery. When @p repeat == 0 the function + * computes and caches the sparsity/seed information (stored on the tape). + * Subsequent calls (repeat != 0) reuse cached seeds to compute numerical + * Jacobian values more efficiently. + * + * @tparam SM Default SparseMethod used for pattern detection (IndexDomains). + * @tparam CM Default CompressionMode used for recovery (Column). + * @tparam CFM Default ControlFlowMode for propagation (Safe). + * @tparam BPPD Default bit-propagation direction (Auto). + * + * @param tag Tape identifier. + * @param depen Number of dependent variables (rows). + * @param indep Number of independent variables (columns). + * @param repeat If 0: compute sparsity and prepare seed (no numeric + * recovery). If >0: perform numeric recovery using cached + * seed. + * @param basepoint Array of independent values (required if tight + * control-flow). + * @param[in,out] nnz + * Input/Output: when repeat==0 set by the routine to the + * number of nonzeros discovered. When repeat!=0 must equal + * the earlier computed number of nonzeros (consistency + * check). + * @param[out] rind Pointer-to-pointer receiving row indices (coordinate + * format). + * @param[out] cind Pointer-to-pointer receiving column indices (coordinate + * format). + * @param[out] values Pointer-to-pointer receiving numerical nonzero values. + * + * @return 0 or positive status on success, negative error code on failure. + * + * @note Memory ownership semantics for rind/cind/values mirror the underlying + * recovery routines: if user provides non-null pointers they will be used + * (usermem variants), otherwise unmanaged (allocator) variants are used + * and the caller is responsible for freeing the allocated memory. + */ +template < + SparseMethod SM = SparseMethod::IndexDomains, + CompressionMode CM = CompressionMode::Column, + ControlFlowMode CFM = ControlFlowMode::Safe, + BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto> +ADOLC_API int sparse_jac(short tag, int depen, int indep, int repeat, + const double *basepoint, int *nnz, unsigned int **rind, + unsigned int **cind, double **values) { + using namespace detail; + int ret_val = 0; + if (repeat == 0) + ret_val = buildJacPatternAndSeed(tag, depen, indep, + basepoint, nnz); + if (ret_val < 0) { + printf(" ADOL-C error in sparse_jac() \n"); + return ret_val; + } + return computeSparseJac(tag, depen, indep, basepoint, nnz, rind, cind, + values); +} + +/** + * @brief Computes the Jacobian sparsity pattern for an abs-normal form. + * + * This driver handles sparsity detection for functions represented in + * *abs-normal form* — i.e., functions involving piecewise smooth structures + * that depend on switching variables. + * + * Internally, it performs an index-domain propagation through the recorded + * abs-normal tape. + * + * @param tag ADOL-C tape identifier. + * @param depen Number of dependent (output) variables. + * @param indep Number of independent (input) variables. + * @param numsw Number of switching variables in the abs-normal representation. + * @param basepoint Pointer to basepoint values (independent variable array). + * @param[out] compressedRowStorage + * Span over an array of row pointers, each representing one Jacobian + * row in Compressed Row Storage (CRS) format. + * + * @return Return code from `indopro_forward_absnormal()`: + * - `0` → success, + * - nonzero → warning or error code from ADOL-C. + */ +ADOLC_API int absnormal_jac_pat(short tag, int depen, int indep, int numsw, + const double *basepoint, + std::span &compressedRowStorage); + +/** + * @brief Compute Hessian sparsity pattern (dispatch by control-flow mode). + * + * This function computes the sparsity pattern of the Hessian by dispatching to + * the appropriate internal non-linear independent-index propagation driver + * based on the compile-time ControlFlowMode template parameter. + * + * @tparam CFM ControlFlowMode controlling how control-flow is handled: + * - ControlFlowMode::OldTight : legacy tight control-flow handling + * - ControlFlowMode::OldSafe : legacy safe control-flow handling + * - ControlFlowMode::Tight : tight control-flow (evaluate + * branches at basepoint) + * - ControlFlowMode::Safe : safe / conservative control-flow + * + * @param tag ADOL-C tape identifier. + * @param indep Number of independent variables (Hessian + * dimension). + * @param basepoint Pointer to an array of length `indep` with the + * evaluation point used when tight control-flow is + * required. + * @param[out] compressedRowStorage + * Span of length `indep` where each element is a + * pointer to a Compressed Row Storage (CRS) row + * for the Hessian pattern. For row i: + * - entry [i][0] is the count of indices stored, + * - entry [i][1..] are the independents the output + * i depend on. + * + * @return + * - >= 0 : Driver-specific non-negative code. + * - < 0 : Error code forwarded from the underlying driver. + * + * @note + * - The routine resets `compressedRowStorage` pointers before use. + * - Memory allocated for each row will be owned by the caller/tape. + */ +template +ADOLC_API int hess_pat(short tag, int indep, const double *basepoint, + std::span &compressedRowStorage) { + detail::resetInput(compressedRowStorage); + if constexpr (CFM == ControlFlowMode::OldTight) + return nonl_ind_old_forward_tight(tag, 1, indep, basepoint, + compressedRowStorage.data()); + else if (CFM == ControlFlowMode::OldSafe) + return nonl_ind_old_forward_safe(tag, 1, indep, basepoint, + compressedRowStorage.data()); + else if (CFM == ControlFlowMode::Tight) + return nonl_ind_forward_tight(tag, 1, indep, basepoint, + compressedRowStorage.data()); + else if (CFM == ControlFlowMode::Safe) + return nonl_ind_forward_safe(tag, 1, indep, basepoint, + compressedRowStorage.data()); +} + +/** + * @brief Generate a seed matrix for compressed Hessian recovery. + * + * Produces a seed matrix for Hessian recovery using ColPack and + * returns the unmanaged pointer to the seed via @p Seed and the compressed + * dimension via @p p. The function uses ColPack's unmanaged API, so the + * returned @p Seed memory is owned by ColPack and will be freed when the + * ColPack graph object is destroyed (or as per ColPack's unmanaged semantics). + * + * @tparam RCM RecoveryMethod selecting the recovery strategy: + * - RecoveryMethod::Indirect : indirect recovery (acyclic strategy) + * - RecoveryMethod::Direct : direct recovery (star strategy) + * + * @param n Number of variables (dimension of Hessian). + * @param HP Span over CRS row pointers representing Hessian sparsity (HP). + * @param[out] Seed Pointer to the returned 2D seed matrix (ColPack-managed). + * @param[out] p Compressed dimension (number of seed columns/rows depending + * on layout). + * + * @note The caller must respect ColPack unmanaged memory semantics for the + * returned Seed. + */ +template +ADOLC_API void generate_seed_hess(int n, std::span HP, double ***Seed, + int *p) { + if constexpr (RCM == RecoveryMethod::Indirect) + generateSeedHess(n, HP, Seed, p, "ACYCLIC_FOR_INDIRECT_RECOVERY"); + else if (RCM == RecoveryMethod::Direct) + generateSeedHess(n, HP, Seed, p, "STAR"); +} + +namespace detail { + +/** + * @brief Build Hessian sparsity pattern and generate ColPack seed for recovery. + * + * This routine computes the sparsity pattern of the Hessian and constructs a + * seed matrix for Hessian recovery using ColPack. It stores + * pattern and seed/recovery related data inside the tape's `sHInfos()` + * structure for reuse by subsequent numerical recovery calls. + * + * @tparam CFM ControlFlowMode used when computing the Hessian pattern: + * - ControlFlowMode::Safe : conservative control-flow handling + * - ControlFlowMode::Tight : evaluate control-flow at @p basepoint + * @tparam RCM RecoveryMethod selecting Hessian recovery strategy: + * - RecoveryMethod::Indirect : indirect (via two-phase recover) + * - RecoveryMethod::Direct : direct recovery (star pattern) + * + * @param tag ADOL-C tape identifier. + * @param indep Number of independent variables (dimension of Hessian). + * @param basepoint Pointer to an array of length `indep` with the point at + * which control-flow (if any) should be tested. May be + * nullptr if not required by the chosen CFM (Safe vs Tight). + * @param[out] nnz Pointer to an integer where the number of nonzero Hessian + * entries will be stored. + * + * @return + * - >= 0 : Success value forwarded from `hess_pat` + * - < 0 : Error code forwarded from `hess_pat`. + * + * @note + * - The computed sparsity pattern (HP) is stored in `tape.sHInfos().HP_`. + * - The generated seed matrix is stored (via ColPack) and its contents are + * copied into `tape.sHInfos().Xppp_`. The ColPack-owned seed memory will be + * freed when the ColPack graph object (`g_`) is destroyed. + * - Additional workspace arrays (Hcomp_, Xppp_, Yppp_, Zppp_, Upp_) are + * allocated on the tape for use by the numeric recovery routine. + */ +template +int buildHessPatternAndSeed(short tag, int indep, const double *basepoint, + int *nnz) { + ValueTape &tape = findTape(tag); + int ret_val = -1; + // Generate sparsity pattern, determine nnz, allocate memory + tape.sHInfos().setHP(indep, std::vector(indep)); + std::span HPSpan = tape.sHInfos().HP_; + // generate sparsity pattern + ret_val = hess_pat(tag, indep, basepoint, HPSpan); + + if (ret_val < 0) { + printf(" ADOL-C error in sparse_hess() \n"); + return ret_val; + } + + tape.sHInfos().indep_ = indep; + tape.sHInfos().nnzIn_ = 0; + + for (int i = 0; i < indep; i++) { + for (int j = 1; j <= tape.sHInfos().HP_[i][0]; j++) + if (tape.sHInfos().HP_[i][j] >= i) + tape.sHInfos().nnzIn_++; + } + + *nnz = tape.sHInfos().nnzIn_; + + // compute seed matrix => ColPack library + + double **Seed = nullptr; + tape.sHInfos().initColoring(indep); + if constexpr (RCM == RecoveryMethod::Indirect) + tape.sHInfos().generateSeedHess(&Seed, "ACYCLIC_FOR_INDIRECT_RECOVERY"); + else if (RCM == RecoveryMethod::Direct) + tape.sHInfos().generateSeedHess(&Seed, "STAR"); + + tape.sHInfos().Hcomp_ = myalloc2(indep, tape.sHInfos().p_); + tape.sHInfos().Xppp_ = myalloc3(indep, tape.sHInfos().p_, 1); + + for (int i = 0; i < indep; i++) + for (int l = 0; l < tape.sHInfos().p_; l++) + tape.sHInfos().Xppp_[i][l][0] = Seed[i][l]; + + // Seed will be freed by ColPack when g is freed + Seed = nullptr; + + tape.sHInfos().Yppp_ = myalloc3(1, tape.sHInfos().p_, 1); + tape.sHInfos().Zppp_ = myalloc3(tape.sHInfos().p_, indep, 2); + tape.sHInfos().Upp_ = myalloc2(1, 2); + tape.sHInfos().Upp_[0][0] = 1; + tape.sHInfos().Upp_[0][1] = 0; + return ret_val; +} + +/** + * @brief Compute sparse Hessian using precomputed seed/recovery info. + * + * Executes compressed forward/reverse evaluations followed by ColPack recovery + * to produce the Hessian in coordinate (rind, cind, values) form. Supports + * both indirect and direct recovery methods (template parameter RCM). + * + * @tparam RCM RecoveryMethod controlling the chosen ColPack recovery routine. + * + * @param tag ADOL-C tape identifier. + * @param indep Number of independent variables (dimension of Hessian). + * @param basepoint Pointer to an array of length `indep` with evaluation point. + * @param nnz Expected number of nonzeros (must equal the value + * computed and stored by `buildHessPatternAndSeed`). If not, + * the function returns error `-3`. + * @param[out] rind Pointer-to-pointer that will receive the row indices of the + * coordinate-format Hessian (allocated or user-provided). + * @param[out] cind Pointer-to-pointer that will receive the column indices. + * @param[out] values Pointer-to-pointer that will receive the numerical values. + * + * @return + * - >= 0 : Success (driver-specific non-negative codes). + * - < 0 : Error code from internal ADOL-C calls. + * + * @pre + * - `buildHessPatternAndSeed<...>` must have been called previously to + * populate `tape.sHInfos()` (seed, graph `g_`, hr_, workspace arrays). + * + * @note + * - The routine uses `hov_wk_forward` + `hos_ov_reverse` to compute compressed + * Hessian blocks and then invokes ColPack's recover routines. Both usermem + * (if rind/cind/values are non-null) and unmanaged variants are supported. + * - On success, caller is responsible for freeing memory returned by the + * unmanaged recovery functions (if used). + */ +template +int computeSparseHess(short tag, int indep, const double *basepoint, int *nnz, + unsigned int **rind, unsigned int **cind, + double **values) { + ValueTape &tape = findTape(tag); + if (tape.sHInfos().Upp_ == nullptr) { + printf(" ADOL-C error in sparse_hess():" + " First call with repeat = 0 \n"); + return -3; + } + + if (tape.sHInfos().nnzIn_ != *nnz) { + printf(" ADOL-C error in sparse_hess():" + " Number of nonzeros not consistent," + " new call with repeat = 0 \n"); + return -3; + } + + // this is the most efficient variant. However, there was somewhere a + // bug in hos_ov_reverse + double y = 0.0; + int ret_val = + hov_wk_forward(tag, 1, indep, 1, 2, tape.sHInfos().p_, basepoint, + tape.sHInfos().Xppp_, &y, tape.sHInfos().Yppp_); + MINDEC(ret_val, hos_ov_reverse(tag, 1, indep, 1, tape.sHInfos().p_, + tape.sHInfos().Upp_, tape.sHInfos().Zppp_)); + + for (int i = 0; i < tape.sHInfos().p_; ++i) + for (int l = 0; l < indep; ++l) + tape.sHInfos().Hcomp_[l][i] = tape.sHInfos().Zppp_[i][l][1]; + + if (*values != nullptr && *rind != nullptr && *cind != nullptr) { + // everything is preallocated, we assume correctly + // call usermem versions + if constexpr (RCM == RecoveryMethod::Indirect) + tape.sHInfos().indirectRecoverUserMem(rind, cind, values); + else if (RCM == RecoveryMethod::Direct) + tape.sHInfos().directRecoverUserMem(rind, cind, values); + } else { + // at least one of rind cind values is not allocated, deallocate others + // and call unmanaged versions + if (*values != nullptr) + free(*values); + if (*rind != nullptr) + free(*rind); + if (*cind != nullptr) + free(*cind); + if constexpr (RCM == RecoveryMethod::Indirect) + tape.sHInfos().indirectRecover(rind, cind, values); + else if (RCM == RecoveryMethod::Direct) + tape.sHInfos().directRecover(rind, cind, values); + } + return ret_val; +} + +} // namespace detail -ADOLC_API int bit_vector_propagation(short, int, int, const double *, - unsigned int **, int *); +/** + * @brief High-level API to compute a sparse Hessian. + * + * Coordinates Hessian sparsity pattern generation (when `repeat == 0`), seed + * generation and caching, and the numeric recovery step. For `repeat == 0` + * the function computes and caches pattern/seed info; for subsequent calls + * it performs the numeric recovery using cached state. + * + * @tparam CFM ControlFlowMode used when computing the Hessian pattern. + * @tparam RCM RecoveryMethod used by ColPack to perform recovery. + * + * @param tag ADOL-C tape identifier. + * @param indep Number of independent variables (Hessian dimension). + * @param repeat If 0: compute pattern and seed (cache them). If >0: perform + * numeric recovery using cached seed. + * @param basepoint Basepoint for tight control-flow (may be nullptr if not + * required). + * @param[in,out] nnz + * On entry: expected number of nonzeros when repeat != 0. + * On exit: set to the number of nonzeros when repeat == 0. + * @param[out] rind Pointer-to-pointer to receive row indices (coordinate + * format). + * @param[out] cind Pointer-to-pointer to receive column indices. + * @param[out] values Pointer-to-pointer to receive numerical values. + * + * @return + * - >= 0 : Success (driver-specific non-negative codes). + * - < 0 : Error code (forwarded from pattern/compute routines). + * + * @note Memory ownership semantics for rind/cind/values match the underlying + * recovery calls: if user provided non-null pointers, user memory is used + * (usermem variants), otherwise the unmanaged variant allocates memory and the + * caller must free it. + */ +template +ADOLC_API int +sparse_hess(short tag, int indep, int repeat, const double *basepoint, int *nnz, + unsigned int **rind, unsigned int **cind, double **values) { + using namespace detail; + int ret_val = 0; + if (repeat == 0) + ret_val = buildHessPatternAndSeed(tag, indep, basepoint, nnz); + if (ret_val < 0) { + printf(" ADOL-C error in sparse_hess() \n"); + return ret_val; + } + return computeSparseHess(tag, indep, basepoint, nnz, rind, cind, values); +} /****************************************************************************/ -END_C_DECLS -#endif +} // namespace ADOLC::Sparse +#endif // ADOLC_SPARSE_DRIVERS_H diff --git a/ADOL-C/include/adolc/tape_interface.h b/ADOL-C/include/adolc/tape_interface.h index 80e114c4..2270b885 100644 --- a/ADOL-C/include/adolc/tape_interface.h +++ b/ADOL-C/include/adolc/tape_interface.h @@ -288,19 +288,20 @@ ADOLC_API inline size_t get_num_param(short tapeId) { return findTape(tapeId).get_num_param(); } -#ifdef SPARSE - +#ifdef ADOLC_SPARSE +namespace ADOLC::Sparse { /** * @brief Sets the Sparse Jacbian Information of the tape * * @param tapeId ID of the tape to store the jacobian information - * @param sJinfos the information to store + * @param sJInfos the information to store * * @throws ADOLCError::ErrorType::NO_TAPE_ID if the specified tape does not * exist. */ -ADOLC_API void setTapeInfoJacSparse(short tapeId, SparseJacInfos sJinfos) { - findTape(tapeId).setTapeInfoJacSparse(sJinfos); +ADOLC_API inline void setTapeInfoJacSparse(short tapeId, + SparseJacInfos &&sJInfos) { + findTape(tapeId).setTapeInfoJacSparse(std::move(sJInfos)); } /** @@ -312,9 +313,77 @@ ADOLC_API void setTapeInfoJacSparse(short tapeId, SparseJacInfos sJinfos) { * @throws ADOLCError::ErrorType::NO_TAPE_ID if the specified tape does not * exist. */ -ADOLC_API void setTapeInfoHessSparse(short tapeId, SparseHessInfos sHInfos) { - findTape(tapeId).setTageInfoHessSparse(sHInfos); +ADOLC_API inline void setTapeInfoHessSparse(short tapeId, + SparseHessInfos &&sHInfos) { + findTape(tapeId).setTapeInfoHessSparse(std::move(sHInfos)); +} + +/** + * @brief Sets the Hessian sparsity pattern for a given tape. + * + * This function stores a user-provided Hessian sparsity pattern + * (`HPIn`) for the taped function identified by `tapeId`. + * It allows advanced users to bypass the internal sparsity pattern + * computation (`hess_pat()`) when the pattern is already known. + * + * @param tapeId + * Identifier of the taped function. + * @param indep + * Number of independent variables for the taped function. + * @param HPIn + * Pointer to the compressed row storage (CRS) representation + * of the upper-triangular part of the Hessian sparsity pattern: + * - `HPIn[i][0]` → number of nonzero entries in row `i` + * - `HPIn[i][1..]` → column indices (1-based) of the nonzeros. + * + * @note + * - Only the upper-triangular part of the Hessian needs to be provided + * in `HPIn`. + * - This function transfers ownership of the CRS structure into the + * tape's internal storage. + * - Use this function only if you already know the exact sparsity + * pattern and want to skip its recomputation. + * + * @see getHP + * @see hess_pat + */ +ADOLC_API inline void setHP(short tapeId, int indep, + std::vector &&HPIn) { + SparseHessInfos sHinfos; + sHinfos.setHP(indep, std::move(HPIn)); + findTape(tapeId).setTapeInfoHessSparse(std::move(sHinfos)); } -#endif + +/** + * @brief Retrieves the Hessian sparsity pattern from a given tape. + * + * This function returns the Hessian sparsity pattern currently stored + * in the tape identified by `tapeId`. The sparsity pattern must have + * been either: + * - Computed internally by `hess_pat()`, or + * - Set explicitly via `setHP()`. + * + * @param tapeId + * Identifier of the taped function. + * @param HPOut + * Output: Pointer to the CRS representation of the + * upper-triangular Hessian sparsity pattern: + * - `(*HPOut)[i][0]` → number of nonzero entries in row `i` + * - `(*HPOut)[i][1..]` → column indices (1-based) of the nonzeros. + * + * @note + * - The returned CRS structure is owned by the tape; do not free it. + * - Only the upper-triangular part of the Hessian is stored and + * returned. + * + * @see setHP + * @see hess_pat + */ +ADOLC_API inline std::vector getHP(short tapeId) { + return findTape(tapeId).sHInfos().getHP(); +} + +} // namespace ADOLC::Sparse +#endif // ADOLC_SPARSE #endif // ADOLC_TAPE_INTERFACE_H \ No newline at end of file diff --git a/ADOL-C/include/adolc/valuetape/CMakeLists.txt b/ADOL-C/include/adolc/valuetape/CMakeLists.txt index c0e01930..50f08f29 100644 --- a/ADOL-C/include/adolc/valuetape/CMakeLists.txt +++ b/ADOL-C/include/adolc/valuetape/CMakeLists.txt @@ -2,6 +2,7 @@ install(FILES globaltapevarscl.h tapeinfos.h valuetape.h + sparseinfos.h tapeinfos.h persistanttapeinfos.h DESTINATION "include/adolc/valuetape") diff --git a/ADOL-C/include/adolc/valuetape/sparseinfos.h b/ADOL-C/include/adolc/valuetape/sparseinfos.h new file mode 100644 index 00000000..174aba24 --- /dev/null +++ b/ADOL-C/include/adolc/valuetape/sparseinfos.h @@ -0,0 +1,101 @@ +#ifndef ADOLC_SPARSE_INFOS_H +#define ADOLC_SPARSE_INFOS_H +#include +#include +#include +#include +#include + +namespace ADOLC::Sparse { +void generateSeedJac(int dimOut, int dimIn, const std::span JP, + double ***Seed, int *p, + const std::string &coloringVariant); +// stores everything we need to know to compute the sparse jacobian with +// fov_reverse + +struct ADOLC_API SparseJacInfos { + struct Impl; + std::unique_ptr pimpl_; + double *y_{nullptr}; + + // Seed is memory managed by ColPack and will be deleted + double **Seed_{nullptr}; + double **B_{nullptr}; + + // type is dictated by ColPack + std::vector JP_; + int depen_{0}; + int nnzIn_{0}; + int seedClms_{0}; + int seedRows_{0}; + + ~SparseJacInfos(); + SparseJacInfos(); + SparseJacInfos(const SparseJacInfos &) = delete; + SparseJacInfos(SparseJacInfos &&other) noexcept; + SparseJacInfos &operator=(const SparseJacInfos &) = delete; + SparseJacInfos &operator=(SparseJacInfos &&other) noexcept; + + void setJP(std::vector &&JPIn) { + JP_ = std::move(JPIn); + depen_ = JP_.size(); + } + std::vector &getJP() { return JP_; } + void initColoring(int dimOut, int dimIn); + void generateSeedJac(const std::string &coloringVariant); + void recoverRowFormatUserMem(unsigned int **rind, unsigned int **cind, + double **values); + void recoverColFormatUserMem(unsigned int **rind, unsigned int **cind, + double **values); + void recoverRowFormat(unsigned int **rind, unsigned int **cind, + double **values); + void recoverColFormat(unsigned int **rind, unsigned int **cind, + double **values); +}; + +void generateSeedHess(int dimIn, const std::span HP, double ***Seed, + int *p, const std::string &coloringVariant); +// stores everything we have to know to compute the sparse hessian via +// reverse-over-forward +struct ADOLC_API SparseHessInfos { + struct Impl; + std::unique_ptr pimpl_; + + double **Hcomp_{nullptr}; + double ***Xppp_{nullptr}; + double ***Yppp_{nullptr}; + double ***Zppp_{nullptr}; + double **Upp_{nullptr}; + + // type is dictated by ColPack + std::vector HP_; + + int nnzIn_{0}; + int indep_{0}; + int p_{0}; + +public: + ~SparseHessInfos(); + SparseHessInfos(); + SparseHessInfos(const SparseHessInfos &) = delete; + SparseHessInfos(SparseHessInfos &&other) noexcept; + SparseHessInfos &operator=(const SparseHessInfos &) = delete; + SparseHessInfos &operator=(SparseHessInfos &&other) noexcept; + + std::vector &getHP() { return HP_; } + void setHP(int indep, std::vector &&HPIn) { + indep_ = indep; + HP_ = std::move(HPIn); + } + void initColoring(int dimIn); + void generateSeedHess(double ***Seed, const std::string &coloringVariant); + void directRecoverUserMem(unsigned int **rind, unsigned int **cind, + double **values); + void indirectRecoverUserMem(unsigned int **rind, unsigned int **cind, + double **values); + void directRecover(unsigned int **rind, unsigned int **cind, double **values); + void indirectRecover(unsigned int **rind, unsigned int **cind, + double **values); +}; +} // namespace ADOLC::Sparse +#endif // ADOLC_SPARSE_INFOS_H \ No newline at end of file diff --git a/ADOL-C/include/adolc/valuetape/tapeinfos.h b/ADOL-C/include/adolc/valuetape/tapeinfos.h index 129d523b..3f45583d 100644 --- a/ADOL-C/include/adolc/valuetape/tapeinfos.h +++ b/ADOL-C/include/adolc/valuetape/tapeinfos.h @@ -334,162 +334,4 @@ struct TapeInfos { size_t get_val_space(const char *op_fileName, const char *val_fileName); }; -#ifdef SPARSE -// stores everything we need to know to compute the sparse jacobian with -// fov_reverse - -struct ADOLC_API SparseJacInfos { - ~SparseJacInfos(); - SparseJacInfos() = default; - - SparseJacInfos(const SparseJacInfos &) = delete; - SparseJacInfos &operator=(const SparseJacInfos &) = delete; - - SparseJacInfos(SparseJacInfos &&other) noexcept - : g(other.g), jr1d(other.jr1d), y(other.y), Seed(other.Seed), B(other.B), - JP(other.JP), depen(other.depen), nnz_in(other.nnz_in), - seed_clms(other.seed_clms), seed_rows(other.seed_rows) { - - // Null out source object's pointers to prevent double deletion - other.g = nullptr; - other.jr1d = nullptr; - other.y = nullptr; - other.Seed = nullptr; - other.B = nullptr; - other.JP = nullptr; - } - - SparseJacInfos &operator=(SparseJacInfos &&other) noexcept { - if (this != &other) { - // Free existing resources - -#ifdef HAVE_LIBCOLPACK - delete dynamic_cast(g); - delete dynamic_cast(jr1d); -#endif - myfree1(y); - myfree2(B); - for (int i = 0; i < depen; i++) { - free(JP[i]); - } - free(JP); - - // Move resources - g = other.g; - jr1d = other.jr1d; - y = other.y; - Seed = other.Seed; - B = other.B; - JP = other.JP; - depen = other.depen; - nnz_in = other.nnz_in; - seed_clms = other.seed_clms; - seed_rows = other.seed_rows; - - // Null out moved-from object’s pointers - other.g = nullptr; - other.jr1d = nullptr; - other.y = nullptr; - other.Seed = nullptr; - other.B = nullptr; - other.JP = nullptr; - } - return *this; - } - - void *g{nullptr}; - void *jr1d{nullptr}; - double *y{nullptr}; - - // Seed is memory managed by ColPack and will be deleted - double **Seed{nullptr}; - double **B{nullptr}; - unsigned int **JP{nullptr}; - int depen = 0, nnz_in = 0, seed_clms = 0, seed_rows = 0; -}; - -// stores everything we have to know to compute the sparse hessian via -// reverse-over-forward -struct ADOLC_API SparseHessInfos { - ~SparseHessInfos(); - SparseHessInfos() = default; - SparseHessInfos(const SparseHessInfos &) = delete; - SparseHessInfos &operator=(const SparseHessInfos &) = delete; - - SparseHessInfos(SparseHessInfos &&other) noexcept - : g(other.g), hr(other.hr), Hcomp(other.Hcomp), Xppp(other.Xppp), - Yppp(other.Yppp), Zppp(other.Zppp), Upp(other.Upp), HP(other.HP), - nnz_in(other.nnz_in), indep(other.indep), p(other.p) { - - // Null out moved-from object's pointers - other.g = nullptr; - other.hr = nullptr; - other.Hcomp = nullptr; - other.Xppp = nullptr; - other.Yppp = nullptr; - other.Zppp = nullptr; - other.Upp = nullptr; - other.HP = nullptr; - } - - SparseHessInfos &operator=(SparseHessInfos &&other) noexcept { - if (this != &other) { - // Free existing resources - -#ifdef HAVE_LIBCOLPACK - delete dynamic_cast g; - delete dynamic_cast hr; -#endif - myfree2(Hcomp); - myfree3(Xppp); - myfree3(Yppp); - myfree3(Zppp); - myfree2(Upp); - if (HP) { - for (int i = 0; i < indep; i++) - free(HP[i]); - free(HP); - } - - // Move resources - g = other.g; - hr = other.hr; - Hcomp = other.Hcomp; - Xppp = other.Xppp; - Yppp = other.Yppp; - Zppp = other.Zppp; - Upp = other.Upp; - HP = other.HP; - nnz_in = other.nnz_in; - indep = other.indep; - p = other.p; - - // Null out moved-from object's pointers - other.g = nullptr; - other.hr = nullptr; - other.Hcomp = nullptr; - other.Xppp = nullptr; - other.Yppp = nullptr; - other.Zppp = nullptr; - other.Upp = nullptr; - other.HP = nullptr; - } - return *this; - } - - void *g{nullptr}; - void *hr{nullptr}; - - double **Hcomp{nullptr}; - double ***Xppp{nullptr}; - double ***Yppp{nullptr}; - double ***Zppp{nullptr}; - double **Upp{nullptr}; - - unsigned int **HP{nullptr}; - - int nnz_in = 0, indep = 0, p = 0; -}; -#endif - #endif // ADOLC_GLOBALTAPEVARSCL_H \ No newline at end of file diff --git a/ADOL-C/include/adolc/valuetape/valuetape.h b/ADOL-C/include/adolc/valuetape/valuetape.h index 44a029da..5123728d 100644 --- a/ADOL-C/include/adolc/valuetape/valuetape.h +++ b/ADOL-C/include/adolc/valuetape/valuetape.h @@ -11,8 +11,8 @@ #include #include #include +#include #include -#include #include #include @@ -48,23 +48,17 @@ class ADOLC_API ValueTape { using StackElement = double **; std::stack cp_stack_; -#ifdef SPARSE - // updates the tape infos on sparse Jac for the given ID - void setTapeInfoJacSparse(SparseJacInfos sJinfos); - // updates the tape infos n sparse Hess for the given ID - void setTapeInfoHessSparse(SparseHessInfos sHinfos); +#ifdef ADOLC_SPARSE + ADOLC::Sparse::SparseJacInfos sJInfos_; + ADOLC::Sparse::SparseHessInfos sHInfos_; + void initSparse() { sJInfos_.~SparseJacInfos(); - new (&sJInfos_) SparseJacInfos(); + new (&sJInfos_) ADOLC::Sparse::SparseJacInfos(); sHInfos_.~SparseHessInfos(); - new (&sHInfos_) SparseHessInfos(); + new (&sHInfos_) ADOLC::Sparse::SparseHessInfos(); } - // sparse Jacobian matrices - SparseJacInfos sJinfos_; - // sparse Hessian matrices - SparseHessInfos &sHinfos() { return sHinfos_ }; - SparseJacInfos &sJonfos() { return sJInfos_ }; #endif // the compiler can not distinguish between the fct ptrs for ext_diff_fct @@ -96,10 +90,9 @@ class ADOLC_API ValueTape { ext2_buffer_(std::move(other.ext2_buffer_)), cp_buffer_(std::move(other.cp_buffer_)), cp_stack_(std::move(other.cp_stack_)) -#ifdef SPARSE +#ifdef ADOLC_SPARSE , - sJinfos_(std::move(other.sJInfos_)), - sHinfos_(std::move(other.sHInfos_)), + sJInfos_(std::move(other.sJInfos_)), sHInfos_(std::move(other.sHInfos_)) #endif { } @@ -113,14 +106,26 @@ class ADOLC_API ValueTape { ext2_buffer_ = std::move(other.ext2_buffer_); cp_buffer_ = std::move(other.cp_buffer_); cp_stack_ = std::move(other.cp_stack_); -#ifdef SPARSE - sJinfos_ = std::move(other.sJInfos_); - sHinfos_ = std::move(other.sHInfos_); +#ifdef ADOLC_SPARSE + sJInfos_ = std::move(other.sJInfos_); + sHInfos_ = std::move(other.sHInfos_); #endif } return *this; } +#ifdef ADOLC_SPARSE + // updates the tape infos on sparse Jac or Hess for the given ID + void setTapeInfoJacSparse(ADOLC::Sparse::SparseJacInfos &&sJInfos) { + sJInfos_ = std::move(sJInfos); + } + void setTapeInfoHessSparse(ADOLC::Sparse::SparseHessInfos &&sHInfos) { + sHInfos_ = std::move(sHInfos); + } + ADOLC::Sparse::SparseJacInfos &sJInfos() { return sJInfos_; } + ADOLC::Sparse::SparseHessInfos &sHInfos() { return sHInfos_; } +#endif + // Interface to PersistentTapeInfos void tapeBaseNames(size_t loc, const std::string &baseName) { perTapeInfos_.tapeBaseNames_[loc] = baseName; @@ -384,7 +389,7 @@ class ADOLC_API ValueTape { /* fill the constants buffer and write it to disk */ void put_vals_notWriteBlock(double *reals, size_t numReals) { return tapeInfos_.put_vals_notWriteBlock(reals, numReals); - }; + } /* write some constants to the buffer without disk access */ void put_val_block(double *lastValP1) { return tapeInfos_.put_val_block(lastValP1, val_fileName()); @@ -750,21 +755,6 @@ class ADOLC_API ValueTape { void cp_restore(CpInfos *cpInfos); void cp_release(CpInfos *cpInfos); CpInfos *get_cp_fct(int index) { return cp_buffer_.getElement(index); } - -#ifdef SPARSE - // updates the tape infos on sparse Jac for the given ID - void setTapeInfoJacSparse(SparseJacInfos sJinfos); - // updates the tape infos n sparse Hess for the given ID - void setTapeInfoHessSparse(SparseHessInfos sHinfos); - void initSparse() { - sJInfos_.~SparseJacInfos(); - new (&sJInfos_) SparseJacInfos(); - - sHInfos_.~SparseHessInfos(); - new (&sHInfos_) SparseHessInfos(); - } - sHInfos() -#endif }; #endif // ADOLC_VALUETAPE_H \ No newline at end of file diff --git a/ADOL-C/src/adolcerror.cpp b/ADOL-C/src/adolcerror.cpp index 5bb691a3..70ec922c 100644 --- a/ADOL-C/src/adolcerror.cpp +++ b/ADOL-C/src/adolcerror.cpp @@ -655,6 +655,9 @@ void fail(ErrorType error, source_location LocInfo, const FailInfo &failinfo) { throw ADOLCError("To many input directions", LocInfo); break; + case to_underlying(ErrorType::NOT_IMPLEMENTED): + throw ADOLCError("Method or Option is not implemented!", LocInfo); + default: throw ADOLCError("ADOL-C error => unknown error type!\n", LocInfo); } diff --git a/ADOL-C/src/adouble_tl_indo.cpp b/ADOL-C/src/adouble_tl_indo.cpp index e5acd83d..48aaf6db 100644 --- a/ADOL-C/src/adouble_tl_indo.cpp +++ b/ADOL-C/src/adouble_tl_indo.cpp @@ -17,13 +17,12 @@ sparse patterns. #include #include -#include #include -#include +#include using std::cout; -namespace adtl_indo { +namespace ADOLC::Sparse::adtl_indo { /******************* i/o operations ***************************************/ ostream &operator<<(ostream &out, const adouble &a) { @@ -37,7 +36,7 @@ istream &operator>>(istream &in, adouble &a) { } /**************** ADOLC_TRACELESS_SPARSE_PATTERN ****************************/ -int ADOLC_Init_sparse_pattern(adouble *a, int n, unsigned int start_cnt) { +int ADOLC_Init_sparse_pattern(adouble *a, int n, uint start_cnt) { for (unsigned int i = 0; i < n; i++) { a[i].delete_pattern(); a[i].pattern.push_back(i + start_cnt); @@ -45,27 +44,23 @@ int ADOLC_Init_sparse_pattern(adouble *a, int n, unsigned int start_cnt) { return 3; } -int ADOLC_get_sparse_pattern(const adouble *const b, int m, - unsigned int **&pat) { - pat = (unsigned int **)malloc(m * sizeof(unsigned int *)); +int ADOLC_get_sparse_pattern(const adouble *b, int m, + std::vector &pat) { + pat = std::vector(m); for (int i = 0; i < m; i++) { - // const_cast(b[i]).pattern.sort(); - // const_cast(b[i]).pattern.unique(); if (b[i].get_pattern_size() > 0) { - pat[i] = (unsigned int *)malloc(sizeof(unsigned int) * - (b[i].get_pattern_size() + 1)); + pat[i] = new uint[b[i].get_pattern_size() + 1]; pat[i][0] = b[i].get_pattern_size(); const list &tmp_set = b[i].get_pattern(); - list::const_iterator it; unsigned int l = 1; - for (it = tmp_set.begin(); it != tmp_set.end(); it++, l++) + for (auto it = tmp_set.begin(); it != tmp_set.end(); it++, l++) pat[i][l] = *it; } else { - pat[i] = (unsigned int *)malloc(sizeof(unsigned int)); + pat[i] = new uint[1]; pat[i][0] = 0; } } return 3; } -} // namespace adtl_indo +} // namespace ADOLC::Sparse::adtl_indo diff --git a/ADOL-C/src/fo_rev.cpp b/ADOL-C/src/fo_rev.cpp index eeb5caa8..5474a05d 100644 --- a/ADOL-C/src/fo_rev.cpp +++ b/ADOL-C/src/fo_rev.cpp @@ -258,24 +258,24 @@ int fov_reverse(short tnum, /* tape id */ /* First Order Vector version of the reverse mode for bit patterns, tight */ /****************************************************************************/ int int_reverse_tight( - short tnum, /* tape id */ - int depen, /* consistency chk on # of deps */ - int indep, /* consistency chk on # of indeps */ - int nrows, /* # of Jacobian rows being calculated */ - const size_t *const *lagrange, /* domain weight vector[var][row](in)*/ - size_t **results) /* matrix of coeff. vectors[var][row]*/ + short tnum, /* tape id */ + int depen, /* consistency chk on # of deps */ + int indep, /* consistency chk on # of indeps */ + int nrows, /* # of Jacobian rows being calculated */ + const bitword_t *const *lagrange, /* domain weight vector[var][row](in)*/ + bitword_t **results) /* matrix of coeff. vectors[var][row]*/ #elif defined(_NTIGHT_) /****************************************************************************/ /* First Order Vector version of the reverse mode, bit pattern, safe */ /****************************************************************************/ int int_reverse_safe( - short tnum, /* tape id */ - int depen, /* consistency chk on # of deps */ - int indep, /* consistency chk on # of indeps */ - int nrows, /* # of Jacobian rows being calculated */ - const size_t *const *lagrange, /* domain weight vector[var][row](in)*/ - size_t **results) /* matrix of coeff. vectors[var][row]*/ + short tnum, /* tape id */ + int depen, /* consistency chk on # of deps */ + int indep, /* consistency chk on # of indeps */ + int nrows, /* # of Jacobian rows being calculated */ + const bitword_t *const *lagrange, /* domain weight vector[var][row](in)*/ + bitword_t **results) /* matrix of coeff. vectors[var][row]*/ #else #error Neither _TIGHT_ nor _NTIGHT_ defined #endif diff --git a/ADOL-C/src/revolve.cpp b/ADOL-C/src/revolve.cpp index fcc7d755..95535825 100644 --- a/ADOL-C/src/revolve.cpp +++ b/ADOL-C/src/revolve.cpp @@ -184,6 +184,7 @@ *--------------------------------------------------------------------*/ #include +#include #define MAXINT 2147483647 diff --git a/ADOL-C/src/sparse/sparse_fo_rev.cpp b/ADOL-C/src/sparse/sparse_fo_rev.cpp index b04073b2..c4a3e78b 100644 --- a/ADOL-C/src/sparse/sparse_fo_rev.cpp +++ b/ADOL-C/src/sparse/sparse_fo_rev.cpp @@ -16,25 +16,27 @@ #include #include #include +#include -#if defined(__cplusplus) +namespace ADOLC::Sparse { /****************************************************************************/ /* Bit pattern propagation; general call */ /* */ -int forward(short tag, int m, int n, int p, double *x, size_t **X, double *y, - size_t **Y, char mode) +int forward(short tag, int m, int n, int p, double *x, + std::vector &X, double *y, std::vector &Y, + char mode) /* forward(tag, m, n, p, x[n], X[n][p], y[m], Y[m][p], mode) */ { int rc = -1; if (mode == 1) // tight version if (x != NULL) - rc = int_forward_tight(tag, m, n, p, x, X, y, Y); + rc = int_forward_tight(tag, m, n, p, x, X.data(), y, Y.data()); else ADOLCError::fail(ADOLCError::ErrorType::SPARSE_NO_BP, CURRENT_LOCATION); else if (mode == 0) // safe version - rc = int_forward_safe(tag, m, n, p, X, Y); + rc = int_forward_safe(tag, m, n, p, X.data(), Y.data()); else ADOLCError::fail(ADOLCError::ErrorType::SPARSE_BAD_MODE, CURRENT_LOCATION); @@ -44,20 +46,22 @@ int forward(short tag, int m, int n, int p, double *x, size_t **X, double *y, /****************************************************************************/ /* Bit pattern propagation; no basepoint */ /* */ -int forward(short tag, int m, int n, int p, size_t **X, size_t **Y, char mode) +int forward(short tag, int m, int n, int p, std::vector &X, + std::vector &Y, char mode) /* forward(tag, m, n, p, X[n][p], Y[m][p], mode) */ { if (mode != 0) // not safe ADOLCError::fail(ADOLCError::ErrorType::SPARSE_BAD_MODE, CURRENT_LOCATION); - return int_forward_safe(tag, m, n, p, X, Y); + return int_forward_safe(tag, m, n, p, X.data(), Y.data()); } /****************************************************************************/ /* */ /* Bit pattern propagation, general call */ /* */ -int reverse(short tag, int m, int n, int q, size_t **U, size_t **Z, char mode) +int reverse(short tag, int m, int n, int q, std::vector &U, + std::vector &Z, char mode) /* reverse(tag, m, n, q, U[q][m], Z[q][n]) */ { int rc = -1; @@ -65,9 +69,9 @@ int reverse(short tag, int m, int n, int q, size_t **U, size_t **Z, char mode) /* ! use better the tight version, the safe version supports no subscripts*/ if (mode == 0) // safe version - rc = int_reverse_safe(tag, m, n, q, U, Z); + rc = int_reverse_safe(tag, m, n, q, U.data(), Z.data()); else if (mode == 1) - rc = int_reverse_tight(tag, m, n, q, U, Z); + rc = int_reverse_tight(tag, m, n, q, U.data(), Z.data()); else ADOLCError::fail(ADOLCError::ErrorType::SPARSE_BAD_MODE, CURRENT_LOCATION); @@ -76,4 +80,4 @@ int reverse(short tag, int m, int n, int q, size_t **U, size_t **Z, char mode) /****************************************************************************/ -#endif +} // namespace ADOLC::Sparse \ No newline at end of file diff --git a/ADOL-C/src/sparse/sparsedrivers.cpp b/ADOL-C/src/sparse/sparsedrivers.cpp index 667640fb..e8576318 100644 --- a/ADOL-C/src/sparse/sparsedrivers.cpp +++ b/ADOL-C/src/sparse/sparsedrivers.cpp @@ -2,7 +2,7 @@ ADOL-C -- Automatic Differentiation by Overloading in C++ File: sparse/sparsedrivers.cpp Revision: $Id$ - Contents: "Easy To Use" C++ interfaces of SPARSE package + Contents: "Easy To Use" C++ interfaces of ADOLC_SPARSE package Copyright (c) Andrea Walther, Benjamin Letschert, Kshitij Kulshreshtha @@ -14,971 +14,106 @@ ----------------------------------------------------------------------------*/ #include #include +#include #include #include #include #include - -#if defined(ADOLC_INTERNAL) -#if HAVE_CONFIG_H -#include "config.h" -#endif -#endif - -#if HAVE_LIBCOLPACK -#include -#endif - +#include +#include +#include +#include +#include #include -#include +#include -#if HAVE_LIBCOLPACK -using namespace ColPack; -#endif - -/****************************************************************************/ -/******* sparse Jacobains, separate drivers ***************/ -/****************************************************************************/ +namespace ADOLC::Sparse { +namespace detail { -/*--------------------------------------------------------------------------*/ -/* sparsity pattern Jacobian */ -/*--------------------------------------------------------------------------*/ -/* */ - -int jac_pat( - short tag, /* tape identification */ - int depen, /* number of dependent variables */ - int indep, /* number of independent variables */ - const double *basepoint, /* independent variable values */ - unsigned int **crs, - /* returned compressed row block-index storage */ - int *options - /* control options - options[0] : way of sparsity pattern computation - 0 - propagation of index domains (default) - 1 - propagation of bit pattern - options[1] : test the computational graph control flow - 0 - safe mode (default) - 1 - tight mode - options[2] : way of bit pattern propagation - 0 - automatic detection (default) - 1 - forward mode - 2 - reverse mode */ -) { - int rc = -1; - int i, ctrl_options[2]; - - if (!crs) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_CRS, CURRENT_LOCATION); - - else - for (i = 0; i < depen; i++) - crs[i] = NULL; - - if ((options[0] < 0) || (options[0] > 1)) - options[0] = 0; /* default */ - if ((options[1] < 0) || (options[1] > 1)) - options[1] = 0; /* default */ - if ((options[2] < -1) || (options[2] > 2)) - options[2] = 0; /* default */ - - if (options[0] == 0) { - if (options[1] == 1) - rc = indopro_forward_tight(tag, depen, indep, basepoint, crs); - else { - rc = indopro_forward_safe(tag, depen, indep, basepoint, crs); - } - } else { - ctrl_options[0] = options[1]; - ctrl_options[1] = options[2]; - rc = - bit_vector_propagation(tag, depen, indep, basepoint, crs, ctrl_options); - } - - return (rc); -} - -int absnormal_jac_pat(short tag, /* tape identification */ - int depen, /* number of dependent variables */ - int indep, /* number of independent variables */ - int numsw, /* number of switches */ - const double *basepoint, /* independent variable values */ - unsigned int **crs - /* returned compressed row block-index storage */ -) { - - if (!crs) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_CRS, CURRENT_LOCATION); - else - for (int i = 0; i < depen + numsw; i++) - crs[i] = NULL; - return indopro_forward_absnormal(tag, depen, indep, numsw, basepoint, crs); -} -/*--------------------------------------------------------------------------*/ -/* seed matrix for Jacobian */ -/*--------------------------------------------------------------------------*/ - -void generate_seed_jac(int m, int n, unsigned int **JP, double ***Seed, int *p, - int option - /* control options - option : way of compression - 0 - column compression - (default) 1 - row compression */ -) -#if HAVE_LIBCOLPACK -{ - int dummy; - - BipartiteGraphPartialColoringInterface *g = - new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, JP, m, n); - - if (option == 1) - g->GenerateSeedJacobian_unmanaged(Seed, p, &dummy, "SMALLEST_LAST", - "ROW_PARTIAL_DISTANCE_TWO"); - else - g->GenerateSeedJacobian_unmanaged(Seed, &dummy, p, "SMALLEST_LAST", - "COLUMN_PARTIAL_DISTANCE_TWO"); - delete g; +/// Specialization of input consistency check for Tight control-flow mode. +template <> +void checkBVPInput(const double *basepoint) { + if (basepoint == nullptr) + ADOLCError::fail(ADOLCError::ErrorType::SPARSE_JAC_NO_BP, CURRENT_LOCATION); } -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); -} -#endif -/****************************************************************************/ -/******* sparse Hessians, separate drivers ***************/ -/****************************************************************************/ - -/*---------------------------------------------------------------------------*/ -/* sparsity pattern Hessian */ -/* */ - -int hess_pat(short tag, /* tape identification */ - int indep, /* number of independent variables */ - const double *basepoint, /* independent variable values */ - unsigned int **crs, - /* returned compressed row block-index storage */ - int option - /* control option - option : test the computational graph control flow - 0 - safe mode (default) - 1 - tight mode - 2 - old safe mode - 3 - old tight mode */ - -) { - int rc = -1; - int i; - - if (!crs) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_CRS, CURRENT_LOCATION); - else - for (i = 0; i < indep; i++) - crs[i] = nullptr; - - if ((option < 0) || (option > 3)) - option = 0; /* default */ - - if (option == 3) - rc = nonl_ind_old_forward_tight(tag, 1, indep, basepoint, crs); - else if (option == 2) - rc = nonl_ind_old_forward_safe(tag, 1, indep, basepoint, crs); - else if (option == 1) - rc = nonl_ind_forward_tight(tag, 1, indep, basepoint, crs); - else - rc = nonl_ind_forward_safe(tag, 1, indep, basepoint, crs); - - return (rc); -} - -/*--------------------------------------------------------------------------*/ -/* seed matrix for Hessian */ -/*--------------------------------------------------------------------------*/ - -void generate_seed_hess(int n, unsigned int **HP, double ***Seed, int *p, - int option - /* control options - option : way of compression - 0 - indirect recovery - (default) 1 - direct recovery */ -) -#if HAVE_LIBCOLPACK -{ - int seed_rows; - - GraphColoringInterface *g = new GraphColoringInterface(SRC_MEM_ADOLC, HP, n); - - if (option == 0) - g->GenerateSeedHessian_unmanaged(Seed, &seed_rows, p, "SMALLEST_LAST", - "ACYCLIC_FOR_INDIRECT_RECOVERY"); - else - g->GenerateSeedHessian_unmanaged(Seed, &seed_rows, p, "SMALLEST_LAST", - "STAR"); - delete g; -} -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); -} -#endif - -static void deepcopy_HP(unsigned int ***HPnew, unsigned int **HP, int indep) { - int i, j, s; - *HPnew = (unsigned int **)malloc(indep * sizeof(unsigned int *)); - for (i = 0; i < indep; i++) { - s = HP[i][0]; - (*HPnew)[i] = (unsigned int *)malloc((s + 1) * (sizeof(unsigned int))); - for (j = 0; j <= s; j++) - (*HPnew)[i][j] = HP[i][j]; - } -} - -/****************************************************************************/ -/******* sparse Jacobians, complete driver ***************/ -/****************************************************************************/ - -int sparse_jac(short tag, /* tape identification */ - int depen, /* number of dependent variables */ - int indep, /* number of independent variables */ - int repeat, /* indicated repeated call with same seed */ - const double *basepoint, /* independent variable values */ - int *nnz, /* number of nonzeros */ - unsigned int **rind, /* row index */ - unsigned int **cind, /* column index */ - double **values, /* non-zero values */ - int *options - /* control options - options[0] : way of sparsity pattern computation - 0 - propagation of index domains - (default) 1 - propagation of bit pattern options[1] : test the - computational graph control flow 0 - safe mode (default) 1 - - tight mode options[2] : way of bit pattern propagation 0 - - automatic detection (default) 1 - forward mode 2 - reverse - mode options[3] : way of compression 0 - column compression - (default) 1 - row compression */ -) -#if HAVE_LIBCOLPACK -{ - std::shared_ptr tape = getTape(tag); - int i; - unsigned int j; - SparseJacInfos sJinfos; - int ret_val = 0; - BipartiteGraphPartialColoringInterface *g; - JacobianRecovery1D *jr1d; - JacobianRecovery1D jr1d_loc; +void extract(size_t wordIdx, size_t stripIdx, size_t currentStripWordIdx, + BvpData &data, + std::span &compressedRowStorage) { - if (repeat == 0) { - if ((options[0] < 0) || (options[0] > 1)) - options[0] = 0; /* default */ - if ((options[1] < 0) || (options[1] > 1)) - options[1] = 0; /* default */ - if ((options[2] < -1) || (options[2] > 2)) - options[2] = 0; /* default */ - if ((options[3] < 0) || (options[3] > 1)) - options[3] = 0; /* default */ + size_t numNonzeros = std::accumulate(data.indepWordHasNonzero_.begin(), + data.indepWordHasNonzero_.end(), 0UL); - sJinfos.JP = (unsigned int **)malloc(depen * sizeof(unsigned int *)); - ret_val = jac_pat(tag, depen, indep, basepoint, sJinfos.JP, options); - - if (ret_val < 0) { - printf(" ADOL-C error in sparse_jac() \n"); - return ret_val; - } - - sJinfos.depen = depen; - sJinfos.nnz_in = depen; - sJinfos.nnz_in = 0; - for (i = 0; i < depen; i++) { - for (j = 1; j <= sJinfos.JP[i][0]; j++) - sJinfos.nnz_in++; - } - - *nnz = sJinfos.nnz_in; - - if (options[2] == -1) { - (*rind) = (unsigned int *)calloc(*nnz, sizeof(unsigned int)); - (*cind) = (unsigned int *)calloc(*nnz, sizeof(unsigned int)); - unsigned int index = 0; - for (i = 0; i < depen; i++) - for (j = 1; j <= sJinfos.JP[i][0]; j++) { - (*rind)[index] = i; - (*cind)[index++] = sJinfos.JP[i][j]; - } - } - - /* sJinfos.Seed is memory managed by ColPack and will be deleted - * along with g. We only keep it in sJinfos for the repeat != 0 case */ - - g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, sJinfos.JP, - depen, indep); - jr1d = new JacobianRecovery1D; - - if (options[3] == 1) { - g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), - &(sJinfos.seed_clms), "SMALLEST_LAST", - "ROW_PARTIAL_DISTANCE_TWO"); - sJinfos.seed_clms = indep; - ret_val = sJinfos.seed_rows; + if (numNonzeros > 0) { + uint *row = nullptr; + if (stripIdx == 0) { + row = new uint[numNonzeros + 1]; + row[0] = 0; } else { - g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), - &(sJinfos.seed_clms), "SMALLEST_LAST", - "COLUMN_PARTIAL_DISTANCE_TWO"); - sJinfos.seed_rows = depen; - ret_val = sJinfos.seed_clms; + row = new uint[compressedRowStorage[wordIdx][0] + numNonzeros + 1]; + std::copy(compressedRowStorage[wordIdx], + compressedRowStorage[wordIdx] + + compressedRowStorage[wordIdx][0] + 1, + row); } - sJinfos.B = myalloc2(sJinfos.seed_rows, sJinfos.seed_clms); - sJinfos.y = myalloc1(depen); - - sJinfos.g = (void *)g; - sJinfos.jr1d = (void *)jr1d; - setTapeInfoJacSparse(tag, sJinfos); - } else { - sJinfos.depen = tape.sJinfos.depen; - sJinfos.nnz_in = tape.sJinfos.nnz_in; - sJinfos.JP = tape.sJinfos.JP; - sJinfos.B = tape.sJinfos.B; - sJinfos.y = tape.sJinfos.y; - sJinfos.Seed = tape.sJinfos.Seed; - sJinfos.seed_rows = tape.sJinfos.seed_rows; - sJinfos.seed_clms = tape.sJinfos.seed_clms; - g = (BipartiteGraphPartialColoringInterface *)tape.sJinfos.g; - jr1d = (JacobianRecovery1D *)tape.sJinfos.jr1d; - } - - if (sJinfos.nnz_in != *nnz) { - printf(" ADOL-C error in sparse_jac():" - " Number of nonzeros not consistent," - " repeat call with repeat = 0 \n"); - return -3; - } - - if (options[2] == -1) - return ret_val; - - /* compute jacobian times matrix product */ - - if (options[3] == 1) { - ret_val = zos_forward(tag, depen, indep, 1, basepoint, sJinfos.y); - if (ret_val < 0) - return ret_val; - MINDEC(ret_val, fov_reverse(tag, depen, indep, sJinfos.seed_rows, - sJinfos.Seed, sJinfos.B)); - } else - ret_val = fov_forward(tag, depen, indep, sJinfos.seed_clms, basepoint, - sJinfos.Seed, sJinfos.y, sJinfos.B); - - /* recover compressed Jacobian => ColPack library */ - - if (*values != NULL && *rind != NULL && *cind != NULL) { - // everything is preallocated, we assume correctly - // call usermem versions - if (options[3] == 1) - jr1d->RecoverD2Row_CoordinateFormat_usermem(g, sJinfos.B, sJinfos.JP, - rind, cind, values); - else - jr1d->RecoverD2Cln_CoordinateFormat_usermem(g, sJinfos.B, sJinfos.JP, - rind, cind, values); - } else { - // at least one of rind cind values is not allocated, deallocate others - // and call unmanaged versions - if (*values != NULL) - free(*values); - if (*rind != NULL) - free(*rind); - if (*cind != NULL) - free(*cind); - if (options[3] == 1) - jr1d->RecoverD2Row_CoordinateFormat_unmanaged(g, sJinfos.B, sJinfos.JP, - rind, cind, values); - else - jr1d->RecoverD2Cln_CoordinateFormat_unmanaged(g, sJinfos.B, sJinfos.JP, - rind, cind, values); - } - - return ret_val; -} -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); - return -1; -} -#endif - -/****************************************************************************/ -/******* sparse Hessians, complete driver ***************/ -/****************************************************************************/ - -int sparse_hess(short tag, /* tape identification */ - int indep, /* number of independent variables */ - int repeat, /* indicated repeated call with same seed */ - const double *basepoint, /* independent variable values */ - int *nnz, /* number of nonzeros */ - unsigned int **rind, /* row index */ - unsigned int **cind, /* column index */ - double **values, /* non-zero values */ - int *options - /* control options - options[0] :test the computational graph control - flow 0 - safe mode (default) 1 - tight mode 2 - old safe mode - 3 - old tight mode - options[1] : way of recovery - 0 - indirect recovery - 1 - direct recovery */ -) -#if HAVE_LIBCOLPACK -{ - ValueTape &tape = findTape(tag); - int i, l; - unsigned int j; - SparseHessInfos sHinfos; - double **Seed; - int dummy; - double y; - int ret_val = -1; - GraphColoringInterface *g; - HessianRecovery *hr; - - /* Generate sparsity pattern, determine nnz, allocate memory */ - if (repeat <= 0) { - if ((options[0] < 0) || (options[0] > 3)) - options[0] = 0; /* default */ - if ((options[1] < 0) || (options[1] > 1)) - options[1] = 0; /* default */ - - if (repeat == 0) { - sHinfos.HP = (unsigned int **)malloc(indep * sizeof(unsigned int *)); - - /* generate sparsity pattern */ - ret_val = hess_pat(tag, indep, basepoint, sHinfos.HP, options[0]); - - if (ret_val < 0) { - printf(" ADOL-C error in sparse_hess() \n"); - return ret_val; + uint idx = row[0] + 1; + for (size_t i = 0; i < data.indepWordHasNonzero_.size(); i++) { + if (data.indepWordHasNonzero_[i]) { + row[idx++] = currentStripWordIdx + i; + data.indepWordHasNonzero_[i] = 0; } - } else { - if (indep != tape.sHinfos.indep) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_HESS_IND, - CURRENT_LOCATION); - deepcopy_HP(&sHinfos.HP, tape.sHinfos.HP, indep); - } - - sHinfos.indep = indep; - sHinfos.nnz_in = 0; - - for (i = 0; i < indep; i++) { - for (j = 1; j <= sHinfos.HP[i][0]; j++) - if ((int)sHinfos.HP[i][j] >= i) - sHinfos.nnz_in++; } - - *nnz = sHinfos.nnz_in; - - /* compute seed matrix => ColPack library */ - - Seed = NULL; - - g = new GraphColoringInterface(SRC_MEM_ADOLC, sHinfos.HP, indep); - hr = new HessianRecovery; - - if (options[1] == 0) - g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, "SMALLEST_LAST", - "ACYCLIC_FOR_INDIRECT_RECOVERY"); - else - g->GenerateSeedHessian(&Seed, &dummy, &sHinfos.p, "SMALLEST_LAST", - "STAR"); - - sHinfos.Hcomp = myalloc2(indep, sHinfos.p); - sHinfos.Xppp = myalloc3(indep, sHinfos.p, 1); - - for (i = 0; i < indep; i++) - for (l = 0; l < sHinfos.p; l++) - sHinfos.Xppp[i][l][0] = Seed[i][l]; - - /* Seed will be freed by ColPack when g is freed */ - Seed = NULL; - - sHinfos.Yppp = myalloc3(1, sHinfos.p, 1); - - sHinfos.Zppp = myalloc3(sHinfos.p, indep, 2); - - sHinfos.Upp = myalloc2(1, 2); - sHinfos.Upp[0][0] = 1; - sHinfos.Upp[0][1] = 0; - - sHinfos.g = (void *)g; - sHinfos.hr = (void *)hr; - - setTapeInfoHessSparse(tag, sHinfos); - - } else { - sHinfos.nnz_in = tape.sHinfos.nnz_in; - sHinfos.HP = tape.sHinfos.HP; - sHinfos.Hcomp = tape.sHinfos.Hcomp; - sHinfos.Xppp = tape.sHinfos.Xppp; - sHinfos.Yppp = tape.sHinfos.Yppp; - sHinfos.Zppp = tape.sHinfos.Zppp; - sHinfos.Upp = tape.sHinfos.Upp; - sHinfos.p = tape.sHinfos.p; - g = (GraphColoringInterface *)tape.sHinfos.g; - hr = (HessianRecovery *)tape.sHinfos.hr; - } - - if (sHinfos.Upp == NULL) { - printf(" ADOL-C error in sparse_hess():" - " First call with repeat = 0 \n"); - return -3; - } - - if (sHinfos.nnz_in != *nnz) { - printf(" ADOL-C error in sparse_hess():" - " Number of nonzeros not consistent," - " new call with repeat = 0 \n"); - return -3; - } - - if (repeat == -1) - return ret_val; - - // this is the most efficient variant. However, there was somewhere a bug - // in hos_ov_reverse - ret_val = hov_wk_forward(tag, 1, indep, 1, 2, sHinfos.p, basepoint, - sHinfos.Xppp, &y, sHinfos.Yppp); - MINDEC(ret_val, hos_ov_reverse(tag, 1, indep, 1, sHinfos.p, sHinfos.Upp, - sHinfos.Zppp)); - - for (i = 0; i < sHinfos.p; ++i) - for (l = 0; l < indep; ++l) - sHinfos.Hcomp[l][i] = sHinfos.Zppp[i][l][1]; - - if (*values != NULL && *rind != NULL && *cind != NULL) { - // everything is preallocated, we assume correctly - // call usermem versions - if (options[1] == 0) - hr->IndirectRecover_CoordinateFormat_usermem(g, sHinfos.Hcomp, sHinfos.HP, - rind, cind, values); - else - hr->DirectRecover_CoordinateFormat_usermem(g, sHinfos.Hcomp, sHinfos.HP, - rind, cind, values); - } else { - // at least one of rind cind values is not allocated, deallocate others - // and call unmanaged versions - if (*values != NULL) - free(*values); - if (*rind != NULL) - free(*rind); - if (*cind != NULL) - free(*cind); - if (options[1] == 0) - hr->IndirectRecover_CoordinateFormat_unmanaged( - g, sHinfos.Hcomp, sHinfos.HP, rind, cind, values); - else - hr->DirectRecover_CoordinateFormat_unmanaged(g, sHinfos.Hcomp, sHinfos.HP, - rind, cind, values); + row[0] = idx - 1; + delete[] compressedRowStorage[wordIdx]; + compressedRowStorage[wordIdx] = row; } - return ret_val; -} -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); - return -1; -} -#endif - -/****************************************************************************/ -/******* sparse Hessians, set and get sparsity pattern ***************/ -/****************************************************************************/ - -void set_HP(short tag, /* tape identification */ - int indep, /* number of independent variables */ - unsigned int **HP) -#ifdef SPARSE -{ - std::shared_ptr tape = getTape(tag); - SparseHessInfos sHinfos; - - sHinfos.nnz_in = 0; - deepcopy_HP(&sHinfos.HP, HP, indep); - sHinfos.Hcomp = NULL; - sHinfos.Xppp = NULL; - sHinfos.Yppp = NULL; - sHinfos.Zppp = NULL; - sHinfos.Upp = NULL; - sHinfos.p = 0; - sHinfos.g = NULL; - sHinfos.hr = NULL; - sHinfos.indep = indep; - setTapeInfoHessSparse(tag, sHinfos); } -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); -} -#endif - -void get_HP(short tag, /* tape identification */ - int indep, /* number of independent variables */ - unsigned int ***HP) -#ifdef SPARSE -{ - deepcopy_HP(HP, findTape(tag).sHinfos.HP, indep); -} -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); -} -#endif - -/*****************************************************************************/ -/* JACOBIAN BLOCK PATTERN */ - -/* ------------------------------------------------------------------------- */ -int bit_vector_propagation( - short tag, /* tape identification */ - int depen, /* number of dependent variables */ - int indep, /* number of independent variables */ - const double *basepoint, /* independent variable values */ - unsigned int **crs, - /* compressed block row storage */ - int *options /* control options */ - /* options[0] : way of bit pattern propagation - 0 - automatic detection (default) - 1 - forward mode - 2 - reverse mode - options[1] : test the computational graph control flow - 0 - safe variant (default) - 1 - tight variant */ -) { - - int rc = 3; - char forward_mode, tight_mode; - int i, ii, j, jj, k, k_old, bits_per_long, i_blocks_per_strip, - d_blocks_per_strip; - int this_strip_i_bl_idx, next_strip_i_bl_idx, this_strip_d_bl_idx, - next_strip_d_bl_idx; - int stripmined_calls, strip_idx; - int p_stripmine, q_stripmine, p_ind_bl_bp, q_dep_bl_bp, i_bl_idx, d_bl_idx; - size_t value1, v; - size_t **seed = NULL, *s, **jac_bit_pat = NULL, *jac; - unsigned char *indep_blocks_flags = NULL, *i_b_flags; - double *valuepoint = NULL; - - if (options[1] == 0) { - if (depen >= indep / 2) - options[1] = 1; /* forward */ - else - options[1] = 2; /* reverse */ - } - - if (options[1] == 1) - forward_mode = 1; - else - forward_mode = 0; - - if (options[0] == 1) - tight_mode = 1; - else - tight_mode = 0; - - if (!forward_mode) - valuepoint = myalloc1(depen); - - /* bit pattern parameters */ - - /* number of bits in an size_t variable */ - bits_per_long = 8 * sizeof(size_t); - /* olvo 20000214 nl: inserted explicit cast to size_t */ - value1 = (size_t)1 << (bits_per_long - 1); /* 10000....0 */ - - /* =================================================== forward propagation */ - if (forward_mode) { - - if ((tight_mode) && !basepoint) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_JAC_NO_BP, - CURRENT_LOCATION); - - /* indep partial derivatives for the whole Jacobian */ - - /* number of size_ts to store the whole seed / Jacobian matrice */ - p_ind_bl_bp = indep / bits_per_long + ((indep % bits_per_long) != 0); - - /* number of size_ts to store the seed / Jacobian strips */ - if (p_ind_bl_bp <= PQ_STRIPMINE_MAX) { - p_stripmine = p_ind_bl_bp; - stripmined_calls = 1; - } else { - p_stripmine = PQ_STRIPMINE_MAX; - stripmined_calls = p_ind_bl_bp / PQ_STRIPMINE_MAX + - ((p_ind_bl_bp % PQ_STRIPMINE_MAX) != 0); - } - - /* number of independent blocks per seed / Jacobian strip */ - i_blocks_per_strip = p_stripmine * bits_per_long; - - /* allocate memory --------------------------------------------------- */ - - if (!(indep_blocks_flags = - (unsigned char *)calloc(i_blocks_per_strip, sizeof(char)))) - ADOLCError::fail( - ADOLCError::ErrorType::SPARSE_JAC_MALLOC, - ADOLCError::FailInfo{.info2 = i_blocks_per_strip * sizeof(char)}, - CURRENT_LOCATION); - seed = myalloc2_ulong(indep, p_stripmine); - jac_bit_pat = myalloc2_ulong(depen, p_stripmine); +void extract(size_t wordIdx, size_t stripIdx, + BvpData &data, + std::span &compressedRowStorage) { - /* strip-mining : repeated forward calls ----------------------------- - */ + size_t numNonzeros = std::accumulate(data.indepWordHasNonzero_.begin(), + data.indepWordHasNonzero_.end(), 0UL); + uint *row = new uint[numNonzeros + 1]; - for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) { - /* build a partition of the seed matrix (indep x indep_blocks) --- */ - /* (indep x i_blocks_per_strip) as a bit pattern */ - s = seed[0]; - for (i = 0; i < indep; i++) - for (ii = 0; ii < p_stripmine; ii++) /* 2 loops if short -> int !!! */ - *s++ = 0; /* set old seed matrix to 0 */ + row[0] = numNonzeros; - this_strip_i_bl_idx = strip_idx * i_blocks_per_strip; - next_strip_i_bl_idx = (strip_idx + 1) * i_blocks_per_strip; - if (next_strip_i_bl_idx > indep) - next_strip_i_bl_idx = indep; - v = value1; /* 10000....0 */ - - for (i = 0; i < indep; i++) - if ((this_strip_i_bl_idx <= i) && (i < next_strip_i_bl_idx)) { - ii = (i - this_strip_i_bl_idx) / bits_per_long; - seed[i][ii] = v >> ((i - this_strip_i_bl_idx) % bits_per_long); - } - - /* bit pattern propagation by forward ---------------------------- */ - - if (tight_mode) { - rc = int_forward_tight(tag, depen, indep, p_stripmine, basepoint, seed, - valuepoint, jac_bit_pat); - } else { - rc = - int_forward_safe(tag, depen, indep, p_stripmine, seed, jac_bit_pat); - } - - /* extract pattern from bit patterns --------------------- */ - - for (j = 0; j < depen; j++) { - ii = -1; - v = 0; - - jac = jac_bit_pat[j]; - i_b_flags = indep_blocks_flags; - for (i_bl_idx = 0; i_bl_idx < i_blocks_per_strip; i_bl_idx++) { - if (!v) { - v = value1; /* 10000....0 */ - ii++; - } - if (v & jac[ii]) - *i_b_flags = 1; - i_b_flags++; - - v = v >> 1; - } - - if (strip_idx == 0) - k_old = 0; - else - k_old = crs[j][0]; - k = 0; - i_b_flags = indep_blocks_flags; - for (i = 0; i < i_blocks_per_strip; i++) - k += *i_b_flags++; - - if ((k > 0) || (strip_idx == 0)) { - if (!(crs[j] = (unsigned int *)realloc( - crs[j], (k_old + k + 1) * sizeof(unsigned int)))) - ADOLCError::fail( - ADOLCError::ErrorType::SPARSE_JAC_MALLOC, - ADOLCError::FailInfo{.info2 = (k_old + k + 1) * - sizeof(unsigned int)}, - CURRENT_LOCATION); - - if (strip_idx == 0) - crs[j][0] = 0; - if (k > 0) { - k = crs[j][0] + 1; - i_b_flags = indep_blocks_flags; - for (i = 0; i < i_blocks_per_strip; i++) { - if (*i_b_flags) { - crs[j][k++] = this_strip_i_bl_idx + i; - *i_b_flags = 0; - } - i_b_flags++; - } - /* current/total number of non-zero blocks of indep. vars. */ - crs[j][0] = k - 1; - } - } - } - } /* strip_idx */ - - } /* forward */ - - /* =================================================== reverse propagation */ - else { - /* depen weight vectors for the whole Jacobian */ - - /* number of size_ts to store the whole seed / Jacobian matrice */ - q_dep_bl_bp = depen / bits_per_long + ((depen % bits_per_long) != 0); - - /* number of size_ts to store the seed / Jacobian strips */ - if (q_dep_bl_bp <= PQ_STRIPMINE_MAX) { - q_stripmine = q_dep_bl_bp; - stripmined_calls = 1; - } else { - q_stripmine = PQ_STRIPMINE_MAX; - stripmined_calls = q_dep_bl_bp / PQ_STRIPMINE_MAX + - ((q_dep_bl_bp % PQ_STRIPMINE_MAX) != 0); - } - - /* number of dependent blocks per seed / Jacobian strip */ - d_blocks_per_strip = q_stripmine * bits_per_long; - - /* allocate memory --------------------------------------------------- - */ - if (!(indep_blocks_flags = - (unsigned char *)calloc(indep, sizeof(unsigned char)))) - ADOLCError::fail( - ADOLCError::ErrorType::SPARSE_JAC_MALLOC, - ADOLCError::FailInfo{.info2 = indep * sizeof(unsigned char)}, - CURRENT_LOCATION); - seed = myalloc2_ulong(q_stripmine, depen); - jac_bit_pat = myalloc2_ulong(q_stripmine, indep); - - /* olvo 20000214: call to forward required in tight mode only, - in safe mode no basepoint available! */ - if (tight_mode) { - if (!basepoint) - ADOLCError::fail(ADOLCError::ErrorType::SPARSE_JAC_NO_BP, - CURRENT_LOCATION); - - rc = zos_forward(tag, depen, indep, 1, basepoint, valuepoint); + uint idx = 1; + for (size_t i = 0; i < data.indepWordHasNonzero_.size(); ++i) { + if (data.indepWordHasNonzero_[i]) { + row[idx++] = i; + data.indepWordHasNonzero_[i] = 0; } - - /* strip-mining : repeated reverse calls ----------------------------- - */ - - for (strip_idx = 0; strip_idx < stripmined_calls; strip_idx++) { - /* build a partition of the seed matrix (depen_blocks x depen) */ - /* (d_blocks_per_strip x depen) as a bit pattern */ - s = seed[0]; - for (jj = 0; jj < q_stripmine; jj++) /* 2 loops if short -> int !!! */ - for (j = 0; j < depen; j++) - *s++ = 0; /* set old seed matrix to 0 */ - - this_strip_d_bl_idx = strip_idx * d_blocks_per_strip; - next_strip_d_bl_idx = (strip_idx + 1) * d_blocks_per_strip; - if (next_strip_d_bl_idx > depen) - next_strip_d_bl_idx = depen; - v = value1; /* 10000....0 */ - - for (j = 0; j < depen; j++) - if ((this_strip_d_bl_idx <= j) && (j < next_strip_d_bl_idx)) { - jj = (j - this_strip_d_bl_idx) / bits_per_long; - seed[jj][j] = v >> ((j - this_strip_d_bl_idx) % bits_per_long); - } - - /* bit pattern propagation by reverse ---------------------------- */ - - if (tight_mode) - rc = int_reverse_tight(tag, depen, indep, q_stripmine, seed, - jac_bit_pat); - else - rc = - int_reverse_safe(tag, depen, indep, q_stripmine, seed, jac_bit_pat); - - /* extract pattern from bit patterns --------------------- */ - - jj = -1; - v = 0; - for (d_bl_idx = this_strip_d_bl_idx; d_bl_idx < next_strip_d_bl_idx; - d_bl_idx++) { - if (!v) { - v = value1; /* 10000....0 */ - jj++; - } - jac = jac_bit_pat[jj]; - for (i = 0; i < indep; i++) { - if (v & *jac++) { - indep_blocks_flags[i] = 1; - } - } - - v = v >> 1; - - k = 0; - i_b_flags = indep_blocks_flags; - for (i = 0; i < indep; i++) - k += *i_b_flags++; - - if (!(crs[d_bl_idx] = - (unsigned int *)malloc((k + 1) * sizeof(unsigned int)))) - ADOLCError::fail( - ADOLCError::ErrorType::SPARSE_JAC_MALLOC, - ADOLCError::FailInfo{.info2 = (k + 1) * sizeof(unsigned int)}, - CURRENT_LOCATION) - - crs[d_bl_idx][0] = k; /* number of non-zero indep. blocks */ - k = 1; - i_b_flags = indep_blocks_flags; - for (i = 0; i < indep; i++) { - if (*i_b_flags) { - crs[d_bl_idx][k++] = i; - *i_b_flags = 0; - } - i_b_flags++; - } - } - - } /* strip_idx */ - - } /* reverse */ - - if (!forward_mode) { - free((char *)valuepoint); - valuepoint = NULL; } - free((char *)*seed); - free((char *)seed); - seed = NULL; - free((char *)*jac_bit_pat); - free((char *)jac_bit_pat); - jac_bit_pat = NULL; - free((char *)indep_blocks_flags); - indep_blocks_flags = NULL; - - return (rc); + delete[] compressedRowStorage[wordIdx]; + compressedRowStorage[wordIdx] = row; } -#include +} // namespace detail -// namespace adtl { +int absnormal_jac_pat(short tag, int depen, int indep, int numsw, + const double *basepoint, + std::span &compressedRowStorage) { + detail::resetInput(compressedRowStorage); + return indopro_forward_absnormal(tag, depen, indep, numsw, basepoint, + compressedRowStorage.data()); +} -#ifdef SPARSE -SparseJacInfos sJinfos = {NULL, NULL, NULL, NULL, NULL, NULL, 0, 0, 0, 0}; -#endif +#include +SparseJacInfos sJInfos; +namespace adtl_indo { -int ADOLC_get_sparse_jacobian(func_ad *const fun, +int ADOLC_get_sparse_jacobian(func_ad<::adtl::adouble> *const fun, func_ad *const fun_indo, int n, int m, int repeat, double *basepoints, int *nnz, unsigned int **rind, - unsigned int **cind, double **values) -#if HAVE_LIBCOLPACK -{ + unsigned int **cind, double **values) { int i; unsigned int j; int ret_val = -1; if (!repeat) { - freeSparseJacInfos(sJinfos.y, sJinfos.B, sJinfos.JP, sJinfos.g, - sJinfos.jr1d, sJinfos.seed_rows, sJinfos.seed_clms, - sJinfos.depen); // setNumDir(n); // setMode(ADTL_INDO); { @@ -998,96 +133,72 @@ int ADOLC_get_sparse_jacobian(func_ad *const fun, return ret_val; } - ret_val = adtl_indo::ADOLC_get_sparse_pattern(y, m, sJinfos.JP); + ret_val = adtl_indo::ADOLC_get_sparse_pattern(y, m, sJInfos.JP_); delete[] x; delete[] y; } - sJinfos.depen = m; - sJinfos.nnz_in = 0; + sJInfos.depen_ = m; + sJInfos.nnzIn_ = 0; for (i = 0; i < m; i++) { - for (j = 1; j <= sJinfos.JP[i][0]; j++) - sJinfos.nnz_in++; + for (j = 1; j <= sJInfos.JP_[i][0]; j++) + sJInfos.nnzIn_++; } - *nnz = sJinfos.nnz_in; - /* sJinfos.Seed is memory managed by ColPack and will be deleted - * along with g. We only keep it in sJinfos for the repeat != 0 case */ - BipartiteGraphPartialColoringInterface *g; - JacobianRecovery1D *jr1d; + *nnz = sJInfos.nnzIn_; + /* sJInfos.Seed is memory managed by ColPack and will be deleted + * along with g. We only keep it in sJInfos for the repeat != 0 case */ - g = new BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, sJinfos.JP, m, - n); - jr1d = new JacobianRecovery1D; + sJInfos.initColoring(m, n); + sJInfos.generateSeedJac("COLUMN_PARTIAL_DISTANCE_TWO"); + sJInfos.seedRows_ = m; + sJInfos.B_ = myalloc2(sJInfos.seedRows_, sJInfos.seedClms_); + sJInfos.y_ = myalloc1(m); - g->GenerateSeedJacobian(&(sJinfos.Seed), &(sJinfos.seed_rows), - &(sJinfos.seed_clms), "SMALLEST_LAST", - "COLUMN_PARTIAL_DISTANCE_TWO"); - sJinfos.seed_rows = m; - - sJinfos.B = myalloc2(sJinfos.seed_rows, sJinfos.seed_clms); - sJinfos.y = myalloc1(m); - - sJinfos.g = (void *)g; - sJinfos.jr1d = (void *)jr1d; - - if (sJinfos.nnz_in != *nnz) { + if (sJInfos.nnzIn_ != *nnz) { printf(" ADOL-C error in sparse_jac():" " Number of nonzeros not consistent," " repeat call with repeat = 0 \n"); return -3; } } - // ret_val = fov_forward(tag, depen, indep, sJinfos.seed_clms, basepoint, - // sJinfos.Seed, sJinfos.y, sJinfos.B); - adtl::setNumDir(sJinfos.seed_clms); + // ret_val = fov_forward(tag, depen, indep, sJInfos.seed_clms, basepoint, + // sJInfos.Seed, sJInfos.y, sJInfos.B); + ::adtl::setNumDir(sJInfos.seedClms_); // setMode(ADTL_FOV); { - adtl::adouble *x, *y; - x = new adtl::adouble[n]; - y = new adtl::adouble[m]; + ::adtl::adouble *x, *y; + x = new ::adtl::adouble[n]; + y = new ::adtl::adouble[m]; for (i = 0; i < n; i++) { x[i] = basepoints[i]; - for (int jj = 0; jj < sJinfos.seed_clms; jj++) - x[i].setADValue(jj, sJinfos.Seed[i][jj]); + for (int jj = 0; jj < sJInfos.seedClms_; jj++) + x[i].setADValue(jj, sJInfos.Seed_[i][jj]); } ret_val = (*fun)(n, x, m, y); for (i = 0; i < m; i++) - for (int jj = 0; jj < sJinfos.seed_clms; jj++) - sJinfos.B[i][jj] = y[i].getADValue(jj); + for (int jj = 0; jj < sJInfos.seedClms_; jj++) + sJInfos.B_[i][jj] = y[i].getADValue(jj); delete[] x; delete[] y; } /* recover compressed Jacobian => ColPack library */ - if (*values != NULL) + if (*values != nullptr) free(*values); - if (*rind != NULL) + if (*rind != nullptr) free(*rind); - if (*cind != NULL) + if (*cind != nullptr) free(*cind); - BipartiteGraphPartialColoringInterface *g; - JacobianRecovery1D *jr1d; - g = (BipartiteGraphPartialColoringInterface *)sJinfos.g; - jr1d = (JacobianRecovery1D *)sJinfos.jr1d; - jr1d->RecoverD2Cln_CoordinateFormat_unmanaged(g, sJinfos.B, sJinfos.JP, rind, - cind, values); - - // delete g; - // delete jr1d; + sJInfos.recoverColFormat(rind, cind, values); return ret_val; } -#else -{ - ADOLCError::fail(ADOLCError::ErrorType::NO_COLPACK, CURRENT_LOCATION); - return -1; -} -#endif - -//} /****************************************************************************/ -/* THAT'S ALL */ +/* THAT'S ALL + */ +} // namespace adtl_indo +} // namespace ADOLC::Sparse diff --git a/ADOL-C/src/uni5_for.cpp b/ADOL-C/src/uni5_for.cpp index 039b7f7a..b22a2753 100644 --- a/ADOL-C/src/uni5_for.cpp +++ b/ADOL-C/src/uni5_for.cpp @@ -558,14 +558,14 @@ int fos_forward_nk( /* First Order Vector version of the forward mode for bit patterns, tight */ /****************************************************************************/ int int_forward_tight( - short tnum, /* tape id */ - int depcheck, /* consistency chk on # of dependents */ - int indcheck, /* consistency chk on # of independents */ - int p, /* # of taylor series, bit pattern */ - const double *basepoint, /* independent variable values (in)*/ - const size_t *const *argument, /* Taylor coeff. (in)*/ - double *valuepoint, /* dependent variable values (out)*/ - size_t **taylors) /* matrix of coefficient vectors(out)*/ + short tnum, /* tape id */ + int depcheck, /* consistency chk on # of dependents */ + int indcheck, /* consistency chk on # of independents */ + int p, /* # of taylor series, bit pattern */ + const double *basepoint, /* independent variable values (in)*/ + const bitword_t *const *argument, /* Taylor coeff. (in)*/ + double *valuepoint, /* dependent variable values (out)*/ + bitword_t **taylors) /* matrix of coefficient vectors(out)*/ /* int_forward_tight( tag, m, n, p, x[n], X[n][p], y[m], Y[m][p]), @@ -585,12 +585,13 @@ int int_forward_tight( /****************************************************************************/ /* First Order Vector version of the forward mode, bit pattern, safe */ /****************************************************************************/ -int int_forward_safe(short tnum, /* tape id */ - int depcheck, /* consistency chk on # of dependents */ - int indcheck, /* consistency chk on # of independents */ - int p, /* # of taylor series, bit pattern */ - const size_t *const *argument, /* Taylor coeff. (in)*/ - size_t **taylors) /* matrix of coefficient vectors (out)*/ +int int_forward_safe( + short tnum, /* tape id */ + int depcheck, /* consistency chk on # of dependents */ + int indcheck, /* consistency chk on # of independents */ + int p, /* # of taylor series, bit pattern */ + const bitword_t *const *argument, /* Taylor coeff. (in)*/ + bitword_t **taylors) /* matrix of coefficient vectors (out)*/ /* int_forward_safe( tag, m, n, p, X[n][p], Y[m][p]), @@ -616,7 +617,7 @@ int indopro_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ /* indopro_forward_tight( tag, m, n, x[n], *crs[m]), @@ -630,7 +631,7 @@ int indopro_forward_tight( int indcheck, /* consistency chk on # of independents */ int swcheck, /* consistency chk on # of switches */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ /* indopro_forward_absnormal( tag, m, n, s, x[n], *crs[s+m]), */ @@ -643,7 +644,7 @@ int indopro_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ /* indopro_forward_safe( tag, m, n, x[n], *crs[m]), @@ -661,7 +662,7 @@ int nonl_ind_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ #endif #if defined(_NTIGHT_) @@ -673,7 +674,7 @@ int nonl_ind_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ /* indopro_forward_safe( tag, m, n, x[n], *crs[m]), @@ -690,7 +691,7 @@ int nonl_ind_old_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ #endif #if defined(_NTIGHT_) @@ -702,7 +703,7 @@ int nonl_ind_old_forward_tight( int depcheck, /* consistency chk on # of dependents */ int indcheck, /* consistency chk on # of independents */ const double *basepoint, /* independent variable values (in) */ - unsigned int **crs) /* returned row index storage (out) */ + uint **crs) /* returned row index storage (out) */ /* indopro_forward_safe( tag, m, n, x[n], *crs[m]), @@ -1174,7 +1175,6 @@ int hov_forward( ind_dom[i][0] = 0; ind_dom[i][1] = NUMNNZ; } - #endif #if defined(_NONLIND_) maxopind = tape.tapestats(TapeInfos::NUM_OPERATIONS) + @@ -1635,13 +1635,13 @@ int hov_forward( #if defined(_INDO_) #if defined(_INDOPRO_) && !defined(_NONLIND_OLD_) if (ind_dom[res][0] != 0) { - crs[indexd] = new unsigned int[ind_dom[res][0] + 1]; + crs[indexd] = new uint[ind_dom[res][0] + 1]; crs[indexd][0] = ind_dom[res][0]; for (l = 1; l <= crs[indexd][0]; l++) { crs[indexd][l] = ind_dom[res][l + 1]; } } else { - crs[indexd] = new unsigned int[1]; + crs[indexd] = new uint[1]; crs[indexd][0] = 0; } #endif @@ -3898,16 +3898,19 @@ int hov_forward( #if defined(_INDO_) #if defined(_INDOPRO_) #if defined(_ABS_NORM_) + // set index domain of switching variables that of the input independent + // which creates the switching variable. note: the switching variables are + // stored in the Compressed Row Storage (crs) AFTER the dependent + // variables. Therefore we have to use depcheck + switchnum as index if (ind_dom[arg][0] != 0) { - crs[depcheck + switchnum] = (unsigned int *)malloc( - sizeof(unsigned int) * (ind_dom[arg][0] + 1)); + crs[depcheck + switchnum] = new uint[ind_dom[arg][0] + 1]; crs[depcheck + switchnum][0] = ind_dom[arg][0]; + // switching variable switchnum gets index-domain of the arg for (l = 1; l <= crs[depcheck + switchnum][0]; l++) { crs[depcheck + switchnum][l] = ind_dom[arg][l + 1]; } } else { - crs[depcheck + switchnum] = - (unsigned int *)malloc(sizeof(unsigned int)); + crs[depcheck + switchnum] = new uint[1]; crs[depcheck + switchnum][0] = 0; } ind_dom[res][0] = 1; diff --git a/ADOL-C/src/valuetape/CMakeLists.txt b/ADOL-C/src/valuetape/CMakeLists.txt index 0b69741d..3798689a 100644 --- a/ADOL-C/src/valuetape/CMakeLists.txt +++ b/ADOL-C/src/valuetape/CMakeLists.txt @@ -1,5 +1,6 @@ target_sources(adolc PRIVATE valuetape.cpp + sparseinfos.cpp tapeinfos.cpp globaltapevarscl.cpp persistanttapeinfos.cpp diff --git a/ADOL-C/src/valuetape/sparseinfos.cpp b/ADOL-C/src/valuetape/sparseinfos.cpp new file mode 100644 index 00000000..1ca186a7 --- /dev/null +++ b/ADOL-C/src/valuetape/sparseinfos.cpp @@ -0,0 +1,245 @@ +#include +#include +#include +#include + +namespace ADOLC::Sparse { + +void generateSeedJac(int dimOut, int dimIn, const std::span JP, + double ***Seed, int *p, + const std::string &coloringVariant) { + int dummy = 0; + ColPack::BipartiteGraphPartialColoringInterface(SRC_MEM_ADOLC, JP.data(), + dimOut, dimIn) + .GenerateSeedJacobian_unmanaged(Seed, p, &dummy, "SMALLEST_LAST", + coloringVariant); +} +struct SparseJacInfos::Impl { + ColPack::BipartiteGraphPartialColoringInterface g_{SRC_WAIT}; + ColPack::JacobianRecovery1D jr1d_{}; +}; + +SparseJacInfos::~SparseJacInfos() { + if (y_) + myfree1(y_); + y_ = nullptr; + if (B_) + myfree2(B_); + B_ = nullptr; + + for (int i = 0; i < depen_; i++) { + delete[] JP_[i]; + JP_[i] = nullptr; + } +} + +SparseJacInfos::SparseJacInfos() + : pimpl_(std::make_unique()) {} + +SparseJacInfos::SparseJacInfos(SparseJacInfos &&other) noexcept + : pimpl_(std::move(other.pimpl_)), y_(other.y_), Seed_(other.Seed_), + B_(other.B_), JP_(std::move(other.JP_)), depen_(other.depen_), + nnzIn_(other.nnzIn_), seedClms_(other.seedClms_), + seedRows_(other.seedRows_) { + + // Null out source object's pointers to prevent double deletion + other.y_ = nullptr; + other.Seed_ = nullptr; + other.B_ = nullptr; +} + +SparseJacInfos &SparseJacInfos::operator=(SparseJacInfos &&other) noexcept { + if (this != &other) { + // Free existing resources + + myfree1(y_); + myfree2(B_); + for (int i = 0; i < depen_; i++) { + delete[] JP_[i]; + JP_[i] = nullptr; + } + + // Move resources + pimpl_ = std::move(other.pimpl_); + y_ = other.y_; + Seed_ = other.Seed_; + B_ = other.B_; + JP_ = other.JP_; + depen_ = other.depen_; + nnzIn_ = other.nnzIn_; + seedClms_ = other.seedClms_; + seedRows_ = other.seedRows_; + + // Null out moved-from object’s pointers + other.y_ = nullptr; + other.Seed_ = nullptr; + other.B_ = nullptr; + } + return *this; +} + +void SparseJacInfos::initColoring(int dimOut, int dimIn) { + pimpl_->g_ = ColPack::BipartiteGraphPartialColoringInterface( + SRC_MEM_ADOLC, JP_.data(), dimOut, dimIn); + pimpl_->jr1d_ = ColPack::JacobianRecovery1D(); +} + +void SparseJacInfos::generateSeedJac(const std::string &coloringVariant) { + pimpl_->g_.GenerateSeedJacobian(&Seed_, &seedRows_, &seedClms_, + "SMALLEST_LAST", coloringVariant); +} + +void SparseJacInfos::recoverRowFormatUserMem(unsigned int **rind, + unsigned int **cind, + double **values) { + pimpl_->jr1d_.RecoverD2Row_CoordinateFormat_usermem( + &pimpl_->g_, B_, JP_.data(), rind, cind, values); +} + +void SparseJacInfos::recoverColFormatUserMem(unsigned int **rind, + unsigned int **cind, + double **values) { + pimpl_->jr1d_.RecoverD2Cln_CoordinateFormat_usermem( + &pimpl_->g_, B_, JP_.data(), rind, cind, values); +} + +void SparseJacInfos::recoverRowFormat(unsigned int **rind, unsigned int **cind, + double **values) { + pimpl_->jr1d_.RecoverD2Row_CoordinateFormat_unmanaged( + &pimpl_->g_, B_, JP_.data(), rind, cind, values); +} + +void SparseJacInfos::recoverColFormat(unsigned int **rind, unsigned int **cind, + double **values) { + pimpl_->jr1d_.RecoverD2Cln_CoordinateFormat_unmanaged( + &pimpl_->g_, B_, JP_.data(), rind, cind, values); +} + +void generateSeedHess(int dimIn, const std::span HP, double ***Seed, + int *p, const std::string &coloringVariant) { + int seed_rows = 0; + ColPack::GraphColoringInterface(SRC_MEM_ADOLC, HP.data(), dimIn) + .GenerateSeedHessian_unmanaged(Seed, &seed_rows, p, "SMALLEST_LAST", + coloringVariant); +} + +struct SparseHessInfos::Impl { + ColPack::GraphColoringInterface g_{SRC_WAIT}; + ColPack::HessianRecovery hr_; +}; + +SparseHessInfos::~SparseHessInfos() { + myfree2(Hcomp_); + Hcomp_ = nullptr; + + myfree3(Xppp_); + Xppp_ = nullptr; + + myfree3(Yppp_); + Yppp_ = nullptr; + + myfree3(Zppp_); + Zppp_ = nullptr; + + myfree2(Upp_); + Upp_ = nullptr; + + for (int i = 0; i < indep_; i++) { + delete[] HP_[i]; + HP_[i] = nullptr; + } +} + +SparseHessInfos::SparseHessInfos() + : pimpl_(std::make_unique()) {}; + +SparseHessInfos::SparseHessInfos(SparseHessInfos &&other) noexcept + : pimpl_(std::move(other.pimpl_)), Hcomp_(other.Hcomp_), Xppp_(other.Xppp_), + Yppp_(other.Yppp_), Zppp_(other.Zppp_), Upp_(other.Upp_), + HP_(std::move(other.HP_)), nnzIn_(other.nnzIn_), indep_(other.indep_), + p_(other.p_) { + + // Null out moved-from object's pointers + other.Hcomp_ = nullptr; + other.Xppp_ = nullptr; + other.Yppp_ = nullptr; + other.Zppp_ = nullptr; + other.Upp_ = nullptr; +} + +SparseHessInfos &SparseHessInfos::operator=(SparseHessInfos &&other) noexcept { + if (this != &other) { + // Free existing resources + myfree2(Hcomp_); + myfree3(Xppp_); + myfree3(Yppp_); + myfree3(Zppp_); + myfree2(Upp_); + + for (int i = 0; i < indep_; i++) { + delete[] HP_[i]; + HP_[i] = nullptr; + } + + // Move resources + pimpl_ = std::move(other.pimpl_); + Hcomp_ = other.Hcomp_; + Xppp_ = other.Xppp_; + Yppp_ = other.Yppp_; + Zppp_ = other.Zppp_; + Upp_ = other.Upp_; + HP_ = other.HP_; + nnzIn_ = other.nnzIn_; + indep_ = other.indep_; + p_ = other.p_; + + // Null out moved-from object's pointers + other.Hcomp_ = nullptr; + other.Xppp_ = nullptr; + other.Yppp_ = nullptr; + other.Zppp_ = nullptr; + other.Upp_ = nullptr; + } + return *this; +} + +void SparseHessInfos::initColoring(int dimIn) { + pimpl_->g_ = + ColPack::GraphColoringInterface(SRC_MEM_ADOLC, HP_.data(), dimIn); + pimpl_->hr_ = ColPack::HessianRecovery(); +} + +void SparseHessInfos::generateSeedHess(double ***Seed, + const std::string &coloringVariant) { + int dummy = 0; + pimpl_->g_.GenerateSeedHessian(Seed, &dummy, &p_, "SMALLEST_LAST", + coloringVariant); +} + +void SparseHessInfos::indirectRecoverUserMem(unsigned int **rind, + unsigned int **cind, + double **values) { + pimpl_->hr_.IndirectRecover_CoordinateFormat_usermem( + &pimpl_->g_, Hcomp_, HP_.data(), rind, cind, values); +} + +void SparseHessInfos::directRecoverUserMem(unsigned int **rind, + unsigned int **cind, + double **values) { + pimpl_->hr_.DirectRecover_CoordinateFormat_usermem( + &pimpl_->g_, Hcomp_, HP_.data(), rind, cind, values); +} + +void SparseHessInfos::indirectRecover(unsigned int **rind, unsigned int **cind, + double **values) { + pimpl_->hr_.IndirectRecover_CoordinateFormat_unmanaged( + &pimpl_->g_, Hcomp_, HP_.data(), rind, cind, values); +} + +void SparseHessInfos::directRecover(unsigned int **rind, unsigned int **cind, + double **values) { + pimpl_->hr_.DirectRecover_CoordinateFormat_unmanaged( + &pimpl_->g_, Hcomp_, HP_.data(), rind, cind, values); +} + +} // namespace ADOLC::Sparse \ No newline at end of file diff --git a/ADOL-C/src/valuetape/tapeinfos.cpp b/ADOL-C/src/valuetape/tapeinfos.cpp index 9ea92e2e..a2cf417f 100644 --- a/ADOL-C/src/valuetape/tapeinfos.cpp +++ b/ADOL-C/src/valuetape/tapeinfos.cpp @@ -1,7 +1,6 @@ #include #include #include // for memset -#include TapeInfos::TapeInfos(short tapeId) { tapeId_ = tapeId; } @@ -268,110 +267,6 @@ void TapeInfos::freeTapeResources() { } } -#ifdef SPARSE -SparseJacInfos::~SparseJacInfos() { - myfree1(y); - y = nullptr; - myfree2(B); - B = nullptr; - - for (int i = 0; i < depen; i++) { - free(JP[i]); - } - free(JP); - JP = nullptr; - -#ifdef HAVE_LIBCOLPACK - if (g) - delete (BipartiteGraphPartialColoringInterface *)g; - if (jr1d) - delete (JacobianRecovery1D *)jr1d; -#endif -} - -SparseHessInfos::~SparseHessInfos() { - myfree2(Hcomp); - Hcomp = nullptr; - - myfree3(Xppp); - Xppp = nullptr; - - myfree3(Yppp); - Yppp = nullptr; - - myfree3(Zppp); - Zppp = nullptr; - - myfree2(Upp); - Upp = nullptr; - - if (HP) { - for (int i = 0; i < indep; i++) - free(HP[i]); - free(HP); - HP = nullptr; - } - -#ifdef HAVE_LIBCOLPACK - delete (GraphColoringInterface *)g; - g = nullptr; - delete (HessianRecovery *)hr; - hr = nullptr; -#endif -} - -void freeSparseJacInfos(double *y, double **B, unsigned int **JP, void *g, - void *jr1d, int seed_rows, int seed_clms, int depen) { - if (y) - myfree1(y); - - if (B) - myfree2(B); - - for (int i = 0; i < depen; i++) { - free(JP[i]); - } - - free(JP); - -#ifdef HAVE_LIBCOLPACK - delete (BipartiteGraphPartialColoringInterface *)g; - delete (JacobianRecovery1D *)jr1d; -#endif -} - -void freeSparseHessInfos(double **Hcomp, double ***Xppp, double ***Yppp, - double ***Zppp, double **Upp, unsigned int **HP, - void *g, void *hr, int p, int indep) { - if (Hcomp) - myfree2(Hcomp); - - if (Xppp) - myfree3(Xppp); - - if (Yppp) - myfree3(Yppp); - - if (Zppp) - myfree3(Zppp); - - if (Upp) - myfree2(Upp); - - if (HP) { - for (int i = 0; i < indep; i++) - free(HP[i]); - free(HP); - } - -#ifdef HAVE_LIBCOLPACK - delete (GraphColoringInterface *)g; - delete (HessianRecovery *)hr; -#endif -} - -#endif - /****************************************************************************/ /* Writes the block of size depth of taylor coefficients from point loc to */ /* the taylor buffer. If the buffer is filled, then it is written to the */ diff --git a/ADOL-C/src/valuetape/valuetape.cpp b/ADOL-C/src/valuetape/valuetape.cpp index 720cba23..ff55fd4e 100644 --- a/ADOL-C/src/valuetape/valuetape.cpp +++ b/ADOL-C/src/valuetape/valuetape.cpp @@ -82,7 +82,7 @@ int ValueTape::initNewTape() { // creates new tapeInfos object with old buffers // thus, we dont allocate the buffers again if they are already existent initTapeInfos_keep(); -#ifdef SPARSE +#ifdef ADOLC_SPARSE initSparse(); #endif @@ -1204,49 +1204,3 @@ std::array ValueTape::readConfigFile() { } return tapeBaseNames; } - -#ifdef SPARSE -/* updates the tape infos on sparse Jac for the given ID */ -void ValueTape::setTapeInfoJacSparse(dSparseJacInfos sJinfos) { - // free memory of tape entry that had been used previously - freeSparseJacInfos( - tapeInfos_.ptapeInfos_.sJinfos.y, tapeInfos_.ptapeInfos_.sJinfos.B, - tapeInfos_.ptapeInfos_.sJinfos.JP, tapeInfos_.ptapeInfos_.sJinfos.g, - tapeInfos_.ptapeInfos_.sJinfos.jr1d, - tapeInfos_.ptapeInfos_.sJinfos.seed_rows, - tapeInfos_.ptapeInfos_.sJinfos.seed_clms, - tapeInfos_.ptapeInfos_.sJinfos.depen); - tapeInfos_.ptapeInfos_.sJinfos.y = sJinfos.y; - tapeInfos_.ptapeInfos_.sJinfos.Seed = sJinfos.Seed; - tapeInfos_.ptapeInfos_.sJinfos.B = sJinfos.B; - tapeInfos_.ptapeInfos_.sJinfos.JP = sJinfos.JP; - tapeInfos_.ptapeInfos_.sJinfos.depen = sJinfos.depen; - tapeInfos_.ptapeInfos_.sJinfos.nnz_in = sJinfos.nnz_in; - tapeInfos_.ptapeInfos_.sJinfos.seed_clms = sJinfos.seed_clms; - tapeInfos_.ptapeInfos_.sJinfos.seed_rows = sJinfos.seed_rows; - tapeInfos_.ptapeInfos_.sJinfos.g = sJinfos.g; - tapeInfos_.ptapeInfos_.sJinfos.jr1d = sJinfos.jr1d; -} - -/* updates the tape infos on sparse Hess for the given ID */ -void ValueTape::setTapeInfoHessSparse(SparseHessInfos sHinfos) { - // free memory of tape entry that had been used previously - freeSparseHessInfos( - tapeInfos_.ptapeInfos_.sHinfos.Hcomp, tapeInfos_.ptapeInfos_.sHinfos.Xppp, - tapeInfos_.ptapeInfos_.sHinfos.Yppp, tapeInfos_.ptapeInfos_.sHinfos.Zppp, - tapeInfos_.ptapeInfos_.sHinfos.Upp, tapeInfos_.ptapeInfos_.sHinfos.HP, - tapeInfos_.ptapeInfos_.sHinfos.g, tapeInfos_.ptapeInfos_.sHinfos.hr, - tapeInfos_.ptapeInfos_.sHinfos.p, tapeInfos_.ptapeInfos_.sHinfos.indep); - tapeInfos_.ptapeInfos_.sHinfos.Hcomp = sHinfos.Hcomp; - tapeInfos_.ptapeInfos_.sHinfos.Xppp = sHinfos.Xppp; - tapeInfos_.ptapeInfos_.sHinfos.Yppp = sHinfos.Yppp; - tapeInfos_.ptapeInfos_.sHinfos.Zppp = sHinfos.Zppp; - tapeInfos_.ptapeInfos_.sHinfos.Upp = sHinfos.Upp; - tapeInfos_.ptapeInfos_.sHinfos.HP = sHinfos.HP; - tapeInfos_.ptapeInfos_.sHinfos.indep = sHinfos.indep; - tapeInfos_.ptapeInfos_.sHinfos.nnz_in = sHinfos.nnz_in; - tapeInfos_.ptapeInfos_.sHinfos.p = sHinfos.p; - tapeInfos_.ptapeInfos_.sHinfos.g = sHinfos.g; - tapeInfos_.ptapeInfos_.sHinfos.hr = sHinfos.hr; -} -#endif diff --git a/CMakeLists.txt b/CMakeLists.txt index fa023696..3a63de99 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,6 +11,7 @@ set(CMAKE_CXX_STANDARD 20) set(CMAKE_CXX_STANDARD_REQUIRED ON) set(CMAKE_CXX_EXTENSIONS OFF) # only standard c++ + # Check for minimum compiler versions that support most of C++20 if (CMAKE_CXX_COMPILER_ID STREQUAL "GNU") if (CMAKE_CXX_COMPILER_VERSION VERSION_LESS 11) @@ -168,9 +169,40 @@ if (ENABLE_BOOST_POOL) set(USE_BOOST_POOL "#define USE_BOOST_POOL 1") endif() +# choose the number of bits for the bit vector propagation +set(ADOLC_BITWORD_BITS 32) +target_compile_definitions(adolc PUBLIC + ADOLC_BITWORD_BITS=${ADOLC_BITWORD_BITS}) if (ENABLE_SPARSE) - find_package(ColPack REQUIRED) # TODO add colpack installation - # target_link_libraries(...) + if (WIN32) + message(FATAL_ERROR "Sparse drivers require OpenMP, which is not supported on Windows/MSVC. ") + endif() + + find_package(ColPack) + if (ColPack_FOUND) + message(STATUS "Using system-installed ColPack at ${ColPack_DIR}") + + # Prefer shared if it exists (usually the case) + if (TARGET ColPack_static) + target_link_libraries(adolc PRIVATE ColPack_static) + else() + message(FATAL_ERROR "ColPack was found but provides no usable CMake targets") + endif() + endif() + + if (NOT ColPack_FOUND) + message(STATUS "System ColPack not found! Pulling from https://github.com/TimSiebert1/ColPack") + include(FetchContent) + FetchContent_declare( + ColPack + GIT_REPOSITORY https://github.com/TimSiebert1/ColPack + GIT_TAG 96eca73724eb517cfb2cc66d8be448e91a30483f # latest commit + CMAKE_ARGS -DENABLE_OPENMP=ON # colpack needs openmp to build + ) + FetchContent_MakeAvailable(ColPack) + target_link_libraries(adolc PRIVATE ColPack::ColPack_static) + endif() + set(ADOLC_SPARSE "#define ADOLC_SPARSE 1") endif() if(ENABLE_PAREXA) @@ -181,16 +213,14 @@ endif() # Set the macros - -target_compile_definitions(adolc PRIVATE +target_compile_definitions(adolc PUBLIC $<$:ADOLC_TRACK_ACTIVITY=1> $<$:ADOLC_ADVANCED_BRANCHING=1> $<$:USE_ADTL_REFCOUNTING=1> $<$:ENABLE_BOOST_POOL=1> $<$:ADOLC_ADOUBLE_STDCZERO=1> - $<$:SPARSE_DRIVERS=1> - $<$:SPARSE=1> - $<$:ADOLC_TAPE_DOC_VALUES=1>) + $<$:ADOLC_TAPE_DOC_VALUES=1> + $<$:ADOLC_SPARSE=1>) # include subdirectories for handling of includes and source files # ---------------------------------------------------------------- @@ -245,12 +275,16 @@ endif() include(CMakePackageConfigHelpers) include(GNUInstallDirs) -install(TARGETS adolc EXPORT adolcTargets) +install(TARGETS adolc + EXPORT adolcTargets + DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/adolc) + install(EXPORT adolcTargets FILE adolc-targets.cmake NAMESPACE adolc:: DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/adolc) + write_basic_package_version_file(adolc-config-version.cmake VERSION ${adol-c_VERSION} COMPATIBILITY SameMinorVersion) # cmake >= 3.11 @@ -264,7 +298,7 @@ install(FILES ${CMAKE_CURRENT_BINARY_DIR}/adolc-config-version.cmake DESTINATION ${CMAKE_INSTALL_LIBDIR}/cmake/adolc) -export(TARGETS adolc NAMESPACE adolc:: FILE adolc-targets.cmake) + # print a summary of found packages and set options # -------------------------------------------------