comparison: src/linsolve.rs
src/linsolve.rs
- changeset 90
- b3c35d16affe
- parent 55
- 7b2ee3e84c5f
equal
deleted
inserted
replaced
| 25:d14c877e14b7 | 90:b3c35d16affe |
|---|---|
| 1 /*! | 1 /*! |
| 2 Linear equation solvers for small problems stored in Rust arrays. | 2 Linear equation solvers for small problems stored in Rust arrays. |
| 3 */ | 3 */ |
| 4 | 4 |
| 5 use crate::types::Float; | 5 use crate::types::Float; |
| 6 #[cfg(feature = "nightly")] | |
| 6 use std::mem::MaybeUninit; | 7 use std::mem::MaybeUninit; |
| 7 | 8 |
| 8 /// Gaussian elimination for $AX=B$, where $A$ and $B$ are both stored in `ab`, | 9 /// Gaussian elimination for $AX=B$, where $A$ and $B$ are both stored in `ab`, |
| 9 /// $A \in \mathbb{R}^{M \times M}$ and $X, B \in \mathbb{R}^{M \times K}$. | 10 /// $A \in \mathbb{R}^{M \times M}$ and $X, B \in \mathbb{R}^{M \times K}$. |
| 10 pub fn linsolve0<F : Float, const M : usize, const N : usize, const K : usize>( | 11 pub fn linsolve0<F : Float, const M : usize, const N : usize, const K : usize>( |
| 43 } | 44 } |
| 44 } | 45 } |
| 45 | 46 |
| 46 // Solve UAX=UB for X where UA with U presenting the transformations above an | 47 // Solve UAX=UB for X where UA with U presenting the transformations above an |
| 47 // upper triangular matrix. | 48 // upper triangular matrix. |
| 49 // | |
| 50 // If the "nightly" feature is enabled, we will use an uninitialised array for a | |
| 51 // little bit of efficiency. | |
| 48 // This use of MaybeUninit assumes F : Copy. Otherwise undefined behaviour may occur. | 52 // This use of MaybeUninit assumes F : Copy. Otherwise undefined behaviour may occur. |
| 49 let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() ); | 53 #[cfg(feature = "nightly")] |
| 50 //unsafe { std::mem::MaybeUninit::uninit().assume_init() }; | 54 { |
| 51 for i in (0..M).rev() { | 55 let mut x : [[MaybeUninit<F>; K]; M] = core::array::from_fn(|_| MaybeUninit::uninit_array::<K>() ); |
| 52 for 𝓁 in 0..K { | 56 //unsafe { std::mem::MaybeUninit::uninit().assume_init() }; |
| 53 let mut tmp = ab[i][M+𝓁]; | 57 for i in (0..M).rev() { |
| 54 for j in (i+1)..M { | 58 for 𝓁 in 0..K { |
| 55 tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) }; | 59 let mut tmp = ab[i][M+𝓁]; |
| 60 for j in (i+1)..M { | |
| 61 tmp -= ab[i][j] * unsafe { *(x[j][𝓁].assume_init_ref()) }; | |
| 62 } | |
| 63 tmp /= ab[i][i]; | |
| 64 x[i][𝓁].write(tmp); | |
| 56 } | 65 } |
| 57 tmp /= ab[i][i]; | 66 } |
| 58 x[i][𝓁].write(tmp); | 67 unsafe { |
| 68 //core::intrinsics::assert_inhabited::<[[F; K]; M]>(); | |
| 69 (&x as *const _ as *const [[F; K]; M]).read() | |
| 59 } | 70 } |
| 60 } | 71 } |
| 61 //unsafe { MaybeUninit::array_assume_init(x) }; | 72 #[cfg(not(feature = "nightly"))] |
| 62 let xinit = unsafe { | 73 { |
| 63 //core::intrinsics::assert_inhabited::<[[F; K]; M]>(); | 74 let mut x : [[F; K]; M] = [[F::ZERO; K]; M]; |
| 64 (&x as *const _ as *const [[F; K]; M]).read() | 75 for i in (0..M).rev() { |
| 65 }; | 76 for 𝓁 in 0..K { |
| 66 | 77 let mut tmp = ab[i][M+𝓁]; |
| 67 //std::mem::forget(x); | 78 for j in (i+1)..M { |
| 68 xinit | 79 tmp -= ab[i][j] * x[j][𝓁]; |
| 80 } | |
| 81 tmp /= ab[i][i]; | |
| 82 x[i][𝓁] = tmp; | |
| 83 } | |
| 84 } | |
| 85 x | |
| 86 } | |
| 69 } | 87 } |
| 70 | 88 |
| 71 /// Gaussian elimination for $Ax=b$, where $A$ and $b$ are both stored in `ab`, | 89 /// Gaussian elimination for $Ax=b$, where $A$ and $b$ are both stored in `ab`, |
| 72 /// $A \in \mathbb{R}^{M \times M}$ and $x, b \in \mathbb{R}^M$. | 90 /// $A \in \mathbb{R}^{M \times M}$ and $x, b \in \mathbb{R}^M$. |
| 73 #[inline] | 91 #[inline] |
