printsorted() implementation

This commit is contained in:
Jiri Pittner 2025-03-04 16:06:16 +01:00
parent ecdd2fefa5
commit ad1bee99a5
4 changed files with 91 additions and 1 deletions

14
mat.h
View File

@ -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);
}

15
smat.h
View File

@ -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

33
t.cc
View File

@ -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);
}
}

30
vec.h
View File

@ -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