Skip to content

Commit d5c73a9

Browse files
committed
Added non-working tridiag_solve
1 parent a03451a commit d5c73a9

File tree

1 file changed

+49
-0
lines changed

1 file changed

+49
-0
lines changed

ALFI/ALFI/util/linalg.h

Lines changed: 49 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -63,4 +63,53 @@ namespace alfi::util::linalg {
6363

6464
return X;
6565
}
66+
67+
// WRONG CODE
68+
// Tried to implement this for "Not-a-knot" spline, because it had the third coefficient in the first equation,
69+
// but decided that it would be better to derive and equation with only two coefficients.
70+
// See https://github.com/ALFI-lib/alfi-lib.github.io/pull/8
71+
template <typename Number = DefaultNumber, template <typename> class Container = DefaultContainer, typename I1, typename I2, typename I3, typename I4>
72+
Container<Number> tridiag_solve(
73+
Number a11, Number a12, Number a13, Number r1,
74+
Number ann2, Number ann1, Number ann, Number rn,
75+
I1& lower_begin, const I1& lower_end,
76+
I2& diag_begin, const I2& diag_end,
77+
I3& upper_begin, const I3& upper_end,
78+
I4& right_begin, const I4& right_end
79+
) {
80+
const auto n = right_end - right_begin + 2;
81+
assert(n == (lower_end - lower_begin + 2));
82+
assert(n == (diag_end - diag_begin + 2));
83+
assert(n == (upper_end - upper_begin + 2));
84+
auto lower = lower_begin;
85+
auto diag = diag_begin;
86+
auto upper = upper_begin;
87+
auto right = right_begin;
88+
const auto m_1 = *lower / a11;
89+
*diag -= m_1 * a12;
90+
*upper -= m_1 * a13;
91+
*right -= m_1 * r1;
92+
++lower, ++diag, ++upper, ++right;
93+
for (; lower < lower_end; ++lower, ++diag, ++upper, ++right) {
94+
const auto m = *lower / *(diag - 1);
95+
*diag -= m * *(upper - 1);
96+
*right -= m * *(right - 1);
97+
}
98+
if (ann2 != 0) {
99+
const auto m_n = ann2 / *(lower - 1);
100+
ann1 -= m_n * *(diag - 1);
101+
ann -= m_n * *(upper - 1);
102+
rn -= m_n * *(right - 1);
103+
}
104+
const auto m_n = ann1 / *(diag - 1);
105+
ann -= m_n * *(upper - 1);
106+
rn -= m_n * *(right - 1);
107+
Container<Number> X(n);
108+
X[n-1] = rn / ann;
109+
for (auto i = n - 2; i > 0; --i) {
110+
X[i] = (*--right - *--upper * X[i+1]) / *--diag;
111+
}
112+
X[0] = (r1 - a12 * X[1] - a13 * X[2]) / a11;
113+
return X;
114+
}
66115
}

0 commit comments

Comments
 (0)