Impl vec diff and scaling

russell_lab/src/internal/array_minus_op.rs

@@ -0,0 +1,148 @@
+use crate::Complex64;
+use crate::StrError;
+
+/// Subtracts two arrays
+///
+/// ```text
+/// w := u - v
+/// ```
+#[inline]
+pub(crate) fn array_minus_op(w: &mut [f64], u: &[f64], v: &[f64]) -> Result<(), StrError> {
+    let n = w.len();
+    if u.len() != n || v.len() != n {
+        return Err("arrays are incompatible");
+    }
+    if n == 0 {
+        return Ok(());
+    }
+    if n == 1 {
+        w[0] = u[0] - v[0];
+    } else if n == 2 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+    } else if n == 3 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+    } else if n == 4 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+    } else if n == 5 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+    } else if n == 6 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+    } else if n == 7 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+        w[6] = u[6] - v[6];
+    } else if n == 8 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+        w[6] = u[6] - v[6];
+        w[7] = u[7] - v[7];
+    } else {
+        let m = n % 4;
+        for i in 0..m {
+            w[i] = u[i] - v[i];
+        }
+        for i in (m..n).step_by(4) {
+            w[i + 0] = u[i + 0] - v[i + 0];
+            w[i + 1] = u[i + 1] - v[i + 1];
+            w[i + 2] = u[i + 2] - v[i + 2];
+            w[i + 3] = u[i + 3] - v[i + 3];
+        }
+    }
+    Ok(())
+}
+
+/// Subtracts two arrays
+///
+/// ```text
+/// w := u - v
+/// ```
+#[inline]
+pub(crate) fn array_minus_op_complex(w: &mut [Complex64], u: &[Complex64], v: &[Complex64]) -> Result<(), StrError> {
+    let n = w.len();
+    if u.len() != n || v.len() != n {
+        return Err("arrays are incompatible");
+    }
+    if n == 0 {
+        return Ok(());
+    }
+    if n == 1 {
+        w[0] = u[0] - v[0];
+    } else if n == 2 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+    } else if n == 3 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+    } else if n == 4 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+    } else if n == 5 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+    } else if n == 6 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+    } else if n == 7 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+        w[6] = u[6] - v[6];
+    } else if n == 8 {
+        w[0] = u[0] - v[0];
+        w[1] = u[1] - v[1];
+        w[2] = u[2] - v[2];
+        w[3] = u[3] - v[3];
+        w[4] = u[4] - v[4];
+        w[5] = u[5] - v[5];
+        w[6] = u[6] - v[6];
+        w[7] = u[7] - v[7];
+    } else {
+        let m = n % 4;
+        for i in 0..m {
+            w[i] = u[i] - v[i];
+        }
+        for i in (m..n).step_by(4) {
+            w[i + 0] = u[i + 0] - v[i + 0];
+            w[i + 1] = u[i + 1] - v[i + 1];
+            w[i + 2] = u[i + 2] - v[i + 2];
+            w[i + 3] = u[i + 3] - v[i + 3];
+        }
+    }
+    Ok(())
+}

russell_lab/src/internal/array_plus_op.rs

@@ -0,0 +1,148 @@
+use crate::Complex64;
+use crate::StrError;
+
+/// Adds two arrays
+///
+/// ```text
+/// w := u + v
+/// ```
+#[inline]
+pub(crate) fn array_plus_op(w: &mut [f64], u: &[f64], v: &[f64]) -> Result<(), StrError> {
+    let n = w.len();
+    if u.len() != n || v.len() != n {
+        return Err("arrays are incompatible");
+    }
+    if n == 0 {
+        return Ok(());
+    }
+    if n == 1 {
+        w[0] = u[0] + v[0];
+    } else if n == 2 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+    } else if n == 3 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+    } else if n == 4 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+    } else if n == 5 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+    } else if n == 6 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+    } else if n == 7 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+        w[6] = u[6] + v[6];
+    } else if n == 8 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+        w[6] = u[6] + v[6];
+        w[7] = u[7] + v[7];
+    } else {
+        let m = n % 4;
+        for i in 0..m {
+            w[i] = u[i] + v[i];
+        }
+        for i in (m..n).step_by(4) {
+            w[i + 0] = u[i + 0] + v[i + 0];
+            w[i + 1] = u[i + 1] + v[i + 1];
+            w[i + 2] = u[i + 2] + v[i + 2];
+            w[i + 3] = u[i + 3] + v[i + 3];
+        }
+    }
+    Ok(())
+}
+
+/// Adds two arrays
+///
+/// ```text
+/// w := u + v
+/// ```
+#[inline]
+pub(crate) fn array_plus_op_complex(w: &mut [Complex64], u: &[Complex64], v: &[Complex64]) -> Result<(), StrError> {
+    let n = w.len();
+    if u.len() != n || v.len() != n {
+        return Err("arrays are incompatible");
+    }
+    if n == 0 {
+        return Ok(());
+    }
+    if n == 1 {
+        w[0] = u[0] + v[0];
+    } else if n == 2 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+    } else if n == 3 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+    } else if n == 4 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+    } else if n == 5 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+    } else if n == 6 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+    } else if n == 7 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+        w[6] = u[6] + v[6];
+    } else if n == 8 {
+        w[0] = u[0] + v[0];
+        w[1] = u[1] + v[1];
+        w[2] = u[2] + v[2];
+        w[3] = u[3] + v[3];
+        w[4] = u[4] + v[4];
+        w[5] = u[5] + v[5];
+        w[6] = u[6] + v[6];
+        w[7] = u[7] + v[7];
+    } else {
+        let m = n % 4;
+        for i in 0..m {
+            w[i] = u[i] + v[i];
+        }
+        for i in (m..n).step_by(4) {
+            w[i + 0] = u[i + 0] + v[i + 0];
+            w[i + 1] = u[i + 1] + v[i + 1];
+            w[i + 2] = u[i + 2] + v[i + 2];
+            w[i + 3] = u[i + 3] + v[i + 3];
+        }
+    }
+    Ok(())
+}

