diff --git a/la_traits.h b/la_traits.h
index fe240fc..81e4b9c 100644
--- a/la_traits.h
+++ b/la_traits.h
@@ -145,16 +145,30 @@ template<typename T, typename I, int type> struct LA_sort_traits;
 template<typename T, typename I>
 struct LA_sort_traits<T,I,0>
         {
-        static inline bool compare(T object, I i, I j) {return object.bigger(i,j);};
+        static inline bool compare(const T &object, I i, I j) {return object.bigger(i,j);};
+	static inline bool comparestable(const T &object, I i, I j, I *auxkey)
+                {
+                bool r= object.bigger(i,j);
+                if(!r && object[i]==object[j]) r= auxkey[i]>auxkey[j];
+                return r;
+                };
+
         };
 
 template<typename T, typename I>
 struct LA_sort_traits<T,I,1>
         {
-        static inline bool compare(T object, I i, I j) {return object.smaller(i,j);};
+        static inline bool compare(const T &object, I i, I j) {return object.smaller(i,j);};
+	static inline bool comparestable(const T &object, I i, I j, I *auxkey)
+                {
+                bool r= object.smaller(i,j);
+                if(!r && object[i]==object[j]) r= auxkey[i]<auxkey[j];
+                return r;
+                };
         };
 
 
+
 //we will need to treat char and unsigned char as numbers in << and >> I/O operators
 template<typename C>
 struct LA_traits_io
diff --git a/qsort.h b/qsort.h
index c7f4cc1..b42f914 100644
--- a/qsort.h
+++ b/qsort.h
@@ -104,6 +104,66 @@ return parity;
 }
 
 
+//stabilized version of quicksort employing auxiliary key of the original position
+
+template<int type, typename SORTABLE, typename INDEX, typename PERMINDEX>
+int memqsortaux(SORTABLE &object, PERMINDEX *perm, INDEX l, INDEX r, INDEX *aux)
+{
+INDEX i,j,piv;
+int parity=0;
+
+if(r<=l) return parity; //1 element
+if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,l,r,aux)) {parity^=1; object.swap(l,r); {INDEX tmp=aux[l]; aux[l]=aux[r]; aux[r]=tmp;} if(perm) {PERMINDEX tmp=perm[l]; perm[l]=perm[r]; perm[r]=tmp;}}
+if(r-l==1) return parity; //2 elements and preparation for median
+piv= l+(r-l)/2; //pivoting by median of 3 - safer 
+if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,l,piv,aux)) {parity^=1; object.swap(l,piv); {INDEX tmp=aux[l]; aux[l]=aux[piv]; aux[piv]=tmp;} if(perm) {PERMINDEX tmp=perm[l]; perm[l]=perm[piv]; perm[piv]=tmp;}} //and change the pivot element implicitly
+if(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,piv,r,aux)) {parity^=1; object.swap(r,piv); {INDEX tmp=aux[r]; aux[r]=aux[piv]; aux[piv]=tmp;} if(perm) {PERMINDEX tmp=perm[r]; perm[r]=perm[piv]; perm[piv]=tmp;}} //and change the pivot element implicitly
+if(r-l==2) return parity; //in the case of 3 elements we are finished too
+
+//general case , l-th r-th already processed
+i=l+1; j=r-1;
+do{
+  //important sharp inequality - stops at sentinel element for efficiency
+  // this is inefficient if all keys are equal - unnecessary n log n swaps are done, but we assume that it is atypical input
+  while(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,piv,i++,aux));
+  i--;
+  while(LA_sort_traits<SORTABLE,INDEX,type>::comparestable(object,j--,piv,aux));
+  j++;
+  if(i<j)
+	{
+	// swap and keep track of position of pivoting element 
+	parity^=1; 
+	object.swap(i,j);  
+	{INDEX tmp =aux[i]; aux[i]=aux[j]; aux[j]=tmp;}
+	if(perm) {PERMINDEX tmp=perm[i]; perm[i]=perm[j]; perm[j]=tmp;}
+	if(i==piv) piv=j; else if(j==piv) piv=i;
+	}
+  if(i<=j) {i++; j--;}
+  }while(i<=j);
+
+if(j-l < r-i)   //because of the stack in bad case process first the shorter subarray
+	{if(l<j) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j,aux); if(i<r) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r,aux);}
+else
+	{if(i<r) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,i,r,aux); if(l<j) parity ^=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,j,aux);}
+return parity;
+}
+
+
+template<int type, typename SORTABLE, typename INDEX, typename PERMINDEX>
+int memqsortstable(SORTABLE &object, PERMINDEX *perm, INDEX l, INDEX r)
+{
+if(r<=l) return 0; //1 element
+INDEX *aux= new INDEX[r-l+1];
+if(aux==NULL) laerror("allocation failed in memqsortstable");
+INDEX *auxshifted = aux-l;
+for(INDEX i=l; i<=r; ++i) auxshifted[i]=i;
+int parity=memqsortaux<type,SORTABLE,INDEX,PERMINDEX>(object,perm,l,r,auxshifted);
+delete[] aux;
+return parity;
+}
+
+
+
 template<typename S, typename PERMINDEX>
 int ptrqsortup(S *l, S *r, PERMINDEX *perm=NULL)
 {
diff --git a/t.cc b/t.cc
index 7e7a640..9fd614b 100644
--- a/t.cc
+++ b/t.cc
@@ -2758,7 +2758,7 @@ cout <<"Check adjacency of the clique: "<<adj.submatrix(clique)<<endl;
 NRVec<int> cover = cliquecover(adj);
 cout <<"Clique cover is "<<cover<<endl;
 NRPerm<int> p(cover.size());
-cover.sort(0,p);
+cover.sort(0,p,true);
 cout<<"permutation to disentabgle the cliques = "<<p<<endl;
 NRSMat<char> adjperm = adj.permuted(p);
 cout <<"resorted graph = "<<adjperm<<endl;
diff --git a/vec.h b/vec.h
index ae3af15..489d29c 100644
--- a/vec.h
+++ b/vec.h
@@ -451,8 +451,8 @@ public:
 	};
 
 	//! sort by default in ascending order and return the parity of corresponding permutation resulting to this order
