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

some more fourindex-tensor conversions

Browse Source
This commit is contained in:
Jiri Pittner 2025-10-24 15:21:51 +02:00
parent f1060245a0
commit f2d392c2a7
3 changed files with 154 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
fourindex.h
View File

@ -1221,9 +1221,11 @@ protected:
unsigned int nnbas;
friend class explicit_t2;
public:
unsigned int nbas() const {return nnbas;};
fourindex_dense(): NRSMat<T>() {nnbas=0;};
void resize(const int n) {nnbas=n; (*this).NRSMat<T>::resize(nnbas*nnbas);};
explicit fourindex_dense(const int n): NRSMat<T>(n*n) {nnbas=n;};
explicit fourindex_dense(const int n, const NRSMat<T>&data): NRSMat<T>(data) {nnbas=n; if(data.nrows()!=n*n) laerror("data size mismatch in fourindex_dense antisymtwoelectronrealdiracAB constructor");};
//here i,a are alpha j,b beta
inline T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b)
@ -1256,11 +1258,13 @@ return (*this).NRSMat<T>::operator() ((j-1)*nnbas+i-1,(b-1)*nnbas+a-1);
template<class T, class I>
class fourindex_dense<T2ijab_aces,T,I> : public NRMat<T> {
protected:
unsigned int nocc,nvrt,ntri;
unsigned int ntri;
friend class explicit_t2;
public:
unsigned int nocc,nvrt;
fourindex_dense(): NRMat<T>() {nocc=nvrt=ntri=0;};
explicit fourindex_dense(const int noc, const int nvr): NRMat<T>(noc*(noc-1)/2,nvr*(nvr-1)/2) {nocc=noc; nvrt=nvr; ntri=nvr*(nvr-1)/2;};
explicit fourindex_dense(const int noc, const int nvr, const NRMat<T> &data): NRMat<T>(data) {nocc=noc; nvrt=nvr; ntri=nvr*(nvr-1)/2; if(data.nrows()!=noc*(noc-1)/2 || data.ncols()!=nvr*(nvr-1)/2) laerror("data size mismatch in T2ijab_aces constructor"); };
void resize(const int noc, const int nvr) {(*this).NRMat<T>::resize(noc*(noc-1)/2,nvr*(nvr-1)/2); nocc=noc; nvrt=nvr; ntri=nvr*(nvr-1)/2;};
//we cannot return reference due to the possible sign change
@ -1322,8 +1326,10 @@ class fourindex_dense<antisymtwoelectronrealdirac,T,I> : public NRSMat<T> {
private:
int nnbas;
public:
int nbas() const {return nnbas;};
fourindex_dense(): NRSMat<T>() {};
explicit fourindex_dense(const int n): nnbas(n), NRSMat<T>(n*(n-1)/2) {};
explicit fourindex_dense(const int n, const NRSMat<T> &data): nnbas(n), NRSMat<T>(data) {if(data.nrows()!=n*(n-1)/2) laerror("data size mismatch in fourindex_dense antisymtwoelectronrealdirac constructor");};
fourindex_dense(const T &a, const int n): nnbas(n), NRSMat<T>(a,n*(n-1)/2) {};
fourindex_dense(const T *a, const int n): nnbas(n), NRSMat<T>(a,n*(n-1)/2) {};
//and also construct it from sparse and externally stored fourindex classes

45
t.cc
View File

@ -3743,9 +3743,8 @@ tensor2fourindex(t,ff);
cout <<"Error = "<<(f-ff).norm()<<endl;
}
if(1)
if(0)
{
//check also index order
int n;
cin>>n;
fourindex_dense<T2IjAb_aces,double,int> f(n+1,n,n+3,n+2);
@ -3757,7 +3756,49 @@ tensor2fourindex(t,ff);
cout <<"Error = "<<(f-ff).norm()<<endl;
}
if(0)
{
int n;
cin>>n;
fourindex_dense<T2ijab_aces,double,int> f(n,n+2);
f.randomize(1.);
Tensor<double> t=fourindex2tensor(f);
cout <<t;
fourindex_dense<T2ijab_aces,double,int> ff;
tensor2fourindex(t,ff);
cout <<"Error = "<<(f-ff).norm()<<endl;
}
if(0)
{
int n;
cin>>n;
fourindex_dense<antisymtwoelectronrealdiracAB,double,int> f(n);
f.randomize(1.);
Tensor<double> t=fourindex2tensor(f);
cout <<t;
fourindex_dense<antisymtwoelectronrealdiracAB,double,int> ff;
tensor2fourindex(t,ff);
cout <<"Error = "<<(f-ff).norm()<<endl;
}
if(1)
{
int n;
cin>>n;
fourindex_dense<antisymtwoelectronrealdirac,double,int> f(n);
f.randomize(1.);
Tensor<double> t=fourindex2tensor(f);
cout <<t;
fourindex_dense<antisymtwoelectronrealdirac,double,int> ff;
tensor2fourindex(t,ff);
cout <<"Error = "<<(f-ff).norm()<<endl;
}
//@@@check operator() for all fourindex-tensor conversions
}//main

