Source file owl_lapacke.ml

# 1 "src/owl/lapacke/owl_lapacke.ml"
(*
 * OWL - OCaml Scientific and Engineering Computing
 * Copyright (c) 2016-2019 Liang Wang <liang.wang@cl.cam.ac.uk>
 *)

(** LAPACKE interface: high-level interface between Owl and LAPACKE *)

(** Please refer to the documentation of Intel Math Kernel Library on the
  LAPACKE Interface. The interface implemented here is compatible the those
  documented on their website.
  url: https://software.intel.com/en-us/mkl-developer-reference-c
 *)

open Ctypes

open Bigarray


module L = Owl_lapacke_generated


type ('a, 'b) t = ('a, 'b, c_layout) Genarray.t


type lapacke_layout = RowMajor | ColMajor
let lapacke_layout : type a. a layout -> int = function
  | C_layout       -> 101
  | Fortran_layout -> 102

type lapacke_transpose = NoTrans | Trans | ConjTrans
let lapacke_transpose = function
    NoTrans -> 'N' | Trans -> 'T' | ConjTrans -> 'C'

type lapacke_uplo = Upper | Lower
let _lapacke_uplo = function Upper -> 121 | Lower -> 122

type lapacke_diag = NonUnit | Unit
let _lapacke_diag = function NonUnit -> 131 | Unit -> 132

type lapacke_side = Left | Right
let _lapacke_side = function Left -> 141 | Right -> 142

let check_lapack_error ret =
  if ret = 0 then ()
  else if ret < 0 then
    raise (Invalid_argument (string_of_int ret))
  else
    failwith (Printf.sprintf "LAPACKE: %i" ret)

(* calculate the leading dimension of a ~matrix *)
let _stride : type a b c. (a, b, c) Genarray.t -> int = fun x ->
  match (Genarray.layout x) with
  | C_layout       -> (Genarray.dims x).(1)
  | Fortran_layout -> (Genarray.dims x).(0)

let _row_num x = (Genarray.dims x).(0)

let _col_num x = (Genarray.dims x).(1)


