zgetPartialInv.c

Go to the documentation of this file.

/*======================================================================
     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 <unistd.h>
#include "mex.h"
#ifdef __cplusplus
extern "C" 
#endif

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

#ifdef MIC
__declspec( target(mic) ) void getpartialinv_mp_zgetpartialinv_(int *sizeA, int *nonzerosA, comp *a, int *ia, int *ja, int *sizeInv, double *real_invA, double *imag_invA, int *error);
#else
void getpartialinv_mp_zgetpartialinv_(int *sizeA, int *nonzerosA, comp *a, int *ia, int *ja, int *sizeInv, 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)  */
      mwSize sizeInv; /*Size of the Schur complement to be calculated*/
      int sizeAn_fort, nonzerosA_fort, sizeInv_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 != 2) {
         mexErrMsgIdAndTxt ("EQuUs:zgetPartialInv:nInput", "Two inputs required.");
      }
      else if (nlhs > 1) {
         mexErrMsgIdAndTxt ("EQuUs:zgetPartialInv:nOutput", "Too many output arguments.");
      }

/*     Validate inputs */
      
/*     Check that the input is a double.*/
      if ((mxIsDouble(prhs[0]) == 0)) {
         mexErrMsgIdAndTxt ("EQuUs:zgetPartialInv:NonDouble", "Input must be a double.");
      }
      
/*     Check if inputs are complex.*/
      if ((mxIsComplex(prhs[0]) == 0) ) {
         mexErrMsgIdAndTxt ("EQuUs:zgetPartialInv: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
      }

/*     getting the size of the partial inverse */
      sizeInv = *mxGetPr(prhs[1]);


      plhs[0] = mxCreateDoubleMatrix(sizeInv, sizeInv, 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;
      sizeInv_fort = (int)sizeInv;

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

#ifdef __INTEL_COMPILER
      getpartialinv_mp_zgetpartialinv_(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, &sizeInv_fort, real, imag, &error);
#else
      __getpartialinv_MOD_zgetpartialinv(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, &sizeInv_fort, real, imag, &error);
#endif

#ifdef MIC
      }
#endif

      mxFree(ia_fort);
      mxFree(ja_fort);
      mxFree(a);

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

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

}