diff --git a/Cargo.toml b/Cargo.toml
index 47ecfe65..5d9e44bb 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -32,11 +32,11 @@ features = ["blas"]
 default-features = false
 
 [dependencies.blas-src]
-version = "0.2"
+version = "0.3"
 default-features = false
 optional = true
 
 [dependencies.lapack-src]
-version = "0.2"
+version = "0.3"
 default-features = false
 optional = true
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
index a6b1f035..b5f1de0c 100644
--- a/azure-pipelines.yml
+++ b/azure-pipelines.yml
@@ -58,3 +58,16 @@ jobs:
           cargo test -v --features=intel-mkl --no-default-features
           cargo test -v --features=intel-mkl,serde-1 --no-default-features
         displayName: run test
+
+  - job: WindowsIntelMKL
+    pool:
+      vmImage: 'windows-2019'
+    steps:
+      - script: |
+          curl -sSf -o rustup-init.exe https://win.rustup.rs
+          rustup-init.exe -y 2>&1
+          set PATH=%PATH%;%USERPROFILE%\.cargo\bin
+          echo '##vso[task.setvariable variable=PATH;]%PATH%;%USERPROFILE%\.cargo\bin'
+        displayName: install rustup on Windows
+      - script: cargo test -v --features=intel-mkl --no-default-features 2>&1
+        displayName: run test
diff --git a/src/eigh.rs b/src/eigh.rs
index a3ce912f..90115d14 100644
--- a/src/eigh.rs
+++ b/src/eigh.rs
@@ -63,6 +63,12 @@ where
     type EigVal = Array1<A::Real>;
 
     fn eigh_inplace(&mut self, uplo: UPLO) -> Result<(Self::EigVal, &mut Self)> {
+        let layout = self.square_layout()?;
+        // XXX Force layout to be Fortran (see #146)
+        match layout {
+            MatrixLayout::C(_) => self.swap_axes(0, 1),
+            MatrixLayout::F(_) => {}
+        }
         let s = unsafe { A::eigh(true, self.square_layout()?, uplo, self.as_allocated_mut()?)? };
         Ok((ArrayBase::from_vec(s), self))
     }
diff --git a/src/lib.rs b/src/lib.rs
index 90f7c003..834c6f26 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -36,6 +36,9 @@
 //!  - [Random matrix generators](generate/index.html)
 //!  - [Scalar trait](types/trait.Scalar.html)
 
+extern crate blas_src;
+extern crate lapack_src;
+
 pub mod assert;
 pub mod cholesky;
 pub mod convert;
diff --git a/tests/eigh.rs b/tests/eigh.rs
index a8b088d5..77be699b 100644
--- a/tests/eigh.rs
+++ b/tests/eigh.rs
@@ -2,10 +2,15 @@ use ndarray::*;
 use ndarray_linalg::*;
 
 #[test]
-fn eigen_vector_manual() {
+fn fixed() {
     let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
     let (e, vecs): (Array1<_>, Array2<_>) = (&a).eigh(UPLO::Upper).unwrap();
     assert_close_l2!(&e, &arr1(&[2.0, 2.0, 5.0]), 1.0e-7);
+
+    // Check eigenvectors are orthogonalized
+    let s = vecs.t().dot(&vecs);
+    assert_close_l2!(&s, &Array::eye(3), 1.0e-7);
+
     for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
         let av = a.dot(&v);
         let ev = v.mapv(|x| e[i] * x);
@@ -14,12 +19,53 @@ fn eigen_vector_manual() {
 }
 
 #[test]
-fn diagonalize() {
-    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
+fn fixed_t() {
+    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]).reversed_axes();
     let (e, vecs): (Array1<_>, Array2<_>) = (&a).eigh(UPLO::Upper).unwrap();
-    let s = vecs.t().dot(&a).dot(&vecs);
-    for i in 0..3 {
-        assert_rclose!(e[i], s[(i, i)], 1e-7);
+    assert_close_l2!(&e, &arr1(&[2.0, 2.0, 5.0]), 1.0e-7);
+
+    // Check eigenvectors are orthogonalized
+    let s = vecs.t().dot(&vecs);
+    assert_close_l2!(&s, &Array::eye(3), 1.0e-7);
+
+    for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
+        let av = a.dot(&v);
+        let ev = v.mapv(|x| e[i] * x);
+        assert_close_l2!(&av, &ev, 1.0e-7);
+    }
+}
+
+#[test]
+fn fixed_lower() {
+    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]);
+    let (e, vecs): (Array1<_>, Array2<_>) = (&a).eigh(UPLO::Lower).unwrap();
+    assert_close_l2!(&e, &arr1(&[2.0, 2.0, 5.0]), 1.0e-7);
+
+    // Check eigenvectors are orthogonalized
+    let s = vecs.t().dot(&vecs);
+    assert_close_l2!(&s, &Array::eye(3), 1.0e-7);
+
+    for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
+        let av = a.dot(&v);
+        let ev = v.mapv(|x| e[i] * x);
+        assert_close_l2!(&av, &ev, 1.0e-7);
+    }
+}
+
+#[test]
+fn fixed_t_lower() {
+    let a = arr2(&[[3.0, 1.0, 1.0], [1.0, 3.0, 1.0], [1.0, 1.0, 3.0]]).reversed_axes();
+    let (e, vecs): (Array1<_>, Array2<_>) = (&a).eigh(UPLO::Lower).unwrap();
+    assert_close_l2!(&e, &arr1(&[2.0, 2.0, 5.0]), 1.0e-7);
+
+    // Check eigenvectors are orthogonalized
+    let s = vecs.t().dot(&vecs);
+    assert_close_l2!(&s, &Array::eye(3), 1.0e-7);
+
+    for (i, v) in vecs.axis_iter(Axis(1)).enumerate() {
+        let av = a.dot(&v);
+        let ev = v.mapv(|x| e[i] * x);
+        assert_close_l2!(&av, &ev, 1.0e-7);
     }
 }
 
@@ -48,3 +94,29 @@ fn ssqrt_t() {
     println!("ss = {:?}", &ss);
     assert_close_l2!(&ss, &ans, 1e-7);
 }
+
+#[test]
+fn ssqrt_lower() {
+    let a: Array2<f64> = random_hpd(3);
+    let ans = a.clone();
+    let s = a.ssqrt(UPLO::Lower).unwrap();
+    println!("a = {:?}", &ans);
+    println!("s = {:?}", &s);
+    assert_close_l2!(&s.t(), &s, 1e-7);
+    let ss = s.dot(&s);
+    println!("ss = {:?}", &ss);
+    assert_close_l2!(&ss, &ans, 1e-7);
+}
+
+#[test]
+fn ssqrt_t_lower() {
+    let a: Array2<f64> = random_hpd(3).reversed_axes();
+    let ans = a.clone();
+    let s = a.ssqrt(UPLO::Lower).unwrap();
+    println!("a = {:?}", &ans);
+    println!("s = {:?}", &s);
+    assert_close_l2!(&s.t(), &s, 1e-7);
+    let ss = s.dot(&s);
+    println!("ss = {:?}", &ss);
+    assert_close_l2!(&ss, &ans, 1e-7);
+}