feat(linalg): make SVD sweep budget configurable

crates/multicalc/src/linear_algebra/svd.rs

@@ -12,6 +12,54 @@ use crate::scalar::Numeric;
 ///
 /// `u` has orthonormal columns, `singular_values` holds the σ in descending order (all ≥ 0), and
 /// `v` has orthonormal columns.
+
+
+/// Settings for the one-sided Jacobi SVD.
+#[derive(Debug, Clone, Copy)]
+pub struct SvdSettings {
+    /// Maximum number of Jacobi sweeps before giving up on convergence.
+    max_sweeps: usize,
+}
+
+impl Default for SvdSettings {
+    /// Default budget: 60 sweeps, matching the previous hardcoded behaviour.
+    fn default() -> Self {
+        Self { max_sweeps: 60 }
+    }
+}
+
+impl SvdSettings {
+    /// Sets the maximum number of Jacobi sweeps.
+    #[must_use]
+    pub const fn with_max_sweeps(mut self, max_sweeps: usize) -> Self {
+        self.max_sweeps = max_sweeps;
+        self
+    }
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 #[derive(Debug, Clone, Copy)]
 #[must_use]
 pub struct Svd<const M: usize, const N: usize, T = f64> {
@@ -68,23 +116,32 @@ impl<const M: usize, const N: usize, T: Numeric> Matrix<M, N, T> {
     ///     }
     /// }
     /// ```
-    pub fn svd(self) -> Result<Svd<M, N, T>, LinalgError> {
-        if M < N {
-            return Err(LinalgError::Underdetermined);
-        }
-        for r in 0..M {
-            for c in 0..N {
-                if !self[(r, c)].is_finite() {
-                    return Err(LinalgError::NonFinite);
-                }
+    /// Decomposes `self` as `U · diag(σ) · Vᵀ` by one-sided Jacobi (thin form, `M ≥ N`), using
+/// the default settings (60 sweeps).
+///
+/// [...existing doc comments yahan rakho jaise the...]
+pub fn svd(self) -> Result<Svd<M, N, T>, LinalgError> {
+    self.svd_with_settings(SvdSettings::default())
+}
+
+/// Same as [`Matrix::svd`] but with configurable settings, such as the sweep budget.
+pub fn svd_with_settings(self, settings: SvdSettings) -> Result<Svd<M, N, T>, LinalgError> {
+    if M < N {
+        return Err(LinalgError::Underdetermined);
+    }
+    for r in 0..M {
+        for c in 0..N {
+            if !self[(r, c)].is_finite() {
+                return Err(LinalgError::NonFinite);
             }
         }
+    }
 
-        let mut u = self;
-        let mut v = Matrix::<N, N, T>::identity();
+    let mut u = self;
+    let mut v = Matrix::<N, N, T>::identity();
 
-        // One-sided Jacobi: rotate column pairs of U until its columns are orthogonal.
-        let max_sweeps = 60;
+    let max_sweeps = settings.max_sweeps;   
+
         for _ in 0..max_sweeps {
             let mut off_max = T::ZERO;
             for p in 0..N {

crates/multicalc/src/vector_field/flux_integral.rs

@@ -72,3 +72,43 @@ pub fn get_2d_custom<T: Numeric>(
         1,
     )?)
 }
+
+pub fn get_3d<T: Numeric>(
+    vector_field: &[&dyn Fn(&[T; 3]) -> T; 3],
+    transformations: &[&dyn Fn(T) -> T; 3],
+    integration_limit: &[T; 2],
+) -> Result<T, IntegrateError> {
+    get_3d_custom(
+        vector_field,
+        transformations,
+        integration_limit,
+        DEFAULT_TOTAL_ITERATIONS,
+    )
+}
+
+pub fn get_3d_custom<T: Numeric>(
+    vector_field: &[&dyn Fn(&[T; 3]) -> T; 3],
+    transformations: &[&dyn Fn(T) -> T; 3],
+    integration_limit: &[T; 2],
+    total_iterations: u64,
+) -> Result<T, IntegrateError> {
+    Ok(line_integral::get_partial_3d(
+        vector_field,
+        transformations,
+        integration_limit,
+        total_iterations,
+        0,
+    )? - line_integral::get_partial_3d(
+        vector_field,
+        transformations,
+        integration_limit,
+        total_iterations,
+        1,
+    )? - line_integral::get_partial_3d(
+        vector_field,
+        transformations,
+        integration_limit,
+        total_iterations,
+        2,
+    )?)
+}

crates/multicalc/tests/suite/vector_field.rs

@@ -127,6 +127,55 @@ fn test_flux_integral_1() {
     assert!(f64::abs(val + 0.0) < 0.01);
 }
 
+#[test]
+fn test_flux_integral_3d_helix() {
+    // Curve r(t) = (cos(t), sin(t), t)
+    // Vector field = (0, 0, z)
+    //
+    // Flux calculation:
+    // partial_x = 0
+    // partial_y = 0
+    // partial_z = ∫ z dz = ∫ t dt = 2*pi^2
+    //
+    // Flux = partial_x - partial_y - partial_z = -2*pi^2
+
+    let vector_field_matrix: [&dyn Fn(&[f64; 3]) -> f64; 3] = [
+        &(|_: &[f64; 3]| -> f64 { 0.0 }),
+        &(|_: &[f64; 3]| -> f64 { 0.0 }),
+        &(|args: &[f64; 3]| -> f64 { args[2] }),
+    ];
+
+    let transformation_matrix: [&dyn Fn(f64) -> f64; 3] = [
+        &(|t: f64| -> f64 { t.cos() }),
+        &(|t: f64| -> f64 { t.sin() }),
+        &(|t: f64| -> f64 { t }),
+    ];
+
+    let two_pi = 2.0 * core::f64::consts::PI;
+    let integration_limit = [0.0, two_pi];
+
+    let val = flux_integral::get_3d_custom(
+        &vector_field_matrix,
+        &transformation_matrix,
+        &integration_limit,
+        100,
+    )
+    .unwrap();
+
+    let expected = -2.0 * core::f64::consts::PI * core::f64::consts::PI;
+
+    assert!(f64::abs(val - expected) < 1e-6);
+}
+
+
+
+
+
+
+
+
+
+
 #[test]
 fn test_curl_2d_1() {
     //vector field is (2*x*y, 3*cos(y)); curl is known to be -2*x, so -2.0 at x = 1