/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meqh.h"
# ifdef SINGLE
# define FCT   LAPACKE_ssysv
# define SYTRF ssytrf
# define SYSV ssysv
# else
# define FCT   LAPACKE_dsysv
# define SYTRF dsytrf
# define SYSV dsysv
# endif
lapack_int FCT(int matrix_layout, char uplo, lapack_int n, lapack_int nrhs, FLT *a, 
               lapack_int lda, lapack_int *ipiv, FLT *b, lapack_int ldb)
{
 ptrdiff_t n_, nrhs_, lda_, ldb_;
 ptrdiff_t info;
 ptrdiff_t lwork = -1;
 ptrdiff_t *ipiv_;
 FLT *work;
 FLT work_query;
 n_ = n; nrhs_ = nrhs; lda_ = lda; ldb_ = ldb;
 /* Do not use provided pointer for ipiv which is of type int * and 
 * not compatible with the expected types for SYTRF and SYSV.
 * This bit is quite bad, but needed for now in order to maintain compatibility
 * with other BLAS/LAPACKE alternatives */
 ipiv_ = malloc(n*sizeof(ptrdiff_t));
 /* Query optimal work array size */
 SYTRF(&uplo, &n_, a, &lda_, ipiv_, &work_query, &lwork, &info);
 if (info < 0) return (lapack_int) info-1;
 lwork = (lapack_int) work_query;
 work = malloc(sizeof(FLT)*lwork);
 if (work == NULL) return -1;
 /* Actual call */
 SYSV(&uplo, &n_, &nrhs_, a, &lda_, ipiv_, b, &ldb_, work, &lwork, &info);
 free(work);
 free(ipiv_);
 return (lapack_int) (info<0 ? info-1 : info);
}