let _gbtrf
  : type a b. kl:int -> ku:int -> m:int -> ab:(a, b) t -> (a, b) t * (int32, int32_elt) t
  = fun ~kl ~ku ~m ~ab ->
  let n = Owl_dense_matrix_generic.col_num ab in
  let minmn = Stdlib.min m n in
  let _kind = Genarray.kind ab in
  let _layout = Genarray.layout ab in
  let layout = lapacke_layout _layout in

  assert (kl >= 0 && ku >=0 && m >= 0 && n >= 0);

  let ipiv = Genarray.create int32 _layout [|minmn|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _ab = bigarray_start Ctypes_static.Genarray ab in
  let ldab = _stride ab in

  let ret = match _kind with
    | Float32   -> L.sgbtrf ~layout ~m ~n ~kl ~ku ~ab:_ab ~ldab ~ipiv:_ipiv
    | Float64   -> L.dgbtrf ~layout ~m ~n ~kl ~ku ~ab:_ab ~ldab ~ipiv:_ipiv
    | Complex32 -> L.cgbtrf ~layout ~m ~n ~kl ~ku ~ab:_ab ~ldab ~ipiv:_ipiv
    | Complex64 -> L.zgbtrf ~layout ~m ~n ~kl ~ku ~ab:_ab ~ldab ~ipiv:_ipiv
    | _         -> failwith "lapacke:gbtrf"
  in
  check_lapack_error ret;
  ab, ipiv


let gbtrs
  : type a b. trans:lapacke_transpose -> kl:int -> ku:int -> n:int
  -> ab:(a, b) t -> ipiv:(int32, int32_elt) t -> b:(a, b) t -> unit
  = fun ~trans ~kl ~ku ~n ~ab ~ipiv ~b ->
    let m = Owl_dense_matrix_generic.col_num ab in
    assert (n = m && n = Owl_dense_matrix_generic.row_num b);
    let nrhs = Owl_dense_matrix_generic.col_num b in
    let _kind = Genarray.kind ab in
    let _layout = Genarray.layout ab in
    let layout = lapacke_layout _layout in
    let trans = lapacke_transpose trans in

    let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
    let _ab = bigarray_start Ctypes_static.Genarray ab in
    let _b = bigarray_start Ctypes_static.Genarray b in
    let ldab = _stride ab in
    let ldb = _stride b in

    let ret = match _kind with
      | Float32   -> L.sgbtrs ~layout ~trans ~n ~kl ~ku ~nrhs ~ab:_ab ~ldab
                       ~ipiv:_ipiv ~b:_b ~ldb
      | Float64   -> L.dgbtrs ~layout ~trans ~n ~kl ~ku ~nrhs ~ab:_ab ~ldab
                       ~ipiv:_ipiv ~b:_b ~ldb
      | Complex32 -> L.cgbtrs ~layout ~trans ~n ~kl ~ku ~nrhs ~ab:_ab ~ldab
                       ~ipiv:_ipiv ~b:_b ~ldb
      | Complex64 -> L.zgbtrs ~layout ~trans ~n ~kl ~ku ~nrhs ~ab:_ab ~ldab
                       ~ipiv:_ipiv ~b:_b ~ldb

      | _         -> failwith "lapacke:gbtrs"
    in
    check_lapack_error ret


let gebal
  : type a b. ?job:char -> a:(a, b) t -> int * int * (a, b) t
  = fun ?(job='B') ~a ->
  assert (job = 'N' || job = 'P' || job = 'S' || job = 'B');
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ilo = Ctypes.(allocate int32_t 0l) in
  let _ihi = Ctypes.(allocate int32_t 0l) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let scale = ref (Genarray.create _kind _layout [|0;0|]) in

  let ret = match _kind with
    | Float32   -> (
        let scale' = Genarray.create float32 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let r = L.sgebal ~layout ~job ~n ~a:_a ~lda ~ilo:_ilo ~ihi:_ihi
            ~scale:_scale in
        scale := scale';
        r
      )
    | Float64   -> (
        let scale' = Genarray.create float64 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let r = L.dgebal ~layout ~job ~n ~a:_a ~lda ~ilo:_ilo ~ihi:_ihi
            ~scale:_scale in
        scale := scale';
        r
      )
    | Complex32 -> (
        let scale' = Genarray.create float32 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let r = L.cgebal ~layout ~job ~n ~a:_a ~lda ~ilo:_ilo ~ihi:_ihi
            ~scale:_scale in
        scale := Owl_dense_matrix_generic.cast_s2c scale';
        r
      )
    | Complex64 -> (
        let scale' = Genarray.create float64 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let r = L.zgebal ~layout ~job ~n ~a:_a ~lda ~ilo:_ilo ~ihi:_ihi
            ~scale:_scale in
        scale := Owl_dense_matrix_generic.cast_d2z scale';
        r
      )
    | _         -> failwith "lapacke:gebal"
  in
  check_lapack_error ret;
  let ilo = Int32.to_int !@_ilo in
  let ihi = Int32.to_int !@_ihi in
  ilo, ihi, !scale


(* TODO: need a solution for scale parameter *)
let gebak
  : type a b. job:char -> side:char
  -> ilo:int -> ihi:int -> scale:(float ptr) -> v:(a, b) t -> unit
  = fun ~job ~side ~ilo ~ihi ~scale ~v ->
  assert (side = 'L' || side = 'R');
  assert (job = 'N' || job = 'P' || job = 'S' || job = 'B');
  let m = Owl_dense_matrix_generic.row_num v in
  let n = Owl_dense_matrix_generic.col_num v in
  assert (m = n);
  let _kind = Genarray.kind v in
  let _layout = Genarray.layout v in
  let layout = lapacke_layout _layout in

  let _v = bigarray_start Ctypes_static.Genarray v in
  let ldv = _stride v in

  let ret = match _kind with
    | Float32   -> L.sgebak ~layout ~job ~side ~n ~ilo ~ihi ~scale ~m
                     ~v:_v ~ldv
    | Float64   -> L.dgebak ~layout ~job ~side ~n ~ilo ~ihi ~scale ~m
                     ~v:_v ~ldv
    | Complex32 -> L.cgebak ~layout ~job ~side ~n ~ilo ~ihi ~scale ~m
                     ~v:_v ~ldv
    | Complex64 -> L.zgebak ~layout ~job ~side ~n ~ilo ~ihi ~scale ~m
                     ~v:_v ~ldv
    | _         -> failwith "lapacke:gebak"
  in
  check_lapack_error ret


let gebrd
  : type a b.
    a:(a, b) t -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let d = ref (Genarray.create _kind _layout [|0;k|]) in
  let e = ref (Genarray.create _kind _layout [|0;k|]) in
  let tauq = Genarray.create _kind _layout [|1;k|] in
  let taup = Genarray.create _kind _layout [|1;k|] in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tauq = bigarray_start Ctypes_static.Genarray tauq in
  let _taup = bigarray_start Ctypes_static.Genarray taup in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> (
        let d' = Genarray.create float32 _layout [|1;k|] in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        let e' = Genarray.create float32 _layout [|1;k|] in
        let _e = bigarray_start Ctypes_static.Genarray e' in
        let r = L.sgebrd ~layout ~m ~n ~a:_a ~lda ~d:_d ~e:_e ~tauq:_tauq ~taup:_taup in
        d := d';
        e := e';
        r
      )
    | Float64   -> (
        let d' = Genarray.create float64 _layout [|1;k|] in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        let e' = Genarray.create float64 _layout [|1;k|] in
        let _e = bigarray_start Ctypes_static.Genarray e' in
        let r = L.dgebrd ~layout ~m ~n ~a:_a ~lda ~d:_d ~e:_e ~tauq:_tauq ~taup:_taup in
        d := d';
        e := e';
        r
      )
    | Complex32 -> (
        let d' = Genarray.create float32 _layout [|1;k|] in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        let e' = Genarray.create float32 _layout [|1;k|] in
        let _e = bigarray_start Ctypes_static.Genarray e' in
        let r = L.cgebrd ~layout ~m ~n ~a:_a ~lda ~d:_d ~e:_e ~tauq:_tauq ~taup:_taup in
        d := Owl_dense_matrix_generic.cast_s2c d';
        e := Owl_dense_matrix_generic.cast_s2c e';
        r
      )
    | Complex64 -> (
        let d' = Genarray.create float64 _layout [|1;k|] in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        let e' = Genarray.create float64 _layout [|1;k|] in
        let _e = bigarray_start Ctypes_static.Genarray e' in
        let r = L.zgebrd ~layout ~m ~n ~a:_a ~lda ~d:_d ~e:_e ~tauq:_tauq ~taup:_taup in
        d := Owl_dense_matrix_generic.cast_d2z d';
        e := Owl_dense_matrix_generic.cast_d2z e';
        r
      )
    | _         -> failwith "lapacke:gebrd"
  in
  check_lapack_error ret;
  a, !d, !e, tauq, taup


let gelqf
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1;k|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> L.sgelqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dgelqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.cgelqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zgelqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:gelqf"
  in
  check_lapack_error ret;
  a, tau


let geqlf
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1;k|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> L.sgeqlf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dgeqlf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.cgeqlf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zgeqlf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:geqlf"
  in
  check_lapack_error ret;
  a, tau


let geqrf
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1;k|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> L.sgeqrf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dgeqrf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.cgeqrf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zgeqrf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:geqrf"
  in
  check_lapack_error ret;
  a, tau


let gerqf
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1;k|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> L.sgerqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dgerqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.cgerqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zgerqf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:gerqf"
  in
  check_lapack_error ret;
  a, tau


let geqp3
  : type a b. ?jpvt:(int32, int32_elt) t -> a:(a, b) t
  -> (a, b) t * (a, b) t * (int32, int32_elt) t
  = fun ?jpvt ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let jpvt = match jpvt with
    | Some jpvt -> jpvt
    | None      -> (
        let jpvt = Genarray.create int32 _layout [|1;n|] in
        Genarray.fill jpvt 0l;
        jpvt
      )
  in
  assert (n = Owl_dense_matrix_generic.col_num jpvt);

  let tau = Genarray.create _kind _layout [|1;k|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _jpvt = bigarray_start Ctypes_static.Genarray jpvt in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = _stride a in

  let ret = match _kind with
    | Float32   -> L.sgeqp3 ~layout ~m ~n ~a:_a ~lda ~jpvt:_jpvt ~tau:_tau
    | Float64   -> L.dgeqp3 ~layout ~m ~n ~a:_a ~lda ~jpvt:_jpvt ~tau:_tau
    | Complex32 -> L.cgeqp3 ~layout ~m ~n ~a:_a ~lda ~jpvt:_jpvt ~tau:_tau
    | Complex64 -> L.zgeqp3 ~layout ~m ~n ~a:_a ~lda ~jpvt:_jpvt ~tau:_tau
    | _         -> failwith "lapacke:geqp3"
  in
  check_lapack_error ret;
  a, tau, jpvt


let geqrt
  : type a b. nb:int -> a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~nb ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  assert (nb <= minmn);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _a = bigarray_start Ctypes_static.Genarray a in
  (* FIXME: there ~might be something wrong with the lapacke interface. The
    behaviour of this function is not consistent with what has been documented
    on Intel's MKL website. I.e., if we allocate [nb x minmn] space for t, it
    is likely there will be ~memory fault. The lapacke code turns out to use
    [minmn x minmn] space actually.
  *)
  let t = Genarray.create _kind _layout [|minmn; minmn|] in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let lda = Stdlib.max 1 (_stride a) in
  let ldt = Stdlib.max 1 (_stride t) in

  let ret = match _kind with
    | Float32   -> L.sgeqrt ~layout ~m ~n ~nb ~a:_a ~lda ~t:_t ~ldt
    | Float64   -> (
        Genarray.fill t 0.7;
        L.dgeqrt ~layout ~m ~n ~nb ~a:_a ~lda ~t:_t ~ldt )
    | Complex32 -> L.cgeqrt ~layout ~m ~n ~nb ~a:_a ~lda ~t:_t ~ldt
    | Complex64 -> L.zgeqrt ~layout ~m ~n ~nb ~a:_a ~lda ~t:_t ~ldt
    | _         -> failwith "lapacke:geqrt"
  in
  check_lapack_error ret;
  (* resize to the shape of [t] to that is supposed to be *)
  let t = Owl_dense_matrix_generic.resize t [|nb; minmn|] in
  a, t


let geqrt3
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m >= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _a = bigarray_start Ctypes_static.Genarray a in
  let t = Genarray.create _kind _layout [|n;n|] in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let lda = Stdlib.max 1 (_stride a) in
  let ldt = Stdlib.max 1 (_stride t) in

  let ret = match _kind with
    | Float32   -> L.sgeqrt3 ~layout ~m ~n ~a:_a ~lda ~t:_t ~ldt
    | Float64   -> L.dgeqrt3 ~layout ~m ~n ~a:_a ~lda ~t:_t ~ldt
    | Complex32 -> L.cgeqrt3 ~layout ~m ~n ~a:_a ~lda ~t:_t ~ldt
    | Complex64 -> L.zgeqrt3 ~layout ~m ~n ~a:_a ~lda ~t:_t ~ldt
    | _         -> failwith "lapacke:geqrt3"
  in
  check_lapack_error ret;
  a, t


let getrf
  : type a b. a:(a, b) t -> (a, b) t * (int32, int32_elt) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1; minmn|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.sgetrf ~layout ~m ~n ~a:_a ~lda ~ipiv:_ipiv
    | Float64   -> L.dgetrf ~layout ~m ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex32 -> L.cgetrf ~layout ~m ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zgetrf ~layout ~m ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:getrf"
  in
  check_lapack_error ret;
  a, ipiv


let tzrzf
  : type a b. a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m <= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1; m|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.stzrzf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dtzrzf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.ctzrzf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.ztzrzf ~layout ~m ~n ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:tzrzf"
  in
  check_lapack_error ret;
  a, tau


let ormrz
  : type a. side:char -> trans:char -> a:(float, a) t -> tau:(float, a) t
  -> c:(float, a) t -> (float, a) t
  = fun ~side ~trans ~a ~tau ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let k = Owl_dense_matrix_generic.((row_num tau) * (col_num tau)) in
  let l = Owl_dense_matrix_generic.((col_num a) - (row_num a)) in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _a = bigarray_start Ctypes_static.Genarray a in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let lda = Stdlib.max 1 (_stride a) in
  let ldc = Stdlib.max 1 (_stride c) in

  let ret = match _kind with
    | Float32   -> L.sormrz ~layout ~side ~trans ~m ~n ~k ~l ~a:_a ~lda
                     ~tau:_tau ~c:_c ~ldc
    | Float64   -> L.dormrz ~layout ~side ~trans ~m ~n ~k ~l ~a:_a ~lda
                     ~tau:_tau ~c:_c ~ldc
  in
  check_lapack_error ret;
  c


let gels
  : type a b. trans:char -> a:(a, b) t -> b:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t
  = fun ~trans ~a ~b ->
  assert (trans = 'N' || trans = 'T' || trans = 'C');
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  if trans = 'N' then
    assert (mb = m)
  else
    assert (mb = n);

  let l = Stdlib.max m n in
  let b = match mb < l with
    | true  -> Owl_dense_matrix_generic.resize b [|l; nb|]
    | false -> b
  in

  let nrhs = nb in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> L.sgels ~layout ~trans ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Float64   -> L.dgels ~layout ~trans ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex32 -> L.cgels ~layout ~trans ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex64 -> L.zgels ~layout ~trans ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | _         -> failwith "lapacke:gels"
  in
  check_lapack_error ret;

  let k = Stdlib.min m n in
  let a' = Owl_dense_matrix_generic.get_slice [[0;k-1]; [0;k-1]] a in
  let f = match m < n with
    | true  -> Owl_dense_matrix_generic.tril a'
    | false -> Owl_dense_matrix_generic.triu a'
  in
  let sol = match trans = 'N' with
    | true  -> Owl_dense_matrix_generic.resize b [|n; nb|]
    | false -> Owl_dense_matrix_generic.resize b [|m; nb|]
  in
  let ssr = match trans = 'N' with
    | true  ->
        if mb > n then
          Owl_dense_matrix_generic.resize ~head:false b [|mb - n; nb|]
        else
          Genarray.create _kind _layout [|0;0|]
    | false ->
        if mb > m then
          Owl_dense_matrix_generic.resize ~head:false b [|mb - m; nb|]
        else
          Genarray.create _kind _layout [|0;0|]
  in
  f, sol, ssr


let gesv
  : type a b. a:(a, b) t -> b:(a, b) t -> (a, b) t * (a, b) t * (int32, int32_elt) t
  = fun ~a ~b ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (m = n && mb = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let nrhs = nb in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in

  let ret = match _kind with
    | Float32   -> L.sgesv ~layout ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Float64   -> L.dgesv ~layout ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex32 -> L.cgesv ~layout ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex64 -> L.zgesv ~layout ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | _         -> failwith "lapacke:gesv"
  in
  check_lapack_error ret;
  a, b, ipiv


let getrs
  : type a b. trans:char -> a:(a, b) t -> ipiv:(int32, int32_elt) t -> b:(a, b) t -> (a, b) t
  = fun ~trans ~a ~ipiv ~b ->
  assert (trans = 'N' || trans = 'T' || trans = 'C');
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (m = n && mb = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let nrhs = nb in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> L.sgetrs ~layout ~trans ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Float64   -> L.dgetrs ~layout ~trans ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex32 -> L.cgetrs ~layout ~trans ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex64 -> L.zgetrs ~layout ~trans ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | _         -> failwith "lapacke:getrs"
  in
  check_lapack_error ret;
  b


let getri
  : type a b. a:(a, b) t -> ipiv:(int32, int32_elt) t -> (a, b) t
  = fun ~a ~ipiv ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let k = Owl_dense_matrix_generic.((row_num ipiv) * (col_num ipiv)) in
  assert (m = n && n = k);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.sgetri ~layout ~n ~a:_a ~lda ~ipiv:_ipiv
    | Float64   -> L.dgetri ~layout ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex32 -> L.cgetri ~layout ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zgetri ~layout ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:getri"
  in
  check_lapack_error ret;
  a


let gesvx
  : type a b. fact:char -> trans:char -> a:(a, b) t -> af:(a, b) t
  -> ipiv:(int32, int32_elt) t -> equed:char -> r:(a, b) t -> c:(a, b) t -> b:(a, b) t
  -> (a, b) t * char * (a, b) t * (a, b) t * (a, b) t * a * (a, b) t * (a, b) t * a
  = fun ~fact ~trans ~a ~af ~ipiv ~equed ~r ~c ~b ->
  assert (fact = 'E' || fact = 'N' || fact = 'T' || fact = 'C');
  assert (trans = 'N' || trans = 'T' || trans = 'C');
  assert (equed = 'N' || equed = 'R' || equed = 'C' || equed = 'B');
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  assert (Owl_dense_matrix_generic.row_num af = n && Owl_dense_matrix_generic.col_num af = n);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ldaf = Stdlib.max 1 (_stride af) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _af = bigarray_start Ctypes_static.Genarray af in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _equed = Ctypes.(allocate char equed) in
  let x = Genarray.create _kind _layout [|n;nrhs|] in
  let _x = bigarray_start Ctypes_static.Genarray x in
  let ldx = Stdlib.max 1 (_stride x) in

  let r_ref = ref (Genarray.create _kind _layout [|0;0|]) in
  let c_ref = ref (Genarray.create _kind _layout [|0;0|]) in
  let rcond = ref (Genarray.create _kind _layout [|0;0|]) in
  let ferr = ref (Genarray.create _kind _layout [|0;0|]) in
  let berr = ref (Genarray.create _kind _layout [|0;0|]) in
  let rpivot = ref (Genarray.create _kind _layout [|0;0|]) in

  let ret = match _kind with
    | Float32   -> (
        let rcond' = Genarray.create float32 _layout [|1;1|] in
        let ferr' = Genarray.create float32 _layout [|1;nrhs|] in
        let berr' = Genarray.create float32 _layout [|1;nrhs|] in
        let rpivot' = Genarray.create float32 _layout [|1;1|] in
        let _rcond = bigarray_start Ctypes_static.Genarray rcond' in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let _rpivot = bigarray_start Ctypes_static.Genarray rpivot' in

        let _r = bigarray_start Ctypes_static.Genarray r in
        let _c = bigarray_start Ctypes_static.Genarray c in

        let ret = L.sgesvx ~layout ~fact ~trans ~n ~nrhs ~a:_a ~lda ~af:_af ~ldaf ~ipiv:_ipiv ~equed:_equed ~r:_r ~c:_c ~b:_b ~ldb ~x:_x ~ldx ~rcond:_rcond ~ferr:_ferr ~berr:_berr ~rpivot:_rpivot in
        r_ref := r;
        c_ref := c;
        rcond := rcond';
        ferr := ferr';
        berr := berr';
        rpivot := rpivot';
        ret
      )
    | Float64   -> (
        let rcond' = Genarray.create float64 _layout [|1;1|] in
        let ferr' = Genarray.create float64 _layout [|1;nrhs|] in
        let berr' = Genarray.create float64 _layout [|1;nrhs|] in
        let rpivot' = Genarray.create float64 _layout [|1;1|] in
        let _rcond = bigarray_start Ctypes_static.Genarray rcond' in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let _rpivot = bigarray_start Ctypes_static.Genarray rpivot' in

        let _r = bigarray_start Ctypes_static.Genarray r in
        let _c = bigarray_start Ctypes_static.Genarray c in

        let ret = L.dgesvx ~layout ~fact ~trans ~n ~nrhs ~a:_a ~lda ~af:_af ~ldaf ~ipiv:_ipiv ~equed:_equed ~r:_r ~c:_c ~b:_b ~ldb ~x:_x ~ldx ~rcond:_rcond ~ferr:_ferr ~berr:_berr ~rpivot:_rpivot in
        r_ref := r;
        c_ref := c;
        rcond := rcond';
        ferr := ferr';
        berr := berr';
        rpivot := rpivot';
        ret
      )
    | Complex32 -> (
        let rcond' = Genarray.create float32 _layout [|1;1|] in
        let ferr' = Genarray.create float32 _layout [|1;nrhs|] in
        let berr' = Genarray.create float32 _layout [|1;nrhs|] in
        let rpivot' = Genarray.create float32 _layout [|1;1|] in
        let _rcond = bigarray_start Ctypes_static.Genarray rcond' in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let _rpivot = bigarray_start Ctypes_static.Genarray rpivot' in

        let r' = Owl_dense_matrix_c.re r in
        let c' = Owl_dense_matrix_c.re c in
        let _r = bigarray_start Ctypes_static.Genarray r' in
        let _c = bigarray_start Ctypes_static.Genarray c' in

        let ret = L.cgesvx ~layout ~fact ~trans ~n ~nrhs ~a:_a ~lda ~af:_af ~ldaf ~ipiv:_ipiv ~equed:_equed ~r:_r ~c:_c ~b:_b ~ldb ~x:_x ~ldx ~rcond:_rcond ~ferr:_ferr ~berr:_berr ~rpivot:_rpivot in
        r_ref := Owl_dense_matrix_generic.cast_s2c r';
        c_ref := Owl_dense_matrix_generic.cast_s2c c';
        rcond := Owl_dense_matrix_generic.cast_s2c rcond';
        ferr := Owl_dense_matrix_generic.cast_s2c ferr';
        berr := Owl_dense_matrix_generic.cast_s2c berr';
        rpivot := Owl_dense_matrix_generic.cast_s2c rpivot';
        ret
      )
    | Complex64 -> (
        let _rpivot = Ctypes.(allocate double 0.) in
        let rcond' = Genarray.create float64 _layout [|1;1|] in
        let ferr' = Genarray.create float64 _layout [|1;nrhs|] in
        let berr' = Genarray.create float64 _layout [|1;nrhs|] in
        let rpivot' = Genarray.create float64 _layout [|1;1|] in
        let _rcond = bigarray_start Ctypes_static.Genarray rcond' in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let _rpivot = bigarray_start Ctypes_static.Genarray rpivot' in

        let r' = Owl_dense_matrix_z.re r in
        let c' = Owl_dense_matrix_z.re c in
        let _r = bigarray_start Ctypes_static.Genarray r' in
        let _c = bigarray_start Ctypes_static.Genarray c' in

        let ret = L.zgesvx ~layout ~fact ~trans ~n ~nrhs ~a:_a ~lda ~af:_af ~ldaf ~ipiv:_ipiv ~equed:_equed ~r:_r ~c:_c ~b:_b ~ldb ~x:_x ~ldx ~rcond:_rcond ~ferr:_ferr ~berr:_berr ~rpivot:_rpivot in
        r_ref := Owl_dense_matrix_generic.cast_d2z r';
        c_ref := Owl_dense_matrix_generic.cast_d2z c';
        rcond := Owl_dense_matrix_generic.cast_d2z rcond';
        ferr := Owl_dense_matrix_generic.cast_d2z ferr';
        berr := Owl_dense_matrix_generic.cast_d2z berr';
        rpivot := Owl_dense_matrix_generic.cast_d2z rpivot';
        ret
      )
    | _         -> failwith "lapacke:gesvx"
  in
  check_lapack_error ret;
  if ret = n + 1 then Owl_log.warn "matrix is singular to working precision.";
  let rcond = Owl_dense_matrix_generic.get !rcond 0 0 in
  let rpivot = Owl_dense_matrix_generic.get !rpivot 0 0 in
  x, !@_equed, !r_ref, !c_ref, b, rcond, !ferr, !berr, rpivot


let gelsd
  : type a b. a:(a, b) t -> b:(a, b) t -> rcond:float -> (a, b) t * int
  = fun ~a ~b ~rcond ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (mb = m);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let b = match mb < n with
    | true  -> Owl_dense_matrix_generic.resize b [|n; nrhs|]
    | false -> b
  in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _rank = Ctypes.(allocate int32_t 0l) in

  let ret = match _kind with
    | Float32   -> (
        let s = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s in
        L.sgelsd ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~s:_s ~rcond ~rank:_rank
      )
    | Float64   -> (
        let s = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s in
        L.dgelsd ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~s:_s ~rcond ~rank:_rank
      )
    | Complex32 -> (
        let s = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s in
        L.cgelsd ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~s:_s ~rcond ~rank:_rank
      )
    | Complex64 -> (
        let s = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s in
        L.zgelsd ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~s:_s ~rcond ~rank:_rank
      )
    | _         -> failwith "lapacke:gelsd"
  in
  check_lapack_error ret;
  let b = Owl_dense_matrix_generic.resize b [|n; nrhs|] in
  b, Int32.to_int !@_rank


let gelsy
  : type a b. a:(a, b) t -> b:(a, b) t -> rcond:float -> (a, b) t * int
  = fun ~a ~b ~rcond ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (mb = m);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let b = match mb < n with
    | true  -> Owl_dense_matrix_generic.resize b [|n; nrhs|]
    | false -> b
  in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _rank = Ctypes.(allocate int32_t 0l) in

  let jpvt = Genarray.create int32 _layout [|1;n|] in
  Genarray.fill jpvt 0l;
  let _jpvt = bigarray_start Ctypes_static.Genarray jpvt in

  let ret = match _kind with
    | Float32   -> L.sgelsy ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~jpvt:_jpvt ~rcond ~rank:_rank
    | Float64   -> L.dgelsy ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~jpvt:_jpvt ~rcond ~rank:_rank
    | Complex32 -> L.cgelsy ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~jpvt:_jpvt ~rcond ~rank:_rank
    | Complex64 -> L.zgelsy ~layout ~m ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~jpvt:_jpvt ~rcond ~rank:_rank
    | _         -> failwith "lapacke:gelsy"
  in
  check_lapack_error ret;
  let b = Owl_dense_matrix_generic.resize b [|n; nrhs|] in
  b, Int32.to_int !@_rank


let gglse
  : type a b. a:(a, b) t -> b:(a, b) t -> c:(a, b) t -> d:(a, b) t
  -> (a, b) t * a
  = fun ~a ~b ~c ~d ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let p = Owl_dense_matrix_generic.row_num b in
  assert (n = Owl_dense_matrix_generic.col_num b);
  assert (m = Owl_dense_matrix_generic.((row_num c) * (col_num c)));
  assert (p = Owl_dense_matrix_generic.((row_num d) * (col_num d)));
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let x = Genarray.create _kind _layout [|1;n|] in
  let _x = bigarray_start Ctypes_static.Genarray x in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> L.sgglse ~layout ~m ~n ~p ~a:_a ~lda ~b:_b ~ldb ~c:_c ~d:_d ~x:_x
    | Float64   -> L.dgglse ~layout ~m ~n ~p ~a:_a ~lda ~b:_b ~ldb ~c:_c ~d:_d ~x:_x
    | Complex32 -> L.cgglse ~layout ~m ~n ~p ~a:_a ~lda ~b:_b ~ldb ~c:_c ~d:_d ~x:_x
    | Complex64 -> L.zgglse ~layout ~m ~n ~p ~a:_a ~lda ~b:_b ~ldb ~c:_c ~d:_d ~x:_x
    | _         -> failwith "lapacke:gglse"
  in
  check_lapack_error ret;

  let c' = Owl_dense_matrix_generic.resize ~head:false c [|1; m - n + p|] in
  let res = Owl_dense_matrix_generic.(mul c' c' |> sum') in
  x, res


let geev
  : type a b. jobvl:char -> jobvr:char -> a:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~jobvl ~jobvr ~a ->
  assert (jobvl = 'N' || jobvl = 'V');
  assert (jobvr = 'N' || jobvr = 'V');
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let vl = match jobvl with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let vr = match jobvr with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _vl = bigarray_start Ctypes_static.Genarray vl in
  let _vr = bigarray_start Ctypes_static.Genarray vr in
  let lda = Stdlib.max 1 (_stride a) in
  let ldvl = Stdlib.max 1 (_stride vl) in
  let ldvr = Stdlib.max 1 (_stride vr) in

  let wr = ref (Genarray.create _kind _layout [|0;0|]) in
  let wi = ref (Genarray.create _kind _layout [|0;0|]) in

  let ret = match _kind with
    | Float32   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let r = L.sgeev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~wr:_wr ~wi:_wi ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        wr := wr';
        wi := wi';
        r
      )
    | Float64   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let r = L.dgeev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~wr:_wr ~wi:_wi ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        wr := wr';
        wi := wi';
        r
      )
    | Complex32 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.cgeev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~w:_w ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        wr := w';
        wi := w';
        r
      )
    | Complex64 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.cgeev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~w:_w ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        wr := w';
        wi := w';
        r
      )
    | _         -> failwith "lapacke:geev"
  in
  check_lapack_error ret;
  !wr, !wi, vl, vr


let gesdd
  : type a b. ?jobz:char -> a:(a, b) t -> (a, b) t * (a, b) t *  (a, b) t
  = fun ?(jobz='A') ~a ->
  assert (jobz = 'A' || jobz = 'S' || jobz = 'O' || jobz = 'N');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  assert (m > 0 && n > 0);

  let s = ref (Genarray.create _kind _layout [|0;0|]) in
  let u = match jobz with
    | 'A' -> Genarray.create _kind _layout [|m;m|]
    | 'S' -> Genarray.create _kind _layout [|m;minmn|]
    | 'O' -> Genarray.create _kind _layout [|m; if m >=n then 0 else m|]
    | _   -> Genarray.create _kind _layout [|m;0|]
  in
  let vt = match jobz with
    | 'A' -> Genarray.create _kind _layout [|n;n|]
    | 'S' -> Genarray.create _kind _layout [|minmn;n|]
    | 'O' -> Genarray.create _kind _layout [|n; if m >=n then n else 0|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let lda = Stdlib.max 1 (_stride a) in
  let ldu = Stdlib.max 1 (_stride u) in
  let ldvt = Stdlib.max 1 (_stride vt) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _u = bigarray_start Ctypes_static.Genarray u in
  let _vt = bigarray_start Ctypes_static.Genarray vt in

  let ret = match _kind with
    | Float32   -> (
        let s' = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let r = L.sgesdd ~layout ~jobz ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt in
        s := s';
        r
      )
    | Float64   -> (
        let s' = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let r = L.dgesdd ~layout ~jobz ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt in
        s := s';
        r
      )
    | Complex32 -> (
        let s' = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let r = L.cgesdd ~layout ~jobz ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt in
        s := Owl_dense_matrix_generic.cast_s2c s';
        r
      )
    | Complex64 -> (
        let s' = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let r = L.zgesdd ~layout ~jobz ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt in
        s := Owl_dense_matrix_generic.cast_d2z s';
        r
      )
    | _        -> failwith "lapacke:gesdd"
  in
  check_lapack_error ret;

  match jobz, m >= n with
  | 'O', true  -> a, !s, vt
  | 'O', false -> u, !s, a
  | _, _       -> u, !s, vt


let gesvd
  : type a b. ?jobu:char -> ?jobvt:char -> a:(a, b) t -> (a, b) t * (a, b) t *  (a, b) t
  = fun ?(jobu='A') ?(jobvt='A') ~a ->
  assert (jobu = 'A' || jobu = 'S' || jobu = 'O' || jobu = 'N');
  assert (jobvt = 'A' || jobvt = 'S' || jobvt = 'O' || jobvt = 'N');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  assert (jobu <> 'O' || jobvt <> 'O');
  assert (m > 0 && n > 0);

  let s = ref (Genarray.create _kind _layout [|0;0|]) in
  let u = match jobu with
    | 'A' -> Genarray.create _kind _layout [|m;m|]
    | 'S' -> Genarray.create _kind _layout [|m;minmn|]
    | _   -> Genarray.create _kind _layout [|m;0|]
  in
  let vt = match jobvt with
    | 'A' -> Genarray.create _kind _layout [|n;n|]
    | 'S' -> Genarray.create _kind _layout [|minmn;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let lda = Stdlib.max 1 (_stride a) in
  let ldu = Stdlib.max 1 (_stride u) in
  let ldvt = Stdlib.max 1 (_stride vt) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _u = bigarray_start Ctypes_static.Genarray u in
  let _vt = bigarray_start Ctypes_static.Genarray vt in

  let ret = match _kind with
    | Float32   -> (
        let s' = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let superb = Genarray.create float32 _layout [|minmn - 1|]
          |> bigarray_start Ctypes_static.Genarray
        in
        let r = L.sgesvd ~layout ~jobu ~jobvt ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt ~superb in
        s := s';
        r
      )
    | Float64   -> (
        let s' = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let superb = Genarray.create float64 _layout [|minmn - 1|]
          |> bigarray_start Ctypes_static.Genarray
        in
        let r = L.dgesvd ~layout ~jobu ~jobvt ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt ~superb in
        s := s';
        r
      )
    | Complex32 -> (
        let s' = Genarray.create float32 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let superb = Genarray.create float32 _layout [|minmn - 1|]
          |> bigarray_start Ctypes_static.Genarray
        in
        let r = L.cgesvd ~layout ~jobu ~jobvt ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt ~superb in
        s := Owl_dense_matrix_generic.cast_s2c s';
        r
      )
    | Complex64 -> (
        let s' = Genarray.create float64 _layout [|1; minmn|] in
        let _s = bigarray_start Ctypes_static.Genarray s' in
        let superb = Genarray.create float64 _layout [|minmn - 1|]
          |> bigarray_start Ctypes_static.Genarray
        in
        let r = L.zgesvd ~layout ~jobu ~jobvt ~m ~n ~a:_a ~lda ~s:_s ~u:_u ~ldu ~vt:_vt ~ldvt ~superb in
        s := Owl_dense_matrix_generic.cast_d2z s';
        r
      )
    | _        -> failwith "lapacke:gesvd"
  in
  check_lapack_error ret;

  match jobu, jobvt with
  | 'O', _ -> a, !s, vt
  | _, 'O' -> u, !s, a
  | _, _   -> u, !s, vt


let ggsvd3
  : type a b. ?jobu:char -> ?jobv:char -> ?jobq:char -> a:(a, b) t -> b:(a, b) t
    -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * int * int * (a, b) t
  = fun ?(jobu='U') ?(jobv='V') ?(jobq='Q') ~a ~b ->
  assert (jobu = 'U' || jobu = 'N');
  assert (jobv = 'V' || jobu = 'N');
  assert (jobq = 'Q' || jobu = 'N');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let p = Owl_dense_matrix_generic.row_num b in
  assert (n = Owl_dense_matrix_generic.col_num b);
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ldu = Stdlib.max 1 m in
  let ldv = Stdlib.max 1 p in
  let ldq = Stdlib.max 1 n in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let alpha = ref (Genarray.create _kind _layout [|0;0|]) in
  let beta = ref (Genarray.create _kind _layout [|0;0|]) in
  let u = match jobu with
    | 'U' -> Genarray.create _kind _layout [|ldu; m|]
    | _   -> Genarray.create _kind _layout [|0; m|]
  in
  let v = match jobv with
    | 'V' -> Genarray.create _kind _layout [|ldv; p|]
    | _   -> Genarray.create _kind _layout [|0; p|]
  in
  let q = match jobq with
    | 'Q' -> Genarray.create _kind _layout [|ldq; n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let iwork = Genarray.create int32 _layout [|n|] in

  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _u = bigarray_start Ctypes_static.Genarray u in
  let _v = bigarray_start Ctypes_static.Genarray v in
  let _q = bigarray_start Ctypes_static.Genarray q in
  let _iwork = bigarray_start Ctypes_static.Genarray iwork in
  let _k = Ctypes.(allocate int32_t 0l) in
  let _l = Ctypes.(allocate int32_t 0l) in

  let ret = match _kind with
    | Float32   -> (
        let alpha' = Genarray.create float32 _layout [|1;n|] in
        let beta' = Genarray.create float32 _layout [|1;n|] in
        let _alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _beta = bigarray_start Ctypes_static.Genarray beta' in
        let r = L.sggsvd3 ~layout ~jobu ~jobv ~jobq ~m ~n ~p ~k:_k ~l:_l ~a:_a~lda ~b:_b ~ldb ~alpha:_alpha ~beta:_beta ~u:_u ~ldu ~v:_v ~ldv  ~q:_q ~ldq ~iwork:_iwork in
        alpha := alpha';
        beta := beta';
        r
      )
    | Float64   -> (
        let alpha' = Genarray.create float64 _layout [|1;n|] in
        let beta' = Genarray.create float64 _layout [|1;n|] in
        let _alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _beta = bigarray_start Ctypes_static.Genarray beta' in
        let r = L.dggsvd3 ~layout ~jobu ~jobv ~jobq ~m ~n ~p ~k:_k ~l:_l ~a:_a~lda ~b:_b ~ldb ~alpha:_alpha ~beta:_beta ~u:_u ~ldu ~v:_v ~ldv  ~q:_q ~ldq ~iwork:_iwork in
        alpha := alpha';
        beta := beta';
        r
      )
    | Complex32 -> (
        let alpha' = Genarray.create float32 _layout [|1;n|] in
        let beta' = Genarray.create float32 _layout [|1;n|] in
        let _alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _beta = bigarray_start Ctypes_static.Genarray beta' in
        let r = L.cggsvd3 ~layout ~jobu ~jobv ~jobq ~m ~n ~p ~k:_k ~l:_l ~a:_a~lda ~b:_b ~ldb ~alpha:_alpha ~beta:_beta ~u:_u ~ldu ~v:_v ~ldv ~q:_q ~ldq ~iwork:_iwork in
        alpha := Owl_dense_matrix_generic.cast_s2c alpha';
        beta := Owl_dense_matrix_generic.cast_s2c beta';
        r
      )
    | Complex64 -> (
        let alpha' = Genarray.create float64 _layout [|1;n|] in
        let beta' = Genarray.create float64 _layout [|1;n|] in
        let _alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _beta = bigarray_start Ctypes_static.Genarray beta' in
        let r = L.zggsvd3 ~layout ~jobu ~jobv ~jobq ~m ~n ~p ~k:_k ~l:_l ~a:_a~lda ~b:_b ~ldb ~alpha:_alpha ~beta:_beta ~u:_u ~ldu ~v:_v ~ldv ~q:_q ~ldq ~iwork:_iwork in
        alpha := Owl_dense_matrix_generic.cast_d2z alpha';
        beta := Owl_dense_matrix_generic.cast_d2z beta';
        r
      )
    | _         -> failwith "lapacke:ggsvd3"
  in
  check_lapack_error ret;

  (* construct R from a and b *)
  let k = Int32.to_int !@_k in
  let l = Int32.to_int !@_l in
  let r = match m - k -l >= 0 with
    | true  -> (
        let r = Owl_dense_matrix_generic.get_slice [[0; k + l - 1]; [n - k - l; n - 1]] a in
        Owl_dense_matrix_generic.triu r
      )
    | false -> (
        let ra = Owl_dense_matrix_generic.get_slice [[]; [n - k - l; n - 1]] a in
        let rb = Owl_dense_matrix_generic.get_slice [[m - k; l - 1]; [n - k - l; n - 1]] b in
        let r = Owl_dense_matrix_generic.concat_vertical ra rb in
        Owl_dense_matrix_generic.triu r
      )
  in
  u, v, q, !alpha, !beta, k, l, r


let geevx
  : type a b. balanc:char -> jobvl:char -> jobvr:char -> sense:char -> a:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * int * int * (a, b) t * float * (a, b) t * (a, b) t
  = fun ~balanc ~jobvl ~jobvr ~sense ~a ->
  assert (balanc = 'N' || balanc = 'P' || balanc = 'S' || balanc = 'B');
  assert (sense = 'N' || sense = 'E' || sense = 'V' || sense = 'B');
  assert (jobvl = 'N' || jobvl = 'V');
  assert (jobvr = 'N' || jobvr = 'V');
  if sense = 'E' || sense = 'B' then assert (jobvl = 'V' && jobvr = 'V');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let vl = match jobvl with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let vr = match jobvr with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let ldvl = Stdlib.max 1 (_stride vl) in
  let ldvr = Stdlib.max 1 (_stride vr) in

  let _ilo = Ctypes.(allocate int32_t 0l) in
  let _ihi = Ctypes.(allocate int32_t 0l) in
  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _vl = bigarray_start Ctypes_static.Genarray vl in
  let _vr = bigarray_start Ctypes_static.Genarray vr in

  let wr = ref (Genarray.create _kind _layout [|0;0|]) in
  let wi = ref (Genarray.create _kind _layout [|0;0|]) in
  let scale = ref (Genarray.create _kind _layout [|0;0|]) in
  let abnrm = ref 0. in
  let rconde = ref (Genarray.create _kind _layout [|0;0|]) in
  let rcondv = ref (Genarray.create _kind _layout [|0;0|]) in

  let ret = match _kind with
    | Float32   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let scale' = Genarray.create _kind _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let rconde' = Genarray.create _kind _layout [|1;n|] in
        let _rconde = bigarray_start Ctypes_static.Genarray rconde' in
        let rcondv' = Genarray.create _kind _layout [|1;n|] in
        let _rcondv = bigarray_start Ctypes_static.Genarray rcondv' in
        let _abnrm = Ctypes.(allocate float 0.) in
        let r = L.sgeevx ~layout ~balanc ~jobvl ~jobvr ~sense ~n ~a:_a ~lda ~wr:_wr ~wi:_wi ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~ilo:_ilo ~ihi:_ihi ~scale:_scale ~abnrm:_abnrm ~rconde:_rconde ~rcondv:_rcondv in
        wr := wr';
        wi := wi';
        scale := scale';
        abnrm := !@_abnrm;
        rconde := rconde';
        rcondv := rcondv';
        r
      )
    | Float64   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let scale' = Genarray.create _kind _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let rconde' = Genarray.create _kind _layout [|1;n|] in
        let _rconde = bigarray_start Ctypes_static.Genarray rconde' in
        let rcondv' = Genarray.create _kind _layout [|1;n|] in
        let _rcondv = bigarray_start Ctypes_static.Genarray rcondv' in
        let _abnrm = Ctypes.(allocate double 0.) in
        let r = L.dgeevx ~layout ~balanc ~jobvl ~jobvr ~sense ~n ~a:_a ~lda ~wr:_wr ~wi:_wi ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~ilo:_ilo ~ihi:_ihi ~scale:_scale ~abnrm:_abnrm ~rconde:_rconde ~rcondv:_rcondv in
        wr := wr';
        wi := wi';
        scale := scale';
        abnrm := !@_abnrm;
        rconde := rconde';
        rcondv := rcondv';
        r
      )
    | Complex32 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let scale' = Genarray.create float32 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let rconde' = Genarray.create float32 _layout [|1;n|] in
        let _rconde = bigarray_start Ctypes_static.Genarray rconde' in
        let rcondv' = Genarray.create float32 _layout [|1;n|] in
        let _rcondv = bigarray_start Ctypes_static.Genarray rcondv' in
        let _abnrm = Ctypes.(allocate float 0.) in
        let r = L.cgeevx ~layout ~balanc ~jobvl ~jobvr ~sense ~n ~a:_a ~lda ~w:_w ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~ilo:_ilo ~ihi:_ihi ~scale:_scale ~abnrm:_abnrm ~rconde:_rconde ~rcondv:_rcondv in
        wr := w';
        wi := w';
        scale := Owl_dense_matrix_generic.cast_s2c scale';
        abnrm := !@_abnrm;
        rconde := Owl_dense_matrix_generic.cast_s2c rconde';
        rcondv := Owl_dense_matrix_generic.cast_s2c rcondv';
        r
      )
    | Complex64 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let scale' = Genarray.create float64 _layout [|1;n|] in
        let _scale = bigarray_start Ctypes_static.Genarray scale' in
        let rconde' = Genarray.create float64 _layout [|1;n|] in
        let _rconde = bigarray_start Ctypes_static.Genarray rconde' in
        let rcondv' = Genarray.create float64 _layout [|1;n|] in
        let _rcondv = bigarray_start Ctypes_static.Genarray rcondv' in
        let _abnrm = Ctypes.(allocate double 0.) in
        let r = L.zgeevx ~layout ~balanc ~jobvl ~jobvr ~sense ~n ~a:_a ~lda ~w:_w ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~ilo:_ilo ~ihi:_ihi ~scale:_scale ~abnrm:_abnrm ~rconde:_rconde ~rcondv:_rcondv in
       wr := w';
        wi := w';
        scale := Owl_dense_matrix_generic.cast_d2z scale';
        abnrm := !@_abnrm;
        rconde := Owl_dense_matrix_generic.cast_d2z rconde';
        rcondv := Owl_dense_matrix_generic.cast_d2z rcondv';
        r
      )
    | _         -> failwith "lapacke:geevx"
  in
  check_lapack_error ret;
  (* return all the results ~modified in-place *)
  let ilo = Int32.to_int !@_ilo in
  let ihi = Int32.to_int !@_ihi in
  a, !wr, !wi, vl, vr, ilo, ihi, !scale, !abnrm, !rconde, !rcondv


