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