lbfgsbrs 0.1.1

Rust port of L-BFGS-B-C
Documentation
/// Performs the operation y = αx + y where:
/// * n: number of elements to process
/// * da (α): scalar multiplier
/// * dx (x): input vector
/// * incx: stride for vector x
/// * dy (y): input/output vector
/// * incy: stride for vector y
/// 
/// Returns:
///  * 0 on success
///  * -1 if input slices are too small
pub fn daxpy(n: i32, da: f64, dx: &[f64], incx: i32, dy: &mut [f64], incy: i32) -> i32 {
    if n <= 0 {
        return 0;
    }
    if da == 0.0 {
        return 0;
    }

    // Validate slice lengths
    let dx_required = (1 + (n - 1) * incx.abs()) as usize;
    let dy_required = (1 + (n - 1) * incy.abs()) as usize;
    
    if dx.len() < dx_required || dy.len() < dy_required {
        return -1; // Error: slices too small
    }

    if incx == 1 && incy == 1 {
        let m = n % 4;
        if m != 0 {
            for i in 0..m {
                dy[i as usize] += da * dx[i as usize];
            }
            if n < 4 {
                return 0;
            }
        }
        
        let mp1 = m;
        for i in (mp1..n).step_by(4) {
            dy[i as usize] += da * dx[i as usize];
            dy[(i + 1) as usize] += da * dx[(i + 1) as usize];
            dy[(i + 2) as usize] += da * dx[(i + 2) as usize];
            dy[(i + 3) as usize] += da * dx[(i + 3) as usize];
        }
    } else {
        let mut ix = if incx < 0 { (-n + 1) * incx } else { 0 };
        let mut iy = if incy < 0 { (-n + 1) * incy } else { 0 };
        
        for _ in 0..n {
            dy[iy as usize] += da * dx[ix as usize];
            ix += incx;
            iy += incy;
        }
    }
    0
}

/// Copies vector x into vector y where:
/// * n: number of elements to copy
/// * dx (x): source vector
/// * incx: stride for vector x
/// * dy (y): destination vector
/// * incy: stride for vector y
/// 
/// Returns:
///  * 0 on success
///  * -1 if input slices are too small
pub fn dcopy(n: i32, dx: &[f64], incx: i32, dy: &mut [f64], incy: i32) -> i32 {
    if n <= 0 {
        return 0;
    }

    // Validate slice lengths
    let dx_required = (1 + (n - 1) * incx.abs()) as usize;
    let dy_required = (1 + (n - 1) * incy.abs()) as usize;
    
    if dx.len() < dx_required || dy.len() < dy_required {
        return -1; // Error: slices too small
    }

    if incx == 1 && incy == 1 {
        let m = n % 7;
        if m != 0 {
            for i in 0..m {
                dy[i as usize] = dx[i as usize];
            }
            if n < 7 {
                return 0;
            }
        }
        
        let mp1 = m;
        for i in (mp1..n).step_by(7) {
            dy[i as usize] = dx[i as usize];
            dy[(i + 1) as usize] = dx[(i + 1) as usize];
            dy[(i + 2) as usize] = dx[(i + 2) as usize];
            dy[(i + 3) as usize] = dx[(i + 3) as usize];
            dy[(i + 4) as usize] = dx[(i + 4) as usize];
            dy[(i + 5) as usize] = dx[(i + 5) as usize];
            dy[(i + 6) as usize] = dx[(i + 6) as usize];
        }
    } else {
        let mut ix = if incx < 0 { (-n + 1) * incx } else { 0 };
        let mut iy = if incy < 0 { (-n + 1) * incy } else { 0 };
        
        for _ in 0..n {
            dy[iy as usize] = dx[ix as usize];
            ix += incx;
            iy += incy;
        }
    }
    0
}

/// Computes the dot product of two vectors where:
/// * n: number of elements to process
/// * dx (x): first input vector
/// * incx: stride for vector x
/// * dy (y): second input vector
/// * incy: stride for vector y
/// 
/// Returns:
///  * The dot product result
///  * 0.0 if input slices are too small or n <= 0
pub fn ddot(n: i32, dx: &[f64], incx: i32, dy: &[f64], incy: i32) -> f64 {
    if n <= 0 {
        return 0.0;
    }

    // Validate slice lengths
    let dx_required = (1 + (n - 1) * incx.abs()) as usize;
    let dy_required = (1 + (n - 1) * incy.abs()) as usize;
    
    if dx.len() < dx_required || dy.len() < dy_required {
        return 0.0; // Error case
    }

    let mut dtemp = 0.0;
    
    if incx == 1 && incy == 1 {
        let m = n % 5;
        if m != 0 {
            for i in 0..m {
                dtemp += dx[i as usize] * dy[i as usize];
            }
            if n < 5 {
                return dtemp;
            }
        }
        
        let mp1 = m;
        for i in (mp1..n).step_by(5) {
            dtemp += dx[i as usize] * dy[i as usize]
                + dx[(i + 1) as usize] * dy[(i + 1) as usize]
                + dx[(i + 2) as usize] * dy[(i + 2) as usize]
                + dx[(i + 3) as usize] * dy[(i + 3) as usize]
                + dx[(i + 4) as usize] * dy[(i + 4) as usize];
        }
    } else {
        let mut ix = if incx < 0 { (-n + 1) * incx } else { 0 };
        let mut iy = if incy < 0 { (-n + 1) * incy } else { 0 };
        
        for _ in 0..n {
            dtemp += dx[ix as usize] * dy[iy as usize];
            ix += incx;
            iy += incy;
        }
    }
    dtemp
}

/// Scales a vector by a constant where:
/// * n: number of elements to process
/// * da (α): scalar multiplier
/// * dx (x): input/output vector
/// * incx: stride for vector x
/// 
/// Returns:
///  * 0 on success
///  * -1 if input slice is too small
pub fn dscal(n: i32, da: f64, dx: &mut [f64], incx: i32) -> i32 {
    if n <= 0 || incx <= 0 {
        return 0;
    }

    // Validate slice length
    let dx_required = (1 + (n - 1) * incx.abs()) as usize;
    
    if dx.len() < dx_required {
        return -1; // Error: slice too small
    }

    if incx == 1 {
        let m = n % 5;
        if m != 0 {
            for i in 0..m {
                dx[i as usize] = da * dx[i as usize];
            }
            if n < 5 {
                return 0;
            }
        }
        
        let mp1 = m;
        for i in (mp1..n).step_by(5) {
            dx[i as usize] = da * dx[i as usize];
            dx[(i + 1) as usize] = da * dx[(i + 1) as usize];
            dx[(i + 2) as usize] = da * dx[(i + 2) as usize];
            dx[(i + 3) as usize] = da * dx[(i + 3) as usize];
            dx[(i + 4) as usize] = da * dx[(i + 4) as usize];
        }
    } else {
        let mut ix = if incx < 0 { (-n + 1) * incx } else { 0 };
        
        for _ in 0..n {
            dx[ix as usize] = da * dx[ix as usize];
            ix += incx;
        }
    }
    0
}