jiri/LA_library
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

bugfix in permutation of smatrix

Browse Source
This commit is contained in:
Jiri Pittner
2023-07-28 15:38:24 +02:00
parent 9d51ef6fbd
commit ec42999812
2 changed files with 26 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

28
smat.cc
View File

@@ -311,13 +311,35 @@ const NRSMat<T> NRSMat<T>::permuted(const NRPerm<int> &p, const bool inverse) co
if(!p.is_valid()) laerror("invalid permutation of smatrix");
#endif
int n=p.size();
if(n!=(*this).size()) laerror("incompatible permutation and smatrix");
if(n!=(*this).nrows()) laerror("incompatible permutation and smatrix");
#ifdef CUDALA
if(this->getlocation() != cpu || p.getlocation() != cpu ) laerror("permutations can be done only in CPU memory");
#endif
NRSMat<T> r(n);
if(inverse) for(int i=1; i<=n; ++i) {int pi = p[i]-1; r(i-1,i-1) = (*this)(pi,pi);}
else for(int i=1; i<=n; ++i) {int pi = p[i]-1; r(pi,pi) = (*this)(i-1,i-1);}
if(inverse)
{
for(int i=1; i<=n; ++i)
{
int pi = p[i]-1;
for(int j=1; j<=i; ++j)
{
int pj = p[j] - 1;
r(i-1,j-1) = (*this)(pi,pj);
}
}
}
else
{
for(int i=1; i<=n; ++i)
{
int pi = p[i]-1;
for(int j=1; j<=i; ++j)
{
int pj = p[j] - 1;
r(pi,pj) = (*this)(i-1,j-1);
}
}
}
return r;
}