rustyqlib/core/fd_solvers/
psor.rs

1//! Projected SOR (PSOR): iterative solve of the tridiagonal linear
2//! complementarity problem with one- or two-sided obstacles.
3//!
4//! Brennan-Schwartz is exact and O(n) but only handles a one-sided
5//! obstacle reachable by a directional sweep; PSOR is the general tool —
6//! two-sided constraints (e.g. callable/putable structures), and the
7//! smoother of choice inside multi-dimensional splitting schemes.
8
9/// Result of a PSOR solve.
10#[derive(Debug, Clone)]
11pub struct PsorResult {
12    pub x: Vec<f64>,
13    pub iterations: usize,
14    pub converged: bool,
15}
16
17/// Solve `A x = d` projected onto `floor <= x <= cap` by projected SOR,
18/// where `A` is tridiagonal in the same layout as
19/// [`thomas_algorithm`](super::tridiagonal::thomas_algorithm): `a` is the
20/// sub-diagonal (`a[i-1]` multiplies `x[i-1]` in row `i`), `b` the
21/// diagonal and `c` the super-diagonal.
22///
23/// `omega` in `(0, 2)` is the relaxation factor (1 = projected
24/// Gauss-Seidel; ~1.2-1.6 typically accelerates diffusion operators).
25/// Convergence requires the usual SOR conditions (diagonally dominant or
26/// symmetric positive definite `A`), which theta-scheme matrices satisfy.
27/// The iteration stops when the largest update falls below `tol`.
28pub fn psor(
29    a: &[f64],
30    b: &[f64],
31    c: &[f64],
32    d: &[f64],
33    floor: Option<&[f64]>,
34    cap: Option<&[f64]>,
35    omega: f64,
36    tol: f64,
37    max_iter: usize,
38) -> PsorResult {
39    let n = d.len();
40    assert!(b.len() == n && a.len() == n - 1 && c.len() == n - 1);
41    assert!(omega > 0.0 && omega < 2.0, "SOR needs omega in (0, 2)");
42    let project = |i: usize, v: f64| -> f64 {
43        let mut v = v;
44        if let Some(f) = floor {
45            v = v.max(f[i]);
46        }
47        if let Some(cp) = cap {
48            v = v.min(cp[i]);
49        }
50        v
51    };
52
53    // start from the projected diagonal solve
54    let mut x: Vec<f64> = (0..n).map(|i| project(i, d[i] / b[i])).collect();
55    for it in 1..=max_iter {
56        let mut max_update: f64 = 0.0;
57        for i in 0..n {
58            let mut gs = d[i];
59            if i > 0 {
60                gs -= a[i - 1] * x[i - 1];
61            }
62            if i < n - 1 {
63                gs -= c[i] * x[i + 1];
64            }
65            gs /= b[i];
66            let xi = project(i, (1.0 - omega) * x[i] + omega * gs);
67            max_update = max_update.max((xi - x[i]).abs());
68            x[i] = xi;
69        }
70        if max_update <= tol {
71            return PsorResult { x, iterations: it, converged: true };
72        }
73    }
74    PsorResult { x, iterations: max_iter, converged: false }
75}
76
77#[cfg(test)]
78mod tests {
79    use super::super::brennan_schwartz::brennan_schwartz;
80    use super::super::tridiagonal::thomas_algorithm;
81    use super::*;
82
83    const A: [f64; 2] = [1.0, 1.0];
84    const B: [f64; 3] = [3.0, 3.0, 3.0];
85    const C: [f64; 2] = [1.0, 1.0];
86    const D: [f64; 3] = [5.0, 10.0, 11.0];
87
88    #[test]
89    fn unconstrained_psor_matches_thomas() {
90        let free = thomas_algorithm(&A, &B, &C, &D);
91        let r = psor(&A, &B, &C, &D, None, None, 1.2, 1e-13, 10_000);
92        assert!(r.converged);
93        for (x, y) in r.x.iter().zip(&free) {
94            assert!((x - y).abs() < 1e-10, "{:?} vs {free:?}", r.x);
95        }
96    }
97
98    #[test]
99    fn floored_psor_matches_brennan_schwartz_on_a_put_lcp() {
100        // one implicit step of an American-put discretization: obstacle =
101        // convex put payoff, the setting in which Brennan-Schwartz is exact
102        // (Jaillet-Lamberton-Lapeyre), so both solvers must agree. (On
103        // arbitrary obstacles only PSOR solves the true LCP.)
104        let n = 60;
105        let lam = 0.45;
106        let strike = 30.0;
107        let a = vec![-lam; n - 1];
108        let b = vec![1.0 + 2.0 * lam; n];
109        let c = vec![-lam; n - 1];
110        let payoff: Vec<f64> = (0..n).map(|i| (strike - i as f64).max(0.0)).collect();
111
112        let bs = brennan_schwartz(&a, &b, &c, &payoff, &payoff, true);
113        let r = psor(&a, &b, &c, &payoff, Some(&payoff), None, 1.4, 1e-13, 20_000);
114        assert!(r.converged, "psor did not converge");
115        for (x, y) in r.x.iter().zip(&bs) {
116            assert!((x - y).abs() < 1e-7, "psor vs brennan-schwartz mismatch");
117        }
118    }
119
120    #[test]
121    fn two_sided_constraints_are_enforced() {
122        let floor = [1.0, 1.0, 1.0];
123        let cap = [2.0, 2.0, 2.0];
124        let r = psor(&A, &B, &C, &D, Some(&floor), Some(&cap), 1.2, 1e-13, 10_000);
125        assert!(r.converged);
126        assert!(r.x.iter().all(|&v| v >= 1.0 - 1e-12 && v <= 2.0 + 1e-12), "{:?}", r.x);
127        // the cap actually binds somewhere (the free solution exceeds 2)
128        let free = thomas_algorithm(&A, &B, &C, &D);
129        assert!(free.iter().any(|&v| v > 2.0));
130        assert!(r.x.iter().any(|&v| (v - 2.0).abs() < 1e-10));
131    }
132}
rustyqlib/core/fd_solvers/psor.rs

rustyqlib/core/fd_solvers/
psor.rs
