Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix optional build of DG units #267

Merged
merged 1 commit into from
Aug 7, 2024
+220 −174
Merged

Fix optional build of DG units #267

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
204 changes: 103 additions & 101 deletions src/libcadet/AdUtils.cpp
View file
Open in desktop
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
Open in desktop
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
Loading