linear_ip/
lib.rs

1extern crate nalgebra as na;
2
3use na::{OMatrix, OVector, Vector, Matrix, Scalar, Dim, default_allocator::DefaultAllocator, allocator::Allocator, convert, RealField, DimMin, DimSub, Const, storage::Owned};
4
5
6/**
7 * A is Dim(MxN), x,c,s are Dim(N), y,b are Dim(M)
8 * Returns (x (primal), y (dual))
9 */
10#[allow(non_snake_case)]
11pub fn solve<T, M, N>(A: &OMatrix<T,M,N>, b: &OVector<T,M>, c: &OVector<T,N>, eps: T, theta: T, gamma: T, max_it: usize) -> (OVector<T, N>,OVector<T,M>) 
12    where 
13        T: Scalar + RealField + Copy, 
14        M: Dim + DimMin<M, Output = M> + DimSub<Const<1>>, 
15        N: Dim + DimMin<N, Output = N>,
16        DefaultAllocator: Allocator<T, M> + Allocator<T, M, N> + Allocator<T, N, M> + Allocator<T, N, N> + Allocator<T, N> + Allocator<T, M, M> + Allocator<(T, usize), M> + Allocator<(usize, usize), M> + Allocator<T, <M as DimSub<Const<1>>>::Output>  {
17    let one: T = convert(1.0);
18    let max_val = T::max_value().unwrap();
19    let svd_eps: T = convert(1e-20);
20    let A_transpose = A.transpose();
21    let n = A.ncols();
22    let m = A.nrows();
23
24    let mut gamma_mu = Vector::<T,N, Owned<T,N,Const<1>>>::from_element_generic(Dim::from_usize(n), Const::<1>,one);
25    let mut x = Vector::<T,N, Owned<T,N,Const<1>>>::from_element_generic(Dim::from_usize(n), Const::<1>,one);
26    let mut s = Vector::<T,N, Owned<T,N,Const<1>>>::from_element_generic(Dim::from_usize(n), Const::<1>,one);
27    let mut y = Vector::<T,M, Owned<T,M,Const<1>>>::zeros_generic(Dim::from_usize(m), Const::<1>);
28    let mut r_primal =  Vector::<T,M, Owned<T,M,Const<1>>>::from_element_generic(Dim::from_usize(m), Const::<1>,max_val);
29    let mut r_dual = Vector::<T,N, Owned<T,N,Const<1>>>::from_element_generic(Dim::from_usize(n), Const::<1>,max_val);
30    let mut k = 0;
31    let mut S_inv = Matrix::<T,N,N,Owned<T,N,N>>::zeros_generic(Dim::from_usize(n), Dim::from_usize(n));
32    let mut X = Matrix::<T,N,N,Owned<T,N,N>>::zeros_generic(Dim::from_usize(n), Dim::from_usize(n));
33    let n_f64 : T = convert(n as f64);
34
35    while k < max_it && (r_primal.norm() > eps || r_dual.norm() > eps || x.dot(&s) > eps) {
36        r_primal = b - A*(&x);
37        r_dual = c - (&A_transpose)*(&y) - (&s);
38        let mu = x.dot(&s)/n_f64;
39        for i in 0..n {
40            S_inv[(i,i)] = one/s[i];
41            X[(i,i)] = x[i];
42            gamma_mu[i] = gamma*mu;
43        }
44        let M = A*(&S_inv)*(&X)*(&A_transpose);
45        let r = b + A*(&S_inv)*((&X)*(&r_dual) - (&gamma_mu));
46        let d_y = M.svd(true,true).solve(&r, svd_eps).expect("direction y failed");
47        let d_s = (&r_dual) - (&A_transpose)*(&d_y);
48        let d_x = -(&x) + (&S_inv)*((&gamma_mu)-(&X)*(&d_s));
49
50        let mut step_primal = max_val;
51        let mut step_dual = max_val;
52
53        for i in 0..n {
54            let dx_i = d_x[i];
55            let ds_i = d_s[i];
56
57            if dx_i < T::zero() {
58                let v = -x[i]/dx_i;
59                step_primal = step_primal.min(v);
60            }
61
62            if ds_i < T::zero() {
63                let v = -s[i]/ds_i;
64                step_dual = step_dual.min(v);
65            }
66        }
67
68        let step_max = step_primal.min(step_dual);
69        let step = one.min(theta*step_max);
70
71        x = x+d_x.scale(step);
72        y = y+d_y.scale(step);
73        s = s+d_s.scale(step);
74        k = k+1;
75    }
76   
77    (x,y)
78}
linear_ip/lib.rs

linear_ip/
lib.rs
