## Coverage Report

#### Created: 2026-01-25 15:05

next uncovered line (L), next uncovered region (R), next uncovered branch (B)

    /home/noah/src/trueno/src/matrix/ops/linear.rs

    Line
    Count
    Source
    1
    //! Linear algebra operations for Matrix
    2
    //!
    3
    //! This module provides linear operations:
    4
    //! - `transpose()` - Matrix transpose
    5
    //! - `matvec()` - Matrix-vector multiplication
    6
    //! - `vecmat()` - Vector-matrix multiplication
    7
    8
    use crate::{Backend, TruenoError, Vector};
    9
    10
    #[cfg(feature = "tracing")]
    11
    use tracing::instrument;
    12
    13
    /// Backend dispatch macro for dot product - centralizes platform-specific SIMD dispatch
    14
    macro_rules! dispatch_dot {
    15
        ($backend:expr, $a:expr, $b:expr) => {{
    16
            use crate::backends::{scalar::ScalarBackend, VectorBackend};
    17
            #[cfg(target_arch = "x86_64")]
    18
            use crate::backends::{avx2::Avx2Backend, sse2::Sse2Backend};
    19
            // SAFETY: CPU features verified at runtime before backend selection
    20
            unsafe {
    21
                match $backend {
    22
                    Backend::Scalar => ScalarBackend::dot($a, $b),
    23
                    #[cfg(target_arch = "x86_64")]
    24
                    Backend::SSE2 | Backend::AVX => Sse2Backend::dot($a, $b),
    25
                    #[cfg(target_arch = "x86_64")]
    26
                    Backend::AVX2 | Backend::AVX512 => Avx2Backend::dot($a, $b),
    27
                    #[cfg(not(target_arch = "x86_64"))]
    28
                    Backend::SSE2 | Backend::AVX | Backend::AVX2 | Backend::AVX512 => {
    29
                        ScalarBackend::dot($a, $b)
    30
                    }
    31
                    #[cfg(any(target_arch = "aarch64", target_arch = "arm"))]
    32
                    Backend::NEON => {
    33
                        use crate::backends::neon::NeonBackend;
    34
                        NeonBackend::dot($a, $b)
    35
                    }
    36
                    #[cfg(not(any(target_arch = "aarch64", target_arch = "arm")))]
    37
                    Backend::NEON => ScalarBackend::dot($a, $b),
    38
                    #[cfg(target_arch = "wasm32")]
    39
                    Backend::WasmSIMD => {
    40
                        use crate::backends::wasm::WasmBackend;
    41
                        WasmBackend::dot($a, $b)
    42
                    }
    43
                    #[cfg(not(target_arch = "wasm32"))]
    44
                    Backend::WasmSIMD => ScalarBackend::dot($a, $b),
    45
                    Backend::GPU | Backend::Auto => ScalarBackend::dot($a, $b),
    46
                }
    47
            }
    48
        }};
    49
    }
    50
    51
    use super::super::Matrix;
    52
    53
    impl Matrix<f32> {
    54
        /// Transpose this matrix (swap rows and columns)
    55
        ///
    56
        /// Returns a new matrix with dimensions swapped: `self.rows → result.cols`,
    57
        /// `self.cols → result.rows`.
    58
        ///
    59
        /// # Performance
    60
        ///
    61
        /// Uses cache-optimized block-wise transpose with 32x32 blocks.
    62
        /// Sequential writes for output ensure good cache behavior.
    63
        ///
    64
        /// # Example
    65
        ///
    66
        /// ```
    67
        /// use trueno::Matrix;
    68
        ///
    69
        /// let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    70
        /// let t = m.transpose();
    71
        ///
    72
        /// // [[1, 2, 3],     [[1, 4],
    73
        /// //  [4, 5, 6]]  →   [2, 5],
    74
        /// //                  [3, 6]]
    75
        /// assert_eq!(t.rows(), 3);
    76
        /// assert_eq!(t.cols(), 2);
    77
        /// assert_eq!(t.get(0, 0), Some(&1.0));
    78
        /// assert_eq!(t.get(0, 1), Some(&4.0));
    79
        /// assert_eq!(t.get(1, 0), Some(&2.0));
    80
        /// ```
    81
        #[cfg_attr(feature = "tracing", instrument(skip(self), fields(dims = %format!("{}x{}", self.rows, self.cols))))]
    82
    0
        pub fn transpose(&self) -> Matrix<f32> {
    83
    0
            let mut result = Matrix::zeros_with_backend(self.cols, self.rows, self.backend);
    84
    85
            // Use block-wise transpose for better cache locality
    86
            // Block size of 32 balances cache efficiency for both square and non-square matrices
    87
            const BLOCK_SIZE: usize = 32;
    88
    89
            // For non-square matrices, process output rows sequentially for write coalescing
    90
            // This ensures writes are sequential in memory regardless of input shape
    91
            // Fix for issue #65: non-square transpose was slow due to strided writes
    92
    93
            // Process in blocks, iterating output rows first for sequential writes
    94
    0
            for j_block in (0..self.cols).step_by(BLOCK_SIZE) {
    95
    0
                let j_end = (j_block + BLOCK_SIZE).min(self.cols);
    96
    97
    0
                for i_block in (0..self.rows).step_by(BLOCK_SIZE) {
    98
    0
                    let i_end = (i_block + BLOCK_SIZE).min(self.rows);
    99
    100
                    // Within block: iterate output rows (j) in outer loop for sequential writes
    101
    0
                    for j in j_block..j_end {
    102
    0
                        let dst_row_start = j * result.cols;
    103
    0
                        for i in i_block..i_end {
    104
    0
                            // result[j, i] = self[i, j]
    105
    0
                            // Sequential write: dst_row_start + i increments by 1
    106
    0
                            // Strided read: acceptable, CPU prefetch handles this
    107
    0
                            result.data[dst_row_start + i] = self.data[i * self.cols + j];
    108
    0
                        }
    109
                    }
    110
                }
    111
            }
    112
    113
    0
            result
    114
    0
        }
    115
    116
        /// Matrix-vector multiplication (column vector): A × v
    117
        ///
    118
        /// Multiplies this matrix by a column vector, computing `A × v` where the result
    119
        /// is a column vector with length equal to the number of rows in `A`.
    120
        ///
    121
        /// # Mathematical Definition
    122
        ///
    123
        /// For an m×n matrix A and an n-dimensional vector v:
    124
        /// ```text
    125
        /// result[i] = Σ(j=0 to n-1) A[i,j] × v[j]
    126
        /// ```
    127
        ///
    128
        /// # Arguments
    129
        ///
    130
        /// * `v` - Column vector with length equal to `self.cols()`
    131
        ///
    132
        /// # Returns
    133
        ///
    134
        /// A new vector with length `self.rows()`
    135
        ///
    136
        /// # Errors
    137
        ///
    138
        /// Returns `InvalidInput` if `v.len() != self.cols()`
    139
        ///
    140
        /// # Example
    141
        ///
    142
        /// ```
    143
        /// use trueno::{Matrix, Vector};
    144
        ///
    145
        /// let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    146
        /// let v = Vector::from_slice(&[1.0, 2.0, 3.0]);
    147
        /// let result = m.matvec(&v).unwrap();
    148
        ///
    149
        /// // [[1, 2, 3]   [1]   [1×1 + 2×2 + 3×3]   [14]
    150
        /// //  [4, 5, 6]] × [2] = [4×1 + 5×2 + 6×3] = [32]
    151
        /// //               [3]
    152
        /// assert_eq!(result.as_slice(), &[14.0, 32.0]);
    153
        /// ```
    154
    12
        pub fn matvec(&self, v: &Vector<f32>) -> Result<Vector<f32>, TruenoError> {
    155
    12
            if v.len() != self.cols {
    156
    0
                return Err(TruenoError::InvalidInput(format!(
    157
    0
                    "Vector length {} does not match matrix columns {} for matrix-vector multiplication",
    158
    0
                    v.len(),
    159
    0
                    self.cols
    160
    0
                )));
    161
    12
            }
    162
    163
    12
            let v_slice = v.as_slice();
    164
    165
    12
            let mut result_data = vec![0.0; self.rows];
    166
    167
            // Parallel execution for very large matrices (≥4096 rows)
    168
            // Note: Thread overhead dominates for smaller matrices
    169
            #[cfg(feature = "parallel")]
    170
            {
    171
                const PARALLEL_THRESHOLD: usize = 4096;
    172
    173
                if self.rows >= PARALLEL_THRESHOLD {
    174
                    use rayon::prelude::*;
    175
                    use std::sync::atomic::{AtomicPtr, Ordering};
    176
                    use std::sync::Arc;
    177
    178
                    let result_ptr = Arc::new(AtomicPtr::new(result_data.as_mut_ptr()));
    179
    180
                    // Process rows in parallel - each row computes an independent dot product
    181
                    (0..self.rows).into_par_iter().for_each(|i| {
    182
                        let row_start = i * self.cols;
    183
                        let row = &self.data[row_start..(row_start + self.cols)];
    184
    185
                        let dot_result = dispatch_dot!(self.backend, row, v_slice);
    186
    187
                        // Write to non-overlapping memory location (thread-safe)
    188
                        // SAFETY: CPU feature verified at runtime, slices bounds-checked
    189
                        unsafe {
    190
                            let ptr = result_ptr.load(Ordering::Relaxed);
    191
                            *ptr.add(i) = dot_result;
    192
                        }
    193
                    });
    194
    195
                    return Ok(Vector::from_slice(&result_data));
    196
                }
    197
            }
    198
    199
            // SIMD-optimized execution: each row-vector product is a dot product
    200
    2.82k
            for (i, result) in result_data.iter_mut()12.enumerate12() {
    201
    2.82k
                let row_start = i * self.cols;
    202
    2.82k
                let row = &self.data[row_start..(row_start + self.cols)];
    203
    204
                // Use SIMD dot product for each row
    205
    2.82k
                *result = dispatch_dot!0(self.backend, row0, v_slice0);
    206
            }
    207
    208
    12
            Ok(Vector::from_slice(&result_data))
    209
    12
        }
    210
    211
        /// Vector-matrix multiplication (row vector): v^T × A
    212
        ///
    213
        /// Multiplies a row vector by this matrix, computing `v^T × A` where the result
    214
        /// is a row vector with length equal to the number of columns in `A`.
    215
        ///
    216
        /// # Mathematical Definition
    217
        ///
    218
        /// For an m-dimensional vector v and an m×n matrix A:
    219
        /// ```text
    220
        /// result[j] = Σ(i=0 to m-1) v[i] × A[i,j]
    221
        /// ```
    222
        ///
    223
        /// # Arguments
    224
        ///
    225
        /// * `v` - Row vector with length equal to `m.rows()`
    226
        /// * `m` - Matrix to multiply
    227
        ///
    228
        /// # Returns
    229
        ///
    230
        /// A new vector with length `m.cols()`
    231
        ///
    232
        /// # Errors
    233
        ///
    234
        /// Returns `InvalidInput` if `v.len() != m.rows()`
    235
        ///
    236
        /// # Example
    237
        ///
    238
        /// ```
    239
        /// use trueno::{Matrix, Vector};
    240
        ///
    241
        /// let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    242
        /// let v = Vector::from_slice(&[1.0, 2.0]);
    243
        /// let result = Matrix::vecmat(&v, &m).unwrap();
    244
        ///
    245
        /// // [1, 2] × [[1, 2, 3]  = [1×1 + 2×4, 1×2 + 2×5, 1×3 + 2×6]
    246
        /// //           [4, 5, 6]]
    247
        /// //         = [9, 12, 15]
    248
        /// assert_eq!(result.as_slice(), &[9.0, 12.0, 15.0]);
    249
        /// ```
    250
    0
        pub fn vecmat(v: &Vector<f32>, m: &Matrix<f32>) -> Result<Vector<f32>, TruenoError> {
    251
    0
            if v.len() != m.rows {
    252
    0
                return Err(TruenoError::InvalidInput(format!(
    253
    0
                    "Vector length {} does not match matrix rows {} for vector-matrix multiplication",
    254
    0
                    v.len(),
    255
    0
                    m.rows
    256
    0
                )));
    257
    0
            }
    258
    259
            // SIMD-optimized implementation using row-wise accumulation
    260
            // Instead of column-wise access (cache-unfriendly), we compute:
    261
            // result = Σ(i) v[i] * row_i (cache-friendly, vectorizable)
    262
            //
    263
            // This approach:
    264
            // 1. Sequential row access (cache-friendly vs strided column access)
    265
            // 2. Uses SIMD scale and add operations
    266
            // 3. Leverages existing optimized Vector operations
    267
    268
    0
            let mut result = Vector::from_slice(&vec![0.0; m.cols]);
    269
    0
            let v_slice = v.as_slice();
    270
    271
            // Accumulate each scaled row into result
    272
    0
            for (i, &scalar) in v_slice.iter().enumerate().take(m.rows) {
    273
    0
                let row_start = i * m.cols;
    274
    0
                let row = &m.data[row_start..(row_start + m.cols)];
    275
    276
                // Create vector for this row
    277
    0
                let row_vec = Vector::from_slice(row);
    278
    279
                // result += scalar * row (using SIMD scale and add)
    280
    0
                let scaled_row = row_vec.scale(scalar)?;
    281
    0
                result = result.add(&scaled_row)?;
    282
            }
    283
    284
    0
            Ok(result)
    285
    0
        }
    286
    }
    287
    288
    #[cfg(test)]
    289
    mod tests {
    290
        use super::*;
    291
    292
        #[test]
    293
        fn test_transpose_square() {
    294
            let m = Matrix::from_vec(2, 2, vec![1.0, 2.0, 3.0, 4.0]).unwrap();
    295
            let t = m.transpose();
    296
            assert_eq!(t.rows(), 2);
    297
            assert_eq!(t.cols(), 2);
    298
            assert_eq!(t.get(0, 0), Some(&1.0));
    299
            assert_eq!(t.get(0, 1), Some(&3.0));
    300
            assert_eq!(t.get(1, 0), Some(&2.0));
    301
            assert_eq!(t.get(1, 1), Some(&4.0));
    302
        }
    303
    304
        #[test]
    305
        fn test_transpose_rect() {
    306
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    307
            let t = m.transpose();
    308
            assert_eq!(t.rows(), 3);
    309
            assert_eq!(t.cols(), 2);
    310
            assert_eq!(t.get(0, 0), Some(&1.0));
    311
            assert_eq!(t.get(0, 1), Some(&4.0));
    312
            assert_eq!(t.get(1, 0), Some(&2.0));
    313
        }
    314
    315
        #[test]
    316
        fn test_matvec_basic() {
    317
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    318
            let v = Vector::from_slice(&[1.0, 2.0, 3.0]);
    319
            let result = m.matvec(&v).unwrap();
    320
            assert_eq!(result.as_slice(), &[14.0, 32.0]);
    321
        }
    322
    323
        #[test]
    324
        fn test_matvec_dimension_mismatch() {
    325
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    326
            let v = Vector::from_slice(&[1.0, 2.0]); // Wrong size
    327
            assert!(m.matvec(&v).is_err());
    328
        }
    329
    330
        #[test]
    331
        fn test_vecmat_basic() {
    332
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    333
            let v = Vector::from_slice(&[1.0, 2.0]);
    334
            let result = Matrix::vecmat(&v, &m).unwrap();
    335
            assert_eq!(result.as_slice(), &[9.0, 12.0, 15.0]);
    336
        }
    337
    338
        #[test]
    339
        fn test_vecmat_dimension_mismatch() {
    340
            let m = Matrix::from_vec(2, 3, vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0]).unwrap();
    341
            let v = Vector::from_slice(&[1.0, 2.0, 3.0]); // Wrong size
    342
            assert!(Matrix::vecmat(&v, &m).is_err());
    343
        }
    344
    }