-	int sort(int direction = 0, int from = 0, int to = -1, int *perm = NULL);
-	int sort(int direction, NRPerm<int> &perm);
+	int sort(int direction = 0, int from = 0, int to = -1, int *perm = NULL, bool stable=false);
+	int sort(int direction, NRPerm<int> &perm, bool stable=false);
 
 	//! apply given function to each element
 	NRVec& call_on_me(T (*_F)(const T &) ){
@@ -518,21 +518,29 @@ namespace LA {
 
 
 template<typename T>
-int NRVec<T>::sort(int direction, int from, int to, int *perm) {
+int NRVec<T>::sort(int direction, int from, int to, int *perm, bool stable) {
 	NOT_GPU(*this);
 
 	copyonwrite();
 	if(to == -1) to = nn - 1;
-	if(direction) return memqsort<1, NRVec<T>, int, int>(*this, perm, from, to);
-	else return memqsort<0, NRVec<T>, int, int>(*this, perm, from, to);
+	if(stable)
+		{
+		if(direction) return memqsortstable<1, NRVec<T>, int, int>(*this, perm, from, to);
+                else return memqsortstable<0, NRVec<T>, int, int>(*this, perm, from, to);
+		}
+	else
+		{	
+		if(direction) return memqsort<1, NRVec<T>, int, int>(*this, perm, from, to);
+		else return memqsort<0, NRVec<T>, int, int>(*this, perm, from, to);
+		}
 }
 
 template<typename T>
-int NRVec<T>::sort(int direction, NRPerm<int> &perm)
+int NRVec<T>::sort(int direction, NRPerm<int> &perm, bool stable)
 {
 if(nn!=perm.size()) laerror("incompatible vector and permutation");
 perm.identity();
-int r=sort(direction,0,nn-1,&perm[1]);
+int r=sort(direction,0,nn-1,&perm[1],stable);
 return r;
 }