diff --git a/mat.h b/mat.h
index 78ead3f..d50c384 100644
--- a/mat.h
+++ b/mat.h
@@ -409,6 +409,10 @@ public:
 	const typename  LA_traits<T>::normtype  nonsymmetry() const;
 	const typename  LA_traits<T>::normtype  nonhermiticity() const;
 
+	//!print sorted list of matrix elements
+        void printsorted(std::ostream &s, int direction=1, bool absvalue=false, bool stable=false, typename LA_traits<T>::normtype thr=0, int topmax=0) const;
+
+
 #ifndef NO_STRASSEN
 	//! Strassen's multiplication (better than \f$\mathacal{O}(n^3)\f$, analogous syntax to \see NRMat<T>::gemm() )
 	void strassen(const T beta, const NRMat &a, const char transa, const NRMat &b, const char transb, const T alpha);
@@ -1582,6 +1586,16 @@ for(int i=0; i<b.nrows(); ++i)
 }
 
 
+template<typename T>
+void NRMat<T>::printsorted(std::ostream &s, int direction, bool absvalue, bool stable, typename LA_traits<T>::normtype thr, int topmax) const                                                                                                                                {
+NRMat<int> indexmap(nn*mm,2);
+NRVec tmp(*this);
+int ii=0;
+for(int i=0; i<nn; ++i) for(int j=0; j<mm; ++j) {indexmap[ii][0]=i; indexmap[ii][1]=j; ii++;}
+tmp.printsorted(s,direction,absvalue,stable,thr,topmax,&indexmap);
+}       
+
+
 
 
 
diff --git a/smat.h b/smat.h
index a5d5825..6f563c7 100644
--- a/smat.h
+++ b/smat.h
@@ -164,6 +164,7 @@ public:
 	inline bool transp(const int i, const int j) const {return i>j;} //this can be used for compact storage of matrices, which are actually not symmetric, but one triangle of them is redundant
 	const typename LA_traits<T>::normtype norm(const T scalar = (T)0) const;
 	void axpy(const T alpha, const NRSMat &x); // this+= a*x
+	
 
 	inline const T amax() const;
 	inline const T amin() const;
@@ -186,6 +187,9 @@ public:
         void storesubmatrix(const int from, const NRSMat &rhs);
         void storesubmatrix(const NRVec<int> &selection, const NRSMat &rhs);
 
+	//!print sorted list of matrix elements
+	void printsorted(std::ostream &s, int direction=1, bool absvalue=false, bool stable=false, typename LA_traits<T>::normtype thr=0, int topmax=0) const;
+
 
 	inline operator T*();
 	inline operator const T*() const;
@@ -1245,6 +1249,17 @@ for(int i=0; i<b.nrows(); ++i)
                 a(i,j) = (T) b(i,j);
 }
 
+template<typename T>
+void NRSMat<T>::printsorted(std::ostream &s, int direction, bool absvalue, bool stable, typename LA_traits<T>::normtype thr, int topmax) const
+{
+NRMat<int> indexmap(NN2,2);
+NRVec tmp(*this);
+int ii=0;
+for(int i=0; i<nn; ++i) for(int j=0; j<=i; ++j) {indexmap[ii][1]=i; indexmap[ii][0]=j; ii++;}
+tmp.printsorted(s,direction,absvalue,stable,thr,topmax,&indexmap);
+}
+
+
 
 /***************************************************************************//**
  * generate operators relating NRSMat<T> objects and scalars
diff --git a/t.cc b/t.cc
index 1a69781..5d21575 100644
--- a/t.cc
+++ b/t.cc
@@ -3440,7 +3440,7 @@ cout <<b;
 cout <<c;
 }
 
-if(1)
+if(0)
 {
 NRMat<double> m;
 cin >>m;
@@ -3451,6 +3451,37 @@ cout <<mm;
 cout <<m.getcount()<<" "<<v.getcount()<<" "<<mm.getcount()<<endl;
 }
 
+if(1)
+{
+NRVec<double> v(12);
+v.randomize(1.);
+cout <<v;
+cout<<"sorted\n";
+v.printsorted(cout,1,false);
+cout<<"abssorted\n";
+v.printsorted(cout,1,true);
+}
+
+if(1)
+{
+NRSMat<double> v(4);
+v.randomize(1.);
+cout <<v;
+cout<<"smat sorted\n";
+v.printsorted(cout,1,false);
+}
+
+if(1)
+{
+NRMat<double> v(4,5);
+v.randomize(1.);
+cout <<v;
+cout<<"mat sorted\n";
+v.printsorted(cout,1,false);
+}
+
+
+
 
 
 }
diff --git a/vec.h b/vec.h
index 0081cf5..4f248aa 100644
--- a/vec.h
+++ b/vec.h
@@ -491,6 +491,7 @@ 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, bool stable=false);
 	int sort(int direction, NRPerm<int> &perm, bool stable=false);
+	void printsorted(std::ostream &s, int direction=1, bool absvalue=false, bool stable=false, typename LA_traits<T>::normtype thr=0, int topmax=0, const NRMat<int> *indexmap=NULL) const;
 
 	//! apply given function to each element
 	NRVec& call_on_me(T (*_F)(const T &) ){
@@ -583,6 +584,35 @@ return r;
 }
 
 
+template<typename T>
+void NRVec<T>::printsorted(std::ostream &s, int direction, bool absvalue, bool stable, typename LA_traits<T>::normtype thr, int topmax,  const NRMat<int> *indexmap) const
+{
+NRPerm<int> perm(nn); perm.identity();
+NRVec_from1<T> tmp;
+NRVec_from1<typename LA_traits<T>::normtype> tmp2;
+if(absvalue) 
+	{
+	tmp2.resize(nn);
+	for(int i=1; i<=nn; ++i) tmp2[i] = MYABS((*this)[i-1]);
+	tmp2.sort(direction,perm,stable);
+	}
+else
+	{
+	tmp= *this;
+	tmp.sort(direction,perm,stable);
+	}
+int ii=1;
+for(int i=1; i<=nn; ++i)
+	{
+	if(topmax && ii>topmax) break;
+	if(thr!=0 && (absvalue?tmp2[i]:MYABS(tmp[i]))<thr) continue;
+	if(indexmap) s<<(*indexmap)[perm[i]-1][0]<<" "<<(*indexmap)[perm[i]-1][1]; else s<<perm[i]-1;
+	s<<' ' << (*this)[perm[i]-1]<<std::endl;
+	++ii;
+	}
+}
+
+
 /***************************************************************************//**
  * indexing operator giving the element at given position with range checking in
  * the DEBUG mode