// Lapack.sgesvd_s: The singular values of an m x n row-major matrix (consumed), largest

// bend-blas: the platform CBLAS and LAPACK, opened once per process
// (Accelerate on macOS; OpenBLAS, else libcblas + libblas + liblapack,
// on Linux). Bend's clang command is fixed, so nothing is linked: every
// symbol comes through dlsym.
#ifndef BEND_BLAS_LOADER
#define BEND_BLAS_LOADER 1
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <dlfcn.h>
#include <math.h>

static void* bend_lib_open(const char** names) {
  void* lib = NULL;
  for (int i = 0; names[i] != NULL && lib == NULL; i += 1) {
    lib = dlopen(names[i], RTLD_NOW);
  }
  return lib;
}

// a symbol from the BLAS, the CBLAS wrapper beside a reference BLAS, or
// the LAPACK library, in that order; NULL when none holds it
static void* bend_blas_sym(const char* name) {
  static void* libs[3] = { NULL, NULL, NULL };
  static int tried = 0;
  if (!tried) {
    const char* blas[] = {
      "/System/Library/Frameworks/Accelerate.framework/Accelerate",
      "libopenblas.so.0", "libopenblas.so", "libblas.so.3", "libblas.so", NULL };
    const char* cblas[] = { "libcblas.so.3", "libcblas.so", NULL };
    const char* lapack[] = { "liblapack.so.3", "liblapack.so", NULL };
    tried = 1;
    libs[0] = bend_lib_open(blas);
    libs[1] = bend_lib_open(cblas);
    libs[2] = bend_lib_open(lapack);
  }
  for (int i = 0; i < 3; i += 1) {
    void* sym = libs[i] == NULL ? NULL : dlsym(libs[i], name);
    if (sym != NULL) {
      return sym;
    }
  }
  return NULL;
}

// an F32 argument is the low 32 bits of its word; an F32 answer the same
static float bend_f32_of(Term t) {
  u32 bits = (u32)t;
  float x;
  memcpy(&x, &bits, 4);
  return x;
}

static Term bend_f32_term(float x) {
  u32 bits;
  memcpy(&bits, &x, 4);
  return (Term)bits;
}

// a fresh block of the smallest power-of-two size holding n floats, the
// tail zero. Answers 0 past 2^31 floats, the largest block class; the
// heap itself never answers, a full heap is a fail-stop of the runtime.
static u64 bend_block_floats(Env e, u64 n, u32* cls_out) {
  u32 cls = 0;
  while ((1ull << cls) < n) {
    cls += 1;
  }
  if (cls > 31) {
    return 0;
  }
  u64 loc = heap_alloc(e, buf_wcls(cls));
  float* c = (float*)blk_ptr(e.mem, loc, 0);
  for (u64 i = n; i < (1ull << cls); i += 1) {
    c[i] = 0.0f;
  }
  *cls_out = cls;
  return loc;
}

// a fresh block holding the first n floats of an operand
static u64 bend_block_copy(Env e, const float* src, u64 n, u32* cls_out) {
  u64 loc = bend_block_floats(e, n, cls_out);
  if (loc != 0) {
    memcpy(blk_ptr(e.mem, loc, 0), src, (size_t)n * sizeof(float));
  }
  return loc;
}

// a row-major m x n operand copied to a column-major buffer (malloc)
static float* bend_colmajor(const float* a, u32 m, u32 n) {
  float* c = (float*)malloc((size_t)m * n * sizeof(float) + 4);
  for (u32 i = 0; i < m; i += 1) {
    for (u32 j = 0; j < n; j += 1) {
      c[(size_t)j * m + i] = a[(size_t)i * n + j];
    }
  }
  return c;
}

// a symmetric operand copied as it is: its column-major view is itself
static float* bend_symcopy(const float* a, u32 n) {
  float* c = (float*)malloc((size_t)n * n * sizeof(float) + 4);
  memcpy(c, a, (size_t)n * n * sizeof(float));
  return c;
}

// a column-major m x n buffer (leading dimension ld) into a fresh
// row-major block
static u64 bend_block_rowmajor(Env e, const float* c, u32 m, u32 n, u32 ld, u32* cls_out) {
  u64 loc = bend_block_floats(e, (u64)m * n, cls_out);
  if (loc != 0) {
    float* r = (float*)blk_ptr(e.mem, loc, 0);
    for (u32 i = 0; i < m; i += 1) {
      for (u32 j = 0; j < n; j += 1) {
        r[(size_t)i * n + j] = c[(size_t)j * ld + i];
      }
    }
  }
  return loc;
}

