C++ BigMatrix类代码示例

/* Prepares the a matrix based on random sample of examples for modelling. For
   each continuous variable, copies only the in-sample indices from asave to a.
   Data for categorical variables are not copied, as they are stored in x.
   This function should only be called if there are any continuous variables. */
SEXP moda(SEXP asaveP, SEXP aP, SEXP insampP) {
    // Initialize function arguments.
    BigMatrix *asave = (BigMatrix*)R_ExternalPtrAddr(asaveP);
    BigMatrix *a = (BigMatrix*)R_ExternalPtrAddr(aP);
    MatrixAccessor<int> asaveAcc(*asave);
    MatrixAccessor<int> aAcc(*a);
    int *asaveCol, *aCol;
    int *insamp = INTEGER(insampP);
    
    // Set up working variables.
    index_type nCols = asave->ncol();
    index_type nRows = asave->nrow();
    index_type i, ja, jb;
    
    // For each numerical variable, move all the in-sample data to the top rows
    // of a.
    for (i = 0; i < nCols; i++) {
        asaveCol = asaveAcc[i];
        aCol = aAcc[i];
        for (ja = 0, jb = 0; ja < nRows; ja++) {
            if (insamp[asaveCol[ja] - 1] >= 1) {
                aCol[jb++] = asaveCol[ja];
            }
        }
    }
    return R_NilValue;
}

 SepMatrixAccessor( BigMatrix &bm)
 {
   _ppMat = reinterpret_cast<T**>(bm.matrix());
   _rowOffset = bm.row_offset();
   _colOffset = bm.col_offset();
   _totalRows = bm.nrow();
 }

    SEXP read_partial_bignifti_data(SEXP nim_addr, SEXP big_addr, \
                                    SEXP rowIndices, SEXP colIndices, \
                                    SEXP totalVoxs) {
        // Get nifti object pointer
        nifti_image *pnim = (nifti_image *)R_ExternalPtrAddr(nim_addr);

        // Get nifti data
        if (nifti_image_load(pnim) < 0) {
            nifti_image_free(pnim);
            error("Could not load nifti data");
        }

        // Save data to big matrix object
        BigMatrix *pMat = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(big_addr));

        if (pMat->separated_columns()) {
            switch (pMat->matrix_type()) {
            case 1:
                return Read_Partial_Nifti_To_BigMatrix_Step2<char>(pnim, SepMatrixAccessor<char>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 2:
                return Read_Partial_Nifti_To_BigMatrix_Step2<short>(pnim, SepMatrixAccessor<short>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 4:
                return Read_Partial_Nifti_To_BigMatrix_Step2<int>(pnim, SepMatrixAccessor<int>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 8:
                return Read_Partial_Nifti_To_BigMatrix_Step2<double>(pnim, SepMatrixAccessor<double>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            }
        }
        else {
            switch (pMat->matrix_type()) {
            case 1:
                return Read_Partial_Nifti_To_BigMatrix_Step2<char>(pnim, MatrixAccessor<char>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 2:
                return Read_Partial_Nifti_To_BigMatrix_Step2<short>(pnim, MatrixAccessor<short>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 4:
                return Read_Partial_Nifti_To_BigMatrix_Step2<int>(pnim, MatrixAccessor<int>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            case 8:
                return Read_Partial_Nifti_To_BigMatrix_Step2<double>(pnim, MatrixAccessor<double>(*pMat), \
                        rowIndices, colIndices, totalVoxs);
                break;
            }
        }

        error("failed to identify big matrix type");
    }

 SEXP BigSumMain(SEXP addr, SEXP cols, SEXP rows) {
     SEXP ret = R_NilValue;
     ret = PROTECT(NEW_NUMERIC(1));
     double *pRet = NUMERIC_DATA(ret);
     
     BigMatrix *pMat = (BigMatrix*)R_ExternalPtrAddr(addr);
     
     if (pMat->separated_columns()) {
         switch (pMat->matrix_type()) {
             case 1:
                 BigSum<char, SepMatrixAccessor<char> >(
                 pMat, cols, rows, pRet);
                 break;
             case 2:
                 BigSum<short, SepMatrixAccessor<short> >(
                 pMat, cols, rows, pRet);
                 break;
             case 4:
                 BigSum<int, SepMatrixAccessor<int> >(
                 pMat, cols, rows, pRet);
                 break;
             case 8:
                 BigSum<double, SepMatrixAccessor<double> >(
                 pMat, cols, rows, pRet);
                 break;
         }
     }
     else {
         switch (pMat->matrix_type()) {
             case 1:
                 BigSum<char, MatrixAccessor<char> >(
                 pMat, cols, rows, pRet);
                 break;
             case 2:
                 BigSum<short, MatrixAccessor<short> >(
                 pMat, cols, rows, pRet);
                 break;
             case 4:
                 BigSum<int, MatrixAccessor<int> >(
                 pMat, cols, rows, pRet);
                 break;
             case 8:
                 BigSum<double, MatrixAccessor<double> >(
                 pMat, cols, rows, pRet);
                 break;
         }
     }
     
     UNPROTECT(1);
     return(ret);
 }

  SEXP CIPTMatrix(SEXP inAddr)
  {
    BigMatrix *pInMat = reinterpret_cast<BigMatrix*>(
      R_ExternalPtrAddr(inAddr));
    
    // Not sure if there is a better way to do these function calls
    if (pInMat->separated_columns()) {
      // Need method for separated_columns
      //CALL_IPT_2(SepMatrixAccessor)
    }
    else
    {
      CALL_IPT_2(MatrixAccessor)
    }

    return R_NilValue;
  }

 MatrixAccessor( BigMatrix &bm )
 {
   _pMat = reinterpret_cast<T*>(bm.matrix());
   _totalRows = bm.total_rows();
   _totalCols = bm.total_columns();
   _rowOffset = bm.row_offset();
   _colOffset = bm.col_offset();
   _nrow = bm.nrow();
   _ncol = bm.ncol();
 }

/* Pointer utility, returns a double pointer for either a BigMatrix or a
 * standard R matrix.
 */
double *
make_double_ptr (SEXP matrix, SEXP isBigMatrix)
{
  double *matrix_ptr;

  if (LOGICAL_VALUE (isBigMatrix) == (Rboolean) TRUE)   // Big Matrix
    {
      SEXP address = GET_SLOT (matrix, install ("address"));
      BigMatrix *pbm =
        reinterpret_cast < BigMatrix * >(R_ExternalPtrAddr (address));
      if (!pbm)
        return (NULL);

      // Check that have acceptable big.matrix
      if (pbm->row_offset () > 0 && pbm->ncol () > 1)
        {
          std::string errMsg =
            string ("sub.big.matrix objects cannoth have row ") +
            string
            ("offset greater than zero and number of columns greater than 1");
          Rf_error (errMsg.c_str ());
          return (NULL);
        }

      index_type offset = pbm->nrow () * pbm->col_offset ();
      matrix_ptr = reinterpret_cast < double *>(pbm->matrix ()) + offset;
    }
  else                          // Regular R Matrix
    {
      matrix_ptr = NUMERIC_DATA (matrix);
    }

  return (matrix_ptr);
};

SEXP ComputePvalsMain(SEXP Rinmat, SEXP Routmat, SEXP Routcol) {
    BigMatrix *inMat = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(Rinmat));
    BigMatrix *outMat = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(Routmat));
    double outCol = NUMERIC_DATA(Routcol)[0];
    
    if (inMat->separated_columns() != outMat->separated_columns())
        Rf_error("all big matrices are not the same column separated type");
    if (inMat->matrix_type() != outMat->matrix_type())
        Rf_error("all big matrices are not the same matrix type");
    if (inMat->ncol() != outMat->nrow())
        Rf_error("inMat # of cols must be the same as outMat # of rows");
    
    CALL_BIGFUNCTION_ARGS_THREE(ComputePvals, inMat, outMat, outCol)
    return(ret);
}

template<typename MatrixType> void diagonalmatrices(const MatrixType& m)
{
  typedef typename MatrixType::Index Index;
  typedef typename MatrixType::Scalar Scalar;
  enum { Rows = MatrixType::RowsAtCompileTime, Cols = MatrixType::ColsAtCompileTime };
  typedef Matrix<Scalar, Rows, 1> VectorType;
  typedef Matrix<Scalar, 1, Cols> RowVectorType;
  typedef Matrix<Scalar, Rows, Rows> SquareMatrixType;
  typedef DiagonalMatrix<Scalar, Rows> LeftDiagonalMatrix;
  typedef DiagonalMatrix<Scalar, Cols> RightDiagonalMatrix;
  typedef Matrix<Scalar, Rows==Dynamic?Dynamic:2*Rows, Cols==Dynamic?Dynamic:2*Cols> BigMatrix;
  Index rows = m.rows();
  Index cols = m.cols();

  MatrixType m1 = MatrixType::Random(rows, cols),
             m2 = MatrixType::Random(rows, cols);
  VectorType v1 = VectorType::Random(rows),
             v2 = VectorType::Random(rows);
  RowVectorType rv1 = RowVectorType::Random(cols),
             rv2 = RowVectorType::Random(cols);
  LeftDiagonalMatrix ldm1(v1), ldm2(v2);
  RightDiagonalMatrix rdm1(rv1), rdm2(rv2);
  
  Scalar s1 = internal::random<Scalar>();

  SquareMatrixType sq_m1 (v1.asDiagonal());
  VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
  sq_m1 = v1.asDiagonal();
  VERIFY_IS_APPROX(sq_m1, v1.asDiagonal().toDenseMatrix());
  SquareMatrixType sq_m2 = v1.asDiagonal();
  VERIFY_IS_APPROX(sq_m1, sq_m2);
  
  ldm1 = v1.asDiagonal();
  LeftDiagonalMatrix ldm3(v1);
  VERIFY_IS_APPROX(ldm1.diagonal(), ldm3.diagonal());
  LeftDiagonalMatrix ldm4 = v1.asDiagonal();
  VERIFY_IS_APPROX(ldm1.diagonal(), ldm4.diagonal());
  
  sq_m1.block(0,0,rows,rows) = ldm1;
  VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
  sq_m1.transpose() = ldm1;
  VERIFY_IS_APPROX(sq_m1, ldm1.toDenseMatrix());
  
  Index i = internal::random<Index>(0, rows-1);
  Index j = internal::random<Index>(0, cols-1);
  
  VERIFY_IS_APPROX( ((ldm1 * m1)(i,j))  , ldm1.diagonal()(i) * m1(i,j) );
  VERIFY_IS_APPROX( ((ldm1 * (m1+m2))(i,j))  , ldm1.diagonal()(i) * (m1+m2)(i,j) );
  VERIFY_IS_APPROX( ((m1 * rdm1)(i,j))  , rdm1.diagonal()(j) * m1(i,j) );
  VERIFY_IS_APPROX( ((v1.asDiagonal() * m1)(i,j))  , v1(i) * m1(i,j) );
  VERIFY_IS_APPROX( ((m1 * rv1.asDiagonal())(i,j))  , rv1(j) * m1(i,j) );
  VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * m1)(i,j))  , (v1+v2)(i) * m1(i,j) );
  VERIFY_IS_APPROX( (((v1+v2).asDiagonal() * (m1+m2))(i,j))  , (v1+v2)(i) * (m1+m2)(i,j) );
  VERIFY_IS_APPROX( ((m1 * (rv1+rv2).asDiagonal())(i,j))  , (rv1+rv2)(j) * m1(i,j) );
  VERIFY_IS_APPROX( (((m1+m2) * (rv1+rv2).asDiagonal())(i,j))  , (rv1+rv2)(j) * (m1+m2)(i,j) );

  BigMatrix big;
  big.setZero(2*rows, 2*cols);
  
  big.block(i,j,rows,cols) = m1;
  big.block(i,j,rows,cols) = v1.asDiagonal() * big.block(i,j,rows,cols);
  
  VERIFY_IS_APPROX((big.block(i,j,rows,cols)) , v1.asDiagonal() * m1 );
  
  big.block(i,j,rows,cols) = m1;
  big.block(i,j,rows,cols) = big.block(i,j,rows,cols) * rv1.asDiagonal();
  VERIFY_IS_APPROX((big.block(i,j,rows,cols)) , m1 * rv1.asDiagonal() );
  
  
  // scalar multiple
  VERIFY_IS_APPROX(LeftDiagonalMatrix(ldm1*s1).diagonal(), ldm1.diagonal() * s1);
  VERIFY_IS_APPROX(LeftDiagonalMatrix(s1*ldm1).diagonal(), s1 * ldm1.diagonal());
  
  VERIFY_IS_APPROX(m1 * (rdm1 * s1), (m1 * rdm1) * s1);
  VERIFY_IS_APPROX(m1 * (s1 * rdm1), (m1 * rdm1) * s1);
  
  // Diagonal to dense
  sq_m1.setRandom();
  sq_m2 = sq_m1;
  VERIFY_IS_APPROX( (sq_m1 += (s1*v1).asDiagonal()), sq_m2 += (s1*v1).asDiagonal().toDenseMatrix() );
  VERIFY_IS_APPROX( (sq_m1 -= (s1*v1).asDiagonal()), sq_m2 -= (s1*v1).asDiagonal().toDenseMatrix() );
  VERIFY_IS_APPROX( (sq_m1 = (s1*v1).asDiagonal()), (s1*v1).asDiagonal().toDenseMatrix() );
}

  SEXP CtransposeMatrix(SEXP inAddr, SEXP outAddr, SEXP rowInds, SEXP colInds, 
    SEXP typecast_warning)
  {
    BigMatrix *pInMat = reinterpret_cast<BigMatrix*>(
      R_ExternalPtrAddr(inAddr));
    BigMatrix *pOutMat = reinterpret_cast<BigMatrix*>(
      R_ExternalPtrAddr(outAddr));
    
    if ((pOutMat->matrix_type() < pInMat->matrix_type()) & 
      (Rcpp::as<bool>(typecast_warning) == (Rboolean)TRUE))
    {
      string type_names[9] = {
        "", "char", "short", "", "integer", "", "", "", "double"};
      
      std::string warnMsg = string("Assignment will down cast from ") + 
        type_names[pInMat->matrix_type()] + string(" to ") + 
        type_names[pOutMat->matrix_type()] + string("\n") + 
        string("Hint: To remove this warning type: ") + 
        string("options(bigmemory.typecast.warning=FALSE)");
      Rf_warning(warnMsg.c_str());
    }
    
    // Not sure if there is a better way to do these function calls
    if (pInMat->separated_columns() && pOutMat->separated_columns()) {
      CALL_transpose_1(SepMatrixAccessor, SepMatrixAccessor)
    }
    else if(pInMat->separated_columns() && !(pOutMat->separated_columns()))
    {
      CALL_transpose_1(SepMatrixAccessor, MatrixAccessor)
    }
    else if(!(pInMat->separated_columns()) && pOutMat->separated_columns())
    {
      CALL_transpose_1(MatrixAccessor, SepMatrixAccessor)
    }
    else
    {
      CALL_transpose_1(MatrixAccessor, MatrixAccessor)
    }

    return R_NilValue;
  }

SEXP binit1BigMatrix(SEXP x, SEXP col, SEXP breaks)
{
  BigMatrix *pMat =  reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(x));
  if (pMat->separated_columns())
  {
    switch (pMat->matrix_type())
    {
      case 1:
        return CBinIt1<char>(SepMatrixAccessor<char>(*pMat),
          pMat->nrow(), col, breaks);
      case 2:
        return CBinIt1<short>(SepMatrixAccessor<short>(*pMat),
          pMat->nrow(), col, breaks);
      case 4:
        return CBinIt1<int>(SepMatrixAccessor<int>(*pMat),
          pMat->nrow(), col, breaks);
      case 8:
        return CBinIt1<double>(SepMatrixAccessor<double>(*pMat),
          pMat->nrow(), col, breaks);
    }
  }
  else
  {
    switch (pMat->matrix_type())
    {
      case 1:
        return CBinIt1<char>(MatrixAccessor<char>(*pMat),
          pMat->nrow(), col, breaks);
      case 2:
        return CBinIt1<short>(MatrixAccessor<short>(*pMat),
          pMat->nrow(), col, breaks);
      case 4:
        return CBinIt1<int>(MatrixAccessor<int>(*pMat),
          pMat->nrow(), col, breaks);
      case 8:
        return CBinIt1<double>(MatrixAccessor<double>(*pMat),
          pMat->nrow(), col, breaks);
    }
  }
  return R_NilValue;
}

    SEXP write_bignifti(SEXP header, SEXP big_addr, SEXP indices, SEXP outfile) {
        SEXP Rdim, Rdatatype;

        // Load big matrix
        BigMatrix *pMat = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(big_addr));

        // Get dim
        PROTECT(Rdim = GET_LIST_ELEMENT(header, "dim"));
        if (Rdim == R_NilValue)
            error("header must have a proper dim (dimension) attribute");
        if (GET_LENGTH(Rdim) != 4)
            error("header must have a 4D dim (dimension) attribute");

        // Get datatype
        PROTECT(Rdatatype = GET_LIST_ELEMENT(header, "datatype"));
        if (Rdatatype == R_NilValue) {
            int datatype;
            switch(pMat->matrix_type()) {
            case 1:     // char
                datatype = DT_INT8;
                break;
            case 2:     // short
                datatype = DT_INT16;
                break;
            case 4:     // int
                datatype = DT_INT32;
                break;
            case 8:     // double
                datatype = DT_FLOAT64;
                break;
            default:
                error("unrecognized big matrix data type");
                break;
            }
            Rdatatype = int_to_SEXP(datatype);
        }
        UNPROTECT(1);

        // Load nifti object
        nifti_image *pnim = create_nifti_image(header, Rdim, Rdatatype, outfile);

        if (pMat->separated_columns()) {
            switch (pMat->matrix_type()) {
            case 1:
                Write_BigMatrix_To_Nifti_Step2<char, SepMatrixAccessor<char> >(pnim, pMat, indices);
                break;
            case 2:
                Write_BigMatrix_To_Nifti_Step2<short, SepMatrixAccessor<short> >(pnim, pMat, indices);
                break;
            case 4:
                Write_BigMatrix_To_Nifti_Step2<int, SepMatrixAccessor<int> >(pnim, pMat, indices);
                break;
            case 8:
                Write_BigMatrix_To_Nifti_Step2<double, SepMatrixAccessor<double> >(pnim, pMat, indices);
                break;
            }
        }
        else {
            switch (pMat->matrix_type()) {
            case 1:
                Write_BigMatrix_To_Nifti_Step2<char, MatrixAccessor<char> >(pnim, pMat, indices);
                break;
            case 2:
                Write_BigMatrix_To_Nifti_Step2<short, MatrixAccessor<short> >(pnim, pMat, indices);
                break;
            case 4:
                Write_BigMatrix_To_Nifti_Step2<int, MatrixAccessor<int> >(pnim, pMat, indices);
                break;
            case 8:
                Write_BigMatrix_To_Nifti_Step2<double, MatrixAccessor<double> >(pnim, pMat, indices);
                break;
            }
        }

        if (!nifti_nim_is_valid(pnim, 1))
            error("data seems invalid");

        if (pnim!=NULL)
            nifti_image_write(pnim);
        else
            error("pnim was NULL");

        nifti_image_free(pnim);

        return R_NilValue;
    }