russell_lab/src/internal/add_arrays.rs → russell_lab/src/internal/array_plus_opx.rs

@@ -32,15 +32,13 @@ extern "C" {
     fn cblas_zscal(n: i32, alpha: *const Complex64, x: *mut Complex64, incx: i32);
 }
 
-/// Adds two arrays
-///
-/// **Note:** This is an internal function used by `vec_add` and `mat_add`.
+/// Adds two arrays (with multipliers)
 ///
 /// ```text
 /// w := α⋅u + β⋅v
 /// ```
 #[inline]
-pub(crate) fn add_arrays(w: &mut [f64], alpha: f64, u: &[f64], beta: f64, v: &[f64]) -> Result<(), StrError> {
+pub(crate) fn array_plus_opx(w: &mut [f64], alpha: f64, u: &[f64], beta: f64, v: &[f64]) -> Result<(), StrError> {
     let n = w.len();
     if u.len() != n || v.len() != n {
         return Err("arrays are incompatible");
@@ -119,15 +117,13 @@ pub(crate) fn add_arrays(w: &mut [f64], alpha: f64, u: &[f64], beta: f64, v: &[f
     Ok(())
 }
 
-/// Adds two arrays
-///
-/// **Note:** This is an internal function used by `vec_add` and `mat_add`.
+/// Adds two arrays (with multipliers)
 ///
 /// ```text
 /// w := α⋅u + β⋅v
 /// ```
 #[inline]
-pub fn add_arrays_complex(
+pub(crate) fn array_plus_opx_complex(
     w: &mut [Complex64],
     alpha: Complex64,
     u: &[Complex64],

russell_lab/src/internal/mod.rs

@@ -1,10 +1,14 @@
 //! This module implements some functions for internal use
 
-mod add_arrays;
+mod array_minus_op;
+mod array_plus_op;
+mod array_plus_opx;
 mod constants;
 mod dgeev_data;
 mod to_i32;
-pub(crate) use crate::internal::add_arrays::*;
-pub(crate) use crate::internal::constants::*;
-pub(crate) use crate::internal::dgeev_data::*;
-pub(crate) use crate::internal::to_i32::*;
+pub(crate) use array_minus_op::*;
+pub(crate) use array_plus_op::*;
+pub(crate) use array_plus_opx::*;
+pub(crate) use constants::*;
+pub(crate) use dgeev_data::*;
+pub(crate) use to_i32::*;

russell_lab/src/matrix/complex_mat_add.rs

@@ -1,5 +1,5 @@
 use super::ComplexMatrix;
-use crate::{add_arrays_complex, Complex64, StrError};
+use crate::{array_plus_opx_complex, Complex64, StrError};
 
 /// Performs the addition of two matrices
 ///
@@ -44,7 +44,7 @@ pub fn complex_mat_add(
     if a.nrow() != m || a.ncol() != n || b.nrow() != m || b.ncol() != n {
         return Err("matrices are incompatible");
     }
-    add_arrays_complex(c.as_mut_data(), alpha, a.as_data(), beta, b.as_data())
+    array_plus_opx_complex(c.as_mut_data(), alpha, a.as_data(), beta, b.as_data())
 }
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -116,7 +116,7 @@ mod tests {
     }
 
     #[test]
-    fn complex_add_matrix_oblas_works() {
+    fn complex_add_matrix_5x5_works() {
         const NOISE: Complex64 = cpx!(1234.567, 3456.789);
         let a = ComplexMatrix::from(&[
             [1.0, 2.0, 3.0, 4.0, 5.0],