RustedSciThe 0.4.6

Rust framework for symbolic and numerical computing:BVP ( Newton-Raphson frozen/damped/with collocations ), IVP( BDF, Radau, Backward Euler, LSODE, LSODA, RK45, DoPri), nonlinear equations ( Levenberg, Gavin) and more
Documentation
#[cfg(test)]
mod tests {
    use crate::numerical::BVP_sci::BVP_sci_bordered_banded::{
        BvpSciBandedRoute, BvpSciBandedRoutePolicy, profile_bordered_banded_global_jacobian,
    };
    use crate::numerical::BVP_sci::BVP_sci_faer::{construct_global_jac, faer_col};
    use faer::sparse::{SparseColMat, Triplet};

    fn sparse_from_triplets(
        nrows: usize,
        ncols: usize,
        triplets: Vec<Triplet<usize, usize, f64>>,
    ) -> SparseColMat<usize, f64> {
        SparseColMat::try_new_from_triplets(nrows, ncols, &triplets).unwrap()
    }

    fn dense_local_jacobian(n: usize, diag: f64, off_diag: f64) -> SparseColMat<usize, f64> {
        let mut triplets = Vec::with_capacity(n * n);
        for row in 0..n {
            for col in 0..n {
                let value = if row == col { diag } else { off_diag };
                triplets.push(Triplet::new(row, col, value));
            }
        }
        sparse_from_triplets(n, n, triplets)
    }

    fn endpoint_bc_jacobian(n: usize) -> SparseColMat<usize, f64> {
        let triplets = (0..n)
            .map(|idx| Triplet::new(idx, idx, 1.0))
            .collect::<Vec<_>>();
        sparse_from_triplets(n, n, triplets)
    }

    fn empty_bc_jacobian(n: usize) -> SparseColMat<usize, f64> {
        sparse_from_triplets(n, n, Vec::new())
    }

    fn global_jacobian_for_dense_block_bvp(
        n: usize,
        m: usize,
        left_bc: bool,
        right_bc: bool,
    ) -> SparseColMat<usize, f64> {
        let k = 0;
        let h = faer_col::from_fn(m - 1, |_| 1.0 / (m - 1) as f64);
        let df_dy = (0..m)
            .map(|_| dense_local_jacobian(n, 2.0, 0.05))
            .collect::<Vec<_>>();
        let df_dy_middle = (0..m - 1)
            .map(|_| dense_local_jacobian(n, 2.1, 0.04))
            .collect::<Vec<_>>();
        let dbc_dya = if left_bc {
            endpoint_bc_jacobian(n)
        } else {
            empty_bc_jacobian(n)
        };
        let dbc_dyb = if right_bc {
            endpoint_bc_jacobian(n)
        } else {
            empty_bc_jacobian(n)
        };

        construct_global_jac(
            n,
            m,
            k,
            &h,
            &df_dy,
            &df_dy_middle,
            None,
            None,
            &dbc_dya,
            &dbc_dyb,
            None,
        )
    }

    #[test]
    fn compact_full_scalar_matrix_recommends_full_banded() {
        let mat = sparse_from_triplets(
            4,
            4,
            vec![
                Triplet::new(0, 0, 4.0),
                Triplet::new(1, 0, -1.0),
                Triplet::new(0, 1, -1.0),
                Triplet::new(1, 1, 4.0),
                Triplet::new(2, 1, -1.0),
                Triplet::new(1, 2, -1.0),
                Triplet::new(2, 2, 4.0),
                Triplet::new(3, 2, -1.0),
                Triplet::new(2, 3, -1.0),
                Triplet::new(3, 3, 4.0),
            ],
        );

        let profile = profile_bordered_banded_global_jacobian(
            &mat,
            1,
            4,
            0,
            BvpSciBandedRoutePolicy::default(),
        );
        assert_eq!(
            profile.recommended_route,
            BvpSciBandedRoute::FullScalarBanded
        );
        assert!(profile.full_amplification.unwrap() <= 8.0);
    }

    #[test]
    fn endpoint_bc_rows_recommend_bordered_banded() {
        let n = 6;
        let m = 200;
        let mat = global_jacobian_for_dense_block_bvp(n, m, true, true);

        let profile = profile_bordered_banded_global_jacobian(
            &mat,
            n,
            m,
            0,
            BvpSciBandedRoutePolicy::default(),
        );

        assert_eq!(profile.recommended_route, BvpSciBandedRoute::BorderedBanded);
        assert_eq!(profile.full_scalar.kl, (m - 1) * n);
        assert!(profile.collocation_scalar.kl <= n);
        assert!(profile.full_amplification.unwrap() > 50.0);
        assert!(profile.collocation_amplification.unwrap() < 2.0);
    }

    #[test]
    fn shape_mismatch_recommends_sparse_fallback() {
        let mat = sparse_from_triplets(2, 3, vec![Triplet::new(0, 0, 1.0)]);

        let profile = profile_bordered_banded_global_jacobian(
            &mat,
            1,
            2,
            0,
            BvpSciBandedRoutePolicy::default(),
        );

        assert_eq!(profile.recommended_route, BvpSciBandedRoute::SparseFallback);
        assert_eq!(
            profile.reason,
            "matrix shape does not match BVP_sci global layout"
        );
    }

    #[test]
    fn wide_sparse_collocation_recommends_sparse_fallback() {
        let n = 8;
        let m = 20;
        let size = n * m;
        let collocation_rows = n * (m - 1);
        let mut triplets = Vec::with_capacity(collocation_rows * 2);
        for row in 0..collocation_rows {
            triplets.push(Triplet::new(row, row, 2.0));
            triplets.push(Triplet::new(row, size - 1, 0.01));
        }
        let mat = sparse_from_triplets(size, size, triplets);

        let profile = profile_bordered_banded_global_jacobian(
            &mat,
            n,
            m,
            0,
            BvpSciBandedRoutePolicy::default(),
        );

        assert_eq!(profile.recommended_route, BvpSciBandedRoute::SparseFallback);
        assert!(profile.full_amplification.unwrap() > 8.0);
        assert!(profile.collocation_amplification.unwrap() > 4.0);
    }
}