SEXP kmeansMatrixEuclid(MatrixType x, index_type n, index_type m,
                  SEXP pcen, SEXP pclust, SEXP pclustsizes,
                  SEXP pwss, SEXP itermax)
{

  index_type j, col, nchange;

  int maxiters = Rf_asInteger(itermax);
  SEXP Riter;
  Rf_protect(Riter = Rf_allocVector(INTSXP, 1));
  int *iter = INTEGER(Riter);
  iter[0] = 0;

  BigMatrix *pcent = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(pcen));
  MatrixAccessor<double> cent(*pcent);
  BigMatrix *Pclust = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(pclust));
  MatrixAccessor<int> clust(*Pclust);
  BigMatrix *Pclustsizes = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(pclustsizes));
  MatrixAccessor<double> clustsizes(*Pclustsizes);
  BigMatrix *Pwss = reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(pwss));
  MatrixAccessor<double> ss(*Pwss);

  int k = (int) pcent->nrow();                // number of clusters
  int cl, bestcl, oldcluster, newcluster;
  int done = 0;

  double temp;
  vector<double> d(k);                        // Vector of distances, internal only.
  vector<double> temp1(k);
  vector<vector<double> > tempcent(m, temp1); // For copy of global centroids k x m

  // At this point I can use [][] to access things, with ss[0][cl]
  // being used for the vectors, for example.
  // Before starting the loop, we only have cent (centers) as passed into the function.
  // Calculate clust and clustsizes, then update cent as centroids.
  
  for (cl=0; cl<k; cl++) clustsizes[0][cl] = 0.0;
  for (j=0; j<n; j++) {
    bestcl = 0;
    for (cl=0; cl<k; cl++) {
      d[cl] = 0.0;
      for (col=0; col<m; col++) {
        temp = (double)x[col][j] - cent[col][cl];
        d[cl] += temp * temp;
      }
      if (d[cl]<d[bestcl]) bestcl = cl;
    }
    clust[0][j] = bestcl + 1;          // Saving the R cluster number, not the C index.
    clustsizes[0][bestcl]++;
    for (col=0; col<m; col++)
      tempcent[col][bestcl] += (double)x[col][j];
  }
  for (cl=0; cl<k; cl++)
    for (col=0; col<m; col++)
      cent[col][cl] = tempcent[col][cl] / clustsizes[0][cl];

  do {

    nchange = 0;
    for (j=0; j<n; j++) { // For each of my points, this is offset from hash position

      oldcluster = clust[0][j] - 1;
      bestcl = 0;
      for (cl=0; cl<k; cl++) {         // Consider each of the clusters
        d[cl] = 0.0;                   // We'll get the distance to this cluster.
        for (col=0; col<m; col++) {    // Loop over the dimension of the data
          temp = (double)x[col][j] - cent[col][cl];
          d[cl] += temp * temp;
        }
        if (d[cl]<d[bestcl]) bestcl = cl;
      } // End of looking over the clusters for this j

      if (d[bestcl] < d[oldcluster]) {           // MADE A CHANGE!
        newcluster = bestcl;
        clust[0][j] = newcluster + 1;
        nchange++;
        clustsizes[0][newcluster]++;
        clustsizes[0][oldcluster]--;
        for (col=0; col<m; col++) {
          cent[col][oldcluster] += ( cent[col][oldcluster] - (double)x[col][j] ) / clustsizes[0][oldcluster];
          cent[col][newcluster] += ( (double)x[col][j] - cent[col][newcluster] ) / clustsizes[0][newcluster];
        }
      }

    } // End of this pass over my points.

    iter[0]++;
    if ( (nchange==0) || (iter[0]>=maxiters) ) done = 1;

  } while (done==0);

  // Collect the sums of squares now that we're done.
  for (cl=0; cl<k; cl++) ss[0][cl] = 0.0;
  for (j=0; j<n; j++) {
    for (col=0; col<m; col++) {
      cl = clust[0][j]-1;
      temp = (double)x[col][j] - cent[col][cl];
      ss[0][cl] += temp * temp;
    }
  }
