zgetDiagInv.c

Go to the documentation of this file.

/*======================================================================
     Copyright (C) 2018 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"
#include <stdio.h>
#include <omp.h>

#ifdef __cplusplus
extern "C" 
#endif

typedef struct dcomplex { 
     double real; 
     double imag; 
} comp;

#ifdef MIC
__declspec( target(mic) ) void getdiaginv_mp_zgetdiaginv_(int *sizeA, int *nonzerosA, comp *a, int *ia, int *ja, double *real_invA, double *imag_invA, int *error);
#else
void getdiaginv_mp_zgetdiaginv_(int *sizeA, int *nonzerosA, comp *a, int *ia, int *ja, double *real_invA, double *imag_invA, int *error);
#endif

  

void mexFunction(int nlhs, mxArray *plhs[],
                 int nrhs, const mxArray *prhs[])
{
/*     Array information: */
      mwSize nonzerosA; /*number of nonzeros elements in A*/
      mwSize sizeAm;
      mwSize sizeAn; /*size of the matrix A (rows, and cols)  */
      int sizeAn_fort, nonzerosA_fort; /*fortran compatible matrix dimensions*/

      mwSize *ia, *ja;
      int *ia_fort, *ja_fort; /* fortran compatible indices */
      int error;
      comp* a;
      double *real, *imag;


/*-----------------------------------------------------------------------*/
/*     Check for proper number of arguments. */
      if(nrhs != 1) {
         mexErrMsgIdAndTxt ("EQuUs:zgetdiaginv:nInput", "One input required.");
      }
      else if (nlhs > 1) {
         mexErrMsgIdAndTxt ("EQuUs:zgetdiaginv:nOutput", "Too many output arguments.");
      }

/*     Validate inputs */
      
/*     Check that the input is a double.*/
      if ((mxIsDouble(prhs[0]) == 0)) {
         mexErrMsgIdAndTxt ("EQuUs:zgetdiaginv:NonDouble", "Input must be a double.");
      }
      
/*     Check if inputs are complex.*/
      if ((mxIsComplex(prhs[0]) == 0) ) {
         mexErrMsgIdAndTxt ("EQuUs:zgetdiaginv:NonComplex", "Inputs must be complex.");
      }

/*     Get the size of the input array. (MATLAB) */
      nonzerosA = mxGetNzmax(prhs[0]);
      sizeAm = mxGetM(prhs[0]);
      sizeAn = mxGetN(prhs[0]);

      ia = mxGetJc(prhs[0]);
      ja = mxGetIr(prhs[0]);
      real = mxGetPr(prhs[0]);
      imag = mxGetPi(prhs[0]);

      ia_fort = (int*)mxMalloc( (sizeAm+1)*sizeof(int) );
      ja_fort = (int*)mxMalloc( nonzerosA*sizeof(int) );
#ifdef matlab
      a = (comp*)mxMalloc( nonzerosA*sizeof(comp) );
#else
      a = (comp*)real;
#endif

      mwSize idx; 
      for( idx=1; idx<=sizeAm+1; idx++ ) {
        ia_fort[idx-1] = (int)ia[idx-1] + 1;
      }

      for( idx=1; idx<=nonzerosA; idx++ ) {
        ja_fort[idx-1] = (int)ja[idx-1] + 1;
#ifdef matlab
        a[idx-1].real = real[idx-1];
        a[idx-1].imag = imag[idx-1];
#endif
      }

      
      plhs[0] = mxCreateDoubleMatrix(sizeAn, 1, mxCOMPLEX); /* create output matrix */
      /* populate the real and imaginary part of the created array */
      real = mxGetPr(plhs[0]);
      imag = mxGetPi(plhs[0]);

      sizeAn_fort = (int)sizeAn;
      nonzerosA_fort = (int)nonzerosA;

#ifdef MIC
#pragma offload target(mic) \
    in(sizeAn_fort, nonzerosA_fort) \
     in(a:length(nonzerosA)) \
     in(ia_fort:length(sizeAn+1)) \
     in(ja_fort:length(nonzerosA)) \
     out(real:length(sizeAn_fort)) \
     out(imag:length(sizeAn_fort)) \
     out(error) 
      {
#endif

#ifdef __INTEL_COMPILER
      getdiaginv_mp_zgetdiaginv_(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, real, imag, &error);
#else
      __getdiaginv_MOD_zgetdiaginv(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, real, imag, &error);
#endif

#ifdef MIC
      }
#endif

      mxFree(ia_fort);
      mxFree(ja_fort);

/*    Check the error flag   .*/
      if ((error != 0) ) {
         char output[100];
         sprintf(output, "PARDISO solver ended with error flag:  %i\n", error);
         mexErrMsgIdAndTxt ("EQuUs:zgetdiaginv:PardisoError", output);
      }

#ifdef DEBUG
      mexPrintf("Pardiso calculations finished for real matrix\n");
#endif

}