DCSBMcluster<-function(adjMat,k,rowNorm,gamma,nrep){

	adj<-1e-15

	pckg = try(require(rARPACK))
	if(!pckg) {
		cat("Installing 'rARPACK' from CRAN\n")
		install.packages("rARPACK", repos = "http://cran.r-project.org")
		require(rARPACK)
	}
	pckg = try(require(flexclust))	
	if(!pckg) {
	cat("Installing 'flexclust' from CRAN\n")
		install.packages("flexclust", repos = "http://cran.r-project.org")
		require(flexclust)
	}
		
	degDist<-rowSums(adjMat);
	D2<-1/sqrt(degDist+max(1,median(degDist)));
	L<-t(adjMat*D2)*D2;
	m2<-sum(degDist);
	modMat<-adjMat-degDist%*%t(degDist)/m2;		

	eigDecomp<-tryCatch(eigs(L,k,which="LM"),error=function(e){return(NA)})
	if(all(is.na(eigDecomp))){
		svdDecomp<-T
	}else if(any(is.complex(eigDecomp[["vectors"]]))){
		svdDecomp<-T	
	}else{
		svdDecomp<-F
	}
	if(svdDecomp){
		tempDecomp<-svd(L)
		eigDecomp<-list(vectors=tempDecomp[["u"]][,1:k],values=tempDecomp[["d"]][1:k])
		rm(list="tempDecomp")
	}
	if(k>1){
		Xadj<-eigDecomp[["vectors"]][,2:k]
	}else{
		Xadj<-eigDecomp[["vectors"]][,1]
	}
	if(is.vector(Xadj)){
		Xadj<-as.matrix(Xadj)
	}
	rm(list=c("eigDecomp","L"))
	
	leverages<-sqrt(rowSums(Xadj^2));
	toKMcluster<-which(leverages>=gamma/sqrt(nrow(Xadj)));
	notKMcluster<-setdiff(1:nrow(Xadj),toKMcluster);
	if(rowNorm){
		Xadj<-Xadj/sqrt(rowSums(Xadj^2)+adj);
	}

	tempK<-max(min(k,length(toKMcluster)-1),1) 

	if(tempK>1){	        	
        		        		   
		tempClustering<-lapply(1:nrep,function(irep){
        	    
			DCSBMclustersTemp<-rep(0,nrow(Xadj));
			tempClust<-tryCatch(cclust(Xadj[toKMcluster,],tempK,control=list(initcent="kmeanspp",iter.max=100)),error=function(e){return(NA)})
			if(mode(tempClust)=="S4"){
				DCSBMclustersTemp[toKMcluster]<-tempClust@cluster
				C<-tempClust@centers
			}else{
				DCSBMclustersTemp[toKMcluster]<-rep(1,length(toKMcluster))
				C<-matrix(1,nrow=tempK,ncol=ncol(Xadj))		    		
			}
			if(length(notKMcluster)>0){
				for(i in 1:length(notKMcluster)){
					DCSBMclustersTemp[notKMcluster[i]]<-which.min(sqrt(rowSums((C-rep(1,nrow(C))%*%t(Xadj[notKMcluster[i],]))^2)));
				}
			}
		 		
			modTemp<-sum(sapply(unique(DCSBMclustersTemp),function(tempCom){
				tempLoc<-DCSBMclustersTemp==tempCom;
				return(sum(modMat[tempLoc,tempLoc]))
			}))/m2
		        
			return(list(mod=modTemp,clusters=DCSBMclustersTemp))
	
		})
		
		tempLoc<-which.max(sapply(tempClustering,function(temp){return(temp[["mod"]])}))
		DCSBMclusters<-tempClustering[[tempLoc]][["clusters"]]
			
	}else{
		modBest<-0;
		DCSBMclusters<-rep(1,nrow(Xadj));  			
	}

	localMod<-sapply(sort(unique(DCSBMclusters)),function(tempNum){
		tempLoc<-DCSBMclusters==tempNum
		return(mean(modMat[tempLoc,tempLoc]))
	})
	ord<-order(localMod,decreasing=T)
	for(i in 1:k){
		DCSBMclusters[DCSBMclusters==ord[i]]<-(-i);
	}
	DCSBMclusters<-(-DCSBMclusters);
    	
	return(DCSBMclusters)
	
}