1 //
2 // BAGEL - Brilliantly Advanced General Electronic Structure Library
3 // Filename: _vrr.h
4 // Copyright (C) 2012 Toru Shiozaki
5 //
6 // Author: Toru Shiozaki <shiozaki@northwestern.edu>
7 // Maintainer: Shiozaki group
8 //
9 // This file is part of the BAGEL package.
10 //
11 // This program is free software: you can redistribute it and/or modify
12 // it under the terms of the GNU General Public License as published by
13 // the Free Software Foundation, either version 3 of the License, or
14 // (at your option) any later version.
15 //
16 // This program is distributed in the hope that it will be useful,
17 // but WITHOUT ANY WARRANTY; without even the implied warranty of
18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU General Public License
22 // along with this program.  If not, see <http://www.gnu.org/licenses/>.
23 //
24 
25 
26 // replaces generated codes _vrr_xxxx.cc etc
27 
28 #ifndef __SRC_RYSINT____VRR_H
29 #define __SRC_RYSINT____VRR_H
30 
31 #include <algorithm>
32 #include <stdexcept>
33 
34 namespace bagel {
35 
36 template<int a_, int c_, int rank_, typename DataType>
vrr(DataType * data_,const DataType * C00,const DataType * D00,const DataType * B00,const DataType * B01,const DataType * B10)37 void vrr(DataType* data_, const DataType* C00, const DataType* D00, const DataType* B00, const DataType* B01, const DataType* B10) {
38 
39   static_assert(a_>=0 && c_>=0 && rank_ >= 1, "parameter(s) wrong in vrr");
40 
41   alignas(32) DataType C00_[rank_];
42   alignas(32) DataType D00_[rank_];
43   alignas(32) DataType B00_[rank_];
44   alignas(32) DataType B01_[rank_];
45   alignas(32) DataType B10_[rank_];
46   std::copy_n(C00, rank_, C00_);
47   std::copy_n(D00, rank_, D00_);
48   std::copy_n(B00, rank_, B00_);
49   std::copy_n(B01, rank_, B01_);
50   std::copy_n(B10, rank_, B10_);
51 
52   if ((a_ > 1 && c_ > 1) || (a_ == 1 && c_ >  1)) {
53     // c == 0
54     for (int t = 0; t != rank_; ++t)
55       data_[rank_*0+t] = 1.0;
56 
57     for (int t = 0; t != rank_; ++t)
58       data_[rank_*1+t] = C00_[t];
59 
60     if (a_ == 2) {
61       for (int t = 0; t != rank_; ++t)
62         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_[t];
63     } else if (a_ > 1) {
64       alignas(32) DataType B10_current[rank_];
65       // a == 2
66       for (int t = 0; t != rank_; ++t)
67         B10_current[t] = B10_[t];
68       for (int t = 0; t != rank_; ++t)
69         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_current[t];
70       // a > 2
71       for (int a = 3; a != a_+1; ++a) {
72         for (int t = 0; t != rank_; ++t)
73           B10_current[t] += B10_[t];
74         for (int t = 0; t != rank_; ++t)
75           data_[rank_*a+t] = C00_[t] * data_[rank_*(a-1)+t] + B10_current[t] * data_[rank_*(a-2)+t];
76       }
77     }
78     // c == 1
79     for (int t = 0; t != rank_; ++t)
80       data_[rank_*(a_+1)+t] = D00_[t];
81 
82     alignas(32) DataType cB00_current[rank_];
83     for (int t = 0; t != rank_; ++t)
84       cB00_current[t] = B00_[t];
85 
86     for (int t = 0; t != rank_; ++t)
87       data_[rank_*(a_+2)+t] = C00_[t] * data_[rank_*(a_+1)+t] + cB00_current[t];
88 
89     if (a_ == 2) {
90       for (int t = 0; t != rank_; ++t)
91         data_[rank_*(a_+3)+t] = C00_[t] * data_[rank_*(a_+2)+t] + B10_[t] * data_[rank_*(a_+1)+t] + cB00_current[t] * data_[rank_+t];
92     } else if (a_ > 1) {
93 
94       alignas(32) DataType B10_current[rank_];
95       for (int t = 0; t != rank_; ++t)
96         B10_current[t] = B10_[t];
97 
98       for (int t = 0; t != rank_; ++t)
99         data_[rank_*(a_+3)+t] = C00_[t] * data_[rank_*(a_+2)+t] + B10_current[t] * data_[rank_*(a_+1)+t] + cB00_current[t] * data_[rank_+t];
100 
101       for (int a = 3; a != a_+1; ++a) {
102         for (int t = 0; t != rank_; ++t)
103           B10_current[t] += B10_[t];
104 
105         for (int t = 0; t != rank_; ++t)
106          data_[rank_*(a_+1+a)+t] = C00_[t] * data_[rank_*(a_+a)+t] + B10_current[t] * data_[rank_*(a_-1+a)+t] + cB00_current[t] * data_[rank_*(a-1)+t];
107       }
108     }
109     // c > 1
110     alignas(32) DataType B01_current[rank_] = {0.0};
111     for (int c = 2; c != c_+1; ++c) {
112       for (int t = 0; t != rank_; ++t)
113         B01_current[t] += B01_[t];
114 
115       for (int t = 0; t != rank_; ++t)
116         data_[rank_*(a_+1)*c+t] = D00_[t] * data_[rank_*(a_+1)*(c-1)+t] + B01_current[t] * data_[rank_*(a_+1)*(c-2)+t];
117 
118       for (int t = 0; t != rank_; ++t)
119         cB00_current[t] += B00_[t];
120 
121       for (int t = 0; t != rank_; ++t)
122         data_[rank_*((a_+1)*c+1)+t] = C00_[t] * data_[rank_*((a_+1)*c)+t] + cB00_current[t] * data_[rank_*((a_+1)*(c-1))+t];
123 
124       if (a_ > 1) {
125         alignas(32) DataType B10_current[rank_];
126         for (int t = 0; t != rank_; ++t)
127           B10_current[t] = B10_[t];
128 
129         for (int t = 0; t != rank_; ++t)
130           data_[rank_*((a_+1)*c+2)+t] = C00_[t] * data_[rank_*((a_+1)*c+1)+t] + B10_current[t] * data_[rank_*(a_+1)*c+t] + cB00_current[t] * data_[rank_*((a_+1)*(c-1)+1)+t];
131 
132         for (int a = 3; a != a_+1; ++a) {
133           for (int t = 0; t != rank_; ++t)
134             B10_current[t] += B10_[t];
135 
136           for (int t = 0; t != rank_; ++t)
137             data_[rank_*((a_+1)*c+a)+t] = C00_[t] * data_[rank_*((a_+1)*c+a-1)+t] + B10_current[t] * data_[rank_*((a_+1)*c+a-2)+t] + cB00_current[t] * data_[rank_*((a_+1)*(c-1)+a-1)+t];
138         }
139       }
140     }
141 
142   } else if (a_ == 0 && c_ >  0) {
143     for (int t = 0; t != rank_; ++t)
144       data_[rank_*0+t] = 1.0;
145 
146     for (int t = 0; t != rank_; ++t)
147       data_[rank_*1+t] = D00_[t];
148 
149     if (c_ == 2) {
150       for (int t = 0; t != rank_; ++t)
151         data_[rank_*2+t] = D00_[t] * data_[rank_*1+t] + B01_[t];
152     } else if (c_ > 2) {
153       alignas(32) DataType B01_current[rank_];
154       // c == 2
155       for (int t = 0; t != rank_; ++t)
156         B01_current[t] = B01_[t];
157       for (int t = 0; t != rank_; ++t)
158         data_[rank_*2+t] = D00_[t] * data_[rank_+t] + B01_current[t];
159       // c > 2
160       for (int c = 3; c != c_+1; ++c) {
161         for (int t = 0; t != rank_; ++t)
162           B01_current[t] += B01_[t];
163 
164         for (int t = 0; t != rank_; ++t)
165           data_[rank_*c+t] = D00_[t] * data_[rank_*(c-1)+t] + B01_current[t] * data_[rank_*(c-2)+t];
166       }
167     }
168   } else if (a_ >  1 && c_ == 1) {
169     for (int t = 0; t != rank_; ++t)
170       data_[rank_*0+t] = 1.0;
171 
172     for (int t = 0; t != rank_; ++t)
173       data_[rank_*1+t] = C00_[t];
174 
175     if (a_ == 2) {
176       for (int t = 0; t != rank_; ++t)
177         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_[t];
178 
179       for (int t = 0; t != rank_; ++t)
180         data_[rank_*(a_+1)+t] = D00_[t];
181 
182       for (int t = 0; t != rank_; ++t)
183         data_[rank_*(a_+2)+t] = C00_[t] * data_[rank_*(a_+1)+t] + B00_[t];
184 
185       for (int t = 0; t != rank_; ++t)
186         data_[rank_*(a_+3)+t] = C00_[t] * data_[rank_*(a_+2)+t] + B10_[t] * data_[rank_*(a_+1)+t] + B00_[t] * data_[rank_+t];
187     } else {
188       alignas(32) DataType B10_current[rank_];
189       // a == 2
190       for (int t = 0; t != rank_; ++t)
191         B10_current[t] = B10_[t];
192       for (int t = 0; t != rank_; ++t)
193         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_current[t];
194       // a > 2
195       for (int a = 3; a != a_+1; ++a) {
196         for (int t = 0; t != rank_; ++t)
197           B10_current[t] += B10_[t];
198         for (int t = 0; t != rank_; ++t)
199           data_[rank_*a+t] = C00_[t] * data_[rank_*(a-1)+t] + B10_current[t] * data_[rank_*(a-2)+t];
200       }
201 
202       for (int t = 0; t != rank_; ++t)
203         data_[rank_*(a_+1)+t] = D00_[t];
204 
205       for (int t = 0; t != rank_; ++t)
206         data_[rank_*(a_+2)+t] = C00_[t] * data_[rank_*(a_+1)+t] + B00_[t];
207 
208       // a == 2
209       for (int t = 0; t != rank_; ++t)
210         B10_current[t] = B10_[t];
211       for (int t = 0; t != rank_; ++t)
212         data_[rank_*(a_+3)+t] = C00_[t] * data_[rank_*(a_+2)+t] + B10_current[t] * data_[rank_*(a_+1)+t] + B00_[t] * data_[rank_+t];
213 
214       for (int a = 3; a != a_+1; ++a) {
215         for (int t = 0; t != rank_; ++t)
216           B10_current[t] += B10_[t];
217         for (int t = 0; t != rank_; ++t)
218           data_[rank_*(a_+1+a)+t] = C00_[t] * data_[rank_*(a_+a)+t] + B10_current[t] * data_[rank_*(a_+a-1)+t] + B00_[t] * data_[rank_*(a-1)+t];
219       }
220     }
221 
222   } else if (a_ == 1 && c_ == 1) {
223     for (int t = 0; t != rank_; ++t)
224       data_[rank_*0+t] = 1.0;
225 
226     for (int t = 0; t != rank_; ++t)
227       data_[rank_*1+t] = C00_[t];
228 
229     for (int t = 0; t != rank_; ++t)
230       data_[rank_*2+t] = D00_[t];
231 
232     for (int t = 0; t != rank_; ++t)
233       data_[rank_*3+t] = C00_[t] * data_[rank_*2+t] + B00_[t];
234 
235   } else if (a_ >  0 && c_ == 0) {
236     for (int t = 0; t != rank_; ++t)
237       data_[rank_*0+t] = 1.0;
238 
239     for (int t = 0; t != rank_; ++t)
240       data_[rank_*1+t] = C00_[t];
241 
242     if (a_ == 2) {
243       for (int t = 0; t != rank_; ++t)
244         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_[t];
245     } else if (a_ > 2) {
246       alignas(32) DataType B10_current[rank_];
247       // a == 2
248       for (int t = 0; t != rank_; ++t)
249         B10_current[t] = B10_[t];
250       for (int t = 0; t != rank_; ++t)
251         data_[rank_*2+t] = C00_[t] * data_[rank_*1+t] + B10_current[t];
252       // a > 2
253       for (int a = 3; a != a_+1; ++a) {
254         for (int t = 0; t != rank_; ++t)
255           B10_current[t] += B10_[t];
256         for (int t = 0; t != rank_; ++t)
257           data_[rank_*a+t] = C00_[t] * data_[rank_*(a-1)+t] + B10_current[t] * data_[rank_*(a-2)+t];
258       }
259     }
260   } else if (a_ == 0 && c_ == 0) {
261     for (int t = 0; t != rank_; ++t)
262       data_[t] = 1.0;
263   } else {
264     throw std::runtime_error("strange");
265   }
266 }
267 
268 
269 }
270 #endif
271