*** empty log message ***

la_traits.h

@@ -31,6 +31,22 @@ template<typename C> class NRMat;
 template<typename C> class NRSMat;
 template<typename C> class SparseMat;
 
+//for general sortable classes
+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);};
+        };
+
+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);};
+        };
+
+
 //we will need to treat char and unsigned char as numbers in << and >> I/O operators
 template<typename C>
 struct LA_traits_io

qsort.h

@@ -43,18 +43,19 @@ return parity;
 }
 
 
-template<typename INDEX>
-int genqsort2(INDEX l, INDEX r, bool (*bigger)(const INDEX, const INDEX), void (*swap)(const INDEX,const INDEX)) 
+//for SORTABLE classes which provide LA_sort_traits<SORTABLE,INDEX,type>::compare and swap member functions
+template<int type, typename SORTABLE, typename INDEX>
+int memqsort(SORTABLE &object, INDEX l, INDEX r)
 {
 INDEX i,j,piv;
 int parity=0;
 
 if(r<=l) return parity; //1 element
-if(bigger(l,r)) {parity^=1; swap(l,r);}
+if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,r)) {parity^=1; object.swap(l,r);}
 if(r-l==1) return parity; //2 elements and preparation for median
 piv= (l+r)/2; //pivoting by median of 3 - safer 
-if(bigger(l,piv)) {parity^=1; swap(l,piv);} //and change the pivot element implicitly
-if(bigger(piv,r)) {parity^=1; swap(r,piv);} //and change the pivot element implicitly
+if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,l,piv)) {parity^=1; object.swap(l,piv);} //and change the pivot element implicitly
+if(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,piv,r)) {parity^=1; object.swap(r,piv);} //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
@@ -62,23 +63,23 @@ 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(bigger(piv,i++));
+  while(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,piv,i++));
   i--;
-  while(bigger(j--,piv));
+  while(LA_sort_traits<SORTABLE,INDEX,type>::compare(object,j--,piv));
   j++;
   if(i<j)
 	{
 	// swap and keep track of position of pivoting element 
-	parity^=1; swap(i,j); 
+	parity^=1; object.swap(i,j); 
 	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 ^=genqsort2(l,j,bigger,swap); if(i<r) parity ^=genqsort2(i,r,bigger,swap);}
+	{if(l<j) parity ^=memqsort<type,SORTABLE,INDEX>(object,l,j); if(i<r) parity ^=memqsort<type,SORTABLE,INDEX>(object,i,r);}
 else
-	{if(i<r) parity ^=genqsort2(i,r,bigger,swap); if(l<j) parity ^=genqsort2(l,j,bigger,swap);}
+	{if(i<r) parity ^=memqsort<type,SORTABLE,INDEX>(object,i,r); if(l<j) parity ^=memqsort<type,SORTABLE,INDEX>(object,l,j);}
 return parity;
 }
 

sparsemat.h

@@ -3,14 +3,6 @@
 
 #include "la_traits.h"
 
-template<typename T>
-inline const T MAX(const T &a, const T &b)
-        {return b > a ? (b) : (a);}
-
-template<typename T>
-inline void SWAP(T &a, T &b)
-        {T dum=a; a=b; b=dum;}
-
 
 //threshold for neglecting elements, if not defined, no tests are done except exact zero test in simplify - might be even faster
 //seems to perform better with a threshold, in spite of abs() tests

vec.cc

@@ -493,36 +493,11 @@ NRVec< complex<double> >::operator|(const NRVec< complex<double> > &b) const
 }
 
 
-//sorting of elements in the vector
 template<typename T>
-static inline void SWAP(T &a, T &b)
-        {T dum=a; a=b; b=dum;}
-
-static void *sortbase; //global for sort !!! is not thread-safe
-
-template<typename T>
-static void myswap(int i, int j)
-{
-SWAP(((T *)sortbase)[i],((T *)sortbase)[j]);
-}
-
-template<typename T>
-static bool mybigger(int i, int j)
-{
-return LA_traits<T>::bigger(((T *)sortbase)[i],((T *)sortbase)[j]);
-}
-
-template<typename T>
-static bool mysmaller(int i, int j)
-{
-return LA_traits<T>::smaller(((T *)sortbase)[i],((T *)sortbase)[j]);
-}
-
-template<typename T>
-int NRVec<T>::sort(int direction)
+int NRVec<T>::sort(int direction, int from, int to)
 {
 copyonwrite();
-sortbase = v;
-if(direction) return genqsort2(0,nn-1,mysmaller<T>,myswap<T>);
-else return genqsort2(0,nn-1,mybigger<T>,myswap<T>);
+if(to == -1) to=nn-1;
+if(direction) return memqsort<1,NRVec<T>,int>(*this,from,to);
+else return memqsort<0,NRVec<T>,int>(*this,from,to);
 }

vec.h

@@ -103,9 +103,13 @@ public:
 //sparse matrix concerning members
 	explicit NRVec(const SparseMat<T> &rhs);                // dense from sparse matrix with one of dimensions =1
 	inline void simplify() {}; //just for compatibility with sparse ones
-	int sort(int direction=0); //sort, ascending by default, returns parity of permutation
+	bool bigger(int i, int j) const {return LA_traits<T>::bigger(v[i],v[j]);};
+	bool smaller(int i, int j) const {return LA_traits<T>::smaller(v[i],v[j]);};
+	void swap(int i, int j) {T tmp; tmp=v[i]; v[i]=v[j]; v[j]=tmp;};
+	int sort(int direction=0, int from=0, int to= -1); //sort, ascending by default, returns parity of permutation
 };
 
+
 //due to mutual includes this has to be after full class declaration
 #include "mat.h"
 #include "smat.h"