dggev.c

Go to the documentation of this file.

/*======================================================================     
     (See: http://www.netlib.org/lapack/explore-html/d9/d52/dggev_8f.html for details)
     Copyright (C) 2016 Peter Rakyta, Ph.D.

     This program is free software: you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
     This program is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.

     You should have received a copy of the GNU General Public License
     along with this program.  If not, see http://www.gnu.org/licenses/.
     
====================================================================== */

#include "mex.h"
#define lapack_int long
#include "mkl.h"
#ifdef __cplusplus
extern "C" 
#endif
 

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
/* variable declarations here */
/*     Array information: */
      mwSize nonzerosA; /*number of nonzeros elements in A*/
      mwSize sizeAm, sizeAn; /*size of the matrix A (rows, and cols)  */
      mwSize sizeBm, sizeBn; /*size of the matrix B (rows, and cols)  */ 
      lapack_int sizeAm_mkl, sizeAn_mkl;
      lapack_int sizeBm_mkl, sizeBn_mkl;
      

      double* a, *b; /* pointers to the matrices */
      double *alphar, *alphai, *beta;
      double *l_eigvec, *r_eigvec;  

      double Work1[1];
      double *Work;
      
      lapack_int Lwork, info;



/* code here */
/*-----------------------------------------------------------------------*/
/*     Check for proper number of arguments. */
      if(nrhs != 2) {
         mexErrMsgIdAndTxt ("MATLAB:dggev:nInput", "Two inputs required.");
      }
      else if (nlhs > 5) {
         mexErrMsgIdAndTxt ("MATLAB:dggev:nOutput", "Too many output arguments.");
      }

/*     Validate inputs */
      
/*     Check that the input is a double.*/
      if ((mxIsDouble(prhs[0]) == 0) && (mxIsDouble(prhs[1]) == 0)) {
         mexErrMsgIdAndTxt ("MATLAB:dggev:NonDouble", "Inputs must be both a double.");
      }
      
/*     Check if inputs are real.*/
      if ((mxIsComplex(prhs[0]) != 0) && (mxIsComplex(prhs[1]) != 0)) {
         mexErrMsgIdAndTxt ("MATLAB:dggev:NonReal", "Inputs must be both real.");
      }

/*     Get the size of the input array. (MATLAB) */
      sizeAm = mxGetM(prhs[0]);
      sizeAn = mxGetN(prhs[0]);
      sizeBm = mxGetM(prhs[1]);
      sizeBn = mxGetN(prhs[1]);
      
      a = mxGetPr(prhs[0]);
      b = mxGetPr(prhs[1]);

/* Creating output arrays */
      plhs[0] = mxCreateDoubleMatrix(sizeAm, 1, mxREAL); /* alphar*/
      plhs[1] = mxCreateDoubleMatrix(sizeAm, 1, mxREAL); /* alphai*/
      plhs[2] = mxCreateDoubleMatrix(sizeAm, 1, mxREAL); /* beta*/
      plhs[3] = mxCreateDoubleMatrix(sizeAm, sizeAn, mxREAL); /* l_eigvec*/
      plhs[4] = mxCreateDoubleMatrix(sizeAm, sizeAn, mxREAL); /* r_eigvec*/


      alphar = mxGetPr(plhs[0]);
      alphai = mxGetPr(plhs[1]);
      beta = mxGetPr(plhs[2]);
      l_eigvec = mxGetPr(plhs[3]);
      r_eigvec = mxGetPr(plhs[4]);



      sizeAm_mkl = (lapack_int)sizeAm;
      sizeAn_mkl = (lapack_int)sizeAn;
      sizeBm_mkl = (lapack_int)sizeBm;
      sizeBn_mkl = (lapack_int)sizeBn;

      if ( sizeAm != sizeAn || sizeAm != sizeBm || sizeAm != sizeBn ) {
         mexErrMsgIdAndTxt ("MATLAB:dggev:MatrixSize", "Inputs must be square matrices and must have the same dimensions.");
      }      

#ifdef DEBUG
      mexPrintf("Size of the matrix: %d\n", sizeAm_mkl);
#endif

      info =  LAPACKE_dggev( LAPACK_COL_MAJOR, 'V', 'V',
                          sizeAm_mkl, a, sizeAn_mkl, b,
                          sizeBm_mkl, alphar, alphai,
                          beta, l_eigvec, sizeAm_mkl, r_eigvec,
                          sizeAm_mkl );

/* TODO_ check info */

}