DCSBMcoCluster<-function(adjMat,kR,kC,rowNorm,gamma,nrep){
  
  adj<-1e-15
  
  pckg = try(require(flexclust))
  if(!pckg) {
    cat("Installing 'flexclust' from CRAN\n")
    install.packages("flexclust", repos = "http://cran.r-project.org")
    require(flexclust)
  }
  pckg = try(require(rARPACK))
  if(!pckg) {
    cat("Installing 'rARPACK' from CRAN\n")
    install.packages("rARPACK", repos = "http://cran.r-project.org")
    require(rARPACK)
  }
  
  degDistR<-rowSums(adjMat);
  degDistC<-colSums(adjMat);
  m2<-sum(adjMat);
  modMat<-adjMat-degDistR%*%t(degDistC)/m2;
  
  NGcoModularity<-function(tempClustersR,tempClustersC){
    localCoModMat<-sapply(sort(unique(tempClustersC)),function(tempC){
      return(sapply(sort(unique(tempClustersR)),function(tempR){
        return(sum(modMat[tempClustersR==tempR,tempClustersC==tempC]))
      }))
    })	
    globalCoMod<-sum(abs(localCoModMat))/(2*m2);	
    localSizeMat<-sapply(sort(unique(tempClustersC)),function(tempC){
      return(sapply(sort(unique(tempClustersR)),function(tempR){
        return(sum(tempClustersR==tempR)*sum(tempClustersC==tempC))
      }))
    })		
    return(list(globalCoMod=globalCoMod,localCoModMat=localCoModMat,localSizeMat))		
  }
  
  nR<-nrow(adjMat);
  nC<-ncol(adjMat);
  DR_2<-((degDistR+max(1,median(degDistR)))^-(1/2))%*%t(rep(1,nC));
  DC_2<-rep(1,nR)%*%t((degDistC+max(1,median(degDistC)))^-(1/2));
  L<-DR_2*adjMat*DC_2;
  tempSVD<-svds(L,k=min(max(c(kR,kC)),ncol(L),nrow(L)),nu=kR,nv=kC)
  if(kR>1){
    U<-tempSVD[["u"]][,2:kR]
  }else{
    U<-tempSVD[["u"]][,1]
  }	
  if(kC>1){
    V<-tempSVD[["v"]][,2:kC]
  }else{
    V<-tempSVD[["v"]][,1]
  }
  if(is.vector(U)){
    U<-as.matrix(U)
  }
  if(is.vector(V)){
    V<-as.matrix(V)
  }
  leveragesR<-sqrt(rowSums(U^2));
  leveragesC<-sqrt(rowSums(V^2));
  toKMclusterR<-which(leveragesR>=gamma/sqrt(nR));
  toKMclusterC<-which(leveragesC>=gamma/sqrt(nC));
  notKMclusterR<-setdiff(1:nR,toKMclusterR);
  notKMclusterC<-setdiff(1:nC,toKMclusterC);
  if(rowNorm){
    U<-U/sqrt(rowSums(U^2)+adj);
    V<-V/sqrt(rowSums(V^2)+adj);
  }
  modBest<-(-Inf);

  for(irep in 1:nrep){
    
    tempKR<-max(min(kR,length(toKMclusterR)-1),1)
    if(tempKR>1){    
      DCSBMclustersTempR<-rep(0,nR);
      tempClustR<-tryCatch(cclust(U[toKMclusterR,],tempKR,control=list(initcent="kmeanspp",iter.max=1000)),error=function(e){return(NA)})
      if(mode(tempClustR)=="S4"){
        DCSBMclustersTempR[toKMclusterR]<-tempClustR@cluster
        CR<-tempClustR@centers
      }else{
        DCSBMclustersTempR[toKMclusterR]<-rep(1,length(toKMclusterR))
        CR<-matrix(1,nrow=tempKR,ncol=ncol(U))
      }
      if(length(notKMclusterR)>0){
        for(i in 1:length(notKMclusterR)){
          DCSBMclustersTempR[notKMclusterR[i]]<-which.min(sqrt(rowSums((CR-rep(1,nrow(CR))%*%t(U[notKMclusterR[i],]))^2)));
        }
      }
    }else{
      DCSBMclustersTempR<-rep(1,nR);   			
    }
    
    tempKC<-max(min(kC,length(toKMclusterC)-1))
    if(tempKC>1){ 		
      DCSBMclustersTempC<-rep(0,nC);        
      tempClustC<-tryCatch(cclust(V[toKMclusterC,],tempKC,control=list(initcent="kmeanspp",iter.max=1000)),error=function(e){return(NA)})
      if(mode(tempClustC)=="S4"){
        DCSBMclustersTempC[toKMclusterC]<-tempClustC@cluster
        CC<-tempClustC@centers
      }else{
        DCSBMclustersTempC[toKMclusterC]<-rep(1,length(toKMclusterC))
        CC<-matrix(1,nrow=tempKC,ncol=ncol(V))		    		
      }
      if(length(notKMclusterC)>0){
        for(i in 1:length(notKMclusterC)){
          DCSBMclustersTempC[notKMclusterC[i]]<-which.min(sqrt(rowSums((CC-rep(1,nrow(CC))%*%t(V[notKMclusterC[i],]))^2)));
        }
      }
    }else{
      DCSBMclustersTempC<-rep(1,nC);   			
    }
    
    modTemp<-NGcoModularity(DCSBMclustersTempR,DCSBMclustersTempC)[["globalCoMod"]];
    if(modTemp>modBest){
      modBest<-modTemp;
      DCSBMclustersR<-DCSBMclustersTempR;
      DCSBMclustersC<-DCSBMclustersTempC;
    }
    
  }
  
  coModMat<-NGcoModularity(DCSBMclustersR,DCSBMclustersC)[["localCoModMat"]];
  if(length(unique(DCSBMclustersC))==1){
    coModMat<-as.matrix(coModMat)
  }
  if(length(unique(DCSBMclustersR))==1&length(unique(DCSBMclustersC))!=1){
    coModMat<-t(coModMat)
  }
  rowCoMod<-rowSums(abs(coModMat));
  colCoMod<-colSums(abs(coModMat));
  rowOrd<-order(rowCoMod,decreasing=T)
  colOrd<-order(colCoMod,decreasing=T)
  for(i in 1:kR){
    DCSBMclustersR[DCSBMclustersR==rowOrd[i]]<-(-i);
  }
  for(j in 1:kC){
    DCSBMclustersC[DCSBMclustersC==colOrd[j]]<-(-j);
  }
  DCSBMclustersR<-(-DCSBMclustersR);
  DCSBMclustersC<-(-DCSBMclustersC);
  
  return(list(clustersR=DCSBMclustersR,clustersC=DCSBMclustersC))
  
}