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

use crate::errors::SGError;
use crate::utilities::float::Float;
use crate::basis::linear::POW2_F64;

use crate::{basis::base::Basis, const_generic::{iterators::grid_iterator::GridIteratorT, storage::SparseGridData}};
#[inline]
#[allow(clippy::too_many_arguments)]
pub(crate) fn eval_boundary<const D: usize, const DIM_OUT: usize, BASIS: Basis, T: Float +std::ops::AddAssign, Iterator: GridIteratorT<D>>(storage: &SparseGridData<D>, basis: &[BASIS; D], x: &[f64; D],
    dim: usize, value: T, iterator: &mut Iterator, alpha: &[[T; DIM_OUT]], result: &mut [T; DIM_OUT]) -> Result<(), SGError>
{
    // Precompute x[dim] as T once per call — eliminates one vcvtsd2ss per node visit.
    let x_t: T = T::from(x[dim]);
    let mut level = 0;
    loop
    {
        let node_index = iterator.index().ok_or(SGError::InvalidIteratorSequence)?;
        let work_index = storage[node_index].index[dim];
        if level > 0
        {
            let new_value = basis[dim].eval_t(level, work_index, x_t);
            if dim == D - 1
            {
                T::accumulate(result, &alpha[node_index], value * new_value);
            }
            else
            {
                eval_boundary(storage, basis, x, dim + 1, value * new_value, iterator, alpha, result)?;
            }
        }
        else
        {
            // handle boundaries if we are on level 0
            // level 0, index 0
            // reset_to_left_level_zero now checks if the node exists - after grid coarsening some boundary nodes are removed.
            if iterator.reset_to_left_level_zero(dim)
            {
                let seq_l = iterator.index().ok_or(SGError::InvalidIteratorSequence)?;
                let new_value_l = basis[dim].eval_t(0, 0, x_t);
                if dim == D - 1
                {
                    T::accumulate(result, &alpha[seq_l], value * new_value_l);
                }
                else
                {
                    eval_boundary(storage, basis, x, dim + 1, value * new_value_l, iterator, alpha, result)?;
                }
            }
            // reset_to_right_level_zero now checks if the node exists - after grid coarsening some boundary nodes are removed.
            if iterator.reset_to_right_level_zero(dim)
            {
                let seq_r = iterator.index().ok_or(SGError::InvalidIteratorSequence)?;
                let new_value_r = basis[dim].eval_t(0, 1, x_t);
                if dim == D - 1
                {
                    T::accumulate(result, &alpha[seq_r], value * new_value_r);
                }
                else
                {
                    eval_boundary(storage, basis, x, dim + 1, value * new_value_r, iterator, alpha, result)?;
                }
            }
        }
        // Can't descend any further.
        if iterator.is_leaf()
        {
            break;
        }
        // this decides in which direction we should descend by evaluating
        // the corresponding bit
        // the bits are coded from left to right starting with level 1
        // being in position max_level
        let x = x[dim];
        if level > 0
        {
            // Use table instead of variable shift + vcvtsi2sd.
            let xh = x * POW2_F64[level as usize];
            let wi = work_index as f64;
            if (xh - wi).abs() < 1e-15
            {
                break;
            }
            if xh > wi && !iterator.right_child(dim)
            {
                break;
            }
            if xh <= wi && !iterator.left_child(dim)
            {                 
                break;
            }
        }
        else 
        {
            if x.abs() < 1e-15 || (x-1.0).abs() < 1e-15
            {
                break;
            }            
            if !iterator.reset_to_level_one(dim)
            {
                break;
            }            
        }
        level += 1;
    }
    iterator.reset_to_left_level_zero(dim); 
    Ok(())   
}