anchor.diff:          [<DllImport("libopenblas", EntryPoint="sgetri_")>]
           extern void sgetri_(int *n, float32 *a, int *lda, int *ipiv, float32 *work, int *lwork, int *info)
   
  +        // Only works for square matrices
           let sgesv(a:float32[,], b:float32[]) =
  -            let m = Array2D.length1 a
  -            let n = Array2D.length2 a
  -            let a' = Array2D.zeroCreate n m
  -            for i = 0 to n - 1 do
  -                for j = 0 to m - 1 do
  -                    a'.[i, j] <- a.[j, i]
  -            let b' = Array.copy b
  +            let n = Array2D.length1 a
  +            BLASExtensions.simatcopyT(1.f, a)
               let ipiv = Array.zeroCreate n
               let mutable arg_n = n
               let mutable arg_nrhs = 1
               let mutable arg_lda = n
               let mutable arg_ldb = n
               let mutable arg_info = 0
  -            use arg_a = new PinnedArray2D<float32>(a')
  +            use arg_a = new PinnedArray2D<float32>(a)
               use arg_ipiv = new PinnedArray<int>(ipiv)
  -            use arg_b = new PinnedArray<float32>(b')
  +            use arg_b = new PinnedArray<float32>(b)
               sgesv_(&&arg_n, &&arg_nrhs, arg_a.Ptr, &&arg_lda, arg_ipiv.Ptr, arg_b.Ptr, &&arg_ldb, &&arg_info)
               if arg_info = 0 then
  -                Some(b')
  +                Some(b)
               else
                   None
   
           let ssysv(a:float32[,], b:float32[]) =
               let n = Array2D.length1 a
  -            let b' = Array.copy b
               let ipiv = Array.zeroCreate n
  -            let work = Array.zeroCreate 1
  +            let work = Array.zeroCreate n
               let mutable arg_uplo = 'U' // Assume upper triangular. TODO: check if LAPACK implementation requires the lower triangle to be zeroed
               let mutable arg_n = n
               let mutable arg_nrhs = 1
               let mutable arg_lda = n
               let mutable arg_ldb = n
  -            let mutable arg_lwork = 1
  +            let mutable arg_lwork = work.Length
               let mutable arg_info = 0
               use arg_a = new PinnedArray2D<float32>(a)
               use arg_ipiv = new PinnedArray<int>(ipiv)
  
positive.code:         [<DllImport("libopenblas", EntryPoint="dgetri_")>]
          extern void dgetri_(int *n, float *a, int *lda, int *ipiv, float *work, int *lwork, int *info)
  
          let dgesv(a:float[,], b:float[]) =
              let m = Array2D.length1 a
              let n = Array2D.length2 a
              let a' = Array2D.zeroCreate n m
              for i = 0 to n - 1 do
                  for j = 0 to m - 1 do
                      a'.[i, j] <- a.[j, i]
              let b' = Array.copy b
              let ipiv = Array.zeroCreate n
              let mutable arg_n = n
              let mutable arg_nrhs = 1
              let mutable arg_lda = n
              let mutable arg_ldb = n
              let mutable arg_info = 0
              use arg_a = new PinnedArray2D<float>(a')
              use arg_ipiv = new PinnedArray<int>(ipiv)
              use arg_b = new PinnedArray<float>(b')
              dgesv_(&&arg_n, &&arg_nrhs, arg_a.Ptr, &&arg_lda, arg_ipiv.Ptr, arg_b.Ptr, &&arg_ldb, &&arg_info)
              if arg_info = 0 then
                  Some(b')
              else
                  None
  
          let dsysv(a:float[,], b:float[]) =
              let n = Array2D.length1 a
              let b' = Array.copy b
              let ipiv = Array.zeroCreate n
              let work = Array.zeroCreate 1
              let mutable arg_uplo = 'U' // Assume upper triangular. TODO: check if LAPACK implementation requires the lower triangle to be zeroed
              let mutable arg_n = n
              let mutable arg_nrhs = 1
              let mutable arg_lda = n
              let mutable arg_ldb = n
              let mutable arg_lwork = 1
              let mutable arg_info = 0
              use arg_a = new PinnedArray2D<float>(a)
              use arg_ipiv = new PinnedArray<int>(ipiv)
  
negative.code:             cblas_somatcopy(arg_ordering, arg_trans, arg_rows, arg_cols, arg_alpha, arg_a.Ptr, arg_lda, arg_b.Ptr, arg_ldb)
  
          [<SuppressUnmanagedCodeSecurity>]
          [<DllImport("libopenblas", EntryPoint="cblas_dimatcopy")>]
          extern void cblas_dimatcopy(int ordering, int trans, int rows, int cols, float alpha, float *a, int lda, int ldb)
  
          [<SuppressUnmanagedCodeSecurity>]
          [<DllImport("libopenblas", EntryPoint="cblas_domatcopy")>]
          extern void cblas_domatcopy(int ordering, int trans, int rows, int cols, float alpha, float *a, int lda, float *b, int ldb)
  
          // A <- alpha * transpose(A)
          // Only works for square matrices, modification of .NET array metadata might be needed for the correct in-memory transposition of non-square matrices
          let dimatcopyT(alpha:float, a:float[,]) =
              let m = Array2D.length1 a
              let n = Array2D.length2 a
              let arg_ordering = 101 // cblas.h: typedef enum {CblasRowMajor=101, CblasColMajor=102} CBLAS_LAYOUT;
              let arg_trans = 112 // cblas.h: typedef enum {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113} CBLAS_TRANSPOSE;
              let arg_rows = m
              let arg_cols = n
              let arg_alpha = alpha
              let arg_a = new PinnedArray2D<float>(a)
              let arg_lda = n
              let arg_ldb = m
              cblas_dimatcopy(arg_ordering, arg_trans, arg_rows, arg_cols, arg_alpha, arg_a.Ptr, arg_lda, arg_ldb)
  
          // B <- alpha * transpose(A)
          let domatcopyT(alpha:float, a:float[,], b:float[,]) =
  