Hi all. I'm trying to create a FFI to an old Fortran script. It's the one used by scipy.optimize.nnls. I think the official webpage is here. Seeing that there's a C version, I tested it before.
// libc = "0.2"
use libc::{c_double, c_int};
use ndarray::{Array1, Array2, ShapeBuilder};
extern "C" {
fn nnls_(
a: *const c_double,
mda: *const c_int,
m: *const c_int,
n: *const c_int,
b: *const c_double,
x: *mut c_double,
r_norm: *mut c_double,
w: *mut c_double,
zz: *mut c_double,
index_bn: *mut c_int,
mode: *mut c_int,
maxiter: *const c_int,
);
}
pub fn nnls(a: &Array2<f64>, b: &Array1<f64>, max_iter: Option<isize>) -> (Array1<f64>, f64) {
let (m, n) = a.dim();
// This is calling a Fortran library so we must use 'f' order. Moreover, nnls_ modifies it
let mut fortran_a = Array2::zeros(a.dim().f());
fortran_a.assign(&a);
// `b` is copied here because the inner Fortran code modifies it.
let mut b = b.to_owned();
let mut x = Array1::<f64>::zeros(n);
let mut r_norm = 0.0;
let mut w = vec![0.0; n];
let mut zz = vec![0.0; m];
let mut index = vec![0i32; n];
let mut mode = 0i32;
let m = m as i32;
let n = n as i32;
let max_iter = max_iter.unwrap_or(-1) as i32;
unsafe {
nnls_(
fortran_a.as_mut_ptr(),
&m,
&m,
&n,
b.as_mut_ptr(),
x.as_mut_ptr(),
&mut r_norm,
w.as_mut_ptr(),
zz.as_mut_ptr(),
index.as_mut_ptr(),
&mut mode,
&max_iter,
);
}
match mode {
2 => panic!("THE DIMENSIONS OF THE PROBLEM ARE BAD"),
3 => panic!("ITERATION COUNT EXCEEDED. MORE THAN MAXITER ITERATIONS."),
_ => { /* 1 and all other codes are ok */ }
}
(x, r_norm)
}
I'm using pointers even for simple i32 because that's what the C declaration is expecting. It does work. It gives the same results as scipy, it's as fast, etc. However, it crashes when I test it with several threads: Segmentation fault (core dumped)
I tried swapping some 'm' for 'n', it creates different memory errors (I'm pretty sure I got them right anyway) Changing *const
to *mut
changes nothing.
fn main() {
// ndarray = "0.15"
// ndarray-rand = "0.14"
use ndarray::{Array1, Array2};
use ndarray_rand::rand_distr::Uniform;
use ndarray_rand::RandomExt;
let mut a = Array2::random((107, 145), Uniform::new(0.0, 1.0));
a.row_mut(0).fill(1.0);
let b = Array1::random(107, Uniform::new(0.0, 1.0));
rayon::scope(|s| {
s.spawn(|_| {
for _ in 0..5 {
let _ = nnls(&a, &b, None).0;
}
});
s.spawn(|_| {
for _ in 0..5 {
let _ = nnls(&a, &b, None).0;
}
});
});
}
At that point, I'm not sure what to test and I gladly accept your help.