let ggev
  : type a b. jobvl:char -> jobvr:char -> a:(a, b) t -> b:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~jobvl ~jobvr ~a ~b ->
  assert (jobvl = 'N' || jobvl = 'V');
  assert (jobvr = 'N' || jobvr = 'V');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (m = n && mb = n && mb = nb);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ldvl = match jobvl with
    | 'V' -> n
    | _   -> 1
  in
  let ldvr = match jobvr with
    | 'V' -> n
    | _   -> 1
  in
  let vl = Genarray.create _kind _layout [|n;ldvl|] in
  let vr = Genarray.create _kind _layout [|n;ldvr|] in

  let alphar = ref (Genarray.create _kind _layout [|0;0|]) in
  let alphai = ref (Genarray.create _kind _layout [|0;0|]) in
  let beta = Genarray.create _kind _layout [|1;n|] in

  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _vl = bigarray_start Ctypes_static.Genarray vl in
  let _vr = bigarray_start Ctypes_static.Genarray vr in
  let _beta = bigarray_start Ctypes_static.Genarray beta in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let r = L.sggev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~b:_b ~ldb ~alphar:_alphar ~alphai:_alphai ~beta:_beta ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Float64   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let r = L.dggev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~b:_b ~ldb ~alphar:_alphar ~alphai:_alphai ~beta:_beta ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Complex32 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let r = L.cggev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~b:_b ~ldb ~alpha ~beta:_beta ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | Complex64 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let r = L.cggev ~layout ~jobvl ~jobvr ~n ~a:_a ~lda ~b:_b ~ldb ~alpha ~beta:_beta ~vl:_vl ~ldvl ~vr:_vr ~ldvr in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | _         -> failwith "lapacke:ggev"
  in
  check_lapack_error ret;
  (* note alphar and alphai are the same for complex flavour *)
  !alphar, !alphai, beta, vl, vr


