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