Skip to content

Commit

Permalink
Fix optional build of DG units (#267)
Browse files Browse the repository at this point in the history
  • Loading branch information
@jbreue16
jbreue16 authored Aug 7, 2024
1 parent c488249 commit fd5ddac
Show file tree
Hide file tree
Showing 11 changed files with 220 additions and 174 deletions.
204 changes: 103 additions & 101 deletions src/libcadet/AdUtils.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -73,65 +73,6 @@ void extractBandedJacobianFromAd(active const* const adVec, int adDirOffset, int
}
}

void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat)
{
const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = 0; eq < mat.rows(); ++eq)
{
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;

// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= 0 && // left block boundary
eq - lowerBandwidth + diag < mat.rows() && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
mat.coeffRef(eq, eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir);

// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
++dir;
}
}
}

#ifdef ENABLE_DG
void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nRows, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
{
const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = blockOffset; eq < blockOffset + nRows; ++eq)
{
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;

// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary
eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
mat.coeffRef(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir);

// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
++dir;
}
}
}
#endif

void prepareAdVectorSeedsForDenseMatrix(active* const adVec, int adDirOffset, int cols)
{
for (int col = 0; col < cols; ++col)
Expand Down Expand Up @@ -285,63 +226,124 @@ double compareDenseJacobianWithBandedAd(active const* const adVec, int row, int
return maxDiff;
}

double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nRows, const Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
void adMatrixVectorMultiply(const linalg::SparseMatrix<active>& mat, double const* x, double* y, double alpha, double beta, int adDir)
{
double maxDiff = 0.0;
const std::vector<int>& rows = mat.rows();
const std::vector<int>& cols = mat.cols();
const std::vector<active>& values = mat.values();
const int numNonZero = mat.numNonZero();

const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = blockOffset; eq < blockOffset + nRows; ++eq)
for (int i = 0; i < numNonZero; ++i)
{
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;
y[rows[i]] = alpha * values[i].getADValue(adDir) * x[cols[i]] + beta * y[rows[i]];
}
}

// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
#ifdef ENABLE_DG

void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat)
{
const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = 0; eq < mat.rows(); ++eq)
{
if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary
eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
{
double baseVal = adVec[eq].getADValue(adDirOffset + dir);
double matVal = mat.coeff(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag);
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;

if (std::isnan(matVal) || std::isnan(baseVal))
return std::numeric_limits<double>::quiet_NaN();
const double diff = std::abs(matVal - baseVal);
// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= 0 && // left block boundary
eq - lowerBandwidth + diag < mat.rows() && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
mat.coeffRef(eq, eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir);

baseVal = std::abs(baseVal);
if (baseVal > 0.0)
maxDiff = std::max(maxDiff, diff / baseVal);
// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
maxDiff = std::max(maxDiff, diff);
++dir;
}

// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
++dir;
}
}
return maxDiff;
}

void adMatrixVectorMultiply(const linalg::SparseMatrix<active>& mat, double const* x, double* y, double alpha, double beta, int adDir)
{
const std::vector<int>& rows = mat.rows();
const std::vector<int>& cols = mat.cols();
const std::vector<active>& values = mat.values();
const int numNonZero = mat.numNonZero();
void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nRows, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
{
const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = blockOffset; eq < blockOffset + nRows; ++eq)
{
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;

for (int i = 0; i < numNonZero; ++i)
// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary
eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
mat.coeffRef(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag) = adVec[eq].getADValue(adDirOffset + dir);

// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
++dir;
}
}
}

double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nRows, const Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset)
{
y[rows[i]] = alpha * values[i].getADValue(adDir) * x[cols[i]] + beta * y[rows[i]];
double maxDiff = 0.0;

const int stride = lowerBandwidth + 1 + upperBandwidth;
for (int eq = blockOffset; eq < blockOffset + nRows; ++eq)
{
// Start with lowest subdiagonal and stay in the range of the columns:
// diagDir might be too big for the matrix and, hence, dir ranges between
// diagDir - lowerBandwidth <= dir <= diagDir + upperBandwidth
int dir = diagDir - lowerBandwidth + eq % stride;

// Loop over diagonals
for (int diag = 0; diag < stride; ++diag)
{
if (eq - lowerBandwidth + diag >= blockOffset && // left block boundary
eq - lowerBandwidth + diag < blockOffset + nRows && // right block boundary
adVec[eq].getADValue(adDirOffset + dir) != 0.0 // keep pattern
)
{
double baseVal = adVec[eq].getADValue(adDirOffset + dir);
double matVal = mat.coeff(matrixOffset + eq, matrixOffset + eq - lowerBandwidth + diag);

if (std::isnan(matVal) || std::isnan(baseVal))
return std::numeric_limits<double>::quiet_NaN();
const double diff = std::abs(matVal - baseVal);

baseVal = std::abs(baseVal);
if (baseVal > 0.0)
maxDiff = std::max(maxDiff, diff / baseVal);
else
maxDiff = std::max(maxDiff, diff);
}

// Wrap around at end of row and jump to lowest subdiagonal
if (dir == diagDir + upperBandwidth)
dir = diagDir - lowerBandwidth;
else
++dir;
}
}
return maxDiff;
}
}

