sparsedrivers.h

Go to the documentation of this file.

/*----------------------------------------------------------------------------
 ADOL-C -- Automatic Differentiation by Overloading in C++
 File:     sparse/sparsedrivers.h
 Revision: $Id$
 Contents: This file contains some "Easy To Use" interfaces of the sparse
package.
 
 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.
 
----------------------------------------------------------------------------*/
#ifndef ADOLC_SPARSE_DRIVERS_H
#define ADOLC_SPARSE_DRIVERS_H
#include <adolc/adalloc.h>
#include <adolc/adolcerror.h>
#include <adolc/interfaces.h>
#include <adolc/internal/common.h>
#include <adolc/sparse/sparse_options.h>
#include <adolc/sparse/sparsematrix.h>
#include <adolc/tape_interface.h>
#include <adolc/valuetape/sparseinfos.h>
#include <adolc/valuetape/valuetape.h>
#include <cmath>
#include <cstddef>
#include <span>
#include <utility>
#include <vector>
 
// 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 {
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 int BITS_PER_WORD = 8 * sizeof(bitword_t);
static constexpr bitword_t MOST_SIGNIFICANT_BIT = static_cast<bitword_t>(1)
                                                  << (BITS_PER_WORD - 1);
 
template <BitPatternPropagationDirection BPPD> struct BvpData {
  int rc_{0};            
  int depen_{-1};        
  int indep_{-1};        
  int wordsPerBatch_{0}; 
  int bitsPerStrip_{0};  
  int numStripmineBatches_{
      0}; 
  std::vector<unsigned char>
      indepWordHasNonzero_;        
  std::vector<double> valuepoint_; 
  std::vector<bitword_t *> seed_;  
  std::vector<bitword_t *>
      jacobianBitPattern_; 
 
  BvpData(const BvpData &) = delete;
  BvpData &operator=(const BvpData &) = delete;
  BvpData(BvpData &&) = delete;
  BvpData &operator=(BvpData &&) = delete;
 
  ~BvpData() {
    for (auto &s : seed_)
      delete[] s;
    for (auto &j : jacobianBitPattern_)
      delete[] j;
  }
 
  BvpData(int depen, int indep) : depen_(depen), indep_(indep) {
    int 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 constexpr (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_);
    }
  }
};
 
template <ControlFlowMode CFM>
void checkBVPInput([[maybe_unused]] const double *basepoint) {};
 
inline void resetInput(std::span<uint *> &compressedRowStorage) {
  for (auto &row : compressedRowStorage) {
    row = nullptr;
  }
}
 
template <BitPatternPropagationDirection BPPD>
void resetOldSeed(BvpData<BPPD> &data, size_t dim2) {
  for (auto &s : data.seed_)
    for (size_t j = 0; j < dim2; ++j)
      s[j] = 0;
}
 
template <BitPatternPropagationDirection BPPD>
void prepareSeedMatrix(size_t stripIdx, BvpData<BPPD> &data) {
  if constexpr (BPPD == BitPatternPropagationDirection::Forward) {
    resetOldSeed(data, data.wordsPerBatch_);
    createSeed(data, stripIdx, data.indep_);
  } else if constexpr (BPPD == BitPatternPropagationDirection::Reverse) {
    resetOldSeed(data, data.depen_);
    createSeed(data, stripIdx, data.depen_);
  }
}
 
template <BitPatternPropagationDirection BPPD>
void setSeed(BvpData<BPPD> &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 constexpr (BPPD == BitPatternPropagationDirection::Reverse)
    data.seed_[(idx - currentStripWordIdx) / BITS_PER_WORD][idx] =
        MOST_SIGNIFICANT_BIT >> ((idx - currentStripWordIdx) % BITS_PER_WORD);
}
 
template <BitPatternPropagationDirection BPPD>
void createSeed(BvpData<BPPD> &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);
}
 
void extract(size_t wordIdx, size_t stripIdx, size_t currentStripWordIdx,
             BvpData<BitPatternPropagationDirection::Forward> &data,
             std::span<uint *> &compressedRowStorage);
 