105
tensor.h
View File

@ -47,7 +47,7 @@
//@@@ will not be particularly efficient
//
//@@@conversions to/from fourindex, optional negative range for beta spin handling in some cases
//@@@use the fact that fourindex_dense is inherited from Mat/SMat and construct tensor from the (unsymmetrized) NRMat sharing data, just rewrite then the shape
//@@@check operator() for all fourindex-tensor conversions
//
//@@@?general permutation of individual indices - check the indices in sym groups remain adjacent, calculate result's shape, loopover the result and permute using unwind_callback
//
@ -414,6 +414,109 @@ f=fourindex_dense<T2IjAb_aces,T,I>(noca,nocb,nvra,nvrb,mat);
template<typename T, typename I>
Tensor<T> fourindex2tensor(const fourindex_dense<T2ijab_aces,T,I> &f)
{
NRVec<indexgroup> shape(2);
for(int i=0; i<2; ++i)
{
shape[i].number=2;
shape[i].symmetry= -1;
shape[i].offset=1;
}
shape[1].range=f.nocc;
shape[0].range=f.nvrt;
NRVec<T> data(f);
return Tensor<T>(shape,data);
}
template<typename T, typename I>
void tensor2fourindex(const Tensor<T> &t, fourindex_dense<T2ijab_aces,T,I> &f)
{
if(t.rank()!=4) laerror("wrong rank in tensor2fourindex");
if(t.shape.size()!=2) laerror("wrong symmetry groups in tensor2fourindex");
int nvrt=t.shape[0].range;
int nocc=t.shape[1].range;
int offset=t.shape[0].offset;
for(int i=0; i<t.shape.size(); ++i)
{
if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
if(-1!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
}
NRMat<T> mat(t.data,nocc*(nocc-1)/2,nvrt*(nvrt-1)/2);
f=fourindex_dense<T2ijab_aces,T,I>(nocc,nvrt,mat);
}
template<typename T, typename I>
Tensor<T> fourindex2tensor(const fourindex_dense<antisymtwoelectronrealdiracAB,T,I> &f)
{
NRVec<indexgroup> shape(4);
int n=f.nbas();
for(int i=0; i<4; ++i)
{
shape[i].number=1;
shape[i].symmetry=0;
shape[i].offset=1;
shape[i].range=n;
}
NRMat<T> data(f); //expand to full matrix
return Tensor<T>(shape,NRVec<T>(data));
}
template<typename T, typename I>
void tensor2fourindex(const Tensor<T> &t, fourindex_dense<antisymtwoelectronrealdiracAB,T,I> &f)
{
if(t.rank()!=4) laerror("wrong rank in tensor2fourindex");
if(t.shape.size()!=4) laerror("wrong symmetry groups in tensor2fourindex");
int range=t.shape[0].range;
int offset=t.shape[0].offset;
for(int i=0; i<t.shape.size(); ++i)
{
if(range!=t.shape[i].range) laerror("range mismatch in tensor2fourindex");
if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
if(0!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
}
NRMat<T> mat(t.data,range*range,range*range);
f=fourindex_dense<antisymtwoelectronrealdiracAB,T,I>(range,NRSMat<T>(mat));
}
template<typename T, typename I>
Tensor<T> fourindex2tensor(const fourindex_dense<antisymtwoelectronrealdirac,T,I> &f)
{
NRVec<indexgroup> shape(2);
int n=f.nbas();
for(int i=0; i<2; ++i)
{
shape[i].number=2;
shape[i].symmetry= -1;
shape[i].offset=1;
shape[i].range=n;
}
NRMat<T> mat(f);
NRVec<T> data(mat); //expand symmetric to full matrix for index group symmetry
return Tensor<T>(shape,data);
}
template<typename T, typename I>
void tensor2fourindex(const Tensor<T> &t, fourindex_dense<antisymtwoelectronrealdirac,T,I> &f)
{
if(t.rank()!=4) laerror("wrong rank in tensor2fourindex");
if(t.shape.size()!=2) laerror("wrong symmetry groups in tensor2fourindex");
int range=t.shape[0].range;
int offset=t.shape[0].offset;
for(int i=0; i<t.shape.size(); ++i)
{
if(range!=t.shape[i].range) laerror("range mismatch in tensor2fourindex");
if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
if(-1!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
}
NRMat<T> mat(t.data,range*(range-1)/2,range*(range-1)/2);
f=fourindex_dense<antisymtwoelectronrealdirac,T,I>(range,NRSMat<T>(mat)); //symmetrize mat
}