/*
 * This mex function is called from MATLAB using the form
 *
 *    [V, D] = dconquer(A);
 *
 * Use the divide and conquer algorithm to compute the
 * eigenvalues and corresponding eigenvectors of a
 * symmetric matrix. The LAPACK function "dsyevd"
 * coded in FORTRAN 77 performs the job.
*/

#include "mex.h"
#include "lapack.h"
#include "string.h"

/* nlhs is the number of left-hand sides (2).
 *
 * plhs is an array of pointers to the variables V and D on the
 * left-hand side
 *
 * nrhs is the number of input arguments (1)
 *
 * prhs[0] points at the matrix A.
*/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
   long m;
   int i, increment;
   /* "V" - compute eigenvalues and eigenvectors.
    * "U" - use the upper triangular portion of A, and
    * assume the triangle above the diagonal is reflected
    * across the diagonal.
   */
   char *jobz = "V", *uplo = "U";
   ptrdiff_t n1, lda1, lwork1, liwork1, info1;
   ptrdiff_t *n, *lda, *lwork, *liwork, *info, *iwork;
   double *A, *w, *work, *eigenvalues;
   mxArray *Awork, *w1, *work1, *iwork1, *A1;
   
   /* Check for the proper number of arguments, and make
    * sure the matrix is square.
   */
   if (nrhs != 1)
       mexErrMsgIdAndTxt("lapack_dsyevd:Error",
                         "The symmetric matrix is the only input argument.");
   else if (nlhs != 2)
      mexErrMsgIdAndTxt("lapack_dsyevd:Error",
                        "There must be two output arguments.");
   n1 = mxGetM(prhs[0]);
   n = &n1;
   m = mxGetN(prhs[0]);
   if (m != n1)
      mexErrMsgIdAndTxt("lapack_dsyevd:Error",
                        "The matrix must be square.");

   /* dsyevd requires that its arugments lwork and liwork have
    * the following minimum values. Optimal values can be computed
    * by calling dsyevd first with lwork = -1 and liwork = -1.
    * In this case, the function returns optimal values in
    * work[0] and iwork[0], respectively. We do not use this
    * feature.
   */
   lwork1 = 1 + 6*n1 + 2*n1*n1;
   lwork = &lwork1;
   liwork1 = 3 + 5*n1;
   liwork = &liwork1;
   info1 = -1;
   info = &info1;

   /* Copy the input matrix to the output matrix. */
   plhs[0] = mxDuplicateArray(prhs[0]);
   A = mxGetPr(plhs[0]);
      
   /* The leading dimension is *n */
   lda = n;

   /* The vector w will contain the *n eigenvalues after
    * calling dsyevd.
   */
   w1 = mxCreateDoubleMatrix(*n,1,mxREAL);
   w = (double *)mxGetPr(w1);
   
   /* Allocate vector work with *lwork entries,
    * and allocate vector iwork with *liwork entries.
   */
   work1 = mxCreateDoubleMatrix(*lwork,1,mxREAL);
   work = (double *)mxGetPr(work1);
   iwork1 = mxCreateNumericMatrix(*liwork, 1, mxINT64_CLASS, mxREAL);
   iwork = ( ptrdiff_t *)mxGetData(iwork1);
   
   /* Call dsyevd. Divide and conquer will be used since we are asking
    * for both the eigenvalues and the eigenvectors.
   */
   dsyevd(jobz, uplo, n, A, lda, w, work, lwork, iwork, liwork, info);
   
   /* The *n eigenvalues are in vector w, but we want to return
    * a diagonal matrix of eigenvalues. Allocate the matrix.
   */
   plhs[1] = mxCreateDoubleMatrix(*n,*n, mxREAL);
   eigenvalues = (double *)mxGetPr(plhs[1]);
   
   /* FORTRAN and MATLAB store matrices in column-major format.
    * Copy the eigenvalues to the diagonal.
   */
   increment = *n + 1;
   for (i = 0;i < *n;i++)
      eigenvalues[i*increment] = w[i];
}