void extract(size_t wordIdx, size_t stripIdx,
             BvpData<BitPatternPropagationDirection::Reverse> &data,
             std::span<uint *> &compressedRowStorage);
 
template <BitPatternPropagationDirection BPPD>
void extractCompressedRowStorage(int stripIdx,
                                 std::span<uint *> &compressedRowStorage,
                                 BvpData<BPPD> &data) {
  int currentStripWordIdx = stripIdx * data.bitsPerStrip_;
  int nextStripWordIdx = (stripIdx + 1) * data.bitsPerStrip_;
  if constexpr (BPPD == BitPatternPropagationDirection::Forward)
    for (int wordIdx = 0; wordIdx < data.depen_; wordIdx++) {
      int idx = 0;
      bitword_t currentBit = MOST_SIGNIFICANT_BIT;
      for (int 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 constexpr (BPPD == BitPatternPropagationDirection::Reverse) {
    if (nextStripWordIdx > data.depen_)
      nextStripWordIdx = data.depen_;
    int idx = 0;
    bitword_t currentBit = MOST_SIGNIFICANT_BIT;
    for (int 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);
    }
  }
}
 
template <BitPatternPropagationDirection BPPD, ControlFlowMode CFM>
void prepareADCalls(short tapeId, const double *basepoint,
                    BvpData<BPPD> &data) {
  if constexpr (BPPD == BitPatternPropagationDirection::Reverse &&
                CFM == ControlFlowMode::Tight)
    data.rc_ = zos_forward(tapeId, data.depen_, data.indep_, 1, basepoint,
                           data.valuepoint_.data());
}
 
template <BitPatternPropagationDirection BPPD, ControlFlowMode CFM>
void ADCalls(short tapeId, const double *basepoint, BvpData<BPPD> &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 constexpr (CFM == ControlFlowMode::Safe)
      data.rc_ = int_forward_safe(tapeId, data.depen_, data.indep_,
                                  data.wordsPerBatch_, data.seed_.data(),
                                  data.jacobianBitPattern_.data());
 
  } else if constexpr (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 constexpr (CFM == ControlFlowMode::Safe)
      data.rc_ = int_reverse_safe(tapeId, data.depen_, data.indep_,
                                  data.wordsPerBatch_, data.seed_.data(),
                                  data.jacobianBitPattern_.data());
  }
}
 
template <ControlFlowMode CFM, BitPatternPropagationDirection BPPD>
int bitVectorPropagation(short tapeId, int depen, int indep,
                         const double *basepoint,
                         std::span<uint *> &compressedRowStorage) {
  using namespace detail;
  checkBVPInput<CFM>(basepoint);
  BvpData<BPPD> data{depen, indep};
  prepareADCalls<BPPD, CFM>(tapeId, basepoint, data);
  for (int stripIdx = 0; stripIdx < data.numStripmineBatches_; stripIdx++) {
    prepareSeedMatrix(stripIdx, data);
    ADCalls<BPPD, CFM>(tapeId, basepoint, data);
    extractCompressedRowStorage(stripIdx, compressedRowStorage, data);
  }
  return data.rc_;
}
 
} // namespace detail
 
template <ControlFlowMode CFM, BitPatternPropagationDirection BPPD =
                                   BitPatternPropagationDirection::Auto>
ADOLC_API int bit_vector_propagation(short tapeId, int depen, int indep,
                                     const double *basepoint,
                                     std::span<uint *> &compressedRowStorage) {
  using namespace detail;
  if constexpr (BPPD == BitPatternPropagationDirection::Auto) {
    // we allow here run-time selection
    if (depen >= indep / 2)
      return bitVectorPropagation<CFM, BitPatternPropagationDirection::Forward>(
          tapeId, depen, indep, basepoint, compressedRowStorage);
 
    else
      return bitVectorPropagation<CFM, BitPatternPropagationDirection::Reverse>(
          tapeId, depen, indep, basepoint, compressedRowStorage);
  } else
    return bitVectorPropagation<CFM, BPPD>(tapeId, depen, indep, basepoint,
                                           compressedRowStorage);
}
 
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<uint *> &compressedRowStorage) {
  detail::resetInput(compressedRowStorage);
  if constexpr (SM == SparseMethod::IndexDomains &&
                CFM == ControlFlowMode::Tight)
    return indopro_forward_tight(tag, depen, indep, basepoint,
                                 compressedRowStorage.data());
  else if constexpr (SM == SparseMethod::IndexDomains &&
                     CFM == ControlFlowMode::Safe)
    return indopro_forward_safe(tag, depen, indep, basepoint,
                                compressedRowStorage.data());
  else if constexpr (SM == SparseMethod::BitPattern)
    return bit_vector_propagation<CFM, BPPD>(tag, depen, indep, basepoint,
                                             compressedRowStorage);
};
 
