#ifndef INCLUDE_GUARD_MLM_H
#define INCLUDE_GUARD_MLM_H

/**
 * Multi Linear Multiplication
 * 
 *  C = alpha A x_modes[0] op(Bs[0]) ... x_modes[nrhs] op(Bs[nrhs]) + beta C
 * 
 * @param trans boolean vector of length `nrhs` indicating if `i`th RHS matrix
 *  is to be transposed. That is `op(Bs[i])` is the transposed of `Bs[i]` iff
 *  `trans[i]` is true, otherwise no-op on `Bs[i]`. Can be `NULL`, then `op` is
 *  always the identity.
 * @param modes integer vector of length `nrhs` specifying the product modes.
 * 
 * @todo TODO: continue doc. !!!
 * 
 * @param alpha scaling factor
 * @param beta scaling factor
 * @param C output memory addr.
 * 
 * @param work_mem NULL or temporary working memory of size
 *  `2 * prod(pmax(dim(A), dim(C))) + ord`
 */
int mlm(
    /*  options */  const int* trans, const int* modes, const int nrhs,
    /*     dims */  const int* dimA, const int* dimC, const int ord,
    /*  scalars */  const double alpha,
    /*   tensor */  const double* A,
    /* matrices */  const double** Bs, const int* ldBs,
    /*   scalar */  const double beta,
    /*   tensor */  double* C,
                    double* work_mem
);

#endif /* INCLUDE_GUARD_MLM_H */