1 2 namespace sc { 3 4 class LocalOSSContribution { 5 private: 6 double * const gmat; 7 double * const gmata; 8 double * const gmatb; 9 10 double * const pmat; 11 double * const pmata; 12 double * const pmatb; 13 public: LocalOSSContribution(double * g,double * p,double * ga,double * pa,double * gb,double * pb)14 LocalOSSContribution(double *g, double *p, double *ga, double *pa, 15 double *gb, double *pb) : 16 gmat(g), gmata(ga), gmatb(gb), pmat(p), pmata(pa), pmatb(pb) {} ~LocalOSSContribution()17 ~LocalOSSContribution() {} 18 set_bound(double,double)19 void set_bound(double,double) {} 20 cont1(int ij,int kl,double val)21 inline void cont1(int ij, int kl, double val) { 22 gmat[ij] += val*pmat[kl]; 23 gmat[kl] += val*pmat[ij]; 24 } 25 cont2(int ij,int kl,double val)26 inline void cont2(int ij, int kl, double val) { 27 val *= 0.25; 28 gmat[ij] -= val*pmat[kl]; 29 gmat[kl] -= val*pmat[ij]; 30 31 gmata[ij] += val*pmata[kl]; 32 gmata[kl] += val*pmata[ij]; 33 34 gmatb[ij] += val*pmatb[kl]; 35 gmatb[kl] += val*pmatb[ij]; 36 37 val *= -3.0; 38 gmatb[ij] += val*pmata[kl]; 39 gmatb[kl] += val*pmata[ij]; 40 41 gmata[ij] += val*pmatb[kl]; 42 gmata[kl] += val*pmatb[ij]; 43 } 44 cont3(int ij,int kl,double val)45 inline void cont3(int ij, int kl, double val) { 46 val *= 0.5; 47 gmat[ij] -= val*pmat[kl]; 48 gmat[kl] -= val*pmat[ij]; 49 50 gmata[ij] += val*pmata[kl]; 51 gmata[kl] += val*pmata[ij]; 52 53 gmatb[ij] += val*pmatb[kl]; 54 gmatb[kl] += val*pmatb[ij]; 55 56 val *= -3.0; 57 gmata[ij] += val*pmatb[kl]; 58 gmata[kl] += val*pmatb[ij]; 59 60 gmatb[ij] += val*pmata[kl]; 61 gmatb[kl] += val*pmata[ij]; 62 } 63 cont4(int ij,int kl,double val)64 inline void cont4(int ij, int kl, double val) { 65 gmat[ij] += 0.75*val*pmat[kl]; 66 gmat[kl] += 0.75*val*pmat[ij]; 67 68 gmata[ij] += 0.25*val*pmata[kl]; 69 gmata[kl] += 0.25*val*pmata[ij]; 70 71 gmatb[ij] += 0.25*val*pmatb[kl]; 72 gmatb[kl] += 0.25*val*pmatb[ij]; 73 74 gmata[ij] -= 0.75*val*pmatb[kl]; 75 gmata[kl] -= 0.75*val*pmatb[ij]; 76 77 gmatb[ij] -= 0.75*val*pmata[kl]; 78 gmatb[kl] -= 0.75*val*pmata[ij]; 79 } 80 cont5(int ij,int kl,double val)81 inline void cont5(int ij, int kl, double val) { 82 val *= 0.5; 83 gmat[ij] += val*pmat[kl]; 84 gmat[kl] += val*pmat[ij]; 85 86 gmata[ij] += val*pmata[kl]; 87 gmata[kl] += val*pmata[ij]; 88 89 gmatb[ij] += val*pmatb[kl]; 90 gmatb[kl] += val*pmatb[ij]; 91 92 val *= -3.0; 93 gmata[ij] += val*pmatb[kl]; 94 gmata[kl] += val*pmatb[ij]; 95 96 gmatb[ij] += val*pmata[kl]; 97 gmatb[kl] += val*pmata[ij]; 98 } 99 }; 100 101 class LocalOSSEnergyContribution { 102 private: 103 double * const pmat; 104 double * const pmata; 105 double * const pmatb; 106 107 public: 108 double ec; 109 double ex; 110 set_bound(double,double)111 void set_bound(double,double) {}; 112 LocalOSSEnergyContribution(double * p,double * pa,double * pb)113 LocalOSSEnergyContribution(double *p, double *pa, double *pb) : 114 pmat(p), pmata(pa), pmatb(pb) { 115 ec=ex=0; 116 } ~LocalOSSEnergyContribution()117 ~LocalOSSEnergyContribution() {} 118 cont1(int ij,int kl,double val)119 inline void cont1(int ij, int kl, double val) { 120 ec += val*pmat[ij]*pmat[kl]; 121 } 122 cont2(int ij,int kl,double val)123 inline void cont2(int ij, int kl, double val) { 124 ex -= 0.25*val*(pmat[ij]*pmat[kl] + pmata[ij]*pmata[kl] + 125 pmatb[ij]*pmatb[kl]); 126 ex += 0.75*val*(pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]); 127 } 128 cont3(int ij,int kl,double val)129 inline void cont3(int ij, int kl, double val) { 130 ex -= 0.5*val*(pmat[ij]*pmat[kl] + pmata[ij]*pmata[kl] + 131 pmatb[ij]*pmatb[kl]); 132 ex += 1.5*val*(pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]); 133 } 134 cont4(int ij,int kl,double val)135 inline void cont4(int ij, int kl, double val) { 136 ec += val*pmat[ij]*pmat[kl]; 137 ex -= 0.25*val*(pmat[ij]*pmat[kl] + pmata[ij]*pmata[kl] + 138 pmatb[ij]*pmatb[kl]); 139 ex += 0.75*val*(pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]); 140 } 141 cont5(int ij,int kl,double val)142 inline void cont5(int ij, int kl, double val) { 143 ec += val*pmat[ij]*pmat[kl]; 144 ex -= 0.5*val*(pmat[ij]*pmat[kl] + pmata[ij]*pmata[kl] + 145 pmatb[ij]*pmatb[kl]); 146 ex += 1.5*val*(pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]); 147 } 148 }; 149 150 class LocalOSSGradContribution { 151 private: 152 double * const pmat; 153 double * const pmata; 154 double * const pmatb; 155 156 public: LocalOSSGradContribution(double * p,double * pa,double * pb)157 LocalOSSGradContribution(double *p, double *pa, double *pb) : 158 pmat(p), pmata(pa), pmatb(pb) {} ~LocalOSSGradContribution()159 ~LocalOSSGradContribution() {} 160 cont1(int ij,int kl)161 inline double cont1(int ij, int kl) { 162 return pmat[ij]*pmat[kl] + 163 0.5*(pmata[ij]*pmat[kl] + pmat[ij]*pmata[kl] + 164 pmatb[ij]*pmat[kl] + pmat[ij]*pmatb[kl]) + 165 0.25*(pmata[ij]*pmata[kl] + pmatb[ij]*pmatb[kl] + 166 pmata[ij]*pmatb[kl] + pmatb[ij]*pmata[kl]); 167 } 168 cont2(int ij,int kl)169 inline double cont2(int ij, int kl) { 170 return pmat[ij]*pmat[kl] + 171 0.5*(pmata[ij]*pmat[kl] + pmat[ij]*pmata[kl] + 172 pmatb[ij]*pmat[kl] + pmat[ij]*pmatb[kl] + 173 pmata[ij]*pmata[kl] + pmatb[ij]*pmatb[kl] - 174 pmata[ij]*pmatb[kl] - pmatb[ij]*pmata[kl]); 175 } 176 }; 177 178 } 179 180