#endif

} // namespace ad

Expand Down
67 changes: 34 additions & 33 deletions src/libcadet/AdUtils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -85,20 +85,41 @@ void extractBandedJacobianFromAd(active const* const adVec, int adDirOffset, int
*/
void extractBandedBlockEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth,
const int blockOffset, const int nCols, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset = 0);
#endif

/**
* @brief Extracts a band matrix (Eigen lib) from band compressed AD seed vectors
* @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to
assemble the Jacobian block which is a band matrix. The block must be on the main diagonal.
* @param [in] adVec Vector of AD datatypes with band compressed seed vectors
* @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions)
* @param [in] diagDir Diagonal direction index
* @param [in] lowerBandwidth lower band width
* @param [in] upperBandwidth upper band width
* @param [out] mat Eigen matrix to be populated with the Jacobian block
*/
void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat);
/**
* @brief Extracts a band matrix (Eigen lib) from band compressed AD seed vectors
* @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to
assemble the Jacobian block which is a band matrix. The block must be on the main diagonal.
* @param [in] adVec Vector of AD datatypes with band compressed seed vectors
* @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions)
* @param [in] diagDir Diagonal direction index
* @param [in] lowerBandwidth lower band width
* @param [in] upperBandwidth upper band width
* @param [out] mat Eigen matrix to be populated with the Jacobian block
*/
void extractBandedEigenJacobianFromAd(active const* const adVec, int adDirOffset, int diagDir, const int lowerBandwidth, const int upperBandwidth, Eigen::SparseMatrix<double, Eigen::RowMajor>& mat);

/**
* @brief Compares a (block-) banded Jacobian in Eigen library row-major format with an AD version derived by band compressed AD seed vectors
* @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to
compare the results with a given banded Jacobian. The AD Jacobian is treated as base and the analytic
Jacobian is compared against it. The relative difference
@f[ \Delta_{ij} = \begin{cases} \left\lvert \frac{ J_{\text{ana},ij} - J_{\text{ad},ij} }{ J_{\text{ad},ij} }\right\rvert, & J_{\text{ad},ij} \neq 0 \\
\left\lvert J_{\text{ana},ij} - J_{\text{ad},ij} \right\rvert, & J_{\text{ad},ij} = 0 \end{cases} @f]
is computed for each matrix entry. The maximum of all @f$ \Delta_{ij} @f$ is returned.
* @param [in] adVec Vector of AD datatypes with band compressed seed vectors
* @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions)
* @param [in] diagDir Diagonal direction index
* @param [in] lowerBandwidth of the block
* @param [in] upperBandwidth of the block
* @param [in] blockOffset of the currently considered block
* @param [in] nRows of the matrix or number of equations to be compared
* @param [in] mat BandMatrix populated with the analytic Jacobian
* @return The maximum absolute relative difference between the matrix elements
*/
double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, const int blockOffset, const int nRows, const Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset);

#endif

/**
* @brief Sets seed vectors on an AD vector for computing a dense Jacobian
Expand Down Expand Up @@ -151,26 +172,6 @@ void extractDenseJacobianFromBandedAd(active const* const adVec, int row, int ad
*/
double compareBandedJacobianWithAd(active const* const adVec, int adDirOffset, int diagDir, const linalg::BandMatrix& mat);

/**
* @brief Compares a (block-) banded Jacobian in Eigen library row-major format with an AD version derived by band compressed AD seed vectors
* @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForBandMatrix() to
compare the results with a given banded Jacobian. The AD Jacobian is treated as base and the analytic
Jacobian is compared against it. The relative difference
@f[ \Delta_{ij} = \begin{cases} \left\lvert \frac{ J_{\text{ana},ij} - J_{\text{ad},ij} }{ J_{\text{ad},ij} }\right\rvert, & J_{\text{ad},ij} \neq 0 \\
\left\lvert J_{\text{ana},ij} - J_{\text{ad},ij} \right\rvert, & J_{\text{ad},ij} = 0 \end{cases} @f]
is computed for each matrix entry. The maximum of all @f$ \Delta_{ij} @f$ is returned.
* @param [in] adVec Vector of AD datatypes with band compressed seed vectors
* @param [in] adDirOffset Offset in the AD directions (can be used to move past parameter sensitivity directions)
* @param [in] diagDir Diagonal direction index
* @param [in] lowerBandwidth of the block
* @param [in] upperBandwidth of the block
* @param [in] blockOffset of the currently considered block
* @param [in] nRows of the matrix or number of equations to be compared
* @param [in] mat BandMatrix populated with the analytic Jacobian
* @return The maximum absolute relative difference between the matrix elements
*/
double compareBandedEigenJacobianWithAd(active const* const adVec, const int adDirOffset, const int diagDir, const int lowerBandwidth, const int upperBandwidth, const int blockOffset, const int nRows, const Eigen::SparseMatrix<double, Eigen::RowMajor>& mat, const int matrixOffset);

/**
* @brief Compares a dense Jacobian with an AD version derived by AD seed vectors
* @details Uses the results of an AD computation with seed vectors set by prepareAdVectorSeedsForDenseMatrix() to
Expand Down
Loading

0 comments on commit fd5ddac

Please sign in to comment.