//.........这里部分代码省略.........

SEXP kmeansBigMatrix(SEXP x, SEXP cen, SEXP clust, SEXP clustsizes,
                     SEXP wss, SEXP itermax, SEXP dist)
{
  BigMatrix *pMat =  reinterpret_cast<BigMatrix*>(R_ExternalPtrAddr(x));
  int dist_calc = INTEGER(dist)[0];
  if (dist_calc == 0)
  {
    if (pMat->separated_columns())
    {
      switch (pMat->matrix_type())
      {
        case 1:
          return kmeansMatrixEuclid<char>(SepMatrixAccessor<char>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 2:
          return kmeansMatrixEuclid<short>(SepMatrixAccessor<short>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 4:
          return kmeansMatrixEuclid<int>(SepMatrixAccessor<int>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 8:
          return kmeansMatrixEuclid<double>(SepMatrixAccessor<double>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
      }
    }
    else
    {
      switch (pMat->matrix_type())
      {
        case 1:
          return kmeansMatrixEuclid<char>(MatrixAccessor<char>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 2:
          return kmeansMatrixEuclid<short>(MatrixAccessor<short>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 4:
          return kmeansMatrixEuclid<int>(MatrixAccessor<int>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 8:
          return kmeansMatrixEuclid<double>(MatrixAccessor<double>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
      }
    }
  }
  else
  {
    if (pMat->separated_columns())
    {
      switch (pMat->matrix_type())
      {
        case 1:
          return kmeansMatrixCosine<char>(SepMatrixAccessor<char>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 2:
          return kmeansMatrixCosine<short>(SepMatrixAccessor<short>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 4:
          return kmeansMatrixCosine<int>(SepMatrixAccessor<int>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 8:
           return kmeansMatrixCosine<double>(SepMatrixAccessor<double>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
      }
    }
    else
    {
      switch (pMat->matrix_type())
      {
        case 1:
          return kmeansMatrixCosine<char>(MatrixAccessor<char>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 2:
          return kmeansMatrixCosine<short>(MatrixAccessor<short>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 4:
          return kmeansMatrixCosine<int>(MatrixAccessor<int>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
        case 8:
           return kmeansMatrixCosine<double>(MatrixAccessor<double>(*pMat),
            pMat->nrow(), pMat->ncol(), cen, clust, clustsizes, wss, itermax);
      }
    }
  }
  return R_NilValue;
}

// [[Rcpp::export]]
SEXP CDeepCopy(SEXP inAddr, SEXP outAddr, SEXP rowInds, SEXP colInds,
               SEXP typecast_warning)
{

#define CALL_DEEP_COPY_2(IN_CTYPE, IN_ACCESSOR, OUT_ACCESSOR) \
    switch(pOutMat->matrix_type()) \
    { \
      case 1: \
        DeepCopy<IN_CTYPE, IN_ACCESSOR<IN_CTYPE>, char, OUT_ACCESSOR<char> >( \
          pInMat, pOutMat, rowInds, colInds); \
        break; \
      case 2: \
        DeepCopy<IN_CTYPE, IN_ACCESSOR<IN_CTYPE>, short, OUT_ACCESSOR<short> >( \
          pInMat, pOutMat, rowInds, colInds); \
        break; \
      case 4: \
        DeepCopy<IN_CTYPE, IN_ACCESSOR<IN_CTYPE>, int, OUT_ACCESSOR<int> >( \
          pInMat, pOutMat, rowInds, colInds); \
        break; \
      case 8: \
        DeepCopy<IN_CTYPE, IN_ACCESSOR<IN_CTYPE>, double, OUT_ACCESSOR<double> >( \
          pInMat, pOutMat, rowInds, colInds); \
        break; \
    }

#define CALL_DEEP_COPY_1(IN_ACCESSOR, OUT_ACCESSOR) \
    switch(pInMat->matrix_type()) \
    { \
      case 1: \
        CALL_DEEP_COPY_2(char, IN_ACCESSOR, OUT_ACCESSOR) \
        break; \
      case 2: \
        CALL_DEEP_COPY_2(short, IN_ACCESSOR, OUT_ACCESSOR) \
        break; \
      case 4: \
        CALL_DEEP_COPY_2(int, IN_ACCESSOR, OUT_ACCESSOR) \
        break; \
      case 8: \
        CALL_DEEP_COPY_2(double, IN_ACCESSOR, OUT_ACCESSOR) \
        break; \
    }

    BigMatrix *pInMat = reinterpret_cast<BigMatrix*>(
                            R_ExternalPtrAddr(inAddr));
    BigMatrix *pOutMat = reinterpret_cast<BigMatrix*>(
                             R_ExternalPtrAddr(outAddr));

    if ((pOutMat->matrix_type() < pInMat->matrix_type()) &
            (LOGICAL_VALUE(typecast_warning) == (Rboolean)TRUE))
    {
        string type_names[9] = {
            "", "char", "short", "", "integer", "", "", "", "double"
        };

        std::string warnMsg = string("Assignment will down cast from ") +
                              type_names[pInMat->matrix_type()] + string(" to ") +
                              type_names[pOutMat->matrix_type()] + string("\n") +
                              string("Hint: To remove this warning type: ") +
                              string("options(bigmemory.typecast.warning=FALSE)");
        Rf_warning(warnMsg.c_str());
    }

    // Not sure if there is a better way to do these function calls
    if (pInMat->separated_columns() && pOutMat->separated_columns()) {
        CALL_DEEP_COPY_1(SepMatrixAccessor, SepMatrixAccessor)
    }
    else if(pInMat->separated_columns() && !(pOutMat->separated_columns()))
    {
        CALL_DEEP_COPY_1(SepMatrixAccessor, MatrixAccessor)
    }
    else if(!(pInMat->separated_columns()) && pOutMat->separated_columns())
    {
        CALL_DEEP_COPY_1(MatrixAccessor, SepMatrixAccessor)
    }
    else
    {
        CALL_DEEP_COPY_1(MatrixAccessor, MatrixAccessor)
    }

    return R_NilValue;
}