zggev.c

Go to the documentation of this file.

/*======================================================================
     (See: http://www.netlib.org/lapack/explore-html/d9/d52/zggev_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;
      

      lapack_complex_double *a, *b, *alpha, *beta, *l_eigvec, *r_eigvec; 
/* pointers to the real and imaginary marts of the MATLAB matrices */
      double *ar, *ai, *br, *bi;
      double *alphar, *alphai, *betar, *betai;
      double *l_eigvec_r, *l_eigvec_i, *r_eigvec_r, *r_eigvec_i;  

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



/* code here */
/*-----------------------------------------------------------------------*/
/*     Check for proper number of arguments. */
      if(nrhs != 2) {
         mexErrMsgIdAndTxt ("MATLAB:zggev:nInput", "Two inputs required.");
      }
      else if (nlhs > 5) {
         mexErrMsgIdAndTxt ("MATLAB:zggev: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:zggev:NonDouble", "Inputs must be both a double.");
      }
      
/*     Check if inputs are complex.*/
      if ((mxIsComplex(prhs[0]) == 0) && (mxIsComplex(prhs[1]) == 0)) {
         mexErrMsgIdAndTxt ("MATLAB:convec:NonReal", "Inputs must be both complex");
      }

/*     Get the size of the input array. (MATLAB) */
      sizeAm = mxGetM(prhs[0]);
      sizeAn = mxGetN(prhs[0]);
      sizeBm = mxGetM(prhs[1]);
      sizeBn = mxGetN(prhs[1]);
      
      ar = mxGetPr(prhs[0]);
      br = mxGetPr(prhs[1]);
      ai = mxGetPi(prhs[0]);
      bi = mxGetPi(prhs[1]);

/* Creating output arrays */
      plhs[0] = mxCreateDoubleMatrix(sizeAm, 1, mxCOMPLEX); /* alpha*/
      plhs[1] = mxCreateDoubleMatrix(sizeAm, 1, mxCOMPLEX); /* beta*/
      plhs[2] = mxCreateDoubleMatrix(sizeAm, sizeAn, mxCOMPLEX); /* l_eigvec*/
      plhs[3] = mxCreateDoubleMatrix(sizeAm, sizeAn, mxCOMPLEX); /* r_eigvec*/


      alphar = mxGetPr(plhs[0]);
      alphai = mxGetPi(plhs[0]);
      betar = mxGetPr(plhs[1]);
      betai = mxGetPi(plhs[1]);
      l_eigvec_r = mxGetPr(plhs[2]);
      l_eigvec_i = mxGetPi(plhs[2]);
      r_eigvec_r = mxGetPr(plhs[3]);
      r_eigvec_i = mxGetPi(plhs[3]);



      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

/* Creating LAPACKE arrays */
      a = (lapack_complex_double*)malloc(sizeAm_mkl*sizeAm_mkl*sizeof(lapack_complex_double));
      if ( a == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix A.");
      }

      b = (lapack_complex_double*)malloc(sizeBm_mkl*sizeBm_mkl*sizeof(lapack_complex_double));
      if ( b == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix B.");
      }

      for( idx=1; idx<=sizeBm_mkl*sizeBm_mkl; idx++) {
        a[idx-1].real = ar[idx-1];
        a[idx-1].imag = ai[idx-1];
        b[idx-1].real = br[idx-1];
        b[idx-1].imag = bi[idx-1];
      } 

      l_eigvec = (lapack_complex_double*)malloc(sizeBm_mkl*sizeBm_mkl*sizeof(lapack_complex_double));
      if ( l_eigvec == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix l_eigvec.");
      }

      r_eigvec = (lapack_complex_double*)malloc(sizeBm_mkl*sizeBm_mkl*sizeof(lapack_complex_double));
      if ( r_eigvec == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix r_eigvec.");
      }


      alpha = (lapack_complex_double*)malloc(sizeBm_mkl*sizeof(lapack_complex_double));
      if ( alpha == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix alpha.");
      }

      beta = (lapack_complex_double*)malloc(sizeBm_mkl*sizeof(lapack_complex_double));
      if ( beta == NULL ) {
        mexErrMsgIdAndTxt ("MATLAB:zggev:malloc", "Failed to allocate space for matrix beta.");
      }
mkl_set_num_threads_local(12);
      info =  LAPACKE_zggev( LAPACK_COL_MAJOR, 'V', 'V',
                          sizeAm_mkl, a, sizeAn_mkl, b,
                          sizeBm_mkl, alpha,
                          beta, l_eigvec, sizeAm_mkl, r_eigvec,
                          sizeAm_mkl );

/* TODO_ check info */

/* Copy output arrays to MATLAB format  */

      for( idx=1; idx<=sizeAm_mkl; idx++) {
        alphar[idx-1] = alpha[idx-1].real;
        alphai[idx-1] = alpha[idx-1].imag;
        betar[idx-1] = beta[idx-1].real;
        betai[idx-1] = beta[idx-1].imag;
      } 


      for( idx=1; idx<=sizeAm_mkl*sizeAm_mkl; idx++) {
        l_eigvec_r[idx-1] = l_eigvec[idx-1].real;
        l_eigvec_i[idx-1] = l_eigvec[idx-1].imag;
        r_eigvec_r[idx-1] = r_eigvec[idx-1].real;
        r_eigvec_i[idx-1] = r_eigvec[idx-1].imag;
      } 

      free( l_eigvec );
      free( r_eigvec );
      free( alpha );
      free( beta );

}