sgrust 0.8.6

A sparse grid library written in Rust.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113

use std::f64::consts::PI;

use rustfft::{FftPlanner, num_complex::Complex};
///
/// Compute Clenshaw-Curtis weights using approach by J. Waldvogel (2013)
/// "Fast construction of the Fejer and Clenshaw-Curtis quadrature rules". 
/// 
pub(crate) fn cc_weights(level: u32) -> Vec<f64> {
    let n = (1_u32 << level) as usize;
    if level == 0
    {
        return vec![1.0];
    }
    let mut n_vals = Vec::new();
    for i in 0..n
    {
        let val = 2 * i + 1;
        if val >= n
        {
            break;
        }
        n_vals.push(val as f64);
    }
    let l = n_vals.len();
    let  m = n - l;    
    let mut v0: Vec<f64> = n_vals.iter()
        .map(|&i| 2.0 / (i * (i - 2.0)))
        .collect();
    v0.push(1.0 / n_vals[n_vals.len()-1]);
    v0.extend(vec![0.0; m]);

    let mut v2 = vec![0.0; v0.len()-1];
    let end = v0.len();
    for i in 0..v2.len()
    {
        v2[i] = -v0[i] - v0[end - i - 1];
    }

    let mut planner = FftPlanner::new();
    let mut weights = vec![Complex::new(-1.0, 0.0); n];
    weights[l] += n as f64;
    weights[m] += n as f64;

    let g_scale = (n * n - 1 + n % 2) as f64;
    // Perform inverse FFT using rustfft
    for (g,w) in weights.iter_mut().zip(v2)
    {
        g.re = g.re / g_scale + w;
    }    
    let fft =   planner.plan_fft_inverse(n);
    fft.process(&mut weights);
    for w in weights.iter_mut()
    {
        w.re *= 0.5 / n as f64;
    }    
    weights.push(weights[0]);
    weights.iter().map(|x|x.re).collect()
}

/// Return Clenshaw-Curtis nodes over a (0,1) interval
pub(crate) fn cc_nodes(level: u32) -> Vec<f64>
{
    let n =  if level == 0 {1} else { 2_u32.pow(level) + 1 } ;
    let mut nodes = vec![0.0; n as usize];    
    if level > 0
    {
        #[allow(clippy::needless_range_loop)]
        for i in 0..n as usize
        {
            nodes[i] = 0.5*(1.0 + f64::cos(PI*(n as usize - 1 - i) as f64 / (n - 1) as f64));
        }
    }
    else 
    {
        nodes[0] = 0.5;
    }
    nodes
}
pub(crate) fn cc_num_nodes(level: u32) -> u32
{
    if level == 0 {1} else { 2_u32.pow(level) + 1 }
}

#[test]
fn check_weights() 
{
    for i in 0..2
    {
        let _ = cc_weights(i);
    }
    let level = 2;
    let cc_weights = cc_weights(level);
    
    println!("Number of Curtis Weights: {}", cc_weights.len());
    println!("Clenshaw-Curtis Weights: {:?}", cc_weights);
    assert!((1.0-cc_weights.iter().sum::<f64>()/1.0).abs() < 1e-15);

   

    // Weights computed from CLENSHAW_CURTIS_RULE by J. Burkardt. These are computed over a (-1,+1) interval, so we multiply
    // these values by 0.5 to correct weighting for our (0,1) interval.
    let weights5 = [0.06666666666666668, 0.5333333333333333, 0.7999999999999999, 0.5333333333333334, 0.06666666666666668];
    let nodes  = cc_nodes(level);
    for i in 0..5
    {
        println!("{},{}",  nodes[i], cc_weights[i]);
        assert!((1.0-cc_weights[i]/(0.5*weights5[i])).abs() < 1e-15);
    }
    
    
}