dgetPartialInv.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 "mex.h"
#include <stdio.h>
#include <omp.h>

#ifdef __cplusplus
extern "C" 
#endif

#ifdef MIC
__declspec( target(mic) ) void getpartialinv_mp_dgetpartialinv_(int *sizeA, int *nonzerosA, double *a, int *ia, int *ja, int *sizeInv, double *invA, int *error);
#else
void getpartialinv_mp_dgetpartialinv_(int *sizeA, int *nonzerosA, double *a, int *ia, int *ja, int *sizeInv, double *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;
      double* a, *invA_fort;


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

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

/*     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]);
      a = mxGetPr(prhs[0]);

      ia_fort = (int*)malloc( (sizeAm+1)*sizeof(int) );
      ja_fort = (int*)malloc( nonzerosA*sizeof(int) );

      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;
      }

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


      plhs[0] = mxCreateDoubleMatrix(sizeInv, sizeInv, mxREAL); /* invA */
      
      invA_fort = (double *)mxGetData(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(invA_fort:length(sizeInv_fort*sizeInv_fort)) \
     out(error) 
      {
#endif

#ifdef __INTEL_COMPILER
      getpartialinv_mp_dgetpartialinv_(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, &sizeInv_fort, invA_fort, &error);
#else
      __getpartialinv_MOD_dgetpartialinv(&sizeAn_fort, &nonzerosA_fort, a, ia_fort, ja_fort, &sizeInv_fort, invA_fort, &error);
#endif

#ifdef MIC
      }
#endif

      free(ia_fort);
      free(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:dgetPartialInv:PardisoError", output);
      }

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

}