00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef _chemistry_qc_intv3_storage_h
00029 #define _chemistry_qc_intv3_storage_h
00030
00031 #ifdef __GNUC__
00032 #pragma interface
00033 #endif
00034
00035 #ifdef __cplusplus
00036
00037 #include <stddef.h>
00038 #include <util/class/class.h>
00039 #include <util/keyval/keyval.h>
00040 #include <util/container/eavlmmap.h>
00041
00042 namespace sc {
00043
00044
00045 #define SH_BITS 15 // the number of bits holding a shell index
00046 #define PE_BITS 1 // the number of bits holding a permutation
00047
00048 #define SH_MASK ((1<<SH_BITS)-1)
00049 #define PE_MASK ((1<<PE_BITS)-1)
00050
00051 #define SH0_SHIFT 0
00052 #define SH1_SHIFT (SH_BITS + SH0_SHIFT)
00053 #define P12_SHIFT (SH_BITS + SH1_SHIFT)
00054 #define P34_SHIFT (PE_BITS + P12_SHIFT)
00055 #define SH2_SHIFT 0
00056 #define SH3_SHIFT (SH_BITS + SH2_SHIFT)
00057 #define P13P24_SHIFT (SH_BITS + SH3_SHIFT)
00058 class IntegralKey {
00059 public:
00060 unsigned int sh0_sh1_p12_p34;
00061 unsigned int sh2_sh3_p13p24;
00062 public:
00063 IntegralKey(int,int,int,int,int,int,int);
00064 IntegralKey(const IntegralKey&);
00065 bool operator == (const IntegralKey &k) const {
00066 return (sh0_sh1_p12_p34 == k.sh0_sh1_p12_p34)
00067 && (sh2_sh3_p13p24 == k.sh2_sh3_p13p24);
00068 }
00069 bool operator < (const IntegralKey &k) const {
00070 if (sh0_sh1_p12_p34 < k.sh0_sh1_p12_p34) return true;
00071 else if (sh0_sh1_p12_p34 > k.sh0_sh1_p12_p34) return false;
00072 return sh2_sh3_p13p24 < k.sh2_sh3_p13p24;
00073 }
00074 int sh0() const { return (sh0_sh1_p12_p34>>SH0_SHIFT) & SH_MASK; }
00075 int sh1() const { return (sh0_sh1_p12_p34>>SH1_SHIFT) & SH_MASK; }
00076 int p12() const { return (sh0_sh1_p12_p34>>P12_SHIFT) & PE_MASK; }
00077 int p34() const { return (sh0_sh1_p12_p34>>P34_SHIFT) & PE_MASK; }
00078 int sh2() const { return (sh2_sh3_p13p24>>SH2_SHIFT) & SH_MASK; }
00079 int sh3() const { return (sh2_sh3_p13p24>>SH3_SHIFT) & SH_MASK; }
00080 int p13p24() const { return (sh2_sh3_p13p24>>P13P24_SHIFT) & PE_MASK; }
00081 };
00082
00083 inline std::ostream&
00084 operator << (std::ostream &o, const IntegralKey &k)
00085 {
00086 o << '[' << k.sh0()
00087 << ' ' << k.sh1()
00088 << ' ' << k.sh2()
00089 << ' ' << k.sh3()
00090 << ' ' << k.p12() << k.p34() << k.p13p24()
00091 << " (" << k.sh0_sh1_p12_p34 << ' ' << k.sh2_sh3_p13p24 << ')'
00092 << ']';
00093 return o;
00094 }
00095
00096 inline
00097 IntegralKey::IntegralKey(int sh1_, int sh2_, int sh3_, int sh4_,
00098 int p12_, int p34_, int p13p24_)
00099 {
00100 sh0_sh1_p12_p34 = (sh1_<<SH0_SHIFT)
00101 |(sh2_<<SH1_SHIFT)
00102 |(p12_<<P12_SHIFT)
00103 |(p34_<<P34_SHIFT);
00104 sh2_sh3_p13p24 = (sh3_<<SH2_SHIFT)
00105 |(sh4_<<SH3_SHIFT)
00106 |(p13p24_<<P13P24_SHIFT);
00107 }
00108
00109 inline
00110 IntegralKey::IntegralKey(const IntegralKey& ik)
00111 {
00112 sh0_sh1_p12_p34 = ik.sh0_sh1_p12_p34;
00113 sh2_sh3_p13p24 = ik.sh2_sh3_p13p24;
00114 }
00115
00116 inline int
00117 compare(const IntegralKey&k1, const IntegralKey&k2)
00118 {
00119 if (k1.sh0_sh1_p12_p34 < k2.sh0_sh1_p12_p34) return -1;
00120 else if (k1.sh0_sh1_p12_p34 > k2.sh0_sh1_p12_p34) return 1;
00121
00122 if (k1.sh2_sh3_p13p24 < k2.sh2_sh3_p13p24) return -1;
00123 else if (k1.sh2_sh3_p13p24 > k2.sh2_sh3_p13p24) return 1;
00124 else return 0;
00125 }
00126
00127 class IntegralLink {
00128 public:
00129 EAVLMMapNode<IntegralKey, IntegralLink> intlist;
00130 EAVLMMapNode<int, IntegralLink> costlist;
00131 int size;
00132 public:
00133 IntegralLink(IntegralKey& key, int cost, int size);
00134 static int size_to_actualsize(int size);
00135 ~IntegralLink();
00136 int actualsize() const;
00137 int hash() const;
00138 static int shells_to_hash(int,int,int,int);
00139 int cost() const { return costlist.key; }
00140 void print();
00141
00142
00143 double* buffer() { return (double*)&this[1]; }
00144 void* operator new(size_t, int);
00145 void operator delete(void*, int);
00146 void operator delete(void*);
00147 };
00148
00149 inline int
00150 IntegralLink::shells_to_hash(int sh1,int sh2,int sh3,int sh4)
00151 {
00152 return sh1 ^ (sh4<<4) ^ (sh2<<8) ^ (sh3<<12);
00153 }
00154
00155 inline int
00156 IntegralLink::hash() const
00157 {
00158 return shells_to_hash(intlist.key.sh0(),
00159 intlist.key.sh1(),
00160 intlist.key.sh2(),
00161 intlist.key.sh3());
00162 }
00163
00164 inline int
00165 IntegralLink::size_to_actualsize(int size)
00166 {
00167 return size*sizeof(double) + sizeof(IntegralLink) + sizeof(void*)*2;
00168 }
00169
00170 inline int
00171 IntegralLink::actualsize() const
00172 {
00173 return size_to_actualsize(size);
00174 }
00175
00176 class IntegralStorer: public DescribedClass {
00177 private:
00178 int table_size_;
00179 EAVLMMap<int,IntegralLink> costlist;
00180 EAVLMMap<IntegralKey,IntegralLink>* table_;
00181 int maxsize_;
00182 int currentsize_;
00183 int n_integrals_;
00184 int n_shellquart_;
00185 public:
00186 IntegralStorer();
00187 IntegralStorer(const Ref<KeyVal>&);
00188 ~IntegralStorer();
00189 void init(int nbytes);
00190 void done();
00191 IntegralLink *find(IntegralKey&);
00192 int should_store(int cost, int actualsize);
00193 void store(IntegralKey& key, const double *buf,
00194 int size, int cost, int actualsize);
00195 void print_stats();
00196 int table_size() const { return table_size_; }
00197 EAVLMMap<IntegralKey,IntegralLink>&table_entry(int i){return table_[i];}
00198 };
00199
00200 }
00201
00202 #endif
00203
00204 #endif
00205
00206
00207
00208
00209