// the capacity of an Array<F32> in floats: its block class is in the term
#define BEND_CAP(t) (1ull << (u32)term_aux(t))
#define BEND_ARR(t) ((const float*)blk_ptr(e.mem, term_loc(t), 0))
// an operand must hold the floats the sizes claim, and a size product
// must fit a block: both are checked before anything is read
#define BEND_NEED(t, n) \
  if ((u64)(n) > 0x7fffffffull) { return io_fail(e, 2, "bend-blas: a size past 2^31 floats"); } \
  if (BEND_CAP(t) < (u64)(n)) { return io_fail(e, 2, "bend-blas: an operand smaller than its sizes claim"); }
#define BEND_SYM(var, type, name) \
  type var = (type)bend_blas_sym(name); \
  if (var == NULL) { return io_fail(e, 2, "bend-blas: no BLAS or LAPACK found (Accelerate, OpenBLAS, libcblas, liblapack): " name); }
#define BEND_OUT(loc, cls, n) \
  u32 cls = 0; \
  u64 loc = bend_block_floats(e, (u64)(n), &cls); \
  if (loc == 0) { return io_fail(e, 2, "bend-blas: a result past 2^31 floats"); }
#define BEND_OUT_COPY(loc, cls, src, n) \
  u32 cls = 0; \
  u64 loc = bend_block_copy(e, (src), (u64)(n), &cls); \
  if (loc == 0) { return io_fail(e, 2, "bend-blas: a result past 2^31 floats"); }
#define BEND_DONE(loc, cls) return io_done(e, term_blk(false, cls, loc))
#endif

static Term lapack_sgesvd_s_body(Env e, Term* f, IoWork* w) {
  u32 m = (u32)f[0];
  u32 n = (u32)f[1];
  BEND_NEED(f[2], (u64)m * n);
  (void)m;
  (void)n;
  typedef void (*fn_t)(const char*, const char*, int*, int*, float*, int*, float*, float*, int*, float*, int*, float*, int*, int*);
  BEND_SYM(sgesvd, fn_t, "sgesvd_");
  u32 k = m < n ? m : n;
  float* A = bend_colmajor(BEND_ARR(f[2]), m, n);
  float* S = (float*)malloc((size_t)k * sizeof(float) + 4);
  float* U = (float*)malloc((0 ? (size_t)m * k : 1) * sizeof(float) + 4);
  float* VT = (float*)malloc((0 ? (size_t)k * n : 1) * sizeof(float) + 4);
  int M = (int)m, N = (int)n, LDU = 0 ? (int)m : 1, LDVT = 0 ? (int)k : 1, info = 0, lwork = -1;
  float wkopt = 0.0f;
  const char* job = 0 ? "S" : "N";
  if (m > 0 && n > 0) {
    sgesvd(job, job, &M, &N, A, &M, S, U, &LDU, VT, &LDVT, &wkopt, &lwork, &info);
    lwork = (int)wkopt + 1;
    float* work = (float*)malloc((size_t)lwork * sizeof(float) + 4);
    sgesvd(job, job, &M, &N, A, &M, S, U, &LDU, VT, &LDVT, work, &lwork, &info);
    free(work);
  }
  free(A);
  if (info != 0) {
    free(S); free(U); free(VT);
    return io_fail(e, (u32)(info > 0 ? info : 1), "Lapack.sgesvd: did not converge");
  }
  u32 cs = 0;
  u64 ls = bend_block_copy(e, S, k, &cs);
  free(S);
  if (ls == 0) { free(U); free(VT); return io_fail(e, 2, "bend-blas: a result past 2^31 floats"); }
  free(U); free(VT);
  BEND_DONE(ls, cs);
}

// the runtime does not drop an effect's arguments: every array operand
// is dropped here, after the body has read it
Term lapack_sgesvd_s_run(Env e, Term* f, IoWork* w) {
  Term r = lapack_sgesvd_s_body(e, f, w);
  term_drop(e, f[2]);
  return r;
}

static void __attribute__((constructor)) lapack_sgesvd_s_use(void) {
  io_eff(CID(Lapack.sgesvd_s), lapack_sgesvd_s_run, 0);
}