Skip to content

Commit

Permalink
Surround oneliners with {}. (#1398)
Browse files Browse the repository at this point in the history
  • Loading branch information
1uc committed Aug 14, 2024
1 parent 32d7e48 commit db16425
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 13 deletions.
21 changes: 14 additions & 7 deletions src/solver/crout/crout.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -50,9 +50,11 @@ EIGEN_DEVICE_FUNC inline int Crout(int n, T* const a, int* const perm, double* c
for (i = 0; i < n; i++) {
perm[i] = i;
k = 0;
for (j = 1; j < n; j++)
if (std::fabs(a[i * n + j]) > std::fabs(a[i * n + k]))
for (j = 1; j < n; j++) {
if (std::fabs(a[i * n + j]) > std::fabs(a[i * n + k])) {
k = j;
}
}
rowmax[i] = a[i * n + k];
}

Expand Down Expand Up @@ -99,8 +101,9 @@ EIGEN_DEVICE_FUNC inline int Crout(int n, T* const a, int* const perm, double* c

/* Check that pivot element is not too small */

if (std::fabs(a[pivot * n + r]) < roundoff)
if (std::fabs(a[pivot * n + r]) < roundoff) {
return -1;
}

/*
* Operate on row in rth position. This produces the upper
Expand Down Expand Up @@ -151,8 +154,9 @@ EIGEN_DEVICE_FUNC inline int solveCrout(int n,
for (i = 0; i < n; i++) {
pivot = perm[i];
sum = 0.0;
for (j = 0; j < i; j++)
for (j = 0; j < i; j++) {
sum += a[pivot * n + j] * (y_(j));
}
y_(i) = (b_(pivot) - sum) / a[pivot * n + i];
}

Expand All @@ -166,16 +170,18 @@ EIGEN_DEVICE_FUNC inline int solveCrout(int n,
for (i = n - 1; i >= 0; i--) {
pivot = perm[i];
sum = 0.0;
for (j = i + 1; j < n; j++)
for (j = i + 1; j < n; j++) {
sum += a[pivot * n + j] * (y_(j));
}
y_(i) -= sum;
}
} else {
for (i = 0; i < n; i++) {
pivot = perm[i];
sum = 0.0;
for (j = 0; j < i; j++)
for (j = 0; j < i; j++) {
sum += a[pivot * n + j] * (p[j]);
}
p[i] = (b_(pivot) - sum) / a[pivot * n + i];
}

Expand All @@ -189,8 +195,9 @@ EIGEN_DEVICE_FUNC inline int solveCrout(int n,
for (i = n - 1; i >= 0; i--) {
pivot = perm[i];
sum = 0.0;
for (j = i + 1; j < n; j++)
for (j = i + 1; j < n; j++) {
sum += a[pivot * n + j] * (p[j]);
}
p[i] -= sum;
}
}
Expand Down
17 changes: 11 additions & 6 deletions src/solver/newton/newton.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -76,13 +76,15 @@ EIGEN_DEVICE_FUNC int newton_solver(Eigen::Matrix<double, N, 1>& X,
// Crout's implementation requires matrices stored in RowMajor order (C-style arrays).
// Therefore, the transposeInPlace is critical such that the data() method to give the rows
// instead of the columns.
if (!J.IsRowMajor)
if (!J.IsRowMajor) {
J.transposeInPlace();
}
Eigen::Matrix<int, N, 1> pivot;
Eigen::Matrix<double, N, 1> rowmax;
// Check if J is singular
if (nmodl::crout::Crout<double>(N, J.data(), pivot.data(), rowmax.data()) < 0)
if (nmodl::crout::Crout<double>(N, J.data(), pivot.data(), rowmax.data()) < 0) {
return -1;
}
Eigen::Matrix<double, N, 1> X_solve;
nmodl::crout::solveCrout<double>(N, J.data(), F.data(), X_solve.data(), pivot.data());
X -= X_solve;
Expand Down Expand Up @@ -156,13 +158,15 @@ EIGEN_DEVICE_FUNC int newton_numerical_diff_solver(Eigen::Matrix<double, N, 1>&
// restore X
X[i] += dX;
}
if (!J.IsRowMajor)
if (!J.IsRowMajor) {
J.transposeInPlace();
}
Eigen::Matrix<int, N, 1> pivot;
Eigen::Matrix<double, N, 1> rowmax;
// Check if J is singular
if (nmodl::crout::Crout<double>(N, J.data(), pivot.data(), rowmax.data()) < 0)
if (nmodl::crout::Crout<double>(N, J.data(), pivot.data(), rowmax.data()) < 0) {
return -1;
}
Eigen::Matrix<double, N, 1> X_solve;
nmodl::crout::solveCrout<double>(N, J.data(), F.data(), X_solve.data(), pivot.data());
X -= X_solve;
Expand Down Expand Up @@ -195,10 +199,11 @@ EIGEN_DEVICE_FUNC int newton_solver_small_N(Eigen::Matrix<double, N, 1>& X,
// The inverse can be called from within OpenACC regions without any issue, as opposed to
// Eigen::PartialPivLU.
J.computeInverseWithCheck(J_inv, invertible);
if (invertible)
if (invertible) {
X -= J_inv * F;
else
} else {
return -1;
}
}
return -1;
}
Expand Down

0 comments on commit db16425

Please sign in to comment.