00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #ifndef _fourindex_included
00020 #define _fourindex_included
00021 #include <iostream>
00022 #include <string.h>
00023 #include <sys/types.h>
00024 #include <sys/vfs.h>
00025 #include <sys/mman.h>
00026 #include <sys/stat.h>
00027 #include <unistd.h>
00028 #include <sys/stat.h>
00029 #include "laerror.h"
00030 #include "vec.h"
00031 #include "smat.h"
00032 #include "mat.h"
00033 #include "nonclass.h"
00034
00035 namespace LA {
00036
00037 static unsigned int hcd0(unsigned int big,unsigned int small)
00038 {
00039 register unsigned int help;
00040
00041
00042 if(big==0)
00043 {
00044 if(small==0) laerror("bad arguments in hcd");
00045 return small;
00046 }
00047 if(small==0) return big;
00048 if(small==1||big==1) return 1;
00049
00050 if(small>big) {help=big; big=small; small=help;}
00051 do {
00052 help=small;
00053 small= big%small;
00054 big=help;
00055 }
00056 while(small != 0);
00057 return big;
00058 }
00059
00060
00061 static inline unsigned int lcm0(unsigned int i,unsigned int j)
00062 {
00063 return (i/hcd0(i,j)*j);
00064 }
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 template<class I>
00075 union packed_index {
00076 I packed[4];
00077 struct {
00078 I i;
00079 I j;
00080 I k;
00081 I l;
00082 } indiv;
00083 };
00084
00085 template<class I, class T>
00086 struct matel4
00087 {
00088 T elem;
00089 matel4 *next;
00090 union packed_index<I> index;
00091 };
00092
00093
00094 template<class I, class T>
00095 struct matel4stored
00096 {
00097 T elem;
00098 union packed_index<I> index;
00099 }
00100 #ifdef __GNUC__
00101 __attribute__((packed))
00102 #endif
00103 ;
00104
00105
00106
00107 typedef enum {undefined_symmetry=-1,nosymmetry=0, twoelectronrealmullikan=1, twoelectronrealdirac=2, T2ijab_aces=3, antisymtwoelectronrealdirac=4, T2IjAb_aces=5} fourindexsymtype;
00108
00109 static const int fourindex_n_symmetrytypes=6;
00110 static const int fourindex_permnumbers[fourindex_n_symmetrytypes]={1,8,8,4,8,8};
00111 static const int fourindex_permutations[fourindex_n_symmetrytypes][8][5]=
00112 {
00113 {{0,1,2,3,1}},
00114 {{0,1,2,3,1}, {1,0,2,3,1}, {0,1,3,2,1}, {1,0,3,2,1}, {2,3,0,1,1}, {3,2,0,1,1}, {2,3,1,0,1}, {3,2,1,0,1}},
00115 {{0,1,2,3,1},{2,1,0,3,1},{0,3,2,1,1},{2,3,0,1,1},{1,0,3,2,1},{3,0,1,2,1},{1,2,3,0,1},{3,2,1,0,1}},
00116 {{0,1,2,3,1},{1,0,2,3,-1},{0,1,3,2,-1},{1,0,3,2,1}},
00117 {{0,1,2,3,1},{1,0,2,3,-1},{0,1,3,2,-1},{1,0,3,2,1},{2,3,0,1,1},{3,2,0,1,-1},{2,3,1,0,-1},{3,2,1,0,1}},
00118 {{0,1,2,3,1}},
00119 };
00120
00121
00122 template <class I, class T>
00123 void symmetry_faktor(const fourindexsymtype symmetry,const union packed_index<I> &index, T &elem)
00124 {
00125 switch(symmetry)
00126 {
00127 case antisymtwoelectronrealdirac:
00128 laerror("not implemented");
00129 case twoelectronrealmullikan:
00130 if(index.indiv.i==index.indiv.j) elem*=.5;
00131 if(index.indiv.k==index.indiv.l) elem*=.5;
00132 if(index.indiv.i==index.indiv.k && index.indiv.j==index.indiv.l
00133 || index.indiv.i==index.indiv.l && index.indiv.j==index.indiv.k ) elem*=.5;
00134 break;
00135 case twoelectronrealdirac:
00136 if(index.indiv.i==index.indiv.k) elem*=.5;
00137 if(index.indiv.j==index.indiv.l) elem*=.5;
00138 if(index.indiv.i==index.indiv.j && index.indiv.k==index.indiv.l
00139 || index.indiv.k==index.indiv.j && index.indiv.i==index.indiv.l) elem*=.5;
00140 break;
00141 case T2ijab_aces: break;
00142 case T2IjAb_aces: break;
00143 case nosymmetry: break;
00144 default: laerror("illegal symmetry");
00145 }
00146 }
00147
00148
00149 template <class I, class T>
00150 class fourindex {
00151 protected:
00152 fourindexsymtype symmetry;
00153 I nn;
00154 int *count;
00155 matel4<I,T> *list;
00156 private:
00157 void deletelist();
00158 void copylist(const matel4<I,T> *l);
00159 public:
00160
00161 class iterator {
00162 private:
00163 matel4<I,T> *p;
00164 public:
00165 iterator() {};
00166 ~iterator() {};
00167 iterator(matel4<I,T> *list): p(list) {};
00168 bool operator==(const iterator &rhs) const {return p==rhs.p;}
00169 bool operator!=(const iterator &rhs) const {return p!=rhs.p;}
00170 iterator &operator++() {p=p->next; return *this;}
00171 iterator operator++(int) {iterator q(p); p=p->next; return q;}
00172 matel4<I,T> & operator*() {return *p;}
00173 matel4<I,T> * operator->() {return p;}
00174 const matel4<I,T> * operator->() const {return p;}
00175 const matel4<I,T> & operator*() const {return *p;}
00176 };
00177 iterator begin() const {return list;}
00178 iterator end() const {return NULL;}
00179
00180
00181
00182 class piterator {
00183 private:
00184 fourindexsymtype symmetry;
00185 matel4<I,T> *p;
00186 matel4<I,T> my;
00187 int permindex;
00188 void setup(void)
00189 {
00190 if(symmetry==undefined_symmetry) laerror("fourindex symmetry has not been set");
00191 if(!p) {permindex=0; memset(&my,0,sizeof(my)); return;}
00192 for(int i=0; i<4; ++i)
00193 my.index.packed[i] = p->index.packed[fourindex_permutations[symmetry][permindex][i]];
00194 my.elem = p->elem * fourindex_permutations[symmetry][permindex][4];
00195
00196
00197 symmetry_faktor(symmetry, p->index, my.elem);
00198 };
00199 public:
00200 piterator() {};
00201 piterator(matel4<I,T> *pp): symmetry(nosymmetry),p(pp),permindex(0){};
00202 ~piterator() {};
00203 piterator(const fourindex &x): symmetry(x.symmetry),p(x.list),permindex(0) {setup();};
00204 piterator& operator++() {if(++permindex>=fourindex_permnumbers[symmetry]) {permindex=0; p=p->next;} setup(); return *this;}
00205 const matel4<I,T> & operator*() const {return my;}
00206 const matel4<I,T> * operator->() const {return &my;}
00207 piterator operator++(int) {laerror("postincrement not possible on permute-iterator");}
00208 bool operator!=(const piterator &rhs) const {return p!=rhs.p;}
00209 bool end(void) {return !p;}
00210 bool notend(void) {return p;}
00211 };
00212 piterator pbegin() const {return piterator(*this);}
00213 piterator pend() const {return piterator(NULL);}
00214
00215
00216 inline fourindex() :symmetry(undefined_symmetry),nn(0),count(NULL),list(NULL) {};
00217 inline fourindex(const I n) :symmetry(undefined_symmetry),nn(n),count(new int(1)),list(NULL) {};
00218 fourindex(const fourindex &rhs);
00219 inline int getcount() const {return count?*count:0;}
00220 fourindex & operator=(const fourindex &rhs);
00221 fourindex & operator+=(const fourindex &rhs);
00222 void setsymmetry(fourindexsymtype s) {symmetry=s;}
00223 fourindexsymtype getsymmetry() const {return symmetry;}
00224 fourindex & join(fourindex &rhs);
00225 inline ~fourindex();
00226 inline matel4<I,T> *getlist() const {return list;}
00227 inline I size() const {return nn;}
00228 void resize(const I n);
00229 void copyonwrite();
00230 unsigned long length() const;
00231 inline void add(const I i, const I j, const I k, const I l, const T elem)
00232 {matel4<I,T> *ltmp= new matel4<I,T>; ltmp->next=list; list=ltmp; list->index.indiv.i=i;list->index.indiv.j=j;list->index.indiv.k=k;list->index.indiv.l=l; list->elem=elem;}
00233
00234 inline void add(const union packed_index<I> &index , const T elem)
00235 {matel4<I,T> *ltmp= new matel4<I,T>; ltmp->next=list; list=ltmp; list->index=index; list->elem=elem;}
00236
00237 inline void add(const I (&index)[4], const T elem)
00238 {matel4<I,T> *ltmp= new matel4<I,T>; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &index, sizeof(union packed_index<I>)); list->elem=elem;}
00239 inline void add(const matel4<I,T> &rhs)
00240 {matel4<I,T> *ltmp= new matel4<I,T>; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &rhs.index, sizeof(union packed_index<I>)); list->elem=rhs.elem;}
00241 inline void add(const matel4stored<I,T> &rhs)
00242 {matel4<I,T> *ltmp= new matel4<I,T>; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &rhs.index, sizeof(union packed_index<I>)); list->elem=rhs.elem;}
00243 unsigned long put(int fd,bool withattr=true) const;
00244 unsigned long get(int fd,bool withattr=true);
00245 };
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 template <class I, class T>
00256 class fourindex_ext {
00257 private:
00258 fourindex_ext();
00259 fourindex_ext(const fourindex_ext &rhs);
00260 fourindex_ext & operator=(const fourindex_ext &rhs);
00261 protected:
00262 matel4stored<I,T> *buffer0;
00263 matel4stored<I,T> *buffer;
00264 matel4stored<I,T> *current;
00265 int fd;
00266 unsigned int bufsize;
00267 unsigned int nread;
00268 fourindexsymtype symmetry;
00269 I nn;
00270
00271
00272 void tryread() const
00273 {
00274 const_cast<fourindex_ext<I,T> *>(this)->current=NULL;
00275 ssize_t r=read(fd,buffer,bufsize*sizeof(matel4stored<I,T>));
00276 if(r<0) {perror("read error"); laerror("read error in fourindex_ext");}
00277 if(r%sizeof(matel4stored<I,T>)) laerror("read inconsistency in fourindex_ext");
00278 const_cast<fourindex_ext<I,T> *>(this)->nread = r/sizeof(matel4stored<I,T>);
00279 if(nread) const_cast<fourindex_ext<I,T> *>(this)->current=buffer;
00280 }
00281 void next() const {
00282 if(current)
00283 {
00284 if( (unsigned int) (++ const_cast<fourindex_ext<I,T> *>(this)->current - buffer) >=nread) tryread();
00285 }
00286 }
00287 bool eof() const {return !current;};
00288
00289
00290 public:
00291 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)
00292 {
00293 struct statfs sfs;
00294 struct stat64 sf;
00295 if(fstat64(fd,&sf)) {perror("cannot fstat");laerror("I/O error");}
00296 if(fstatfs(fd,&sfs)) {perror("cannot fstatfs");laerror("I/O error");}
00297 const unsigned int pagesize=getpagesize();
00298
00299 bufsize=b*sizeof(matel4stored<I,T>);
00300 bufsize=lcm0(bufsize,pagesize);
00301 bufsize=lcm0(bufsize,sfs.f_bsize);
00302 bufsize=lcm0(bufsize,sf.st_blksize);
00303 buffer0 = new matel4stored<I,T>[(bufsize+pagesize)/sizeof(matel4stored<I,T>)+1];
00304 unsigned char *buf= (unsigned char *) buffer0;
00305 buf= buf + pagesize - ((unsigned long)buf % pagesize);
00306 buffer = (matel4stored<I,T> *) buf;
00307 mlock(buf,bufsize);
00308 bufsize /= sizeof(matel4stored<I,T>);
00309 }
00310 ~fourindex_ext() {if(buffer0) delete[] buffer0;}
00311 void setsymmetry(fourindexsymtype s) {symmetry=s;};
00312 fourindexsymtype getsymmetry() const {return symmetry;}
00313 void rewind() const {if(0!=lseek64(fd,0L,SEEK_SET)) {perror("seek error"); laerror("cannot seek in fourindex_ext");} };
00314
00315
00316 void put(const matel4stored<I,T> x)
00317 {
00318 if(!current) current=buffer;
00319 *current++ = x;
00320 if(current-buffer >= bufsize ) flush();
00321 }
00322 void put(I i, I j, I k, I l, const T &elem)
00323 {
00324 if(!current) current=buffer;
00325 current->index.indiv.i=i;
00326 current->index.indiv.j=j;
00327 current->index.indiv.k=k;
00328 current->index.indiv.l=l;
00329 current->elem = elem;
00330 ++current;
00331 if(current-buffer >= bufsize ) flush();
00332 }
00333 void flush()
00334 {
00335 if(current)
00336 {
00337 ssize_t r=write(fd,buffer,(current-buffer)*sizeof(matel4stored<I,T>));
00338 if(r!=(current-buffer)*sizeof(matel4stored<I,T>)) laerror("write error in fourindex_ext");
00339 }
00340 current=NULL;
00341 }
00342
00343 inline I size() const {return nn;}
00344
00345
00346
00347
00348
00349
00350 class iterator {
00351 private:
00352 const fourindex_ext *base;
00353 public:
00354 iterator() {};
00355 iterator(const fourindex_ext *p): base(p) {};
00356 ~iterator() {};
00357 bool operator!=(const iterator &rhs) const {return base!=rhs.base;}
00358 iterator &operator++() {if(base) base->next(); if(base->eof()) base=NULL; return *this;}
00359 iterator operator++(int) {laerror("postincrement not possible");}
00360 const matel4stored<I,T> * operator->() const {return base->current;}
00361 const matel4stored<I,T> & operator*() const {return *base->current;}
00362 bool notNULL() const {return base;}
00363 };
00364 iterator begin() const {rewind(); tryread(); if(!eof()) return this; else return NULL;}
00365 iterator end() const {return iterator(NULL);}
00366
00367
00368
00369 class piterator {
00370 private:
00371 fourindex_ext *base;
00372 matel4<I,T> my;
00373 int permindex;
00374 typename fourindex_ext::iterator it;
00375
00376
00377 void setup(void)
00378 {
00379 if(base->symmetry==undefined_symmetry) laerror("fourindex symmetry has not been set");
00380 if(!it.notNULL()) {permindex=0; memset(&my,0,sizeof(my)); return;}
00381 for(int i=0; i<4; ++i)
00382 my.index.packed[i] = it->index.packed[fourindex_permutations[base->symmetry][permindex][i]];
00383 my.elem = it->elem * fourindex_permutations[base->symmetry][permindex][4];
00384
00385
00386 symmetry_faktor(base->symmetry, it->index, my.elem);
00387 };
00388 public:
00389 piterator() {};
00390 piterator(fourindex_ext *p): base(p),permindex(0) {if(p) {it=p->begin(); setup();}};
00391 piterator(fourindex_ext &x): base(&x),permindex(0) {it= x.begin(); setup();};
00392 ~piterator() {};
00393 bool operator!=(const piterator &rhs) const {return base!=rhs.base;}
00394 piterator &operator++() {if(++permindex>=fourindex_permnumbers[base->symmetry]) {permindex=0; ++it;} if(it.notNULL()) setup(); else base=NULL; return *this;}
00395 piterator operator++(int) {laerror("postincrement not possible");}
00396 const matel4<I,T> * operator->() const {return &my;}
00397 const matel4<I,T> & operator*() const {return my;}
00398 bool end(void) {return !base;}
00399 bool notend(void) {return base;}
00400 };
00401 piterator pbegin() {return piterator(*this);}
00402 piterator pend() const {return piterator(NULL);}
00403
00404
00405 };
00406
00407
00408
00410
00411 template <class I,class T>
00412 unsigned long fourindex<I,T>::put(int fd, bool withattr) const
00413 {
00414 unsigned long n=0;
00415 matel4<I,T> *l=list;
00416 matel4stored<I,T> buf;
00417 if(withattr)
00418 {
00419 union {fourindexsymtype sym; I n; T padding;} u;
00420 u.sym=symmetry;
00421 if(sizeof(u)!=write(fd,&u,sizeof(u))) laerror("write error in fourindex::put");
00422 u.n=nn;
00423 if(sizeof(u)!=write(fd,&u,sizeof(u))) laerror("write error in fourindex::put");
00424 }
00425 while(l)
00426 {
00427 ++n;
00428 buf.elem= l->elem;
00429 buf.index= l->index;
00430 if(sizeof(buf)!=write(fd,&buf,sizeof(buf))) laerror("write error in fourindex::put");
00431 l=l->next;
00432 }
00433 return n;
00434 }
00435
00436
00437 template <class I,class T>
00438 unsigned long fourindex<I,T>::get(int fd,bool withattr)
00439 {
00440 unsigned long n=0;
00441 matel4stored<I,T> buf;
00442 if(withattr)
00443 {
00444 union {fourindexsymtype sym; I n; T padding;} u;
00445 if(sizeof(u)!=read(fd,&u,sizeof(u))) laerror("read inconsistency in fourindex::put");
00446 symmetry=u.sym;
00447 if(sizeof(u)!=read(fd,&u,sizeof(u))) laerror("read inconsistency in fourindex::put");
00448 nn=u.n;
00449
00450 }
00451 while(sizeof(buf)==read(fd,&buf,sizeof(buf))) {++n; add(buf.index,buf.elem);}
00452 return n;
00453 }
00454
00455
00456
00457 template <class I,class T>
00458 fourindex<I,T>::~fourindex()
00459 {
00460 if(!count) return;
00461 if(--(*count)<=0)
00462 {
00463 deletelist();
00464 delete count;
00465 }
00466 }
00467
00468
00469 template <class I, class T>
00470 fourindex<I,T>::fourindex(const fourindex<I,T> &rhs)
00471 {
00472 #ifdef debug
00473 if(! &rhs) laerror("fourindex copy constructor with NULL argument");
00474 #endif
00475 nn=rhs.nn;
00476 if(rhs.list&&!rhs.count) laerror("some inconsistency in fourindex contructors or assignments");
00477 list=rhs.list;
00478 if(list) {count=rhs.count; (*count)++;} else count=new int(1);
00479 }
00480
00481
00482
00483
00484 template <class I, class T>
00485 fourindex<I,T> & fourindex<I,T>::operator=(const fourindex<I,T> &rhs)
00486 {
00487 if (this != &rhs)
00488 {
00489 if(count)
00490 if(--(*count) ==0) {deletelist(); delete count;}
00491 list=rhs.list;
00492 nn=rhs.nn;
00493 if(list) count=rhs.count; else count= new int(0);
00494 if(count) (*count)++;
00495 }
00496 return *this;
00497 }
00498
00499
00500 template <class I, class T>
00501 fourindex<I,T> & fourindex<I,T>::operator+=(const fourindex<I,T> &rhs)
00502 {
00503 if(nn!=rhs.nn) laerror("incompatible dimensions for +=");
00504 if(!count) {count=new int; *count=1; list=NULL;}
00505 else copyonwrite();
00506 register matel4<I,T> *l=rhs.list;
00507 while(l)
00508 {
00509 add( l->index,l->elem);
00510 l=l->next;
00511 }
00512 return *this;
00513 }
00514
00515 template <class I, class T>
00516 fourindex<I,T> & fourindex<I,T>::join(fourindex<I,T> &rhs)
00517 {
00518 if(nn!=rhs.nn) laerror("incompatible dimensions for join");
00519 if(*rhs.count!=1) laerror("shared rhs in join()");
00520 if(!count) {count=new int; *count=1; list=NULL;}
00521 else copyonwrite();
00522 matel4<I,T> **last=&list;
00523 while(*last) last= &((*last)->next);
00524 *last=rhs.list;
00525 rhs.list=NULL;
00526 return *this;
00527 }
00528
00529 template <class I, class T>
00530 void fourindex<I,T>::resize(const I n)
00531 {
00532 if(n<=0 ) laerror("illegal fourindex dimension");
00533 if(count)
00534 {
00535 if(*count > 1) {(*count)--; count=NULL; list=NULL;}
00536 else if(*count==1) deletelist();
00537 }
00538 nn=n;
00539 count=new int(1);
00540 list=NULL;
00541 }
00542
00543
00544 template <class I, class T>
00545 void fourindex<I,T>::deletelist()
00546 {
00547 if(*count >1) laerror("trying to delete shared list");
00548 matel4<I,T> *l=list;
00549 while(l)
00550 {
00551 matel4<I,T> *ltmp=l;
00552 l=l->next;
00553 delete ltmp;
00554 }
00555 list=NULL;
00556 delete count;
00557 count=NULL;
00558 }
00559
00560 template <class I, class T>
00561 void fourindex<I,T>::copylist(const matel4<I,T> *l)
00562 {
00563 list=NULL;
00564 while(l)
00565 {
00566 add(l->index,l->elem);
00567 l=l->next;
00568 }
00569 }
00570
00571 template <class I, class T>
00572 void fourindex<I,T>::copyonwrite()
00573 {
00574 if(!count) laerror("probably an assignment to undefined fourindex");
00575 if(*count > 1)
00576 {
00577 (*count)--;
00578 count = new int; *count=1;
00579 if(!list) laerror("empty list with count>1");
00580 copylist(list);
00581 }
00582 }
00583
00584 template <class I, class T>
00585 unsigned long fourindex<I,T>::length() const
00586 {
00587 unsigned long n=0;
00588 matel4<I,T> *l=list;
00589 while(l)
00590 {
00591 ++n;
00592 l=l->next;
00593 }
00594 return n;
00595 }
00596
00597 template <class I, class T>
00598 std::ostream& operator<<(std::ostream &s, const fourindex_ext<I,T> &x)
00599 {
00600 int n;
00601 n=x.size();
00602 s << n << '\n';
00603 typename fourindex_ext<I,T>::iterator it=x.begin();
00604 while(it!=x.end())
00605 {
00606 s << (typename LA_traits_io<I>::IOtype)it->index.indiv.i << ' ' << (typename LA_traits_io<I>::IOtype)it->index.indiv.j<< ' ' <<(typename LA_traits_io<I>::IOtype)it->index.indiv.k << ' ' << (typename LA_traits_io<I>::IOtype)it->index.indiv.l << ' ' << (typename LA_traits_io<T>::IOtype) it->elem << '\n';
00607 ++it;
00608 }
00609 s << "-1 -1 -1 -1\n";
00610 return s;
00611 }
00612
00613
00614
00615 template <class I, class T>
00616 std::ostream& operator<<(std::ostream &s, const fourindex<I,T> &x)
00617 {
00618 int n;
00619 n=x.size();
00620 s << n << '\n';
00621 typename fourindex<I,T>::iterator it=x.begin(),end=x.end();
00622 while(it!=end)
00623 {
00624 s << (typename LA_traits_io<I>::IOtype)it->index.indiv.i << ' ' << (typename LA_traits_io<I>::IOtype)it->index.indiv.j<< ' ' <<(typename LA_traits_io<I>::IOtype)it->index.indiv.k << ' ' << (typename LA_traits_io<I>::IOtype)it->index.indiv.l << ' ' << (typename LA_traits_io<T>::IOtype) it->elem << '\n';
00625 ++it;
00626 }
00627 s << "-1 -1 -1 -1\n";
00628 return s;
00629 }
00630
00631 template <class I, class T>
00632 std::istream& operator>>(std::istream &s, fourindex<I,T> &x)
00633 {
00634 typename LA_traits_io<I>::IOtype i,j,k,l;
00635 typename LA_traits_io<T>::IOtype elem;
00636 int n;
00637 s >> n ;
00638 x.resize(n);
00639 s >> i >> j >>k >>l;
00640 while(i!= (typename LA_traits_io<I>::IOtype)-1 && j!= (typename LA_traits_io<I>::IOtype)-1 && k != (typename LA_traits_io<I>::IOtype)-1 && l!= (typename LA_traits_io<I>::IOtype)-1)
00641 {
00642 s>>elem;
00643 x.add((I)i,(I)j,(I)k,(I)l,(T)elem);
00644 s >> i >> j >>k >>l;
00645 }
00646 return s;
00647 }
00648
00649
00651
00652
00653
00654
00655
00656
00657
00658
00659 template<fourindexsymtype S, class T, class DUMMY> class fourindex_dense;
00660
00661
00662 template<fourindexsymtype S, class T, class DUMMY>
00663 struct LA_traits<fourindex_dense<S,T,DUMMY> > {
00664 typedef T elementtype;
00665 typedef typename LA_traits<T>::normtype normtype;
00666 };
00667
00668
00669
00670 template<class T, class I>
00671 class fourindex_dense<twoelectronrealmullikan,T,I> : public NRSMat<T> {
00672 public:
00673 fourindex_dense(): NRSMat<T>() {};
00674 explicit fourindex_dense(const int n): NRSMat<T>(n*(n+1)/2) {};
00675 fourindex_dense(const NRSMat<T> &rhs): NRSMat<T>(rhs) {};
00676 fourindex_dense(const T &a, const int n): NRSMat<T>(a,n*(n+1)/2) {};
00677 fourindex_dense(const T *a, const int n): NRSMat<T>(a,n*(n+1)/2) {};
00678
00679
00680 fourindex_dense(const fourindex<I,T> &rhs);
00681 fourindex_dense(const fourindex_ext<I,T> &rhs);
00682
00683 T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l);
00684 const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
00685 void resize(const int n) {(*this).NRSMat<T>::resize(n*(n+1)/2);};
00686 void putext(int f, T thr=1e-15);
00687 int nbas() const {return (int)std::sqrt(2*(*this).nrows());};
00688
00689 };
00690
00691 template<class T, class I>
00692 void fourindex_dense<twoelectronrealmullikan,T,I>::putext(int f, T thr)
00693 {
00694 T y;
00695 for(int i=1; i<=nbas(); ++i) for(int j=1; j<=i; ++j)
00696 for(int k=1; k<=i; ++k) for(int l=1; l<=(i==k?j:k); ++l)
00697 if((y=abs((*this)(i,j,k,l))) > thr)
00698 {
00699 matel4stored<I,T> x;
00700 x.elem= y;
00701 x.index.indiv.i=i;
00702 x.index.indiv.j=j;
00703 x.index.indiv.k=k;
00704 x.index.indiv.l=l;
00705 if(sizeof(matel4stored<I,T>) != write(f,&x,sizeof(matel4stored<I,T>)) )
00706 laerror("write error in putext");
00707 }
00708 }
00709
00710
00711 template<class T, class I>
00712 fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense(const fourindex<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
00713 {
00714 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
00715 typename fourindex<I,T>::iterator p;
00716 #ifdef DEBUG
00717 unsigned int IJ = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
00718 unsigned int KL = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
00719 if (IJ<0 || IJ>=(unsigned int)NRSMat<T>::nn || KL<0 || KL>=(unsigned int)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
00720 #endif
00721 for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j,p->index.indiv.k,p->index.indiv.l) = p->elem;
00722 }
00723
00724 template<class T, class I>
00725 fourindex_dense<twoelectronrealmullikan,T,I>::fourindex_dense(const fourindex_ext<I,T> &rhs) : NRSMat<T>((T)0,rhs.size()*(rhs.size()+1)/2)
00726 {
00727 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
00728 typename fourindex_ext<I,T>::iterator p;
00729 for(p=rhs.begin(); p!= rhs.end(); ++p)
00730 {
00731 #ifdef DEBUG
00732 unsigned int IJ = SMat_index_1(p->index.indiv.i,p->index.indiv.j);
00733 unsigned int KL = SMat_index_1(p->index.indiv.k,p->index.indiv.l);
00734 if (IJ<0 || IJ>=(unsigned int)NRSMat<T>::nn || KL<0 || KL>=(unsigned int)NRSMat<T>::nn) laerror("fourindex_dense index out of range in constructor");
00735 #endif
00736 (*this)(p->index.indiv.i,p->index.indiv.j ,p->index.indiv.k,p->index.indiv.l) = p->elem;
00737 }
00738 }
00739
00740
00741
00742 template<class T, class DUMMY>
00743 T& fourindex_dense<twoelectronrealmullikan,T,DUMMY>::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l)
00744 {
00745 int I = SMat_index_1(i,j);
00746 int J = SMat_index_1(k,l);
00747
00748 #ifdef DEBUG
00749 if (*NRSMat<T>::count != 1) laerror("lval (i,j,k,l) with count > 1 in fourindex_dense");
00750 if (I<0 || I>=NRSMat<T>::nn || J<0 || J>=NRSMat<T>::nn) laerror("fourindex_dense index out of range");
00751 if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
00752 #endif
00753 return NRSMat<T>::v[SMat_index(I,J)];
00754 }
00755
00756
00757 template<class T, class DUMMY>
00758 const T& fourindex_dense<twoelectronrealmullikan,T,DUMMY>::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
00759 {
00760 int I = SMat_index_1(i,j);
00761 int J = SMat_index_1(k,l);
00762
00763 #ifdef DEBUG
00764 if (I<0 || I>=NRSMat<T>::nn || J<0 || J>=NRSMat<T>::nn) laerror("fourindex_dense index out of range");
00765 if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
00766 #endif
00767 return NRSMat<T>::v[SMat_index(I,J)];
00768 }
00769
00770
00771
00772
00773
00774 template<class T, class I>
00775 class fourindex_dense<T2IjAb_aces,T,I> : public NRMat<T> {
00776 protected:
00777 unsigned int noca,nocb,nvra,nvrb;
00778 friend class explicit_t2;
00779 public:
00780 fourindex_dense(): NRMat<T>() {noca=nocb=nvra=nvrb=0;};
00781 void resize(const int nocca, const int noccb, const int nvrta, const int nvrtb) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb; (*this).NRMat<T>::resize(nocca*noccb,nvrta*nvrtb);};
00782 explicit fourindex_dense(const int nocca, const int noccb, const int nvrta, const int nvrtb): NRMat<T>(nocca*noccb,nvrta*nvrtb) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb;};
00783
00784
00785 inline T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b)
00786 {
00787 #ifdef DEBUG
00788 if(i<1||i>noca ||j<1||j>nocb|| a<1||a>nvra||b<1||b>nvrb) laerror("T2IjAb_aces fourindex out of range");
00789 if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
00790 #endif
00791 return (*this).NRMat<T>::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1);
00792 }
00793 inline const T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) const
00794 {
00795 #ifdef DEBUG
00796 if(i<1||i>noca ||j<1||j>nocb|| a<1||a>nvra||b<1||b>nvrb) laerror("T2IjAb_aces fourindex out of range");
00797 if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
00798 #endif
00799 return (*this).NRMat<T>::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1);
00800 }
00801
00802 void print(std::ostream &out) const
00803 {
00804 unsigned int i,j,a,b;
00805 for(i=1; i<=noca; ++i) for(j=1; j<=nocb; ++j) for(a=1; a<=nvra; ++a) for(b=1; b<=nvrb; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
00806 }
00807 };
00808
00809
00810 template<class T, class I>
00811 class fourindex_dense<T2ijab_aces,T,I> : public NRMat<T> {
00812 protected:
00813 unsigned int nocc,nvrt,ntri;
00814 friend class explicit_t2;
00815 public:
00816 fourindex_dense(): NRMat<T>() {nocc=nvrt=ntri=0;};
00817 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;};
00818 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;};
00819
00820
00821
00822 inline T operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) const
00823 {
00824 #ifdef DEBUG
00825 if(i<1||i>nocc ||j<1||j>nocc|| a<1||a>nvrt||b<1||b>nvrt) laerror("T2ijab_aces fourindex out of range");
00826 if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
00827 #endif
00828 int minus=0;
00829 if(i==j||a==b) return (T)0;
00830 if(i<j) {minus++; unsigned int t=i; i=j; j=t;}
00831 if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
00832 T val=(*this).NRMat<T>::operator() ((i-2)*(i-1)/2+j-1, (a-2)*(a-1)/2+b-1);
00833 return (minus&1)? -val:val;
00834 }
00835 inline void set(unsigned int i, unsigned int j, unsigned int a, unsigned int b, T elem)
00836 {
00837 #ifdef DEBUG
00838 if(i<1||i>nocc ||j<1||j>nocc|| a<1||a>nvrt||b<1||b>nvrt) laerror("T2ijab_aces fourindex out of range");
00839 if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
00840 if(i==j||a==b && elem) laerror("antisymmetry violation in fourindex_dense<T2ijab_aces>");
00841 #endif
00842 int minus=0;
00843 if(i<j) {minus++; unsigned int t=i; i=j; j=t;}
00844 if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
00845 (*this).NRMat<T>::operator() ((i-2)*(i-1)/2+j-1, (a-2)*(a-1)/2+b-1) = minus? -elem : elem;
00846 }
00847 inline void add(unsigned int i, unsigned int j, unsigned int a, unsigned int b, T elem)
00848 {
00849 #ifdef DEBUG
00850 if(i<1||i>nocc ||j<1||j>nocc|| a<1||a>nvrt||b<1||b>nvrt) laerror("T2ijab_aces fourindex out of range");
00851 if (!NRMat<T>::v) laerror("access to unallocated fourindex_dense");
00852 if(i==j||a==b && elem) laerror("antisymmetry violation in fourindex_dense<T2ijab_aces>");
00853 #endif
00854 int minus=0;
00855 if(i<j) {minus++; unsigned int t=i; i=j; j=t;}
00856 if(a<b) {minus++; unsigned int t=a; a=b; b=t;}
00857 (*this).NRMat<T>::operator() ((i-2)*(i-1)/2+j-1, (a-2)*(a-1)/2+b-1) += minus? -elem : elem;
00858 }
00859
00860
00861 void print(std::ostream &out) const
00862 {
00863 unsigned int i,j,a,b;
00864 for(i=1; i<=nocc; ++i) for(j=1; j<i; ++j) for(a=1; a<=nvrt; ++a) for(b=1; b<a; ++b) out << i<<" "<<j<<" "<<a<<" "<<b<<" "<<(*this)(i,j,a,b)<<std::endl;
00865 }
00866
00867
00868 };
00869
00870
00871
00872
00873
00874 template<class T, class I>
00875 class fourindex_dense<antisymtwoelectronrealdirac,T,I> : public NRSMat<T> {
00876 private:
00877 int nbas;
00878 public:
00879 fourindex_dense(): NRSMat<T>() {};
00880 explicit fourindex_dense(const int n): nbas(n), NRSMat<T>(n*(n-1)/2) {};
00881 fourindex_dense(const T &a, const int n): nbas(n), NRSMat<T>(a,n*(n-1)/2) {};
00882 fourindex_dense(const T *a, const int n): nbas(n), NRSMat<T>(a,n*(n-1)/2) {};
00883
00884
00885 fourindex_dense(const fourindex<I,T> &rhs);
00886 fourindex_dense(const fourindex_ext<I,T> &rhs);
00887
00888 void set(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
00889 void add(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
00890 void add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem);
00891 const T& operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const;
00892 void resize(const int n) {nbas=n; (*this).NRSMat<T>::resize(n*(n-1)/2);};
00893 void print(std::ostream &out) const
00894 {
00895 unsigned int i,j,k,l;
00896 for(i=1; i<=nbas; ++i)
00897 for(k=1;k<i; ++k)
00898 for(j=1; j<=i; ++j)
00899 for(l=1; l<j && (j==i ? l<=k : 1); ++l)
00900 std::cout << i<<" "<<k<<" "<<j<<" "<<l<<" "<<(*this)(i,k,j,l)<<std::endl;
00901 }
00902
00903
00904 };
00905
00906
00907
00908
00909 template<class T, class DUMMY>
00910 const T& fourindex_dense<antisymtwoelectronrealdirac,T,DUMMY>::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const
00911 {
00912 int I = ASMat_index_1(i,j);
00913 int J = ASMat_index_1(k,l);
00914 if (I<0 || J<0) return 0.;
00915 #ifdef DEBUG
00916 if (I>=(unsigned int)NRSMat<T>::nn || J>=(unsigned int)NRSMat<T>::nn) laerror("index out of range");
00917 if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
00918 #endif
00919 return NRSMat<T>::v[SMat_index(I,J)];
00920 }
00921
00922
00923 template<class T, class DUMMY>
00924 void fourindex_dense<antisymtwoelectronrealdirac,T,DUMMY>::set(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem)
00925 {
00926 if(i<j) elem = -elem;
00927 if(k<l) elem = -elem;
00928 int I = ASMat_index_1(i,j);
00929 int J = ASMat_index_1(k,l);
00930 if (I<0 || J<0) laerror("assignment to nonexisting element");
00931 #ifdef DEBUG
00932 if (I>=NRSMat<T>::nn || J>=NRSMat<T>::nn) laerror("index out of range");
00933 if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
00934 #endif
00935 NRSMat<T>::v[SMat_index(I,J)] = elem;
00936 }
00937
00938
00939 template<class T, class DUMMY>
00940 void fourindex_dense<antisymtwoelectronrealdirac,T,DUMMY>::add(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem)
00941 {
00942 if(i<j) elem = -elem;
00943 if(k<l) elem = -elem;
00944 int I = ASMat_index_1(i,j);
00945 int J = ASMat_index_1(k,l);
00946 if (I<0 || J<0) laerror("assignment to nonexisting element");
00947 #ifdef DEBUG
00948 if (I>=NRSMat<T>::nn || J>=NRSMat<T>::nn) laerror("index out of range");
00949 if (!NRSMat<T>::v) laerror("access to unallocated fourindex_dense");
00950 #endif
00951 NRSMat<T>::v[SMat_index(I,J)] += elem;
00952 }
00953
00954 template<class T, class DUMMY>
00955 void fourindex_dense<antisymtwoelectronrealdirac,T,DUMMY>::add_unique(unsigned int i, unsigned int j, unsigned int k, unsigned int l, T elem)
00956 {
00957 if(i<=j || k<=l) return;
00958 int I = ASMat_index_1(i,j);
00959 int J = ASMat_index_1(k,l);
00960 if (I<0 || J<0 || I<J) return;
00961 NRSMat<T>::v[SMat_index(I,J)] += elem;
00962 }
00963
00964
00965 template<class T, class I>
00966 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)
00967 {
00968 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
00969 typename fourindex_ext<I,T>::piterator p;
00970 for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
00971 {
00972 I i=p->index.indiv.i;
00973 I j=p->index.indiv.j;
00974 I k=p->index.indiv.k;
00975 I l=p->index.indiv.l;
00976 add_unique(i,k,j,l,p->elem);
00977 add_unique(i,k,l,j,-p->elem);
00978 }
00979 }
00980
00981
00982 template<class T, class I>
00983 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)
00984 {
00985 if(rhs.getsymmetry() != twoelectronrealmullikan ) laerror("fourindex_dense symmetry mismatch");
00986 typename fourindex_ext<I,T>::piterator p;
00987 for(p= const_cast<fourindex_ext<I,T> *>(&rhs)->pbegin(); p.notend(); ++p)
00988 {
00989 I i=p->index.indiv.i;
00990 I j=p->index.indiv.j;
00991 I k=p->index.indiv.k;
00992 I l=p->index.indiv.l;
00993 add_unique(i,k,j,l,p->elem);
00994 add_unique(i,k,l,j,-p->elem);
00995 }
00996 }
00997
00998
00999 }
01000
01001 #endif