diff --git a/ALFI/ALFI/spline/cubic.h b/ALFI/ALFI/spline/cubic.h index 8781cda..c290f23 100644 --- a/ALFI/ALFI/spline/cubic.h +++ b/ALFI/ALFI/spline/cubic.h @@ -115,52 +115,29 @@ namespace alfi::spline { coeffs = util::spline::simple_spline(X, Y, 3); return; } - - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); - - Container B(n); - - Container diag(n - 2), right(n - 2); - diag[0] = ((dX[0] + dX[1]) * (dX[0] + 2 * dX[1])) / dX[1]; - for (SizeT i = 1; i < n - 3; ++i) { - diag[i] = 2 * (dX[i] + dX[i+1]); - } - diag[n-3] = ((dX[n-2] + dX[n-3]) * (dX[n-2] + 2 * dX[n-3])) / dX[n-3]; - for (SizeT i = 0; i < n - 2; ++i) { - right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + Container lower(n), diag(n), upper(n), right(n); + for (SizeT i = 1; i < n - 1; ++i) { + lower[i] = dX[i-1]; + diag[i] = 2 * (dX[i-1] + dX[i]); + upper[i] = dX[i]; + right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]); } - - const auto m_1 = dX[1] / diag[0]; - diag[1] -= m_1 * (dX[1] - dX[0]*dX[0]/dX[1]); - right[1] -= m_1 * right[0]; - for (SizeT i = 2; i < n - 3; ++i) { - const auto m = dX[i] / diag[i-1]; - diag[i] -= m * dX[i]; - right[i] -= m * right[i-1]; - } - const auto m_n_3 = (dX[n-3] - dX[n-2]*dX[n-2]/dX[n-3]) / diag[n-4]; - diag[n-3] -= m_n_3 * dX[n-3]; - right[n-3] -= m_n_3 * right[n-4]; - - B[n-2] = right[n-3] / diag[n-3]; - for (SizeT i = n - 3; i >= 2; --i) { - B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1]; + lower[0] = NAN; + diag[0] = dX[0] - dX[1]; + upper[0] = 2*dX[0] + dX[1]; + right[0] = dX[0] / (dX[0]+dX[1]) * right[1]; + lower[n-1] = 2*dX[n-2] + dX[n-3]; + diag[n-1] = dX[n-2] - dX[n-3]; + upper[n-1] = NAN; + right[n-1] = dX[n-2] / (dX[n-2]+dX[n-3]) * right[n-2]; + const auto C = util::linalg::tridiag_solve(lower, diag, upper, right); + for (SizeT i = 0, index = 0; i < n - 1; ++i) { + coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]); + coeffs[index++] = C[i]; + coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3); + coeffs[index++] = Y[i]; } - B[1] = (right[0] - (dX[1] - dX[0]*dX[0]/dX[1]) * B[2]) / diag[0]; - B[0] = B[1] + dX[0]/dX[1] * (B[1] - B[2]); - B[n-1] = B[n-2] + dX[n-2]/dX[n-3] * (B[n-2] - B[n-3]); - - Container A(n - 1); - for (SizeT i = 1; i < A.size() - 1; ++i) { - A[i] = (B[i+1] - B[i]) / (3*dX[i]); - } - A[0] = A[1]; - A[n-2] = A[n-3]; - Container C(n - 1); - for (SizeT i = 0; i < C.size(); ++i) { - C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3); - } - util::spline::merge_coeffs(coeffs, {A, B, C, Y}); }, [&](const typename Types::Periodic&) { if (n <= 2) { @@ -168,13 +145,13 @@ namespace alfi::spline { return; } - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); - Container B(n); + Container C(n); if (n == 3) { - B[0] = 3 * (dY[0]/dX[0] - dY[1]/dX[1]) / (dX[0] + dX[1]); - B[1] = -B[0]; - B[2] = B[0]; + C[0] = 3 * (dY[0]/dX[0] - dY[1]/dX[1]) / (dX[0] + dX[1]); + C[1] = -C[0]; + C[2] = C[0]; } else { Container diag(n - 1), right(n - 1); for (SizeT i = 0; i < n - 2; ++i) { @@ -204,162 +181,114 @@ namespace alfi::spline { diag[n-2] -= last_row[n-3] / diag[n-3] * last_col[n-3]; right[n-2] -= last_row[n-3] / diag[n-3] * right[n-3]; - B[n-1] = right[n-2] / diag[n-2]; - B[n-2] = (right[n-3] - last_col[n-3] * B[n-1]) / diag[n-3]; + C[n-1] = right[n-2] / diag[n-2]; + C[n-2] = (right[n-3] - last_col[n-3] * C[n-1]) / diag[n-3]; for (SizeT i = n - 3; i >= 1; --i) { - B[i] = (right[i-1] - dX[i]*B[i+1] - last_col[i-1]*B[n-1]) / diag[i-1]; + C[i] = (right[i-1] - dX[i]*C[i+1] - last_col[i-1]*C[n-1]) / diag[i-1]; } - B[0] = B[n-1]; + C[0] = C[n-1]; } - Container A(n - 1); - for (SizeT i = 0; i < A.size(); ++i) { - A[i] = (B[i+1] - B[i]) / (3*dX[i]); + for (SizeT i = 0, index = 0; i < n - 1; ++i) { + coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]); + coeffs[index++] = C[i]; + coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3); + coeffs[index++] = Y[i]; } - Container C(n - 1); - for (SizeT i = 0; i < C.size(); ++i) { - C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3); - } - util::spline::merge_coeffs(coeffs, {A, B, C, Y}); }, [&](const typename Types::ParabolicEnds&) { coeffs = compute_coeffs(X, Y, typename Types::FixedThird{0, 0}); }, [&](const typename Types::Clamped& c) { - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); - - Container B(n); - if (n == 2) { - B[0] = (3*dY[0]/dX[0] - c.df_n - 2*c.df_1) / dX[0]; - B[1] = (c.df_1 + 2*c.df_n - 3*dY[0]/dX[0]) / dX[0]; - } else if (n == 3) { - B[1] = (3*(dY[1]/dX[1] - dY[0]/dX[0]) + (c.df_1 - c.df_n)) / (dX[0] + dX[1]); - B[0] = (3*(dY[0]/dX[0] - c.df_1)/dX[0] - B[1]) / 2; - B[2] = (3*(c.df_n - dY[1]/dX[1])/dX[1] - B[1]) / 2; - } else { - Container diag(n - 2), right(n - 2); - diag[0] = 1.5*dX[0] + 2*dX[1]; - right[0] = 3*dY[1]/dX[1] - 4.5*dY[0]/dX[0] + 1.5*c.df_1; - for (SizeT i = 1; i < n - 3; ++i) { - diag[i] = 2 * (dX[i] + dX[i+1]); - right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]); - } - diag[n-3] = 2*dX[n-3] + 1.5*dX[n-2]; - right[n-3] = 4.5*dY[n-2]/dX[n-2] - 3*dY[n-3]/dX[n-3] - 1.5*c.df_n; - - for (SizeT i = 1; i < n - 2; ++i) { - const auto m = dX[i] / diag[i-1]; - diag[i] -= m * dX[i]; - right[i] -= m * right[i-1]; - } - - B[n-2] = right[n-3] / diag[n-3]; - for (SizeT i = n - 3; i >= 1; --i) { - B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1]; - } - B[n-1] = 1.5*(c.df_n/dX[n-2] - dY[n-2]/dX[n-2]/dX[n-2]) - 0.5*B[n-2]; - B[0] = 1.5*(dY[0]/dX[0]/dX[0] - c.df_1/dX[0]) - 0.5*B[1]; - } - - Container A(n - 1); - for (SizeT i = 0; i < A.size(); ++i) { - A[i] = (B[i+1] - B[i]) / (3*dX[i]); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + Container lower(n), diag(n), upper(n), right(n); + for (SizeT i = 1; i < n - 1; ++i) { + lower[i] = dX[i-1]; + diag[i] = 2 * (dX[i-1] + dX[i]); + upper[i] = dX[i]; + right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]); } - Container C(n - 1); - for (SizeT i = 0; i < C.size(); ++i) { - C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3); + lower[0] = NAN; + diag[0] = 2*dX[0]; + upper[0] = dX[0]; + right[0] = 3 * (dY[0]/dX[0] - c.df_1); + lower[n-1] = dX[n-2]; + diag[n-1] = 2*dX[n-2]; + upper[n-1] = NAN; + right[n-1] = 3 * (c.df_n - dY[n-2]/dX[n-2]); + const auto C = util::linalg::tridiag_solve(lower, diag, upper, right); + for (SizeT i = 0, index = 0; i < n - 1; ++i) { + coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]); + coeffs[index++] = C[i]; + coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3); + coeffs[index++] = Y[i]; } - util::spline::merge_coeffs(coeffs, {A, B, C, Y}); }, [&](const typename Types::FixedSecond& s) { - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); - - Container B(n); - B[0] = s.d2f_1 / 2; - B[n-1] = s.d2f_n / 2; - if (n > 2) { - Container diag(n - 2), right(n - 2); - for (SizeT i = 0; i < n - 2; ++i) { - diag[i] = 2 * (dX[i] + dX[i+1]); - right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]); - } - right[0] -= B[0] * dX[0]; - right[n-3] -= B[n-1] * dX[n-2]; - for (SizeT i = 1; i < n - 2; ++i) { - const auto m = dX[i] / diag[i-1]; - diag[i] -= m * dX[i]; - right[i] -= m * right[i-1]; - } - B[n-2] = right[n-3] / diag[n-3]; - for (SizeT i = n - 3; i >= 1; --i) { - B[i] = (right[i-1] - dX[i] * B[i+1]) / diag[i-1]; - } - } - - Container A(n - 1); - for (SizeT i = 0; i < A.size(); ++i) { - A[i] = (B[i+1] - B[i]) / (3*dX[i]); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + Container lower(n), diag(n), upper(n), right(n); + for (SizeT i = 1; i < n - 1; ++i) { + lower[i] = dX[i-1]; + diag[i] = 2 * (dX[i-1] + dX[i]); + upper[i] = dX[i]; + right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]); } - Container C(n - 1); - for (SizeT i = 0; i < C.size(); ++i) { - C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3); + lower[0] = NAN; + diag[0] = 1; + upper[0] = 0; + right[0] = s.d2f_1 / 2; + lower[n-1] = 0; + diag[n-1] = 1; + upper[n-1] = NAN; + right[n-1] = s.d2f_n / 2; + const auto C = util::linalg::tridiag_solve(lower, diag, upper, right); + for (SizeT i = 0, index = 0; i < n - 1; ++i) { + coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]); + coeffs[index++] = C[i]; + coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3); + coeffs[index++] = Y[i]; } - util::spline::merge_coeffs(coeffs, {A, B, C, Y}); }, [&](const typename Types::FixedThird& t) { - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); - - Container B(n); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + Container C(n); if (n == 2) { - B[0] = -dX[0] * (t.d3f_1 + t.d3f_n) / 8; - B[1] = -B[0]; + C[0] = -dX[0] * (t.d3f_1 + t.d3f_n) / 8; + C[1] = -C[0]; } else { - Container diag(n), right(n); - diag[0] = 1; + Container lower(n), diag(n), upper(n), right(n); for (SizeT i = 1; i < n - 1; ++i) { + lower[i] = dX[i-1]; diag[i] = 2 * (dX[i-1] + dX[i]); - } - diag[n-1] = 1; - right[0] = -dX[0] * t.d3f_1 / 2; - for (SizeT i = 1; i < n - 1; ++i) { + upper[i] = dX[i]; right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]); } + lower[0] = NAN; + diag[0] = 1; + upper[0] = -1; + right[0] = -dX[0] * t.d3f_1 / 2; + lower[n-1] = -1; + diag[n-1] = 1; + upper[n-1] = NAN; right[n-1] = dX[n-2] * t.d3f_n / 2; - - const auto m_1 = dX[0] / diag[0]; - diag[1] -= m_1 * (-1); - right[1] -= m_1 * right[0]; - for (SizeT i = 2; i < n - 1; ++i) { - const auto m = dX[i-1] / diag[i-1]; - diag[i] -= m * dX[i-1]; - right[i] -= m * right[i-1]; - } - const auto m_n_1 = (-1) / diag[n-2]; - diag[n-1] -= m_n_1 * dX[n-2]; - right[n-1] -= m_n_1 * right[n-2]; - - B[n-1] = right[n-1] / diag[n-1]; - for (SizeT i = n - 2; i >= 1; --i) { - B[i] = (right[i] - dX[i] * B[i+1]) / diag[i]; - } - B[0] = right[0] + B[1]; + C = util::linalg::tridiag_solve(lower, diag, upper, right); } - - Container A(n - 1); - for (SizeT i = 1; i < A.size() - 1; ++i) { - A[i] = (B[i+1] - B[i]) / (3*dX[i]); + Container D(n - 1); + for (SizeT i = 1; i < D.size() - 1; ++i) { + D[i] = (C[i+1] - C[i]) / (3*dX[i]); } if (n == 2) { - A[0] = (t.d3f_1 + t.d3f_n) / 12; + D[0] = (t.d3f_1 + t.d3f_n) / 12; } else { - A[0] = t.d3f_1 / 6; - A[n-2] = t.d3f_n / 6; + D[0] = t.d3f_1 / 6; + D[n-2] = t.d3f_n / 6; } - Container C(n - 1); - for (SizeT i = 0; i < C.size(); ++i) { - C[i] = dY[i]/dX[i] - dX[i] * ((2*B[i] + B[i+1]) / 3); + for (SizeT i = 0, index = 0; i < n - 1; ++i) { + coeffs[index++] = D[i]; + coeffs[index++] = C[i]; + coeffs[index++] = dY[i]/dX[i] - dX[i] * ((2*C[i] + C[i+1]) / 3); + coeffs[index++] = Y[i]; } - util::spline::merge_coeffs(coeffs, {A, B, C, Y}); }, [&](const typename Types::NotAKnotStart&) { if (n <= 4) { @@ -367,7 +296,7 @@ namespace alfi::spline { return; } - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); { // first four points const auto h1 = dX[0], h2 = dX[1]; @@ -414,7 +343,7 @@ namespace alfi::spline { return; } - const Container dX = util::arrays::diff(X), dY = util::arrays::diff(Y); + const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y); { // last four points const auto h1 = dX[n-4], h2 = dX[n-3]; diff --git a/ALFI/ALFI/util/linalg.h b/ALFI/ALFI/util/linalg.h index 829f6ec..95410c7 100644 --- a/ALFI/ALFI/util/linalg.h +++ b/ALFI/ALFI/util/linalg.h @@ -63,4 +63,110 @@ namespace alfi::util::linalg { return X; } + + /** + @brief Solves a tridiagonal system of linear equations. + + This function solves a system of linear equations \f(AX = B\f), where \f(A\f) is a tridiagonal matrix. + The input vectors are modified in place. + + Unlike \ref tridiag_solve, this algorithm is generally unstable.\n + The sufficient conditions for this algorithm to be stable are + (see this article: https://stsynkov.math.ncsu.edu/book_sample_material/Sections_5.4-5.5.pdf): + - \f(\abs{diag[0]} \ge \abs{upper[0]}\f), + - \f(\abs{diag[k]} \ge \abs{lower[k]} + \abs{upper[k]}, k = 1, 2, ..., n - 2\f), + - \f(\abs{diag[n-1]} \ge \abs{lower[n-1]}\f), + - one of \f(n\f) inequalities is strict, i.e., \f(>\f) rather than \f(\ge\f). + + The first element of `lower` and the last element of `upper` are ignored. + + @param lower the subdiagonal elements of the tridiagonal matrix (first element is ignored) + @param diag the diagonal elements of the tridiagonal matrix + @param upper the superdiagonal elements of the tridiagonal matrix (last element is ignored) + @param right the right-hand side vector of the system + @return the solution vector + */ + template class Container = DefaultContainer> + Container tridiag_solve_unstable(const auto& lower, auto&& diag, const auto& upper, auto&& right) { + const auto n = right.size(); + assert(n == lower.size()); + assert(n == diag.size()); + assert(n == upper.size()); + if (n == 0) { + return {}; + } + for (SizeT i = 0; i < n - 1; ++i) { + const auto m = lower[i+1] / diag[i]; + diag[i+1] -= m * upper[i]; + right[i+1] -= m * right[i]; + } + Container X(n); + X[n-1] = right[n-1] / diag[n-1]; + for (SizeT iter = 2; iter <= n; ++iter) { + const auto i = n - iter; + X[i] = (right[i] - upper[i] * X[i+1]) / diag[i]; + } + return X; + } + + /** + @brief Solves a tridiagonal system of linear equations. + + This function solves a system of linear equations \f(AX = B\f), where \f(A\f) is a tridiagonal matrix. + The input vectors are modified in place. + + Unlike \ref tridiag_solve_unstable, this algorithm is stable. + + The first element of `lower` and the last element of `upper` are ignored. + + @note Based on the https://github.com/snsinfu/cxx-spline/blob/625d4d325cb2/include/spline.hpp#L40-L105 code + from https://github.com/snsinfu, which was licensed under the BSL-1.0 (Boost Software License 1.0). + + @param lower the subdiagonal elements of the tridiagonal matrix (first element is ignored) + @param diag the diagonal elements of the tridiagonal matrix + @param upper the superdiagonal elements of the tridiagonal matrix (last element is ignored) + @param right the right-hand side vector of the system + @return the solution vector + */ + template class Container = DefaultContainer> + Container tridiag_solve(auto&& lower, auto&& diag, auto&& upper, auto&& right) { + const auto n = right.size(); + assert(n == lower.size()); + assert(n == diag.size()); + assert(n == upper.size()); + if (n == 0) { + return {}; + } + if (n == 1) { + return {right[0]/diag[0]}; + } + for (SizeT i = 0; i < n - 1; ++i) { + if (std::abs(diag[i]) >= std::abs(lower[i+1])) { + const auto m = lower[i+1] / diag[i]; + diag[i+1] -= m * upper[i]; + right[i+1] -= m * right[i]; + lower[i+1] = 0; + } else { + // Swap rows i and (i+1). + // Eliminate the lower[i+1] element by reducing row (i+1) by row i. + // Use lower[i+1] as a buffer for the non-tridiagonal element (i,i+2) (hack, used below). + const auto m = diag[i] / lower[i+1]; + diag[i] = lower[i+1]; + lower[i+1] = upper[i+1]; + upper[i+1] *= -m; + std::swap(upper[i], diag[i+1]); + diag[i+1] -= m * upper[i]; + std::swap(right[i], right[i+1]); + right[i+1] -= m * right[i]; + } + } + Container X(n); + X[n-1] = right[n-1] / diag[n-1]; + X[n-2] = (right[n-2] - upper[n-2] * X[n-1]) / diag[n-2]; + for (SizeT iter = 3; iter <= n; ++iter) { + const auto i = n - iter; + X[i] = (right[i] - upper[i] * X[i+1] - lower[i+1] * X[i+2]) / diag[i]; + } + return X; + } } \ No newline at end of file diff --git a/ALFI/ALFI/util/spline.h b/ALFI/ALFI/util/spline.h index 912c6ce..49396ff 100644 --- a/ALFI/ALFI/util/spline.h +++ b/ALFI/ALFI/util/spline.h @@ -100,13 +100,4 @@ namespace alfi::util::spline { << ") is too big (bigger than 4). Returning an empty array..." << std::endl; return {}; } - - template class C1 = DefaultContainer, template class C2 = DefaultContainer> - void merge_coeffs(auto& coeffs, const C1>& separate) { - for (SizeT index = 0, i = 0; index < coeffs.size(); ++i) { - for (SizeT k = 0; k < separate.size(); ++k) { - coeffs[index++] = separate[k][i]; - } - } - } } \ No newline at end of file diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..749c66e --- /dev/null +++ b/LICENSE @@ -0,0 +1,3 @@ +Copyright © 2024 Petr Alexandrovich Sabanov. + +All rights reserved. No rights granted. Unauthorized use prohibited. \ No newline at end of file diff --git a/licenses/BSL-1.0 b/licenses/BSL-1.0 new file mode 100644 index 0000000..36b7cd9 --- /dev/null +++ b/licenses/BSL-1.0 @@ -0,0 +1,23 @@ +Boost Software License - Version 1.0 - August 17th, 2003 + +Permission is hereby granted, free of charge, to any person or organization +obtaining a copy of the software and accompanying documentation covered by +this license (the "Software") to use, reproduce, display, distribute, +execute, and transmit the Software, and to prepare derivative works of the +Software, and to permit third-parties to whom the Software is furnished to +do so, all subject to the following: + +The copyright notices in the Software and this entire statement, including +the above license grant, this restriction and the following disclaimer, +must be included in all copies of the Software, in whole or in part, and +all derivative works of the Software, unless such copies or derivative +works are solely in the form of machine-executable object code generated by +a source language processor. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE, TITLE AND NON-INFRINGEMENT. IN NO EVENT +SHALL THE COPYRIGHT HOLDERS OR ANYONE DISTRIBUTING THE SOFTWARE BE LIABLE +FOR ANY DAMAGES OR OTHER LIABILITY, WHETHER IN CONTRACT, TORT OR OTHERWISE, +ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER +DEALINGS IN THE SOFTWARE.