implement dgetrf and dgetri to compute matrix inverse

README.md

@@ -5,6 +5,8 @@ The package provides an interface to the [Linear Algebra PACKage][1].
 The list of available routines currently includes
 
 * [DGESVD](http://www.netlib.org/lapack/explore-html/d8/d2d/dgesvd_8f.html) and
+* [DGETRF](http://www.netlib.org/lapack/explore-html/d3/d6a/dgetrf_8f.html) and
+* [DGETRI](http://www.netlib.org/lapack/explore-html/df/da4/dgetri_8f.html) and
 * [DSYEV](http://www.netlib.org/lapack/explore-html/dd/d4c/dsyev_8f.html).
 
 The list is extended upon request.

src/metal.rs

@@ -78,3 +78,32 @@ pub fn dgesvd(jobu: Jobu, jobvt: Jobvt, m: usize, n: usize, a: &mut [f64], lda:
  info as *mut _ as *mut _);
  }
 }
+
+#[inline]
+pub fn dgetrf(m: usize, n: usize, a: &mut [f64], lda: usize,
+ ipiv: &mut [i32], info: &mut isize) {
+
+ unsafe {
+ raw::dgetrf_(&(m as c_int) as *const _ as *mut _,
+ &(n as c_int) as *const _ as *mut _,
+ a.as_mut_ptr(),
+ &(lda as c_int) as *const _ as *mut _,
+ ipiv.as_mut_ptr(),
+ info as *mut _ as *mut _);
+ }
+}
+
+#[inline]
+pub fn dgetri(n: usize, a: &mut [f64], lda: usize,
+ ipiv: &mut [i32], work: &mut [f64], lwork: usize, info: &mut isize) {
+
+ unsafe {
+ raw::dgetri_(&(n as c_int) as *const _ as *mut _,
+ a.as_mut_ptr(),
+ &(lda as c_int) as *const _ as *mut _,
+ ipiv.as_mut_ptr(),
+ work.as_mut_ptr(),
+ &(lwork as c_int) as *const _ as *mut _,
+ info as *mut _ as *mut _);
+ }
+}

tests/metal.rs

@@ -78,6 +78,54 @@ fn dgesvd() {
  assert::within(&vt, &expected_vt, 1e-14);
 }
 
+#[test]
+fn dgetrf_and_dgetri() {
+
+ // Compute the inverse of a square matrix.
+
+ use std::iter::repeat;
+
+ let m = 2;
+ let n = 2;
+ let mut a = vec![ // column major order
+ 1.0, 3.0,
+ 2.0, 4.0,
+ ];
+ let lda = m;
+
+ let mut ipiv = vec![0];
+ let mut info = 0;
+
+ metal::dgetrf(m, n, &mut a, lda,
+ &mut ipiv, &mut info);
+
+ if info < 0 {
+ panic!("illegal argument to dgetrf.");
+ }
+ if info > 0 {
+ panic!("singular matrix.");
+ }
+
+ let lwork = n*n;
+ let mut work = repeat(0.0).take(lwork).collect::<Vec<_>>();
+
+ metal::dgetri(n, &mut a, lda,
+ &mut ipiv, &mut work, lwork, &mut info);
+
+ if info < 0 {
+ panic!("illegal argument to dgetri.");
+ }
+ if info > 0 {
+ panic!("singular matrix.");
+ }
+
+ let expected_a = vec![ // column major order
+ -2.0, 1.5,
+ 1.0, -0.5,
+ ];
+ assert::within(&a, &expected_a, 1e-14);
+}
+
 #[test]
 fn dsyev() {
  use std::iter::repeat;