diff --git a/ChangeLog b/ChangeLog index 6b644d4..e950f76 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,6 @@ +12.11.2019 fourindex_dense with nosymetry implemented +12.11.2019 fourindex.cc created for template specializations +12.11.2019 fourindex optional scaling and terminator 11.11.2019 fourindex for twobody integrals with reduced symmetry 08.11.2019 Conversion constructor for NRMat_from1(NRSMat_from1) 31.03.2019 AUTOCONF files adapted for using MKL as first option diff --git a/Makefile.am b/Makefile.am index e27a4e6..7bebaf3 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,6 +1,6 @@ lib_LTLIBRARIES = libla.la include_HEADERS = fortran.h cuda_la.h auxstorage.h davidson.h laerror.h mat.h qsort.h vec.h bisection.h diis.h la.h noncblas.h smat.h bitvector.h fourindex.h la_traits.h nonclass.h sparsemat.h sparsesmat.h csrmat.h conjgrad.h gmres.h matexp.h permutation.h -libla_la_SOURCES = vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc +libla_la_SOURCES = vec.cc mat.cc smat.cc sparsemat.cc sparsesmat.cc csrmat.cc laerror.cc noncblas.cc bitvector.cc strassen.cc nonclass.cc cuda_la.cc fourindex.cc check_PROGRAMS = t test t_SOURCES = t.cc t2.cc test_SOURCES = test.cc @@ -22,12 +22,14 @@ AM_CXXFLAGS += -g AM_CXXFLAGS += $(OPTIMIZEOPT) $(CUDAOPT) $(FORINTOPT) $(DEBUGOPT) $(MATPTROPT) AM_CXXFLAGS += -DNO_STRASSEN -DFORTRAN_ -AM_CXXFLAGS += $(CBLASOPT) $(CLAPACKOPT) +AM_CXXFLAGS += $(CBLASOPT) $(CLAPACKOPT) +AM_CXXFLAGS += $(MKLOPT) AM_CXXFLAGS += $(TRACEBACKOPT) #AM_LDFLAGS += .libs/libla.a AM_LDFLAGS += $(CBLASLIB) AM_LDFLAGS += $(BLASLIB) AM_LDFLAGS += $(ATLASLIB) +AM_LDFLAGS += $(MKLLIB) AM_LDFLAGS += $(CUDALIBS) AM_LDFLAGS += $(TRACEBACKLIB) diff --git a/fourindex.cc b/fourindex.cc new file mode 100644 index 0000000..04a4c22 --- /dev/null +++ b/fourindex.cc @@ -0,0 +1,55 @@ +/* + LA: linear algebra C++ interface library + Copyright (C) 2008-2019 Jiri Pittner or + + This program is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with this program. If not, see . +*/ + +#include "fourindex.h" +namespace LA { + +template <> +void fourindex::fscanf(FILE *ff) +{ +int i,j,k,l; +double elem; +if(5!= ::fscanf(ff,"%lf %d %d %d %d",&elem,&i,&j,&k,&l)) laerror("read error in fourindex::fscanf"); +while(i!=terminator&&j!=terminator&&k!=terminator&&l!=terminator) + { + add(i,j,k,l,elem); + if(5!= ::fscanf(ff,"%lf %d %d %d %d",&elem,&i,&j,&k,&l)) laerror("read error in fourindex::fscanf"); + } +} + +template <> +void fourindex::fprintf(FILE *f, char *format) const +{ +fourindex::iterator it=this->begin(),end=this->end(); +while(it!=end) + { + ::fprintf(f,format,it->elem,it->index.indiv.i,it->index.indiv.j,it->index.indiv.k,it->index.indiv.l); + ++it; + } +::fprintf(f,format,0.,0,0,0,0); +} + + + +/***************************************************************************//** + * forced instantization in the corresponding object file + ******************************************************************************/ +template class fourindex; +template class fourindex_ext; +template class fourindex_dense; +} diff --git a/fourindex.h b/fourindex.h index d133eee..fa962bc 100644 --- a/fourindex.h +++ b/fourindex.h @@ -1,6 +1,6 @@ /* LA: linear algebra C++ interface library - Copyright (C) 2008 Jiri Pittner or + Copyright (C) 2008-2019 Jiri Pittner or This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -19,6 +19,7 @@ #ifndef _fourindex_included #define _fourindex_included #include +#include #include #include #include @@ -108,7 +109,7 @@ __attribute__((packed)) typedef enum { undefined_symmetry=-1, nosymmetry=0, - twoelectronrealmullikan=1, + twoelectronrealmullikan=1, twoelectronrealmullikanAA=1, twoelectronrealdirac=2, T2ijab_aces=3, trdm2AA=3, @@ -147,6 +148,13 @@ void symmetry_faktor(const fourindexsymtype symmetry,const union packed_index { switch(symmetry) { + case twoelectronrealmullikanreducedsymAA: + if(index.indiv.i==index.indiv.j && index.indiv.k==index.indiv.l) elem*=.5; + if(index.indiv.i==index.indiv.k && index.indiv.j==index.indiv.l) elem*=.5; + break; + case twoelectronrealmullikanreducedsymAB: + if(index.indiv.i==index.indiv.j && index.indiv.k==index.indiv.l) elem*=.5; + break; case antisymtwoelectronrealdirac: case antisymtwoelectronrealdiracAB: laerror("not implemented"); @@ -165,7 +173,7 @@ switch(symmetry) case T2ijab_aces: break; //result will automatically vanish due to generated antisymmetry; i!=a from principle case T2IjAb_aces: break; //no actual symmetry case nosymmetry: break; - default: laerror("illegal symmetry"); + default: laerror("illegal symmetry or symmetry-redundant scaling factor not implemented"); } } @@ -177,6 +185,8 @@ protected: I nn; int *count; matel4 *list; + I terminator; + bool doscaling; private: void deletelist(); void copylist(const matel4 *l); @@ -198,10 +208,14 @@ public: const matel4 * operator->() const {return p;} const matel4 & operator*() const {return *p;} }; + void setterminator(const I terminator0) {terminator=terminator0;} + I getterminator() const {return terminator;} + void setscaling(const bool doscaling0) {doscaling=doscaling0;} iterator begin() const {return list;} iterator end() const {return NULL;} - //permiterator ... iterates also over all permutations, with a possibly scaled matrix element or skips permutations yielding equivalent result + //permiterator ... iterates also over all permutations, with a possibly scaled matrix element for redundant results of the permutations + //it is assumed the original fourindex is the symmetry-reduced petite list, then the results of piterator contributions can be accumulated //has to take into account the symmetry type of the fourindex class piterator { private: @@ -209,6 +223,7 @@ public: matel4 *p; matel4 my; int permindex; + bool doscaling; void setup(void) //make a copy of *p to my with scaled element and anti/permuted indices { if(symmetry==undefined_symmetry) laerror("fourindex symmetry has not been set"); @@ -218,17 +233,17 @@ public: my.elem = p->elem * fourindex_permutations[symmetry][permindex][4]; //now treat the redundancy due to possibly equal indices by a scaling factor //if the processing of individual term becomes very costly, an alternative would be to screen permutations yielding identical result - symmetry_faktor(symmetry, p->index, my.elem); + if(doscaling) symmetry_faktor(symmetry, p->index, my.elem); }; public: piterator() {}; - piterator(matel4 *pp): symmetry(nosymmetry),p(pp),permindex(0){}; + piterator(matel4 *pp): symmetry(nosymmetry),p(pp),permindex(0),doscaling(true){}; ~piterator() {}; - piterator(const fourindex &x): symmetry(x.symmetry),p(x.list),permindex(0) {setup();}; + piterator(const fourindex &x): symmetry(x.symmetry),p(x.list),permindex(0),doscaling(x.doscaling) {setup();}; piterator& operator++() {if(++permindex>=fourindex_permnumbers[symmetry]) {permindex=0; p=p->next;} setup(); return *this;} const matel4 & operator*() const {return my;} const matel4 * operator->() const {return &my;} - piterator operator++(int) {laerror("postincrement not possible on permute-iterator");} + piterator operator++(int) {laerror("postincrement not possible on permute-iterator"); return *this;} bool operator!=(const piterator &rhs) const {return p!=rhs.p;} //should only be used for comparison with pend() bool end(void) {return !p;} bool notend(void) {return p;} @@ -237,8 +252,8 @@ public: piterator pend() const {return piterator(NULL);}//inefficient, use end() or notend() instead //constructors etc. - inline fourindex() :symmetry(undefined_symmetry),nn(0),count(NULL),list(NULL) {}; - inline fourindex(const I n) :symmetry(undefined_symmetry),nn(n),count(new int(1)),list(NULL) {}; + inline fourindex() :symmetry(undefined_symmetry),nn(0),count(NULL),list(NULL),terminator(-1),doscaling(true) {}; + inline fourindex(const I n, const fourindexsymtype symmetry0=undefined_symmetry, const I terminator0= -1, const bool doscaling0=true) :nn(n),count(new int(1)),list(NULL),symmetry(symmetry0),terminator(terminator0),doscaling(doscaling0) {}; fourindex(const fourindex &rhs); //copy constructor inline int getcount() const {return count?*count:0;} fourindex & operator=(const fourindex &rhs); @@ -267,6 +282,8 @@ public: {matel4 *ltmp= new matel4; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &rhs.index, sizeof(union packed_index)); list->elem=rhs.elem;} unsigned long put(int fd,bool withattr=true) const; unsigned long get(int fd,bool withattr=true); + void fscanf(FILE *f); //C-style formatted IO + void fprintf(FILE *f, char *format) const; }; @@ -292,6 +309,8 @@ protected: unsigned int nread; fourindexsymtype symmetry; I nn; + bool doscaling; + I terminator; //methods void tryread() const @@ -315,8 +334,11 @@ protected: public: + void setterminator(const I terminator0) {terminator=terminator0;} + I getterminator() const {return terminator;} + void setscaling(const bool doscaling0) {doscaling=doscaling0;} void resize(I n) {nn=n;} - fourindex_ext(const int file, const fourindexsymtype s=undefined_symmetry, const I n=0, const unsigned int b=1024) :current(NULL),fd(file),nread(0),symmetry(s),nn(n) + fourindex_ext(const int file, const fourindexsymtype s=undefined_symmetry, const I n=0, const unsigned int b=1024, const I terminator0= -1, const bool doscaling0= true) :current(NULL),fd(file),nread(0),symmetry(s),nn(n) { struct statfs sfs; struct stat64 sf; @@ -334,6 +356,8 @@ public: buffer = (matel4stored *) buf; mlock(buf,bufsize); //ignore error when not root, hope we will not be paged out anyway bufsize /= sizeof(matel4stored); + terminator=terminator0; + doscaling=doscaling0; } ~fourindex_ext() {if(buffer0) delete[] buffer0;} void setsymmetry(fourindexsymtype s) {symmetry=s;}; @@ -384,7 +408,7 @@ public: ~iterator() {}; bool operator!=(const iterator &rhs) const {return base!=rhs.base;} //should only be used for comparison with end() iterator &operator++() {if(base) base->next(); if(base->eof()) base=NULL; return *this;} - iterator operator++(int) {laerror("postincrement not possible");} + iterator operator++(int) {laerror("postincrement not possible"); return *this;} const matel4stored * operator->() const {return base->current;} const matel4stored & operator*() const {return *base->current;} bool notNULL() const {return base;} @@ -411,7 +435,7 @@ public: my.elem = it->elem * fourindex_permutations[base->symmetry][permindex][4]; //redundancy due to possibly equal indices //if the processing of individual term becomes very costly, an alternative would be to screen permutations yielding identical result - symmetry_faktor(base->symmetry, it->index, my.elem); + if(base->doscaling) symmetry_faktor(base->symmetry, it->index, my.elem); }; public: piterator() {}; @@ -420,7 +444,7 @@ public: ~piterator() {}; bool operator!=(const piterator &rhs) const {return base!=rhs.base;} //should only be used for comparison with end() piterator &operator++() {if(++permindex>=fourindex_permnumbers[base->symmetry]) {permindex=0; ++it;} if(it.notNULL()) setup(); else base=NULL; return *this;} - piterator operator++(int) {laerror("postincrement not possible");} + piterator operator++(int) {laerror("postincrement not possible"); return *this;} const matel4 * operator->() const {return &my;} const matel4 & operator*() const {return my;} bool end(void) {return !base;} @@ -504,6 +528,8 @@ if(! &rhs) laerror("fourindex copy constructor with NULL argument"); if(rhs.list&&!rhs.count) laerror("some inconsistency in fourindex contructors or assignments"); list=rhs.list; if(list) {count=rhs.count; (*count)++;} else count=new int(1); //make the matrix defined, but empty and not shared + terminator=rhs.terminator; + doscaling=rhs.doscaling; } @@ -518,6 +544,8 @@ fourindex & fourindex::operator=(const fourindex &rhs) if(--(*count) ==0) {deletelist(); delete count;} // old stuff obsolete list=rhs.list; nn=rhs.nn; + terminator=rhs.terminator; + doscaling=rhs.doscaling; if(list) count=rhs.count; else count= new int(0); //make the matrix defined, but empty and not shared, count will be incremented below if(count) (*count)++; } @@ -622,6 +650,7 @@ while(l) return n; } + template std::ostream& operator<<(std::ostream &s, const fourindex_ext &x) { @@ -634,7 +663,7 @@ std::ostream& operator<<(std::ostream &s, const fourindex_ext &x) s << (typename LA_traits_io::IOtype)it->index.indiv.i << ' ' << (typename LA_traits_io::IOtype)it->index.indiv.j<< ' ' <<(typename LA_traits_io::IOtype)it->index.indiv.k << ' ' << (typename LA_traits_io::IOtype)it->index.indiv.l << ' ' << (typename LA_traits_io::IOtype) it->elem << '\n'; ++it; } - s << "-1 -1 -1 -1 0.\n"; + s << (typename LA_traits_io::IOtype) x.getterminator() << ' ' << (typename LA_traits_io::IOtype)x.getterminator() << ' ' <<(typename LA_traits_io::IOtype)x.getterminator() << ' ' << (typename LA_traits_io::IOtype)x.getterminator() << ' ' << (typename LA_traits_io::IOtype) 0 << '\n'; return s; } @@ -652,7 +681,7 @@ std::ostream& operator<<(std::ostream &s, const fourindex &x) s << (typename LA_traits_io::IOtype)it->index.indiv.i << ' ' << (typename LA_traits_io::IOtype)it->index.indiv.j<< ' ' <<(typename LA_traits_io::IOtype)it->index.indiv.k << ' ' << (typename LA_traits_io::IOtype)it->index.indiv.l << ' ' << (typename LA_traits_io::IOtype) it->elem << '\n'; ++it; } - s << "-1 -1 -1 -1 0.\n"; + s << (typename LA_traits_io::IOtype)x.getterminator() << ' ' << (typename LA_traits_io::IOtype)x.getterminator() << ' ' <<(typename LA_traits_io::IOtype)x.getterminator() << ' ' << (typename LA_traits_io::IOtype)x.getterminator() << ' ' << (typename LA_traits_io::IOtype) 0 << '\n'; return s; } @@ -665,7 +694,7 @@ std::istream& operator>>(std::istream &s, fourindex &x) s >> n ; x.resize(n); s >> i >> j >>k >>l; - while(i!= (typename LA_traits_io::IOtype)-1 && j!= (typename LA_traits_io::IOtype)-1 && k != (typename LA_traits_io::IOtype)-1 && l!= (typename LA_traits_io::IOtype)-1) + while(i!= (typename LA_traits_io::IOtype)x.getterminator() && j!= (typename LA_traits_io::IOtype)x.getterminator() && k != (typename LA_traits_io::IOtype)x.getterminator() && l!= (typename LA_traits_io::IOtype)x.getterminator()) { s>>elem; x.add((I)i,(I)j,(I)k,(I)l,(T)elem); @@ -684,7 +713,7 @@ std::istream& operator>>(std::istream &s, fourindex_ext &x) typename LA_traits_io::IOtype elem; s >> i >> j >>k >>l; - while(i!= (typename LA_traits_io::IOtype)-1 && j!= (typename LA_traits_io::IOtype)-1 && k != (typename LA_traits_io::IOtype)-1 && l!= (typename LA_traits_io::IOtype)-1) + while(i!= (typename LA_traits_io::IOtype)x.getterminator() && j!= (typename LA_traits_io::IOtype)x.getterminator() && k != (typename LA_traits_io::IOtype)x.getterminator() && l!= (typename LA_traits_io::IOtype)x.getterminator()) { s>>elem; x.put((I)i,(I)j,(I)k,(I)l,(T)elem); @@ -919,6 +948,40 @@ int J = SMat_index_1(k,l); return NRSMat::v[SMat_index(I,J)]; } +template +class fourindex_dense : public NRMat { +protected: + unsigned int nn; +friend class explicit_t2; +public: + fourindex_dense(): NRMat() {nn=0;}; + void resize(const int nnn) {nn=nnn; (*this).NRMat::resize(nn*nn,nn*nn);}; + explicit fourindex_dense(const int nnn): NRMat(nnn*nnn,nnn*nnn) {nn=nnn;}; + + inline T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) +{ +#ifdef DEBUG +if(i<1||i>nn ||j<1||j>nn|| a<1||a>nn||b<1||b>nn) laerror("nosymmetry fourindex out of range"); +if (!NRMat::v) laerror("access to unallocated fourindex_dense"); +#endif +return (*this).NRMat::operator() ((j-1)*nn+i-1,(b-1)*nn+a-1); +} + inline const T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) const +{ +#ifdef DEBUG +if(i<1||i>nn ||j<1||j>nn|| a<1||a>nn||b<1||b>nn) laerror("nosymmetry fourindex out of range"); +if (!NRMat::v) laerror("access to unallocated fourindex_dense"); +#endif +return (*this).NRMat::operator() ((j-1)*nn+i-1,(b-1)*nn+a-1); +} + +void print(std::ostream &out) const + { + unsigned int i,j,a,b; + for(i=1; i<=nn; ++i) for(j=1; j<=nn; ++j) for(a=1; a<=nn; ++a) for(b=1; b<=nn; ++b) out << i<<" "<