Mercurial > repos > alg_tools / file diff

--- a/src/linsolve.rs	Fri Dec 06 15:30:23 2024 -0500
+++ b/src/linsolve.rs	Fri Dec 06 16:14:41 2024 -0500
@@ -3,6 +3,7 @@
 */
 
 use crate::types::Float;
+#[cfg(feature = "nightly")]
 use std::mem::MaybeUninit;
 
 /// Gaussian elimination for $AX=B$, where $A$ and $B$ are both stored in `ab`,
@@ -45,27 +46,44 @@
 
     // Solve UAX=UB for X where UA with U presenting the transformations above an
     // upper triangular matrix.
+    //
+    // If the "nightly" feature is enabled, we will use an uninitialised array for a
+    // little bit of efficiency.
     // This use of MaybeUninit assumes F : Copy. Otherwise undefined behaviour may occur.
-    let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() );
-    //unsafe { std::mem::MaybeUninit::uninit().assume_init() };
-    for i in (0..M).rev() {
-        for 𝓁 in 0..K {
-            let mut tmp  = ab[i][M+𝓁];
-            for j in (i+1)..M {
-                tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) };
+    #[cfg(feature = "nightly")]
+    {
+        let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() );
+        //unsafe { std::mem::MaybeUninit::uninit().assume_init() };
+        for i in (0..M).rev() {
+            for 𝓁 in 0..K {
+                let mut tmp  = ab[i][M+𝓁];
+                for j in (i+1)..M {
+                    tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) };
+                }
+                tmp /= ab[i][i];
+                x[i][𝓁].write(tmp);
             }
-            tmp /= ab[i][i];
-            x[i][𝓁].write(tmp);
+        }
+        unsafe {
+            //core::intrinsics::assert_inhabited::<[[F; K]; M]>();
+            (&x as *const _ as *const [[F; K]; M]).read()
         }
     }
-    //unsafe { MaybeUninit::array_assume_init(x) };
-    let xinit = unsafe {
-        //core::intrinsics::assert_inhabited::<[[F; K]; M]>();
-        (&x as *const _ as *const [[F; K]; M]).read()
-    };
-
-    //std::mem::forget(x);
-    xinit
+    #[cfg(not(feature = "nightly"))]
+    {
+        let mut x : [[F; K]; M] = [[F::ZERO; K]; M];
+        for i in (0..M).rev() {
+            for 𝓁 in 0..K {
+                let mut tmp  = ab[i][M+𝓁];
+                for j in (i+1)..M {
+                    tmp -= ab[i][j] * x[j][𝓁];
+                }
+                tmp /= ab[i][i];
+                x[i][𝓁] = tmp;
+            }
+        }
+        x
+    }
 }
 
 /// Gaussian elimination for $Ax=b$, where $A$ and $b$ are both stored in `ab`,