let gtsv
  : type a b. dl:(a, b) t -> d:(a, b) t -> du:(a, b) t -> b:(a, b) t -> (a, b) t
  = fun ~dl ~d ~du ~b ->
  let n = Owl_dense_matrix_generic.numel d in
  let n_dl = Owl_dense_matrix_generic.numel dl in
  let n_du = Owl_dense_matrix_generic.numel du in
  assert (n_dl = n || n_dl = n - 1);
  assert (n_du = n || n_du = n - 1);
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (mb = n);
  let _kind = Genarray.kind b in
  let _layout = Genarray.layout b in
  let layout = lapacke_layout _layout in

  let ldb = Stdlib.max 1 (_stride b) in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _dl = bigarray_start Ctypes_static.Genarray dl in
  let _du = bigarray_start Ctypes_static.Genarray du in

  let ret = match _kind with
    | Float32   -> L.sgtsv ~layout ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~b:_b ~ldb
    | Float64   -> L.dgtsv ~layout ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~b:_b ~ldb
    | Complex32 -> L.cgtsv  ~layout ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~b:_b ~ldb
    | Complex64 -> L.zgtsv  ~layout ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~b:_b ~ldb
    | _         -> failwith "lapacke:gtsv"
  in
  check_lapack_error ret;
  b


let gttrf
  : type a b. dl:(a, b) t -> d:(a, b) t -> du:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (int32, int32_elt) t
  = fun ~dl ~d ~du ->
  let n = Owl_dense_matrix_generic.numel d in
  let n_dl = Owl_dense_matrix_generic.numel dl in
  let n_du = Owl_dense_matrix_generic.numel du in
  assert (n_dl = n - 1 && n_du = n - 1);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in

  let du2 = Genarray.create _kind _layout [|1; n - 2|] in
  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _dl = bigarray_start Ctypes_static.Genarray dl in
  let _du = bigarray_start Ctypes_static.Genarray du in
  let _du2 = bigarray_start Ctypes_static.Genarray du2 in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in

  let ret = match _kind with
    | Float32   -> L.sgttrf ~n ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv
    | Float64   -> L.dgttrf ~n ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv
    | Complex32 -> L.cgttrf ~n ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv
    | Complex64 -> L.zgttrf ~n ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv

    | _         -> failwith "lapacke:gttrf"
  in
  check_lapack_error ret;
  dl, d, du, du2, ipiv


let gttrs
  : type a b. trans:char -> dl:(a, b) t -> d:(a, b) t -> du:(a, b) t
  -> du2:(a, b) t -> ipiv:(int32, int32_elt) t -> b:(a, b) t -> (a, b) t
  = fun ~trans ~dl ~d ~du ~du2 ~ipiv:_ ~b ->
  assert (trans = 'N' || trans = 'T' || trans = 'C');

  let n = Owl_dense_matrix_generic.numel d in
  let n_dl = Owl_dense_matrix_generic.numel dl in
  let n_du = Owl_dense_matrix_generic.numel du in
  assert (n_dl = n - 1 && n_du = n - 1);
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (mb = n);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let ldb = Stdlib.max 1 (_stride b) in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _dl = bigarray_start Ctypes_static.Genarray dl in
  let _du = bigarray_start Ctypes_static.Genarray du in
  let _du2 = bigarray_start Ctypes_static.Genarray du2 in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in

  let ret = match _kind with
    | Float32 ->
      L.sgttrs ~layout ~trans ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv ~b:_b ~ldb
    | Float64 ->
      L.dgttrs ~layout ~trans ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv ~b:_b ~ldb
    | Complex32 ->
      L.cgttrs ~layout ~trans ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv ~b:_b ~ldb
     | Complex64 ->
      L.zgttrs ~layout ~trans ~n ~nrhs ~dl:_dl ~d:_d ~du:_du ~du2:_du2 ~ipiv:_ipiv ~b:_b ~ldb
     | _         -> failwith "lapacke:gttrs"
  in
  check_lapack_error ret;
  b


let orglq
  : type a. ?k:int -> a:(float, a) t -> tau:(float, a) t -> (float, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sorglq ~layout ~m:minmn ~n ~k ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dorglq ~layout ~m:minmn ~n ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading rows if necessary *)
  match minmn < m with
  | true  -> Owl_dense_matrix_generic.get_slice [[0;minmn-1]; []] a
  | false -> a


let unglq
  : type a. ?k:int -> a:(Complex.t, a) t -> tau:(Complex.t, a) t -> (Complex.t, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Complex32 -> L.cunglq ~layout ~m:minmn ~n ~k ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zunglq ~layout ~m:minmn ~n ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading rows if necessary *)
  match minmn < m with
  | true  -> Owl_dense_matrix_generic.get_slice [[0;minmn-1]; []] a
  | false -> a


let orgqr
  : type a. ?k:int -> a:(float, a) t -> tau:(float, a) t -> (float, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sorgqr ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dorgqr ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading columns if necessary *)
  match minmn < n with
  | true  -> Owl_dense_matrix_generic.get_slice [[]; [0;minmn-1]] a
  | false -> a


let ungqr
  : type a. ?k:int -> a:(Complex.t, a) t -> tau:(Complex.t, a) t -> (Complex.t, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Complex32   -> L.cungqr ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
    | Complex64   -> L.zungqr ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading columns if necessary *)
  match minmn < n with
  | true  -> Owl_dense_matrix_generic.get_slice [[]; [0;minmn-1]] a
  | false -> a


let orgql
  : type a. ?k:int -> a:(float, a) t -> tau:(float, a) t -> (float, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sorgql ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dorgql ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading columns if necessary *)
  match minmn < n with
  | true  -> Owl_dense_matrix_generic.get_slice [[]; [0;minmn-1]] a
  | false -> a


let orgrq
  : type a. ?k:int -> a:(float, a) t -> tau:(float, a) t -> (float, a) t
  = fun ?k ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let minmn = Stdlib.min m n in
  let k = match k with
    | Some k -> k
    | None   -> Owl_dense_matrix_generic.numel tau
  in
  assert (k <= minmn);
  assert (k <= Owl_dense_matrix_generic.numel tau);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sorgrq ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dorgrq ~layout ~m ~n:minmn ~k ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  (* extract the first leading columns if necessary *)
  match minmn < n with
  | true  -> Owl_dense_matrix_generic.get_slice [[]; [0;minmn-1]] a
  | false -> a