template <PiecewiseLinear>
ADOLC_API int jac_pat(short tag, int depen, int indep, int numsw,
                      const double *basepoint,
                      std::span<uint *> &compressedRowStorage) {
  detail::resetInput(compressedRowStorage);
  return indopro_forward_absnormal(tag, depen, indep, numsw, basepoint,
                                   compressedRowStorage.data());
}
template <CompressionMode CM>
ADOLC_API void generate_seed_jac(int m, int n, const std::span<uint *> JP,
                                 double ***Seed, int *p) {
  generateSeedJac<CM>(m, n, JP, Seed, p);
}
 
namespace detail {
 
template <
    SparseMethod SM, CompressionMode CM, ControlFlowMode CFM,
    BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto>
int buildJacPatternAndSeed(short tag, int depen, int indep,
                           const double *basepoint) {
 
  ValueTape &tape = findTape(tag);
  tape.sJInfos().setJP(std::vector<uint *>(depen));
  std::span<uint *> JPSpan = tape.sJInfos().getJP();
  int ret_val = jac_pat<SM, CFM, BPPD>(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 (uint j = 1; j <= tape.sJInfos().JP_[i][0]; j++)
      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 constexpr (CM == CompressionMode::Column) {
    tape.sJInfos().generateSeedJac("COLUMN_PARTIAL_DISTANCE_TWO");
    tape.sJInfos().seedRows_ = depen;
    ret_val = tape.sJInfos().seedClms_;
  }
  return ret_val;
}
 
template <PiecewiseLinear>
int buildJacPatternAndSeed(short tag, int depen, int indep, int numSwitches,
                           const double *basepoint) {
  auto &tape = findTape(tag);
  tape.sJInfos().setJP(std::vector<uint *>(depen + numSwitches));
  std::span<uint *> JPSpan{tape.sJInfos().getJP()};
 
  int ret_val = jac_pat<PiecewiseLinear{}>(tag, depen, indep, numSwitches,
                                           basepoint, JPSpan);
  if (ret_val < 0) {
    printf(" ADOL-C error in sparse_jac() \n");
    return ret_val;
  }
 
  tape.sJInfos().depen_ = depen + numSwitches;
  tape.sJInfos().nnzIn_ = 0;
  for (int i = 0; i < depen + numSwitches; i++) {
    for (uint j = 1; j <= tape.sJInfos().getJP()[i][0]; j++) {
      tape.sJInfos().nnzIn_++;
    }
  }
  tape.sJInfos().initColoring(depen + numSwitches, indep + numSwitches);
 
  tape.sJInfos().generateSeedJac("ROW_PARTIAL_DISTANCE_TWO");
  tape.sJInfos().seedClms_ = indep + numSwitches;
  ret_val = tape.sJInfos().seedRows_;
 
  return ret_val;
}
namespace {
inline int ensureSparseJacCacheInitialized(const ValueTape &tape) {
  if (tape.sJInfos().Seed_ == nullptr || tape.sJInfos().seedRows_ == 0 ||
      tape.sJInfos().seedClms_ == 0) {
    printf(" ADOL-C error in sparse_jac():"
           " First call with repeat = 0 \n");
    return -3;
  }
  return 0;
}
 
inline int validateManualStorage(const ValueTape &tape,
                                 const SparseMatrix &sparseJac) {
  if (sparseJac.size() != to_size_t(tape.sJInfos().nnzIn_)) {
    printf(" ADOL-C error in sparse_jac():"
           " SparseMatrix size mismatch for MemoryHandler::Manual,"
           " expected %d entries but got %zu \n",
           tape.sJInfos().nnzIn_, sparseJac.size());
    return -3;
  }
  return 0;
}
 
inline SparseShape countSparseANFEntries(const SparseJacInfos &sJInfos,
                                         int depen, int indep) {
  using Coordinates = CoordinateFormatTripled::Coordinates;
  using coordinate_type = CoordinateFormatTripled::coordinate_type;
  using detail::classifySparseANFBlock;
 
  SparseShape counts;
  for (coordinate_type rowIndex = 0; rowIndex < sJInfos.JP_.size();
       ++rowIndex) {
    const auto numOfNonZeros = sJInfos.JP_[rowIndex][0];
    for (unsigned int j = 1; j <= numOfNonZeros; ++j) {
      const Coordinates coords{rowIndex, sJInfos.JP_[rowIndex][j]};
      switch (classifySparseANFBlock(coords, depen, indep)) {
      case SparseANFBlock::Y:
        ++counts.y;
        break;
      case SparseANFBlock::J:
        ++counts.j;
        break;
      case SparseANFBlock::Z:
        ++counts.z;
        break;
      case SparseANFBlock::L:
        ++counts.l;
        break;
      }
    }
  }
  return counts;
}
 
inline int validateManualStorage(const ValueTape &tape,
                                 const SparseANF &sparseANF, int depen,
                                 int indep) {
  const auto counts = countSparseANFEntries(tape.sJInfos(), depen, indep);
  if (sparseANF.dims().y != counts.y || sparseANF.dims().j != counts.j ||
      sparseANF.dims().z != counts.z || sparseANF.dims().l != counts.l) {
    printf(" ADOL-C error in sparse_jac():"
           " SparseANF block sizes mismatch for MemoryHandler::Manual"
           " (expected Y=%zu J=%zu Z=%zu L=%zu,"
           " got Y=%zu J=%zu Z=%zu L=%zu) \n",
           counts.y, counts.j, counts.z, counts.l, sparseANF.dims().y,
           sparseANF.dims().j, sparseANF.dims().z, sparseANF.dims().l);
    return -3;
  }
  return 0;
}
} // namespace
template <CompressionMode CM, MemoryHandler MH>
void recover(const ValueTape &tape, SparseMatrix &sparseJac) {
  if constexpr (MH == MemoryHandler::Manual) {
    assert(sparseJac.size() == to_size_t(tape.sJInfos().nnzIn_) &&
           "Mismatch in allocated size and required size!");
    if constexpr (CM == CompressionMode::Row) {
      tape.sJInfos().recoverRowFormatUserMem(sparseJac);
    } else if constexpr (CM == CompressionMode::Column) {
      tape.sJInfos().recoverColFormatUserMem(sparseJac);
    }
  } else {
    if constexpr (CM == CompressionMode::Row) {
      tape.sJInfos().recoverRowFormat(sparseJac);
    } else if constexpr (CM == CompressionMode::Column) {
      tape.sJInfos().recoverColFormat(sparseJac);
    }
  }
}
 
template <MemoryHandler MH>
void recover(const ValueTape &tape, SparseANF &sparseANF) {
  const int numInds = static_cast<int>(tape.numInds());
  const int numSwitches = static_cast<int>(tape.numSwitches());
  if constexpr (MH == MemoryHandler::Manual) {
    tape.sJInfos().recoverANFUserMem(sparseANF, numInds, numSwitches);
  } else {
    tape.sJInfos().recoverANF(sparseANF, numInds, numSwitches);
  }
}
 
template <CompressionMode CM, MemoryHandler MH>
int computeSparseJac(short tag, int depen, int indep, const double *basepoint,
                     SparseMatrix &sparseJac) {
  ValueTape &tape = findTape(tag);
  int ret_val = 0;
  myfree2(tape.sJInfos().B_);
  myfree1(tape.sJInfos().y_);
  tape.sJInfos().B_ =
      myalloc2(tape.sJInfos().seedRows_, tape.sJInfos().seedClms_);
  tape.sJInfos().y_ = myalloc1(depen);
 
  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 constexpr (CM == CompressionMode::Column)
    ret_val =
        fov_forward(tag, depen, indep, tape.sJInfos().seedClms_, basepoint,
                    tape.sJInfos().Seed_, tape.sJInfos().y_, tape.sJInfos().B_);
  recover<CM, MH>(tape, sparseJac);
  return ret_val;
}
 
template <MemoryHandler MH>
int computeSparseJac(short tag, int depen, int indep, int numSwitches,
                     const double *basepoint, SparseANF &sparseANF) {
  ValueTape &tape = findTape(tag);
  assert(tape.sJInfos().seedClms_ == indep + numSwitches &&
         "Row Compression must have a seed with column number equal to indep + "
         "numSwitches!");
  int ret_val = 0;
  myfree2(tape.sJInfos().B_);
  tape.sJInfos().B_ =
      myalloc2(tape.sJInfos().seedRows_, tape.sJInfos().seedClms_);
 
  std::vector<double *> results(tape.sJInfos().seedRows_);
  std::vector<double *> resultsSwitch(tape.sJInfos().seedRows_);
 
  std::vector<double *> seedResults(tape.sJInfos().seedRows_);
  std::vector<double *> seedSwitches(tape.sJInfos().seedRows_);
 
  for (int i = 0; i < tape.sJInfos().seedRows_; i++) {
    results[i] = tape.sJInfos().B_[i];
    resultsSwitch[i] = tape.sJInfos().B_[i] + indep;
    seedResults[i] = tape.sJInfos().Seed_[i];
    seedSwitches[i] = tape.sJInfos().Seed_[i] + depen;
  }
 
  sparseANF.y.resize(depen);
  sparseANF.z.resize(numSwitches);
  ret_val = zos_pl_forward(tag, depen, indep, 1, basepoint, sparseANF.y.data(),
                           sparseANF.z.data());
  if (ret_val < 0)
    return ret_val;
  MINDEC(ret_val, fov_pl_reverse(tag, depen, indep, numSwitches,
                                 tape.sJInfos().seedRows_, seedResults.data(),
                                 seedSwitches.data(), results.data(),
                                 resultsSwitch.data()));
 
  recover<MH>(tape, sparseANF);
 
  // update consts
  sparseANF.updateCy();
  sparseANF.updateCz();
 
  return ret_val;
}
} // namespace detail
 
template <
    SparseMethod SM = SparseMethod::IndexDomains,
    CompressionMode CM = CompressionMode::Column,
    ControlFlowMode CFM = ControlFlowMode::Safe,
    BitPatternPropagationDirection BPPD = BitPatternPropagationDirection::Auto,
    MemoryHandler MH = MemoryHandler::Auto>
ADOLC_API int sparse_jac(short tag, int depen, int indep, int repeat,
                         const double *basepoint, SparseMatrix &sparseJac) {
  using namespace detail;
  int ret_val = 0;
  if (repeat == 0)
    ret_val =
        buildJacPatternAndSeed<SM, CM, CFM, BPPD>(tag, depen, indep, basepoint);
  else
    ret_val = ensureSparseJacCacheInitialized(findTape(tag));
  if (ret_val < 0) {
    printf(" ADOL-C error in sparse_jac() \n");
    return ret_val;
  }
  if constexpr (MH == MemoryHandler::Manual) {
 
    ret_val = validateManualStorage(findTape(tag), sparseJac);
    if (ret_val < 0)
      return ret_val;
  }
  return computeSparseJac<CM, MH>(tag, depen, indep, basepoint, sparseJac);
}
template <PiecewiseLinear PL = PiecewiseLinear{},
          MemoryHandler MH = MemoryHandler::Auto>
ADOLC_API int sparse_jac(short tag, int depen, int indep, int numSwitches,
                         int repeat, const double *basepoint,
                         SparseANF &sparseANF) {
 
  using namespace detail;
  int ret_val = 0;
  if (repeat == 0)
    ret_val =
        buildJacPatternAndSeed<PL>(tag, depen, indep, numSwitches, basepoint);
  else
    ret_val = ensureSparseJacCacheInitialized(findTape(tag));
  if (ret_val < 0) {
    printf(" ADOL-C error in sparse_jac() \n");
    return ret_val;
  }
  if constexpr (MH == MemoryHandler::Manual) {
    ret_val = validateManualStorage(findTape(tag), sparseANF, depen, indep);
    if (ret_val < 0)
      return ret_val;
  }
  return computeSparseJac<MH>(tag, depen, indep, numSwitches, basepoint,
                              sparseANF);
}
 
template <ControlFlowMode CFM>
ADOLC_API int hess_pat(short tag, int indep, const double *basepoint,
                       std::span<uint *> &compressedRowStorage) {
  detail::resetInput(compressedRowStorage);
  if constexpr (CFM == ControlFlowMode::OldTight)
    return nonl_ind_old_forward_tight(tag, 1, indep, basepoint,
                                      compressedRowStorage.data());
  else if constexpr (CFM == ControlFlowMode::OldSafe)
    return nonl_ind_old_forward_safe(tag, 1, indep, basepoint,
                                     compressedRowStorage.data());
  else if constexpr (CFM == ControlFlowMode::Tight)
    return nonl_ind_forward_tight(tag, 1, indep, basepoint,
                                  compressedRowStorage.data());
  else if constexpr (CFM == ControlFlowMode::Safe)
    return nonl_ind_forward_safe(tag, 1, indep, basepoint,
                                 compressedRowStorage.data());
}
 
template <RecoveryMethod RCM>
ADOLC_API void generate_seed_hess(int n, std::span<uint *> HP, double ***Seed,
                                  int *p) {
  if constexpr (RCM == RecoveryMethod::Indirect)
    generateSeedHess(n, HP, Seed, p, "ACYCLIC_FOR_INDIRECT_RECOVERY");
  else if constexpr (RCM == RecoveryMethod::Direct)
    generateSeedHess(n, HP, Seed, p, "STAR");
}
 
namespace detail {
 
template <ControlFlowMode CFM, RecoveryMethod RCM>
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<uint *>(indep));
  std::span<uint *> HPSpan = tape.sHInfos().HP_;
  // generate sparsity pattern
  ret_val = hess_pat<CFM>(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 (uint i = 0; i < static_cast<uint>(indep); i++) {
    for (uint 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
 
  // Seed is handled by ColPack
  double **Seed = nullptr;
  tape.sHInfos().initColoring(indep);
  if constexpr (RCM == RecoveryMethod::Indirect)
    tape.sHInfos().generateSeedHess(&Seed, "ACYCLIC_FOR_INDIRECT_RECOVERY");
  else if constexpr (RCM == RecoveryMethod::Direct)
    tape.sHInfos().generateSeedHess(&Seed, "STAR");
 
  // data might still be allocated, ensure that its not leaked
  myfree2(tape.sHInfos().Hcomp_);
  myfree3(tape.sHInfos().Xppp_);
  myfree3(tape.sHInfos().Yppp_);
  myfree3(tape.sHInfos().Zppp_);
  myfree2(tape.sHInfos().Upp_);
 
  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];
 
  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;
}
 
template <RecoveryMethod RCM>
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 constexpr (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) {
      delete[] *values;
      *values = nullptr;
    }
    if (*rind != nullptr) {
      delete[] *rind;
      *rind = nullptr;
    }
    if (*cind != nullptr) {
      delete[] *cind;
      *cind = nullptr;
    }
    if constexpr (RCM == RecoveryMethod::Indirect)
      tape.sHInfos().indirectRecover(rind, cind, values);
    else if constexpr (RCM == RecoveryMethod::Direct)
      tape.sHInfos().directRecover(rind, cind, values);
  }
  return ret_val;
}
 
} // namespace detail
 
template <ControlFlowMode CFM, RecoveryMethod RCM>
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<CFM, RCM>(tag, indep, basepoint, nnz);
  if (ret_val < 0) {
    printf(" ADOL-C error in sparse_hess() \n");
    return ret_val;
  }
  return computeSparseHess<RCM>(tag, indep, basepoint, nnz, rind, cind, values);
}
 
/****************************************************************************/
 
} // namespace ADOLC::Sparse
#endif // ADOLC_SPARSE_DRIVERS_H