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

bugfixes in smat and related, mainly for complex hermitian interpretation

Browse Source
This commit is contained in:
jiri committed 2026-10-27 20:36:24 +01:00
1 parent 1d4b33b50c
commit 25d0d40bee
6 files changed
+65 -32
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

mat.cc
+6 -6
View File
@@ -1705,8 +1705,8 @@ NRMat< std::complex<double> >::operator+=(const NRSMat< std::complex<double> > &
p += i + 1; p += i + 1;
} }
p = rhs; p++; p = rhs; p++;
for(register int i=1; i<nn; i++){ for(register int i=1; i<nn; i++){ //upper triangle is conjugated (hermitian)
cblas_zaxpy(i, &CONE, p, 1, (*this)[0]+i, nn); for(register int j=0; j<i; ++j) *((*this)[j]+i) += std::conj(p[j]);
p += i+1; p += i+1;
} }
#ifdef CUDALA #ifdef CUDALA
@@ -1752,7 +1752,7 @@ NRMat<T> & NRMat<T>::operator+=(const NRSMat<T> &rhs) {
} }
p = rhs; p++; p = rhs; p++;
for(register int i=1; i<nn; i++) { for(register int i=1; i<nn; i++) {
for(register int j=0; j<i; ++j) *((*this)[j]+i) += p[j]; for(register int j=0; j<i; ++j) *((*this)[j]+i) += LA_traits<T>::conjugate(p[j]); //hermitian for complex
p += i+1; p += i+1;
} }
return *this; return *this;
@@ -1829,8 +1829,8 @@ NRMat<std::complex<double> >::operator-=(const NRSMat<std::complex<double> > &rh
p += i + 1; p += i + 1;
} }
p = rhs; p++; p = rhs; p++;
for(register int i=1; i<nn; i++){ for(register int i=1; i<nn; i++){ //upper triangle is conjugated (hermitian)
cblas_zaxpy(i, &CMONE, p, 1, (*this)[0]+i, nn); for(register int j=0; j<i; ++j) *((*this)[j]+i) -= std::conj(p[j]);
p += i + 1; p += i + 1;
} }
#ifdef CUDALA #ifdef CUDALA
@@ -1872,7 +1872,7 @@ NRMat<T> & NRMat<T>::operator-=(const NRSMat<T> &rhs) {
} }
p = rhs; p++; p = rhs; p++;
for(register int i=1; i<nn; i++){ for(register int i=1; i<nn; i++){
for(register int j=0; j<i; ++j) *((*this)[j]+i) -= p[j]; for(register int j=0; j<i; ++j) *((*this)[j]+i) -= LA_traits<T>::conjugate(p[j]); //hermitian for complex
p += i+1; p += i+1;
} }
return *this; return *this;
mat.h
+5 -3
View File
@@ -675,7 +675,7 @@ NRMat<T>::NRMat(const NRMat &rhs) {
} }
/***************************************************************************//** /***************************************************************************//**
* create matrix from a \c NRSMat object * create matrix from a \c NRSMat object (for complex types hermitian, the upper triangle is conjugated)
* @param[in] rhs \c NRSMat input object to be converted * @param[in] rhs \c NRSMat input object to be converted
* @see count, v, vec.h, NRSMat<T> * @see count, v, vec.h, NRSMat<T>
******************************************************************************/ ******************************************************************************/
@@ -702,13 +702,15 @@ NRMat<T>::NRMat(const NRSMat<T> &rhs) {
#ifdef MATPTR #ifdef MATPTR
for (i=0; i<nn; i++){ for (i=0; i<nn; i++){
for (j=0; j<=i; j++){ for (j=0; j<=i; j++){
v[i][j] = v[j][i] = rhs[k++]; v[i][j] = rhs[k];
v[j][i] = LA_traits<T>::conjugate(rhs[k++]);
} }
} }
#else #else
for (i=0; i<nn; i++){ for (i=0; i<nn; i++){
for (j=0; j<=i; j++){ for (j=0; j<=i; j++){
v[i*(size_t)nn + j] = v[j*(size_t)nn + i] = rhs[k++]; v[i*(size_t)nn + j] = rhs[k];
v[j*(size_t)nn + i] = LA_traits<T>::conjugate(rhs[k++]);
} }
} }
#endif #endif
nonclass.cc
+6 -4
View File
@@ -102,7 +102,8 @@ void lawritemat(FILE *file,const T *a,int r,int c,const char *form0,
ii=j; ii=j;
jj=i; jj=i;
} }
fprintf(file, f, ((std::complex<double>)a[ii*(ii+1)/2+jj]).real(), ((std::complex<double>)a[ii*(ii+1)/2+jj]).imag()); const std::complex<double> x = (std::complex<double>)a[(ii-1)*ii/2+jj-1]; //packed lower triangle, 1-based ii,jj
fprintf(file, f, x.real(), i>=j ? x.imag() : -x.imag()); //hermitian for complex
} else fprintf(file, f, ((std::complex<double>)a[(i-1)*c+j-1]).real(), ((std::complex<double>)a[(i-1)*c+j-1]).imag()); } else fprintf(file, f, ((std::complex<double>)a[(i-1)*c+j-1]).real(), ((std::complex<double>)a[(i-1)*c+j-1]).imag());
if (j < n2) fputc(' ',file); if (j < n2) fputc(' ',file);
} }
@@ -122,7 +123,8 @@ void lawritemat(FILE *file,const T *a,int r,int c,const char *form0,
ii=j; ii=j;
jj=i; jj=i;
} }
fprintf(file, f, ((std::complex<double>)a[ii*(ii+1)/2+jj]).real(), ((std::complex<double>)a[ii*(ii+1)/2+jj]).imag()); const std::complex<double> x = (std::complex<double>)a[(ii-1)*ii/2+jj-1]; //packed lower triangle, 1-based ii,jj
fprintf(file, f, x.real(), i>=j ? x.imag() : -x.imag()); //hermitian for complex
} else fprintf(file,f,((std::complex<double>)a[(i-1)*c+j-1]).real(), ((std::complex<double>)a[(i-1)*c+j-1]).imag()); } else fprintf(file,f,((std::complex<double>)a[(i-1)*c+j-1]).real(), ((std::complex<double>)a[(i-1)*c+j-1]).imag());
putc(j<c?' ':'\n',file); putc(j<c?' ':'\n',file);
} }
@@ -1555,7 +1557,7 @@ const NRSMat<double> realpart<NRSMat< std::complex<double> > >(const NRSMat< std
#ifdef CUDALA #ifdef CUDALA
if(a.location == cpu){ if(a.location == cpu){
#endif #endif
cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a(0,0)), 2, result, 1); cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a[0]), 2, result, 1);
#ifdef CUDALA #ifdef CUDALA
}else{ }else{
laerror("not implemented for cuda yet"); laerror("not implemented for cuda yet");
@@ -1573,7 +1575,7 @@ const NRSMat<double> imagpart<NRSMat< std::complex<double> > >(const NRSMat< std
if(a.location == cpu){ if(a.location == cpu){
#endif #endif
cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a(0,0))+1, 2, result, 1); cblas_dcopy(a.nrows()*(a.ncols()+1)/2, (const double *)(&a[0])+1, 2, result, 1);
#ifdef CUDALA #ifdef CUDALA
}else{ }else{
laerror("not implemented for cuda yet"); laerror("not implemented for cuda yet");
smat.cc
+4 -3
View File
@@ -102,7 +102,8 @@ void NRSMat<T>::get(int fd, bool dim, bool transp, bool orcaformat) {
/***************************************************************************//** /***************************************************************************//**
* constructor symmetrizing given matrix \f$A\f$ of general type <code>T</code> yielding \f$(A+A^\mathrm{T})/2\f$ * constructor symmetrizing given matrix \f$A\f$ of general type <code>T</code> yielding \f$(A+A^\mathrm{T})/2\f$,
* for complex types hermitizing it, i.e. yielding \f$(A+A^\dagger)/2\f$
* @param[in] rhs matrix \f$A\f$ * @param[in] rhs matrix \f$A\f$
******************************************************************************/ ******************************************************************************/
template <typename T> template <typename T>
@@ -122,7 +123,7 @@ NRSMat<T>::NRSMat(const NRMat<T> &rhs) {
int i, j, k(0); int i, j, k(0);
for(i=0; i<nn; i++){ for(i=0; i<nn; i++){
for(j=0; j<=i; j++){ for(j=0; j<=i; j++){
v[k++] = (rhs[i][j] + rhs[j][i])/((T)2); v[k++] = (rhs[i][j] + LA_traits<T>::conjugate(rhs[j][i]))/((T)2);
} }
} }
} }
@@ -759,7 +760,7 @@ void NRSMat<std::complex<double> >::axpy(const std::complex<double> alpha, const
#ifdef CUDALA #ifdef CUDALA
if(location == cpu){ if(location == cpu){
#endif #endif
cblas_zaxpy(nn, &alpha, x.v, 1, v, 1); cblas_zaxpy(NN2, &alpha, x.v, 1, v, 1);
#ifdef CUDALA #ifdef CUDALA
}else{ }else{
const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag()); const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag());
smat.h
+39 -14
View File
@@ -31,6 +31,15 @@ namespace LA {
template<typename T> class NRPerm; template<typename T> class NRPerm;
template<typename T> class NRSMat_from1; template<typename T> class NRSMat_from1;
//return type of the const element access NRSMat<T>::operator(): a reference, except for complex types,
//where the upper triangle has to be conjugated on the fly
template<typename T> struct NRSMat_const_element {typedef const T& type;};
template<typename C> struct NRSMat_const_element<std::complex<C> > {typedef const std::complex<C> type;};
//conjugation of the stored lower triangle element when the upper triangle element is requested
template<typename T> inline const T& NRSMat_conj_upper(const T &x) {return x;}
template<typename C> inline const std::complex<C> NRSMat_conj_upper(const std::complex<C> &x) {return std::conj(x);}
/***************************************************************************//** /***************************************************************************//**
* This class implements a general symmetric or hermitian matrix the elements * This class implements a general symmetric or hermitian matrix the elements
@@ -154,9 +163,11 @@ public:
inline const T& operator[](const size_t ij) const; inline const T& operator[](const size_t ij) const;
inline T& operator[](const size_t ij); inline T& operator[](const size_t ij);
//NOTE: it stores the matrix as symemtric and operator() assumes it is symmetric, does not support complex hermitean as that would require a smart pointer for l-value version //NOTE: complex NRSMat is hermitian, only the lower triangle is stored, the upper one is its complex conjugate
//complex conjugation options are available for BLAS routines to facilitate hermitan matrices //const operator() returns conj(a(j,i)) for i<j (by value for complex types, by reference otherwise)
inline const T& operator()(const int i, const int j) const; //PITFALL: the non-const operator() cannot conjugate (it would need a smart pointer for the l-value), for i<j it returns
//a reference to the stored element a(j,i) WITHOUT conjugation - this also applies to reading from a non-const complex matrix
inline typename NRSMat_const_element<T>::type operator()(const int i, const int j) const;
inline T& operator()(const int i, const int j); inline T& operator()(const int i, const int j);
inline int nrows() const; inline int nrows() const;
@@ -177,7 +188,7 @@ public:
{ {
T s = NRVec<T>(*this).sum(); T s = NRVec<T>(*this).sum();
if(onlytriangle) return s; if(onlytriangle) return s;
s*=2; s += LA_traits<T>::conjugate(s); //upper triangle is conjugated for complex (hermitian)
for(int i=0; i<nn; ++i) s -= (*this)(i,i); for(int i=0; i<nn; ++i) s -= (*this)(i,i);
return s; return s;
} }
@@ -193,7 +204,16 @@ public:
return r; return r;
} }
const NRVec<T> csum() const {return rsum();}; const NRVec<T> csum() const //differs from rsum() by conjugation for complex (hermitian)
{
NRVec<T> r(nn);
for(int j=0; j<nn;++j)
{
r[j]=(T)0;
for(int i=0; i<nn; ++i) r[j] += (*this)(i,j);
}
return r;
}
const T trace() const; const T trace() const;
@@ -306,6 +326,7 @@ inline NRSMat<T>::NRSMat(const T *a, const int n) : nn(n), count(new int(1)) {
location = DEFAULT_LOC; location = DEFAULT_LOC;
if(location == cpu){ if(location == cpu){
#endif #endif
v = new T[NN2];
if(LA_traits<T>::is_plaindata()) memcpy(v, a, NN2*sizeof(T)); if(LA_traits<T>::is_plaindata()) memcpy(v, a, NN2*sizeof(T));
else for( int i=0; i<NN2; i++) v[i] = a[i]; else for( int i=0; i<NN2; i++) v[i] = a[i];
#ifdef CUDALA #ifdef CUDALA
@@ -763,6 +784,9 @@ return (i>j)? (i-2)*(i-1)/2+j-1 : (j-2)*(j-1)/2+i-1;
/***************************************************************************//** /***************************************************************************//**
* determine matrix element of this symmetric matrix of general type <code>T</code> * determine matrix element of this symmetric matrix of general type <code>T</code>
* PITFALL for complex (hermitian) matrices: for \f$i<j\f$ the reference to the stored element \f$a_{ji}\f$
* is returned WITHOUT complex conjugation (an l-value cannot be conjugated); this also
* applies when only reading from a non-const matrix, use the const version or \f$i\geq{}j\f$
* @param[in] i row index running from 0 * @param[in] i row index running from 0
* @param[in] j column index running from 0 * @param[in] j column index running from 0
* @return reference to the corresponding matrix element * @return reference to the corresponding matrix element
@@ -786,14 +810,15 @@ inline T & NRSMat<T>::operator()(const int i, const int j) {
} }
/***************************************************************************//** /***************************************************************************//**
* determine matrix element of this symmetric matrix of general type <code>T</code> * determine matrix element of this symmetric or hermitian matrix of general type <code>T</code>
* @param[in] i row index running from 0 * @param[in] i row index running from 0
* @param[in] j column index running from 0 * @param[in] j column index running from 0
* @return constant reference to the corresponding matrix element * @return constant reference to the corresponding matrix element; for complex types the value,
* complex conjugated for \f$i<j\f$ (hermitian matrix)
* @see count, SMat_index, NRSMat<T>::operator[] * @see count, SMat_index, NRSMat<T>::operator[]
******************************************************************************/ ******************************************************************************/
template <typename T> template <typename T>
inline const T & NRSMat<T>::operator()(const int i, const int j) const { inline typename NRSMat_const_element<T>::type NRSMat<T>::operator()(const int i, const int j) const {
#ifdef DEBUG #ifdef DEBUG
if(i<0||i>=nn||j<0||j>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<" "<<0<<" "<<j<<" "<<nn-1<<std::endl; if(i<0||i>=nn||j<0||j>=nn) std::cout<<"INDEX PROBLEM "<<0<<" "<<i<<" "<<nn-1<<" "<<0<<" "<<j<<" "<<nn-1<<std::endl;
if(i<0) laerror("T & NRSMat<T>::operator()(const int, const int) first index out of range - low"); if(i<0) laerror("T & NRSMat<T>::operator()(const int, const int) first index out of range - low");
@@ -805,7 +830,7 @@ inline const T & NRSMat<T>::operator()(const int i, const int j) const {
#endif #endif
NOT_GPU(*this); NOT_GPU(*this);
return v[SMat_index(i,j)]; return (i>=j) ? v[SMat_index_igej(i,j)] : NRSMat_conj_upper(v[SMat_index_igej(j,i)]);
} }
/***************************************************************************//** /***************************************************************************//**
@@ -1284,7 +1309,7 @@ std::istream& operator>>(std::istream &s, NRSMat<T> &x) {
template <typename T, typename S> template <typename T, typename S>
void NRSMat_convert(NRSMat<T> &a, const NRSMat<S> &b) void NRSMat_convert(NRSMat<T> &a, const NRSMat<S> &b)
{ {
a.resize(b.nrows(),b.ncols()); a.resize(b.nrows());
for(int i=0; i<b.nrows(); ++i) for(int i=0; i<b.nrows(); ++i)
for(int j=0; j<=i; ++j) for(int j=0; j<=i; ++j)
a(i,j) = (T) b(i,j); a(i,j) = (T) b(i,j);
@@ -1307,7 +1332,7 @@ tmp.printsorted(s,direction,absvalue,stable,thr,topmax,&indexmap);
* corresponding macro is defined in vec.h * corresponding macro is defined in vec.h
******************************************************************************/ ******************************************************************************/
NRVECMAT_OPER(SMat,+) NRVECMAT_OPER(SMat,+)
NRVECMAT_OPER(SMat,-) NRVECMAT_OPER_NONCOMMUTATIVE(SMat,-)
NRVECMAT_OPER(SMat,*) NRVECMAT_OPER(SMat,*)
/***************************************************************************//** /***************************************************************************//**
@@ -1332,7 +1357,7 @@ public:
explicit NRSMat_from1(const NRMat<T> &rhs): NRSMat<T>(rhs) {}; explicit NRSMat_from1(const NRMat<T> &rhs): NRSMat<T>(rhs) {};
explicit NRSMat_from1(const NRVec<T> &rhs, const int n): NRSMat<T>(rhs,n) {}; explicit NRSMat_from1(const NRVec<T> &rhs, const int n): NRSMat<T>(rhs,n) {};
inline const T& operator() (const int i, const int j) const { inline typename NRSMat_const_element<T>::type operator() (const int i, const int j) const { //for complex conjugated for i<j, see NRSMat<T>::operator()
#ifdef DEBUG #ifdef DEBUG
if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl; if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl;
if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low"); if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low");
@@ -1340,10 +1365,10 @@ public:
if(j<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - low"); if(j<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - low");
if(j>NRSMat<T>::nn) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - high"); if(j>NRSMat<T>::nn) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) second index out of range - high");
#endif #endif
return NRSMat<T>::v[SMat_index_1(i,j)]; return (i>=j) ? NRSMat<T>::v[SMat_index_1igej(i,j)] : NRSMat_conj_upper(NRSMat<T>::v[SMat_index_1igej(j,i)]);
} }
inline T& operator() (const int i, const int j){ inline T& operator() (const int i, const int j){ //PITFALL: no conjugation for complex i<j, see NRSMat<T>::operator()
#ifdef DEBUG #ifdef DEBUG
if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl; if(i<=0||i>NRSMat<T>::nn||j<=0||j>NRSMat<T>::nn) std::cout<<"INDEX PROBLEM "<<1<<" "<<i<<" "<<NRSMat<T>::nn<<" "<<1<<" "<<j<<" "<<NRSMat<T>::nn<<std::endl;
if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low"); if(i<=0) laerror("index in const T& NRSMat_from1<T>::operator() (const int, const int) first indexout of range - low");
vec.cc
+5 -2
View File
@@ -655,7 +655,7 @@ void NRVec<std::complex<double> >::gemv(const double beta, const NRSMat<double>
* \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f] * \f[\vec{y}\leftarrow \alpha\operatorname{op}(A)\cdot\vec{x}+\beta\vec{y}\f]
* @param[in] beta complex parameter \f$\beta\f$ * @param[in] beta complex parameter \f$\beta\f$
* @param[in] A <b>complex Hermitian</b> matrix \f$A\f$ stored in packed form * @param[in] A <b>complex Hermitian</b> matrix \f$A\f$ stored in packed form
* @param[in] trans not used * @param[in] trans 'n' or 'c' for \f$A\f$, 't' for \f$A^\mathrm{T}=A^*\f$
* @param[in] alpha complex parameter \f$\alpha\f$ * @param[in] alpha complex parameter \f$\alpha\f$
* @param[in] x complex vector \f$\vec{x}\f$ * @param[in] x complex vector \f$\vec{x}\f$
* @see gemm, NRSMat<T> * @see gemm, NRSMat<T>
@@ -674,7 +674,10 @@ void NRVec<std::complex<double> >::gemv(const std::complex<double> beta,
#ifdef CUDALA #ifdef CUDALA
if(location == cpu){ if(location == cpu){
#endif #endif
cblas_zhpmv(CblasRowMajor, CblasLower, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1); if(tolower(trans) == 't') //row-major lower packed H is column-major upper packed H^T
cblas_zhpmv(CblasColMajor, CblasUpper, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1);
else
cblas_zhpmv(CblasRowMajor, CblasLower, A.ncols(), &alpha, A, x.v, 1, &beta, v, 1);
#ifdef CUDALA #ifdef CUDALA
}else{ }else{
const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag()); const cuDoubleComplex _alpha = make_cuDoubleComplex(alpha.real(), alpha.imag());