let ormlq
  : type a. side:char -> trans:char -> a:(float, a) t -> tau:(float, a) t
  -> c:(float, a) t -> (float, a) t
  = fun ~side ~trans ~a ~tau ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let ma = Owl_dense_matrix_generic.row_num a in
  let na = Owl_dense_matrix_generic.col_num a in
  let k = Owl_dense_matrix_generic.numel tau in
  if side = 'L' then
    assert (ma = k && na = m && k <= m)
  else (* if side = 'R' *)
    assert (ma = k && na = n && k <= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sormlq ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
    | Float64   -> L.dormlq ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
  in
  check_lapack_error ret;
  c


let ormqr
  : type a. side:char -> trans:char -> a:(float, a) t -> tau:(float, a) t
  -> c:(float, a) t -> (float, a) t
  = fun ~side ~trans ~a ~tau ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let ma = Owl_dense_matrix_generic.row_num a in
  let na = Owl_dense_matrix_generic.col_num a in
  let k = Owl_dense_matrix_generic.numel tau in
  if side = 'L' then
    assert (ma = m && na = k && k <= m)
  else (* if side = 'R' *)
    assert (ma = n && na = k && k <= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sormqr ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
    | Float64   -> L.dormqr ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
  in
  check_lapack_error ret;
  c


let ormql
  : type a. side:char -> trans:char -> a:(float, a) t -> tau:(float, a) t
  -> c:(float, a) t -> (float, a) t
  = fun ~side ~trans ~a ~tau ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let ma = Owl_dense_matrix_generic.row_num a in
  let na = Owl_dense_matrix_generic.col_num a in
  let k = Owl_dense_matrix_generic.numel tau in
  if side = 'L' then
    assert (ma = m && na = k && k <= m)
  else (* if side = 'R' *)
    assert (ma = n && na = k && k <= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sormql ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
    | Float64   -> L.dormql ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
  in
  check_lapack_error ret;
  c


let ormrq
  : type a. side:char -> trans:char -> a:(float, a) t -> tau:(float, a) t
  -> c:(float, a) t -> (float, a) t
  = fun ~side ~trans ~a ~tau ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let ma = Owl_dense_matrix_generic.row_num a in
  let na = Owl_dense_matrix_generic.col_num a in
  let k = Owl_dense_matrix_generic.numel tau in
  if side = 'L' then
    assert (ma = k && na = m && k <= m)
  else (* if side = 'R' *)
    assert (ma = k && na = n && k <= n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let _tau = bigarray_start Ctypes_static.Genarray tau in

  let ret = match _kind with
    | Float32   -> L.sormrq ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
    | Float64   -> L.dormrq ~layout ~side ~trans ~m ~n ~k ~a:_a ~lda ~tau:_tau ~c:_c ~ldc
  in
  check_lapack_error ret;
  c


let gemqrt
  : type a b. side:char -> trans:char -> v:(a, b) t -> t:(a, b) t -> c:(a, b) t -> (a, b) t
  = fun ~side ~trans ~v ~t ~c ->
  assert (side = 'L' || side = 'R');
  assert (trans = 'N' || trans = 'T' || trans = 'C');

  let m = Owl_dense_matrix_generic.row_num c in
  let n = Owl_dense_matrix_generic.col_num c in
  let nb = Owl_dense_matrix_generic.row_num t in
  let k = Owl_dense_matrix_generic.col_num t in
  let mv = Owl_dense_matrix_generic.row_num v in
  let ldv = Stdlib.max 1 (_stride v) in
  assert (k >= nb);
  if side = 'L' then
    assert (mv = m && ldv = k && k <= m)
  else (* if side = 'R' *)
    assert (mv = n && ldv = k && k <= n);
  let _kind = Genarray.kind c in
  let _layout = Genarray.layout c in
  let layout = lapacke_layout _layout in

  let ldt = Stdlib.max 1 (_stride t) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _v = bigarray_start Ctypes_static.Genarray v in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let _c = bigarray_start Ctypes_static.Genarray c in

  let ret = match _kind with
    | Float32   -> L.sgemqrt ~layout ~side ~trans ~m ~n ~k ~nb  ~v:_v ~ldv  ~t:_t ~ldt  ~c:_c ~ldc
    | Float64   -> L.dgemqrt ~layout ~side ~trans ~m ~n ~k ~nb  ~v:_v ~ldv  ~t:_t ~ldt  ~c:_c ~ldc
    | Complex32 -> L.cgemqrt ~layout ~side ~trans ~m ~n ~k ~nb  ~v:_v ~ldv  ~t:_t ~ldt  ~c:_c ~ldc
    | Complex64 -> L.zgemqrt ~layout ~side ~trans ~m ~n ~k ~nb  ~v:_v ~ldv  ~t:_t ~ldt  ~c:_c ~ldc
    | _         -> failwith "lapacke:gemqrt"
  in
  check_lapack_error ret;
  c


let posv
  : type a b. uplo:char -> a:(a, b) t -> b:(a, b) t -> (a, b) t * (a, b) t
  = fun ~uplo ~a ~b ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let nrhs = _stride b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in

  let ret = match _kind with
    | Float32   -> L.sposv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Float64   -> L.dposv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex32 -> L.cposv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex64 -> L.zposv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | _         -> failwith "lapacke:posv"
  in
  check_lapack_error ret;
  a, b


let potrf
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in

  let ret = match _kind with
    | Float32   -> L.spotrf ~layout ~uplo ~n ~a:_a ~lda
    | Float64   -> L.dpotrf ~layout ~uplo ~n ~a:_a ~lda
    | Complex32 -> L.cpotrf ~layout ~uplo ~n ~a:_a ~lda
    | Complex64 -> L.zpotrf ~layout ~uplo ~n ~a:_a ~lda
    | _         -> failwith "lapacke:potrf"
  in
  check_lapack_error ret;
  a


let potri
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in

  let ret = match _kind with
    | Float32   -> L.spotri ~layout ~uplo ~n ~a:_a ~lda
    | Float64   -> L.dpotri ~layout ~uplo ~n ~a:_a ~lda
    | Complex32 -> L.cpotri ~layout ~uplo ~n ~a:_a ~lda
    | Complex64 -> L.zpotri ~layout ~uplo ~n ~a:_a ~lda
    | _         -> failwith "lapacke:potri"
  in
  check_lapack_error ret;
  a


let potrs
  : type a b. uplo:char -> a:(a, b) t -> b:(a, b) t -> (a, b) t
  = fun ~uplo ~a ~b ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let nrhs = _stride b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in

  let ret = match _kind with
    | Float32   -> L.spotrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Float64   -> L.dpotrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex32 -> L.cpotrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex64 -> L.zpotrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | _         -> failwith "lapacke:potrs"
  in
  check_lapack_error ret;
  b


let pstrf
  : type a b. uplo:char -> a:(a, b) t -> tol:a -> (a, b) t * (int32, int32_elt) t * int * int
  = fun ~uplo ~a ~tol ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let piv = Genarray.create int32 _layout [|1;n|] in
  let _piv = bigarray_start Ctypes_static.Genarray piv in
  let _rank = Ctypes.(allocate int32_t 0l) in

  let ret = match _kind with
    | Float32   -> L.spstrf ~layout ~uplo ~n ~a:_a ~lda ~piv:_piv ~rank:_rank ~tol
    | Float64   -> L.dpstrf ~layout ~uplo ~n ~a:_a ~lda ~piv:_piv ~rank:_rank ~tol
    | Complex32 -> L.cpstrf ~layout ~uplo ~n ~a:_a ~lda ~piv:_piv ~rank:_rank ~tol:Complex.(tol.re)
    | Complex64 -> L.zpstrf ~layout ~uplo ~n ~a:_a ~lda ~piv:_piv ~rank:_rank ~tol:Complex.(tol.re)
    | _         -> failwith "lapacke:pstrf"
  in
  check_lapack_error ret;
  let rank = Int32.to_int !@_rank in
  a, piv, rank, ret


let ptsv
  : type a b. d:(a, b) t -> e:(a, b) t -> b:(a, b) t -> (a, b) t
  = fun ~d ~e ~b ->
  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (n_e = n - 1);
  assert (mb = n);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let ldb = Stdlib.max 1 (_stride b) in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _b = bigarray_start Ctypes_static.Genarray b in

  (* NOTE: only use the real part of d *)
  let ret = match _kind with
    | Float32   -> (
        let _d = bigarray_start Ctypes_static.Genarray d in
        L.sptsv ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Float64   -> (
        let _d = bigarray_start Ctypes_static.Genarray d in
        L.dptsv ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Complex32 -> (
        let d' = Owl_dense_matrix_c.re d in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        L.cptsv ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Complex64 -> (
        let d' = Owl_dense_matrix_z.re d in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        L.zptsv ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | _         -> failwith "lapacke:ptsv"
  in
  check_lapack_error ret;
  b


let pttrf
  : type a b. d:(a, b) t -> e:(a, b) t -> (a, b) t * (a, b) t
  = fun ~d ~e ->
  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let _kind = Genarray.kind d in
  let _e = bigarray_start Ctypes_static.Genarray e in

  (* NOTE: only use the real part of d *)
  let ret = match _kind with
    | Float32   -> (
        let d = bigarray_start Ctypes_static.Genarray d in
        L.spttrf ~n ~d ~e:_e
      )
    | Float64   -> (
        let d = bigarray_start Ctypes_static.Genarray d in
        L.dpttrf ~n ~d ~e:_e
      )
    | Complex32 -> (
        let d' = Owl_dense_matrix_c.re d in
        let d = bigarray_start Ctypes_static.Genarray d' in
        L.cpttrf ~n ~d ~e:_e
      )
    | Complex64 -> (
        let d' = Owl_dense_matrix_z.re d in
        let d = bigarray_start Ctypes_static.Genarray d' in
        L.zpttrf ~n ~d ~e:_e
      )
    | _         -> failwith "lapacke:pttrf"
  in
  check_lapack_error ret;
  d, e


let pttrs
: type a b. ?uplo:char -> d:(a, b) t -> e:(a, b) t -> b:(a, b) t -> (a, b) t
= fun ?uplo ~d ~e ~b ->
  (* NOTE: ~uplo is only for complex flavour *)
  let uplo = match uplo with
    | Some uplo -> uplo
    | None      -> 'U'
  in
  assert (uplo = 'U' || uplo = 'L');

  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (n_e = n - 1);
  assert (mb = n);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let ldb = Stdlib.max 1 (_stride b) in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _b = bigarray_start Ctypes_static.Genarray b in

  (* NOTE: only use the real part of d *)
  let ret = match _kind with
    | Float32   -> (
        let _d = bigarray_start Ctypes_static.Genarray d in
        L.spttrs ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Float64   -> (
        let _d = bigarray_start Ctypes_static.Genarray d in
        L.dpttrs ~layout ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Complex32 -> (
        let d' = Owl_dense_matrix_c.re d in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        L.cpttrs ~layout ~uplo ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | Complex64 -> (
        let d' = Owl_dense_matrix_z.re d in
        let _d = bigarray_start Ctypes_static.Genarray d' in
        L.zpttrs ~layout ~uplo ~n ~nrhs ~d:_d ~e:_e ~b:_b ~ldb
      )
    | _         -> failwith "lapacke:pttrs"
  in
  check_lapack_error ret;
  b


let trtri
  : type a b. uplo:char -> diag:char -> a:(a, b) t -> (a, b) t
  = fun ~uplo ~diag ~a ->
  assert (uplo = 'U' || uplo = 'L');
  assert (diag = 'N' || diag = 'U');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in

  let ret = match _kind with
    | Float32   -> L.strtri ~layout ~uplo ~diag ~n ~a:_a ~lda
    | Float64   -> L.dtrtri ~layout ~uplo ~diag ~n ~a:_a ~lda
    | Complex32 -> L.ctrtri ~layout ~uplo ~diag ~n ~a:_a ~lda
    | Complex64 -> L.ztrtri ~layout ~uplo ~diag ~n ~a:_a ~lda
    | _         -> failwith "lapacke:trtri"
  in
  check_lapack_error ret;
  a


let trtrs
  : type a b. uplo:char -> trans:char -> diag:char -> a:(a, b) t -> b:(a, b) t -> (a, b) t
  = fun ~uplo ~trans ~diag ~a ~b ->
  assert (uplo = 'U' || uplo = 'L');
  assert (diag = 'N' || diag = 'U');
  assert (trans = 'N' || trans = 'T' || trans = 'C');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let mb = Owl_dense_matrix_generic.row_num b in
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (mb = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in

  let ret = match _kind with
    | Float32   -> L.strtrs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Float64   -> L.dtrtrs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex32 -> L.ctrtrs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | Complex64 -> L.ztrtrs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb
    | _         -> failwith "lapacke:trtrs"
  in
  check_lapack_error ret;
  b


let trcon
  : type a b. norm:char -> uplo:char -> diag:char -> a:(a, b) t -> float
  = fun ~norm ~uplo ~diag ~a ->
  assert (uplo = 'U' || uplo = 'L');
  assert (diag = 'N' || diag = 'U');
  assert (norm = '1' || norm = 'O' || norm = 'I');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let rcond = ref 0. in

  let ret = match _kind with
    | Float32   -> (
        let _rcond = Ctypes.(allocate float 0.) in
        let r = L.strcon ~layout ~norm ~uplo ~diag ~n ~a:_a ~lda ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Float64   -> (
        let _rcond = Ctypes.(allocate double 0.) in
        let r = L.dtrcon ~layout ~norm ~uplo ~diag ~n ~a:_a ~lda ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Complex32 -> (
        let _rcond = Ctypes.(allocate float 0.) in
        let r = L.ctrcon ~layout ~norm ~uplo ~diag ~n ~a:_a ~lda ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Complex64 -> (
        let _rcond = Ctypes.(allocate double 0.) in
        let r = L.ztrcon ~layout ~norm ~uplo ~diag ~n ~a:_a ~lda ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | _         -> failwith "lapacke:trcon"
  in
  check_lapack_error ret;
  !rcond


let trevc
  : type a b. side:char -> howmny:char -> select:(int32, int32_elt) t -> t:(a, b) t
  -> (int32, int32_elt) t * (a, b) t * (a, b) t
  = fun ~side ~howmny ~select ~t ->
  assert (side = 'L' || side = 'R' || side = 'B');
  assert (howmny = 'A' || howmny = 'B' || howmny = 'S');

  let mt = Owl_dense_matrix_generic.row_num t in
  let n = Owl_dense_matrix_generic.col_num t in
  assert (mt = n);
  let _kind = Genarray.kind t in
  let _layout = Genarray.layout t in
  let layout = lapacke_layout _layout in

  (* NOTE: I ~might allocate too ~much ~memory for vl and vr, please refer to Intel
    MKL documentation for ~more detailed ~memory allocation strategy. Fix later.
    url: https://software.intel.com/en-us/mkl-developer-reference-c-trevc
  *)
  let vl = Genarray.create _kind _layout [|n;n|] in
  let vr = Genarray.create _kind _layout [|n;n|] in
  let mm = Owl_dense_matrix_generic.col_num vl in
  let ldt = _stride t in
  let ldvl = _stride vl in
  let ldvr = _stride vr in
  let _m = Ctypes.(allocate int32_t 0l) in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let _vl = bigarray_start Ctypes_static.Genarray vl in
  let _vr = bigarray_start Ctypes_static.Genarray vr in
  let _select = bigarray_start Ctypes_static.Genarray select in

  let ret = match _kind with
    | Float32   -> L.strevc ~layout ~side ~howmny ~select:_select ~n ~t:_t ~ldt  ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~mm ~m:_m
    | Float64   -> L.dtrevc ~layout ~side ~howmny ~select:_select ~n ~t:_t ~ldt  ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~mm ~m:_m
    | Complex32 -> L.ctrevc ~layout ~side ~howmny ~select:_select ~n ~t:_t ~ldt  ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~mm ~m:_m
    | Complex64 -> L.ztrevc ~layout ~side ~howmny ~select:_select ~n ~t:_t ~ldt  ~vl:_vl ~ldvl ~vr:_vr ~ldvr ~mm ~m:_m
    | _         -> failwith "lapacke:trevc"
  in
  check_lapack_error ret;

  let m = Int32.to_int !@_m in
  let _empty = Genarray.create _kind _layout [|0;0|] in
  if howmny = 'S' then (      (* return selected eigenvectors *)
    if side = 'L' then        (* left eigenvectors only *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vl, _empty
    else if side = 'R' then   (* right eigenvectors only *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vr, _empty
    else                      (* both eigenvectors *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vl, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vr
  )
  else (                      (* return all eigenvectors *)
    if side = 'L' then        (* left eigenvectors only *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vl, _empty
    else if side = 'R' then   (* right eigenvectors only *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vr, _empty
    else                      (* both eigenvectors *)
      select, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vl, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] vr
  )


let trrfs
  : type a b. uplo:char -> trans:char -> diag:char -> a:(a, b) t -> b:(a, b) t
  -> x:(a, b) t -> (a, b) t * (a, b) t
  = fun ~uplo ~trans ~diag ~a ~b ~x ->
  assert (uplo = 'U' || uplo = 'L');
  assert (diag = 'N' || diag = 'U');
  assert (trans = 'N' || trans = 'T' || trans = 'C');

  let n = Owl_dense_matrix_generic.col_num a in
  assert (Owl_dense_matrix_generic.row_num b = n);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  assert (Owl_dense_matrix_generic.col_num x = nrhs);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ldx = Stdlib.max 1 (_stride x) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _x = bigarray_start Ctypes_static.Genarray x in

  let ferr = ref (Genarray.create _kind _layout [|0;0|]) in
  let berr = ref (Genarray.create _kind _layout [|0;0|]) in

  let ret = match _kind with
    | Float32   -> (
        let ferr' = Genarray.create _kind _layout [|1;nrhs|] in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let berr' = Genarray.create _kind _layout [|1;nrhs|] in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let r = L.strrfs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~x:_x ~ldx ~ferr:_ferr ~berr:_berr in
        ferr := ferr';
        berr := berr';
        r
      )
    | Float64   -> (
        let ferr' = Genarray.create _kind _layout [|1;nrhs|] in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let berr' = Genarray.create _kind _layout [|1;nrhs|] in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let r = L.dtrrfs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~x:_x ~ldx ~ferr:_ferr ~berr:_berr in
        ferr := ferr';
        berr := berr';
        r
      )
    | Complex32 -> (
        let ferr' = Genarray.create float32 _layout [|1;nrhs|] in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let berr' = Genarray.create float32 _layout [|1;nrhs|] in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let r = L.ctrrfs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~x:_x ~ldx ~ferr:_ferr ~berr:_berr in
        ferr := Owl_dense_matrix_generic.cast_s2c ferr';
        berr := Owl_dense_matrix_generic.cast_s2c berr';
        r
      )
    | Complex64 -> (
        let ferr' = Genarray.create float64 _layout [|1;nrhs|] in
        let _ferr = bigarray_start Ctypes_static.Genarray ferr' in
        let berr' = Genarray.create float64 _layout [|1;nrhs|] in
        let _berr = bigarray_start Ctypes_static.Genarray berr' in
        let r = L.ctrrfs ~layout ~uplo ~trans ~diag ~n ~nrhs ~a:_a ~lda ~b:_b ~ldb ~x:_x ~ldx ~ferr:_ferr ~berr:_berr in
        ferr := Owl_dense_matrix_generic.cast_d2z ferr';
        berr := Owl_dense_matrix_generic.cast_d2z berr';
        r
      )
    | _         -> failwith "lapacke:trrfs"
  in
  check_lapack_error ret;
  !ferr, !berr


let stev
  : type a. jobz:char -> d:(float, a) t -> e:(float, a) t
  -> (float, a) t * (float, a) t
  = fun ~jobz ~d ~e ->
  assert (jobz = 'N' && jobz = 'V');

  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let z = match jobz with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let ldz = Stdlib.max 1 (_stride z) in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _z = bigarray_start Ctypes_static.Genarray z in

  let ret = match _kind with
    | Float32   -> L.sstev ~layout ~jobz ~n ~d:_d ~e:_e ~z:_z ~ldz
    | Float64   -> L.dstev ~layout ~jobz ~n ~d:_d ~e:_e ~z:_z ~ldz
  in
  check_lapack_error ret;
  d, z


let stebz
  : type a. range:char -> order:char -> vl:float -> vu:float -> il:int -> iu:int
  -> abstol:float -> d:(float, a) t -> e:(float, a) t
  -> (float, a) t * (int32, int32_elt) t * (int32, int32_elt) t
  = fun ~range ~order ~vl ~vu ~il ~iu ~abstol ~d ~e ->
  assert (range = 'A' || range = 'V' && range = 'I');
  assert (order = 'B' || order = 'E');

  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in

  let _m = Ctypes.(allocate int32_t 0l) in
  let _nsplit = Ctypes.(allocate int32_t 0l) in
  let w = Genarray.create _kind _layout [|1;n|] in
  let iblock = Genarray.create int32 _layout [|1;n|] in
  let isplit = Genarray.create int32 _layout [|1;n|] in

  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _w = bigarray_start Ctypes_static.Genarray w in
  let _iblock = bigarray_start Ctypes_static.Genarray iblock in
  let _isplit = bigarray_start Ctypes_static.Genarray isplit in

  let ret = match _kind with
    | Float32   ->
      L.sstebz ~range ~order ~n ~vl ~vu ~il ~iu ~abstol ~d:_d ~e:_e ~m:_m
        ~nsplit:_nsplit ~w:_w ~iblock:_iblock ~isplit:_isplit
    | Float64   ->
      L.dstebz ~range ~order ~n ~vl ~vu ~il ~iu ~abstol ~d:_d ~e:_e ~m:_m
        ~nsplit:_nsplit ~w:_w ~iblock:_iblock ~isplit:_isplit
  in
  check_lapack_error ret;
  let m = Int32.to_int !@_m in
  let w = Owl_dense_matrix_generic.resize w [|1; m|] in
  let iblock = Owl_dense_matrix_generic.resize iblock [|1; m|] in
  let isplit = Owl_dense_matrix_generic.resize isplit [|1; m|] in
  w, iblock, isplit


let stegr
  : type a b. kind:(a, b) kind -> jobz:char -> range:char -> d:(float, b) t
  -> e:(float, b) t -> vl:float -> vu:float -> il:int -> iu:int -> (a, b) t *(a, b) t
  = fun ~kind ~jobz ~range ~d ~e ~vl ~vu ~il ~iu ->
  assert (jobz = 'N' && jobz = 'V');
  assert (range = 'A' || range = 'V' && range = 'I');

  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let _kind = kind in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let abstol = 1. in  (* note that ~abstol is unused. *)
  let e = Owl_dense_matrix_generic.resize e [|1; n|] in
  let ldz = match jobz with 'N' -> 1 | _ -> n in
  let z = match range with
    | 'I' -> Genarray.create _kind _layout [|(iu - il + 1); ldz|]
    | _   -> Genarray.create _kind _layout [|n; ldz|]
  in
  let isuppz = Genarray.create int32 _layout [|1; (Owl_dense_matrix_generic.row_num z)|] in

  let w = ref (Genarray.create _kind _layout [|0;0|]) in
  let _m = Ctypes.(allocate int32_t 0l) in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _z = bigarray_start Ctypes_static.Genarray z in
  let _isuppz = bigarray_start Ctypes_static.Genarray isuppz in

  let ret = match _kind with
    | Float32   -> (
        let w' = Genarray.create float32 _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.sstegr ~layout ~jobz ~range ~n ~d:_d ~e:_e ~vl ~vu ~il ~iu
            ~abstol ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz in
        w := w';
        r
      )
    | Float64   -> (
        let w' = Genarray.create float64 _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.dstegr ~layout ~jobz ~range ~n ~d:_d ~e:_e ~vl ~vu ~il ~iu
            ~abstol ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz in
        w := w';
        r
      )
    | Complex32 -> (
        let w' = Genarray.create float32 _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r =  L.cstegr ~layout ~jobz ~range ~n ~d:_d ~e:_e ~vl ~vu ~il ~iu
            ~abstol ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz in
        w := Owl_dense_matrix_generic.cast_s2c w';
        r
      )
    | Complex64 -> (
        let w' = Genarray.create float64 _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.zstegr ~layout ~jobz ~range ~n ~d:_d ~e:_e ~vl ~vu ~il ~iu
            ~abstol ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz in
        w := Owl_dense_matrix_generic.cast_d2z w';
        r
      )
    | _         -> failwith "lapacke:stegr"
  in
  check_lapack_error ret;

  let m = Int32.to_int !@_m in
  let w = match m < Owl_dense_matrix_generic.col_num !w with
    | true  -> Owl_dense_matrix_generic.resize !w [|1; m|]
    | false -> !w
  in
  let z = match m < Owl_dense_matrix_generic.col_num z with
    | true  -> Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] z
    | false -> z
  in
  w, z


let stein
  : type a b. kind:(a, b) kind -> d:(float, b) t -> e:(float, b) t
  -> w:(float, b) t -> iblock:(int32, int32_elt) t -> isplit:(int32, int32_elt) t
  -> (a, b) t * (int32, int32_elt) t
  = fun ~kind ~d ~e ~w ~iblock:_ ~isplit:_ ->
  let n = Owl_dense_matrix_generic.numel d in
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let m = Owl_dense_matrix_generic.numel w in
  assert (n <= m);
  let _kind = kind in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let e = Owl_dense_matrix_generic.resize e [|1; n|] in
  let ldz = Stdlib.max 1 m in
  let z = Genarray.create _kind _layout [|n; |] in
  let ifailv = Genarray.create int32 _layout [|1; m|] in
  (* TODO: cases where inputs are invalid, refer to julia implementation *)
  let iblock = Genarray.create int32 _layout [|1;n|] in
  let isplit = Genarray.create int32 _layout [|1;n|] in

  let _w = bigarray_start Ctypes_static.Genarray w in
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _z = bigarray_start Ctypes_static.Genarray z in
  let _iblock = bigarray_start Ctypes_static.Genarray iblock in
  let _isplit = bigarray_start Ctypes_static.Genarray isplit in
  let _ifailv = bigarray_start Ctypes_static.Genarray ifailv in

  let ret = match _kind with
    | Float32   -> L.sstein ~layout ~n ~d:_d ~e:_e ~m ~w:_w ~iblock:_iblock ~isplit:_isplit ~z:_z ~ldz ~ifailv:_ifailv
    | Float64   -> L.dstein ~layout ~n ~d:_d ~e:_e ~m ~w:_w ~iblock:_iblock ~isplit:_isplit ~z:_z ~ldz ~ifailv:_ifailv
    | Complex32 -> L.cstein ~layout ~n ~d:_d ~e:_e ~m ~w:_w ~iblock:_iblock ~isplit:_isplit ~z:_z ~ldz ~ifailv:_ifailv
    | Complex64 -> L.zstein ~layout ~n ~d:_d ~e:_e ~m ~w:_w ~iblock:_iblock ~isplit:_isplit ~z:_z ~ldz ~ifailv:_ifailv
    | _         -> failwith "lapacke:stein"
  in
  check_lapack_error ret;
  z, ifailv


let syconv
  : type a b. uplo:char -> way:char -> a:(a, b) t -> ipiv:(int32, int32_elt) t
  -> (a, b) t
  = fun ~uplo ~way ~a ~ipiv ->
  assert (uplo = 'U' || uplo = 'L');
  assert (way = 'C' || way = 'R');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let e = Genarray.create _kind _layout [|1;n|] in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.ssyconv ~layout ~uplo ~way ~n ~a:_a ~lda ~ipiv:_ipiv ~e:_e
    | Float64   -> L.dsyconv ~layout ~uplo ~way ~n ~a:_a ~lda ~ipiv:_ipiv ~e:_e
    | Complex32 -> L.csyconv ~layout ~uplo ~way ~n ~a:_a ~lda ~ipiv:_ipiv ~e:_e
    | Complex64 -> L.zsyconv ~layout ~uplo ~way ~n ~a:_a ~lda ~ipiv:_ipiv ~e:_e
    | _         -> failwith "lapacke:syconv"
  in
  check_lapack_error ret;
  e


let sysv
  : type a b. uplo:char -> a:(a, b) t -> b:(a, b) t
  -> (a, b) t * (a, b) t * (int32, int32_elt) t
  = fun ~uplo ~a ~b ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  assert (n = Owl_dense_matrix_generic.row_num b);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> L.ssysv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Float64   -> L.dsysv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex32 -> L.csysv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex64 -> L.zsysv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | _         -> failwith "lapacke:sysv"
  in
  check_lapack_error ret;
  b, a, ipiv


let sytrf
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t * (int32, int32_elt) t * int
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.ssytrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Float64   -> L.dsytrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex32 -> L.csytrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zsytrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:sytrf"
  in
  check_lapack_error ret;
  a, ipiv, ret


let sytrf_rook
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t * (int32, int32_elt) t * int
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.ssytrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Float64   -> L.dsytrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex32 -> L.csytrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zsytrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:sytrf_rook"
  in
  check_lapack_error ret;
  a, ipiv, ret


let sytri
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.ssytri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Float64   -> L.dsytri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex32 -> L.csytri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zsytri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:sytri"
  in
  check_lapack_error ret;
  a


let sytrs
  : type a b. uplo:char -> a:(a, b) t -> ipiv:(int32, int32_elt) t -> b:(a, b) t
  -> (a, b) t
  = fun ~uplo ~a ~ipiv ~b ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  assert (n = Owl_dense_matrix_generic.row_num b);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Float32   -> L.ssytrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Float64   -> L.dsytrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex32 -> L.csytrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex64 -> L.zsytrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | _         -> failwith "lapacke:sytrs"
  in
  check_lapack_error ret;
  b


let hesv
  : type a. uplo:char -> a:(Complex.t, a) t -> b:(Complex.t, a) t
  -> (Complex.t, a) t * (Complex.t, a) t * (int32, int32_elt) t
  = fun ~uplo ~a ~b ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  assert (n = Owl_dense_matrix_generic.row_num b);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Complex32 -> L.chesv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv
                     ~b:_b ~ldb
    | Complex64 -> L.zhesv ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv
                     ~b:_b ~ldb
  in
  check_lapack_error ret;
  b, a, ipiv


let hetrf
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t * (int32, int32_elt) t * int
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Complex32 -> L.chetrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zhetrf ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:hetrf"
  in
  check_lapack_error ret;
  a, ipiv, ret


let hetrf_rook
  : type a b. uplo:char -> a:(a, b) t -> (a, b) t * (int32, int32_elt) t * int
  = fun ~uplo ~a ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let ipiv = Genarray.create int32 _layout [|1;n|] in
  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Complex32 -> L.chetrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zhetrf_rook ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | _         -> failwith "lapacke:hetrf_rook"
  in
  check_lapack_error ret;
  a, ipiv, ret


let hetri
  : type a. uplo:char -> a:(Complex.t, a) t -> ipiv:(int32, int32_elt) t
  -> (Complex.t, a) t
  = fun ~uplo ~a ~ipiv ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Complex32 -> L.chetri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
    | Complex64 -> L.zhetri ~layout ~uplo ~n ~a:_a ~lda ~ipiv:_ipiv
  in
  check_lapack_error ret;
  a


let hetrs
  : type a. uplo:char -> a:(Complex.t, a) t -> ipiv:(int32, int32_elt) t
  -> b:(Complex.t, a) t -> (Complex.t, a) t
  = fun ~uplo ~a ~ipiv ~b ->
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  assert (n = Owl_dense_matrix_generic.row_num b);
  let nrhs = Owl_dense_matrix_generic.col_num b in
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ipiv = bigarray_start Ctypes_static.Genarray ipiv in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in

  let ret = match _kind with
    | Complex32 ->
      L.chetrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
    | Complex64 ->
      L.zhetrs ~layout ~uplo ~n ~nrhs ~a:_a ~lda ~ipiv:_ipiv ~b:_b ~ldb
  in
  check_lapack_error ret;
  b


let syev
  : type a. jobz:char -> uplo:char -> a:(float, a) t -> (float, a) t * (float, a) t
  = fun ~jobz ~uplo ~a ->
  assert (jobz = 'N' || jobz = 'V');
  assert (uplo = 'U' || uplo = 'L');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let w = Genarray.create _kind _layout [|1;n|] in
  let _w = bigarray_start Ctypes_static.Genarray w in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.ssyev ~layout ~jobz ~uplo ~n ~a:_a ~lda ~w:_w
    | Float64   -> L.dsyev ~layout ~jobz ~uplo ~n ~a:_a ~lda ~w:_w
  in
  check_lapack_error ret;

  match jobz with
  | 'V' -> w, a
  | _   -> w, Genarray.create _kind _layout [|0;0|]


let syevr
  : type a. jobz:char -> range:char -> uplo:char -> a:(float, a) t -> vl:float
  -> vu:float -> il:int -> iu:int
  -> abstol:float -> (float, a) t * (float, a) t
  = fun ~jobz ~range ~uplo ~a ~vl ~vu ~il ~iu ~abstol ->
  assert (jobz = 'N' || jobz = 'V');
  assert (range = 'A' || range = 'V' && range = 'I');
  assert (uplo = 'U' || uplo = 'L' );

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _ = match range with
    | 'I' -> assert (1 <= il && il <= iu && iu <= n)
    | 'V' -> assert (vu <= vl)
    | _   -> ()
  in
  let lda = Stdlib.max 1 n in
  let ldz = match jobz with
    | 'V' -> Stdlib.max 1 m
    | _   -> 1
  in
  let z = Genarray.create _kind _layout [|n;ldz|] in
  let w = Genarray.create _kind _layout [|1;n|] in
  let isuppz = Genarray.create int32 _layout [|1;(2 * m)|] in

  let _m = Ctypes.(allocate int32_t 0l) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _z = bigarray_start Ctypes_static.Genarray z in
  let _w = bigarray_start Ctypes_static.Genarray w in
  let _isuppz = bigarray_start Ctypes_static.Genarray isuppz in

  let ret = match _kind with
    | Float32   ->
      L.ssyevr ~layout ~jobz ~range ~uplo ~n ~a:_a ~lda ~vl ~vu ~il ~iu ~abstol
        ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz
    | Float64   ->
      L.dsyevr ~layout ~jobz ~range ~uplo ~n ~a:_a ~lda ~vl ~vu ~il ~iu ~abstol
        ~m:_m ~w:_w ~z:_z ~ldz ~isuppz:_isuppz
  in
  check_lapack_error ret;

  let m = Int32.to_int !@_m in
  let w = Owl_dense_matrix_generic.resize w [|1; m|] in
  match jobz with
  | 'V' -> w, Owl_dense_matrix_generic.get_slice [[]; [0;m-1]] z
  | _   -> w, Genarray.create _kind _layout [|0;0|]


let sygvd : type a. ityp:int -> jobz:char ->
  uplo:char -> a:(float, a) t -> b:(float, a) t
  -> (float, a) t * (float, a) t * (float, a) t
  = fun ~ityp ~jobz ~uplo ~a ~b ->
  assert (ityp > 0 && ityp < 4);
  assert (jobz = 'N' || jobz = 'V');
  assert (uplo = 'U' || uplo = 'L' );

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (m = n && n = mb && n = nb);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let w = Genarray.create _kind _layout [|1;n|] in
  let _w = bigarray_start Ctypes_static.Genarray w in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in

  let ret = match _kind with
    | Float32   ->
      L.ssygvd ~layout ~ityp ~jobz ~uplo ~n ~a:_a ~lda ~b:_b ~ldb ~w:_w
    | Float64   ->
      L.dsygvd ~layout ~ityp ~jobz ~uplo ~n ~a:_a ~lda ~b:_b ~ldb ~w:_w
  in
  check_lapack_error ret;
  w, a, b


let bdsqr
  : type a b. uplo:char -> d:(float, b) t -> e:(float, b) t -> vt:(a, b) t
  -> u:(a, b) t -> c:(a, b) t -> (float, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun  ~uplo ~d ~e ~vt ~u ~c ->
  assert (uplo = 'U' || uplo = 'L' );

  let n = Owl_dense_matrix_generic.numel d in
  let ncvt = Owl_dense_matrix_generic.col_num vt in
  let nru = Owl_dense_matrix_generic.row_num u in
  let ncc = Owl_dense_matrix_generic.col_num c in
  assert (Owl_dense_matrix_generic.row_num vt = n);
  assert (Owl_dense_matrix_generic.row_num c = n);
  let n_e = Owl_dense_matrix_generic.numel e in
  assert (n_e = n - 1);
  let _kind = Genarray.kind vt in
  let _layout = Genarray.layout vt in
  let layout = lapacke_layout _layout in

  let ldvt = Stdlib.max 1 (_stride vt) in
  let ldu = Stdlib.max 1 (_stride u) in
  let ldc = Stdlib.max 1 ncc in
  assert (ldvt >= ncvt);
  assert (ldu >= n);
  assert (ldc >= ncc);
  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _vt = bigarray_start Ctypes_static.Genarray vt in
  let _u = bigarray_start Ctypes_static.Genarray u in
  let _c = bigarray_start Ctypes_static.Genarray c in

  let ret = match _kind with
    | Float32   ->
      L.sbdsqr ~layout ~uplo ~n ~ncvt ~nru ~ncc ~d:_d ~e:_e ~vt:_vt ~ldvt
        ~u:_u ~ldu ~c:_c ~ldc
    | Float64   ->
      L.dbdsqr ~layout ~uplo ~n ~ncvt ~nru ~ncc ~d:_d ~e:_e ~vt:_vt ~ldvt
        ~u:_u ~ldu ~c:_c ~ldc
    | Complex32 ->
      L.cbdsqr ~layout ~uplo ~n ~ncvt ~nru ~ncc ~d:_d ~e:_e ~vt:_vt ~ldvt
        ~u:_u ~ldu ~c:_c ~ldc
    | Complex64 ->
      L.zbdsqr ~layout ~uplo ~n ~ncvt ~nru ~ncc ~d:_d ~e:_e ~vt:_vt ~ldvt
        ~u:_u ~ldu ~c:_c ~ldc
    | _         -> failwith "lapacke:bdsqr"
  in
  check_lapack_error ret;
  d, vt, u, c


let bdsdc
  : type a. uplo:char -> compq:char -> d:(float, a) t -> e:(float, a) t
  -> (float, a) t * (float, a) t * (float, a) t *
     (float, a) t * (float, a) t * (int32, int32_elt) t
  = fun ~uplo ~compq ~d ~e ->
  assert (uplo = 'U' || uplo = 'L');
  assert (compq = 'N' || compq = 'P' || compq = 'I');

  let n = Owl_dense_matrix_generic.numel d in
  let _kind = Genarray.kind d in
  let _layout = Genarray.layout d in
  let layout = lapacke_layout _layout in

  let u = match compq with
    | 'I' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let vt = match compq with
    | 'I' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let q = match compq with
    | 'P' -> (
      let smlsiz = 100. in
      let n = float_of_int n in
      let ldq = Owl_maths.(n *. (11. +. 2. *. smlsiz +. 8. *. round (log ((n /. (smlsiz +. 1.))) /. (log 2.)))) in
      Genarray.create _kind _layout [|1; (Stdlib.max 0 (int_of_float ldq))|]
      )
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let iq = match compq with
    | 'P' -> (
      let smlsiz = 100. in
      let n = float_of_int n in
      let ldiq = Owl_maths.(n *. (3. +. 3. *. round (log (n /. (smlsiz +. 1.)) /. (log 2.)))) in
      Genarray.create int32 _layout [|1; (Stdlib.max 0 (int_of_float ldiq))|]
      )
    | _   -> Genarray.create int32 _layout [|0; n|]
  in
  let ldu = Stdlib.max 1 (_stride u) in
  let ldvt = Stdlib.max 1 (_stride vt) in

  let _d = bigarray_start Ctypes_static.Genarray d in
  let _e = bigarray_start Ctypes_static.Genarray e in
  let _u = bigarray_start Ctypes_static.Genarray u in
  let _vt = bigarray_start Ctypes_static.Genarray vt in
  let _q = bigarray_start Ctypes_static.Genarray q in
  let _iq = bigarray_start Ctypes_static.Genarray iq in

  let ret = match _kind with
    | Float32   ->
      L.sbdsdc ~layout ~uplo ~compq ~n ~d:_d ~e:_e ~u:_u ~ldu
        ~vt:_vt ~ldvt ~q:_q ~iq:_iq
    | Float64   ->
      L.dbdsdc ~layout ~uplo ~compq ~n ~d:_d ~e:_e ~u:_u ~ldu
        ~vt:_vt ~ldvt ~q:_q ~iq:_iq
  in
  check_lapack_error ret;
  d, e, u, vt, q, iq


let gecon
  : type a b. norm:char -> a:(a, b) t -> anorm:float -> float
  = fun ~norm ~a ~anorm ->
  assert (norm = '1' || norm = 'O' || norm = 'I');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let rcond = ref 0. in

  let ret = match _kind with
    | Float32   -> (
        let _rcond = Ctypes.(allocate float 0.) in
        let r = L.sgecon ~layout ~norm ~n ~a:_a ~lda ~anorm ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Float64   -> (
        let _rcond = Ctypes.(allocate double 0.) in
        let r = L.dgecon ~layout ~norm ~n ~a:_a ~lda ~anorm ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Complex32 -> (
        let _rcond = Ctypes.(allocate float 0.) in
        let r = L.cgecon ~layout ~norm ~n ~a:_a ~lda ~anorm ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | Complex64 -> (
        let _rcond = Ctypes.(allocate double 0.) in
        let r = L.zgecon ~layout ~norm ~n ~a:_a ~lda ~anorm ~rcond:_rcond in
        rcond := !@_rcond;
        r
      )
    | _         -> failwith "lapacke:gecon"
  in
  check_lapack_error ret;
  !rcond


let gehrd
  : type a b. ilo:int -> ihi:int -> a:(a, b) t -> (a, b) t * (a, b) t
  = fun ~ilo ~ihi ~a ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let tau = Genarray.create _kind _layout [|1; (Stdlib.max 1 (n - 1))|] in
  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.sgehrd ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dgehrd ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | Complex32 -> L.cgehrd ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zgehrd ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | _         -> failwith "lapacke:gehrd"
  in
  check_lapack_error ret;
  a, tau


let orghr : type a. ilo:int -> ihi:int ->
  a:(float, a) t -> tau:(float, a) t -> (float, a) t
  = fun ~ilo ~ihi ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let n_tau = Owl_dense_matrix_generic.numel tau in
  assert (n_tau = n - 1);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> L.sorghr ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | Float64   -> L.dorghr ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  a


let unghr : type a. ilo:int -> ihi:int ->
  a:(Complex.t, a) t -> tau:(Complex.t, a) t -> (Complex.t, a) t
  = fun ~ilo ~ihi ~a ~tau ->
  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let n_tau = Owl_dense_matrix_generic.numel tau in
  assert (n_tau = n - 1);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let _tau = bigarray_start Ctypes_static.Genarray tau in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Complex32 -> L.cunghr ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
    | Complex64 -> L.zunghr ~layout ~n ~ilo ~ihi ~a:_a ~lda ~tau:_tau
  in
  check_lapack_error ret;
  a


let gees : type a b.
  jobvs:char -> a:(a, b) t -> (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~jobvs ~a ->
  let sort = 'N' in
  assert (jobvs = 'N' || jobvs = 'V');
  assert (sort = 'N' || sort = 'S');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (m = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let vs = match jobvs with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let ldvs = Stdlib.max 1 (_stride vs) in
  let wr = ref (Genarray.create _kind _layout [|0;0|]) in
  let wi = ref (Genarray.create _kind _layout [|0;0|]) in

  let _select = Ctypes.null in
  let _sdim = Ctypes.(allocate int32_t 0l) in
  let _vs = bigarray_start Ctypes_static.Genarray vs in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let lda = Stdlib.max 1 (_stride a) in

  let ret = match _kind with
    | Float32   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let r = L.sgees ~layout ~jobvs ~sort ~select:_select ~n ~a:_a ~lda
            ~sdim:_sdim ~wr:_wr ~wi:_wi ~vs:_vs ~ldvs in
        wr := wr';
        wi := wi';
        r
      )
    | Float64   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let r = L.dgees ~layout ~jobvs ~sort ~select:_select ~n ~a:_a ~lda
            ~sdim:_sdim ~wr:_wr ~wi:_wi ~vs:_vs ~ldvs in
        wr := wr';
        wi := wi';
        r
      )
    | Complex32 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.cgees ~layout ~jobvs ~sort ~select:_select ~n
            ~a:_a ~lda ~sdim:_sdim ~w:_w ~vs:_vs ~ldvs in
        wr := w';
        wi := w';
        r
      )
    | Complex64 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let r = L.zgees ~layout ~jobvs ~sort ~select:_select ~n
            ~a:_a ~lda ~sdim:_sdim ~w:_w ~vs:_vs ~ldvs in
        wr := w';
        wi := w';
        r
      )
    | _         -> failwith "lapacke:gees"
  in
  check_lapack_error ret;
  (* NOTE: wr and wi are the same for complex flavour *)
  a, vs, !wr, !wi


let gges
  : type a b. jobvsl:char -> jobvsr:char -> a:(a, b) t -> b:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~jobvsl ~jobvsr ~a ~b ->
  let sort = 'N' in
  assert (jobvsl = 'N' || jobvsl = 'V');
  assert (jobvsr = 'N' || jobvsr = 'V');
  assert (sort = 'N' || sort = 'S');

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (m = n && n = mb && mb = nb);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let vsl = match jobvsl with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let vsr = match jobvsr with
    | 'V' -> Genarray.create _kind _layout [|n;n|]
    | _   -> Genarray.create _kind _layout [|0;n|]
  in
  let ldvsl = Stdlib.max 1 (_stride vsl) in
  let ldvsr = Stdlib.max 1 (_stride vsr) in
  let alphar = ref (Genarray.create _kind _layout [|0;0|]) in
  let alphai = ref (Genarray.create _kind _layout [|0;0|]) in
  let beta = Genarray.create _kind _layout [|1;n|] in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _vsl = bigarray_start Ctypes_static.Genarray vsl in
  let _vsr = bigarray_start Ctypes_static.Genarray vsr in
  let _beta = bigarray_start Ctypes_static.Genarray beta in
  let _selctg = Ctypes.null in
  let _sdim = Ctypes.(allocate int32_t 0l) in

  let ret = match _kind with
    | Float32   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let r = L.sgges
            ~layout ~jobvsl ~jobvsr ~sort ~selctg:_selctg ~n ~a:_a ~lda ~b:_b
            ~ldb ~sdim:_sdim ~alphar:_alphar ~alphai:_alphai ~beta:_beta
            ~vsl:_vsl ~ldvsl ~vsr:_vsr ~ldvsr in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Float64   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let r = L.dgges ~layout ~jobvsl ~jobvsr ~sort ~selctg:_selctg ~n ~a:_a
            ~lda ~b:_b ~ldb ~sdim:_sdim ~alphar:_alphar ~alphai:_alphai
            ~beta:_beta ~vsl:_vsl ~ldvsl ~vsr:_vsr ~ldvsr in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Complex32 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let r = L.cgges ~layout ~jobvsl ~jobvsr ~sort ~selctg:_selctg ~n ~a:_a
            ~lda ~b:_b ~ldb ~sdim:_sdim ~alpha ~beta:_beta ~vsl:_vsl ~ldvsl
            ~vsr:_vsr ~ldvsr in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | Complex64 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let r = L.zgges ~layout ~jobvsl ~jobvsr ~sort ~selctg:_selctg ~n ~a:_a
            ~lda ~b:_b ~ldb ~sdim:_sdim ~alpha ~beta:_beta ~vsl:_vsl ~ldvsl
            ~vsr:_vsr ~ldvsr in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | _         -> failwith "lapacke:gges"
  in
  check_lapack_error ret;
  (* NOTE: alphar and alphai are the same for complex flavour *)
  a, b, !alphar, !alphai, beta, vsl, vsr


let trexc
  : type a b. compq:char -> t:(a, b) t -> q:(a, b) t -> ifst:int -> ilst:int
  -> (a, b) t * (a, b) t
  = fun ~compq ~t ~q ~ifst ~ilst ->
  assert (compq = 'N' || compq = 'V');

  let m = Owl_dense_matrix_generic.row_num t in
  let n = Owl_dense_matrix_generic.col_num t in
  assert (m = n);
  assert (1 <= ifst && ifst <= n);
  assert (1 <= ilst && ilst <= n);
  let _kind = Genarray.kind t in
  let _layout = Genarray.layout t in
  let layout = lapacke_layout _layout in

  let ldt = Stdlib.max 1 (_stride t) in
  let ldq = Stdlib.max 1 (_stride q) in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let _q = bigarray_start Ctypes_static.Genarray q in
  let _ifst = Ctypes.(allocate int32_t (Int32.of_int ifst)) in
  let _ilst = Ctypes.(allocate int32_t (Int32.of_int ilst)) in

  let ret = match _kind with
    | Float32   -> L.strexc ~layout ~compq ~n ~t:_t ~ldt
                     ~q:_q ~ldq ~ifst:_ifst ~ilst:_ilst
    | Float64   -> L.dtrexc ~layout ~compq ~n ~t:_t ~ldt
                     ~q:_q ~ldq ~ifst:_ifst ~ilst:_ilst
    | Complex32 -> L.ctrexc ~layout ~compq ~n ~t:_t ~ldt
                     ~q:_q ~ldq ~ifst ~ilst
    | Complex64 -> L.ztrexc ~layout ~compq ~n ~t:_t ~ldt
                     ~q:_q ~ldq ~ifst ~ilst
    | _         -> failwith "lapacke:trexc"
  in
  check_lapack_error ret;
  t, q


let trsen : type a b.
  job:char -> compq:char -> select:(int32, int32_elt) t -> t:(a, b) t
  -> q:(a, b) t -> (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~job ~compq ~select ~t ~q ->
  assert (job = 'N' || job = 'E' || job = 'V' || job = 'B');
  assert (compq = 'N' || compq = 'V');

  let mt = Owl_dense_matrix_generic.row_num t in
  let n = Owl_dense_matrix_generic.col_num t in
  assert (mt = n);
  let _kind = Genarray.kind t in
  let _layout = Genarray.layout t in
  let layout = lapacke_layout _layout in

  let ldt = Stdlib.max 1 (_stride t) in
  let ldq = Stdlib.max 1 (_stride q) in
  let _t = bigarray_start Ctypes_static.Genarray t in
  let _q = bigarray_start Ctypes_static.Genarray q in
  let wr = ref (Genarray.create _kind _layout [|0;0|]) in
  let wi = ref (Genarray.create _kind _layout [|0;0|]) in
  let _select = bigarray_start Ctypes_static.Genarray select in
  let _m = Ctypes.(allocate int32_t 0l) in

  let ret = match _kind with
    | Float32   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let _s = Ctypes.(allocate float 0.) in
        let _sep = Ctypes.(allocate float 0.) in
        let r = L.strsen ~layout ~job ~compq ~select:_select ~n
            ~t:_t ~ldt ~q:_q ~ldq ~wr:_wr ~wi:_wi ~m:_m ~s:_s ~sep:_sep in
        wr := wr';
        wi := wi';
        r
      )
    | Float64   -> (
        let wr' = Genarray.create _kind _layout [|1;n|] in
        let _wr = bigarray_start Ctypes_static.Genarray wr' in
        let wi' = Genarray.create _kind _layout [|1;n|] in
        let _wi = bigarray_start Ctypes_static.Genarray wi' in
        let _s = Ctypes.(allocate double 0.) in
        let _sep = Ctypes.(allocate double 0.) in
        let r = L.dtrsen ~layout ~job ~compq ~select:_select ~n ~t:_t ~ldt
            ~q:_q ~ldq ~wr:_wr ~wi:_wi ~m:_m ~s:_s ~sep:_sep in
        wr := wr';
        wi := wi';
        r
      )
    | Complex32 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let _s = Ctypes.(allocate float 0.) in
        let _sep = Ctypes.(allocate float 0.) in
        let r = L.ctrsen ~layout ~job ~compq ~select:_select ~n ~t:_t ~ldt
            ~q:_q ~ldq  ~w:_w ~m:_m ~s:_s ~sep:_sep in
        wr := w';
        wi := w';
        r
      )
    | Complex64 -> (
        let w' = Genarray.create _kind _layout [|1;n|] in
        let _w = bigarray_start Ctypes_static.Genarray w' in
        let _s = Ctypes.(allocate float 0.) in
        let _sep = Ctypes.(allocate float 0.) in
        let r = L.ztrsen ~layout ~job ~compq ~select:_select ~n ~t:_t ~ldt
            ~q:_q ~ldq  ~w:_w ~m:_m ~s:_s ~sep:_sep in
        wr := w';
        wi := w';
        r
      )
    | _         -> failwith "lapacke:trsen"
  in
  check_lapack_error ret;
  (* NOTE: wr and wi are the same for complex flavour *)
  t, q, !wr, !wi


let tgsen
  : type a b. select:(int32, int32_elt) t ->
  a:(a, b) t -> b:(a, b) t -> q:(a, b) t -> z:(a, b) t
  -> (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t * (a, b) t
  = fun ~select ~a ~b ~q ~z ->
  (* set these values by default *)
  let ijob = 0 in
  let wantq = 1 in
  let wantz = 1 in
  assert (0 <= ijob && ijob <= 5);
  assert (wantq = 0 || wantq = 1);
  assert (wantz = 0 || wantz = 1);

  let ma = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  assert (ma = n);
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  assert (mb = nb);
  let mq = Owl_dense_matrix_generic.row_num q in
  let nq = Owl_dense_matrix_generic.col_num q in
  assert (mq = nq);
  let mz = Owl_dense_matrix_generic.row_num z in
  let nz = Owl_dense_matrix_generic.col_num z in
  assert (mz = nz);
  assert (n = nb);
  assert (n = nq);
  assert (n = nz);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ldq = Stdlib.max 1 (_stride q) in
  let ldz = Stdlib.max 1 (_stride z) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _q = bigarray_start Ctypes_static.Genarray q in
  let _z = bigarray_start Ctypes_static.Genarray z in
  let alphar = ref (Genarray.create _kind _layout [|0;0|]) in
  let alphai = ref (Genarray.create _kind _layout [|0;0|]) in
  let beta = Genarray.create _kind _layout [|1;n|] in
  let _beta = bigarray_start Ctypes_static.Genarray beta in
  let _select = bigarray_start Ctypes_static.Genarray select in

  (* FIXME: not sure summation is really needed *)
  let m = ref 0l in
  for i = 0 to Owl_dense_matrix_generic.col_num select - 1 do
    m := Int32.add !m (Owl_dense_matrix_generic.get select 0 i)
  done;
  let _m = Ctypes.(allocate int32_t !m) in

  let ret = match _kind with
    | Float32   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let _pl = Ctypes.(allocate float 0.) in
        let _pr = Ctypes.(allocate float 0.) in
        let dif = Genarray.create float32 _layout [|1;2|] in
        let _dif = bigarray_start Ctypes_static.Genarray dif in
        let r = L.stgsen ~layout ~ijob ~wantq ~wantz ~select:_select ~n
            ~a:_a ~lda ~b:_b ~ldb ~alphar:_alphar ~alphai:_alphai ~beta:_beta
            ~q:_q ~ldq ~z:_z ~ldz ~m:_m ~pl:_pl ~pr:_pr ~dif:_dif in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Float64   -> (
        let alphar' = Genarray.create _kind _layout [|1;n|] in
        let _alphar = bigarray_start Ctypes_static.Genarray alphar' in
        let alphai' = Genarray.create _kind _layout [|1;n|] in
        let _alphai = bigarray_start Ctypes_static.Genarray alphai' in
        let _pl = Ctypes.(allocate double 0.) in
        let _pr = Ctypes.(allocate double 0.) in
        let dif = Genarray.create float64 _layout [|1;2|] in
        let _dif = bigarray_start Ctypes_static.Genarray dif in
        let r = L.dtgsen ~layout ~ijob ~wantq ~wantz ~select:_select ~n
            ~a:_a ~lda ~b:_b ~ldb ~alphar:_alphar ~alphai:_alphai ~beta:_beta
            ~q:_q ~ldq ~z:_z ~ldz  ~m:_m ~pl:_pl ~pr:_pr ~dif:_dif in
        alphar := alphar';
        alphai := alphai';
        r
      )
    | Complex32 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _pl = Ctypes.(allocate float 0.) in
        let _pr = Ctypes.(allocate float 0.) in
        let dif = Genarray.create float32 _layout [|1;2|] in
        let _dif = bigarray_start Ctypes_static.Genarray dif in
        let r = L.ctgsen ~layout ~ijob ~wantq ~wantz ~select:_select ~n
            ~a:_a ~lda ~b:_b ~ldb ~alpha ~beta:_beta ~q:_q ~ldq ~z:_z ~ldz
            ~m:_m ~pl:_pl ~pr:_pr ~dif:_dif in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | Complex64 -> (
        let alpha' = Genarray.create _kind _layout [|1;n|] in
        let alpha = bigarray_start Ctypes_static.Genarray alpha' in
        let _pl = Ctypes.(allocate double 0.) in
        let _pr = Ctypes.(allocate double 0.) in
        let dif = Genarray.create float64 _layout [|1;2|] in
        let _dif = bigarray_start Ctypes_static.Genarray dif in
        let r = L.ztgsen ~layout ~ijob ~wantq ~wantz ~select:_select ~n
            ~a:_a ~lda ~b:_b ~ldb ~alpha ~beta:_beta ~q:_q ~ldq ~z:_z ~ldz
            ~m:_m ~pl:_pl ~pr:_pr ~dif:_dif in
        alphar := alpha';
        alphai := alpha';
        r
      )
    | _         -> failwith "lapacke:tgsen"
  in
  check_lapack_error ret;
  (* NOTE: alphar and alphai are the same for complex flavour *)
  a, b, !alphar, !alphai, beta, q, z


let trsyl
  : type a b. trana:char -> tranb:char -> isgn:int -> a:(a, b) t
  -> b:(a, b) t -> c:(a, b) t -> (a, b) t * float
  = fun ~trana ~tranb ~isgn ~a ~b ~c ->
  assert (trana = 'N' || trana = 'T' || trana = 'C');
  assert (tranb = 'N' || tranb = 'T' || tranb = 'C');
  assert (isgn = -1 || isgn = +1);

  let m = Owl_dense_matrix_generic.row_num a in
  let n = Owl_dense_matrix_generic.col_num a in
  let mb = Owl_dense_matrix_generic.row_num b in
  let nb = Owl_dense_matrix_generic.col_num b in
  let mc = Owl_dense_matrix_generic.row_num c in
  let nc = Owl_dense_matrix_generic.col_num c in
  assert (m = n);
  assert (mb = nb);
  assert (mc = m && nc = n);
  let _kind = Genarray.kind a in
  let _layout = Genarray.layout a in
  let layout = lapacke_layout _layout in

  let lda = Stdlib.max 1 (_stride a) in
  let ldb = Stdlib.max 1 (_stride b) in
  let ldc = Stdlib.max 1 (_stride c) in
  let _a = bigarray_start Ctypes_static.Genarray a in
  let _b = bigarray_start Ctypes_static.Genarray b in
  let _c = bigarray_start Ctypes_static.Genarray c in
  let scale = ref 0. in

  let ret = match _kind with
    | Float32   -> (
        let _scale = Ctypes.(allocate float 0.) in
        let r = L.strsyl ~layout ~trana ~tranb ~isgn ~m ~n ~a:_a ~lda
            ~b:_b ~ldb ~c:_c ~ldc ~scale:_scale in
        scale := !@_scale;
        r
      )
    | Float64   -> (
        let _scale = Ctypes.(allocate double 0.) in
        let r = L.dtrsyl ~layout ~trana ~tranb ~isgn ~m ~n ~a:_a ~lda
            ~b:_b ~ldb ~c:_c ~ldc ~scale:_scale in
        scale := !@_scale;
        r
      )
    | Complex32 -> (
        let _scale = Ctypes.(allocate float 0.) in
        let r = L.ctrsyl ~layout ~trana ~tranb ~isgn ~m ~n ~a:_a
            ~lda ~b:_b ~ldb ~c:_c ~ldc ~scale:_scale in
        scale := !@_scale;
        r
      )
    | Complex64 -> (
        let _scale = Ctypes.(allocate double 0.) in
        let r = L.ztrsyl ~layout ~trana ~tranb ~isgn ~m ~n ~a:_a
            ~lda ~b:_b ~ldb ~c:_c ~ldc ~scale:_scale in
        scale := !@_scale;
        r
      )
    | _         -> failwith "lapacke:trsyl"
  in
  check_lapack_error ret;
  c, !scale



(* ends here *)