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

*** empty log message ***

Browse Source
This commit is contained in:
jiri 2008-06-17 14:37:27 +00:00
parent 57e9b3a599
commit f5796068ec
1 changed files with 29 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
fourindex.h
View File

@ -859,6 +859,7 @@ public:
void set(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem); void set(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
void add(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem); void add(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
void add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const; const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);}; void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);};
void print(ostream &out) const void print(ostream &out) const
@ -922,18 +923,31 @@ if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
NRSMat<T>::v[SMat_index(I,J)] += elem; NRSMat<T>::v[SMat_index(I,J)] += elem;
} }
template<class T, class DUMMY>
void fourindex_dense<antisymtwoelectronrealdirac,T,DUMMY>::add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem)
{
if(i<=j || k<=l) return;
int I = ASMat_index_1(i,j);
int J = ASMat_index_1(k,l);
if (I<0 || J<0 || I<J) return;
NRSMat<T>::v[SMat_index(I,J)] += elem;
}
template<class T, class I> template<class T, class I>
fourindex_dense<antisymtwoelectronrealdirac,T,I>::fourindex_dense(const fourindex<I,T> &rhs) : nbas(rhs.size()), NRSMat<T>((T)0,rhs.size()*(rhs.size()-1)/2) fourindex_dense<antisymtwoelectronrealdirac,T,I>::fourindex_dense(const fourindex<I,T> &rhs) : nbas(rhs.size()), NRSMat<T>((T)0,rhs.size()*(rhs.size()-1)/2)
{ {
if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch"); if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
typename fourindex<I,T>::iterator p; typename fourindex_ext<I,T>::piterator p; //we have to run over equivalents in non-canonical order to build the antisymmetrization properly; it could be done less elegantly but more efficiently moving the if's to outer parts of the piterator loop, if needed
for(p=rhs.begin(); p!= rhs.end(); ++p) for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
if(p->index.indiv.i!=p->index.indiv.k && p->index.indiv.j!=p->index.indiv.l) {
{ I i=p->index.indiv.i;
add(p->index.indiv.i,p->index.indiv.k,p->index.indiv.j,p->index.indiv.l,p->elem); I j=p->index.indiv.j;
add(p->index.indiv.i,p->index.indiv.k,p->index.indiv.l,p->index.indiv.j, -p->elem); I k=p->index.indiv.k;
} I l=p->index.indiv.l;
add_unique(i,k,j,l,p->elem);
add_unique(i,k,l,j,-p->elem);
}
} }
@ -941,12 +955,15 @@ template<class T, class I>
fourindex_dense<antisymtwoelectronrealdirac,T,I>::fourindex_dense(const fourindex_ext<I,T> &rhs) : nbas(rhs.size()), NRSMat<T>((T)0,rhs.size()*(rhs.size()-1)/2) fourindex_dense<antisymtwoelectronrealdirac,T,I>::fourindex_dense(const fourindex_ext<I,T> &rhs) : nbas(rhs.size()), NRSMat<T>((T)0,rhs.size()*(rhs.size()-1)/2)
{ {
if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch"); if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
typename fourindex_ext<I,T>::iterator p; typename fourindex_ext<I,T>::piterator p; //we have to run over equivalents in non-canonical order to build the antisymmetrization properly; it could be done less elegantly but more efficiently moving the if's to outer parts of the piterator loop, if needed
for(p=rhs.begin(); p!= rhs.end(); ++p) for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
if(p->index.indiv.i!=p->index.indiv.k && p->index.indiv.j!=p->index.indiv.l)
{ {
add(p->index.indiv.i,p->index.indiv.k,p->index.indiv.j,p->index.indiv.l,p->elem); I i=p->index.indiv.i;
add(p->index.indiv.i,p->index.indiv.k,p->index.indiv.l,p->index.indiv.j, -p->elem); I j=p->index.indiv.j;
I k=p->index.indiv.k;
I l=p->index.indiv.l;
add_unique(i,k,j,l,p->elem);
add_unique(i,k,l,j,-p->elem);
} }
} }