QMCPACK
SmallMatrixDetCalculator.h
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2019 QMCPACK developers.
6 //
7 // File developed by: Peter Doak, doakpw@ornl.gov, Oak Ridge National Lab
8 // ChangMo Yang, nichthierwohne@gmail.com, University of Illinois at Urbana-Champaign
9 // Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
10 // Miguel Morales, moralessilva2@llnl.gov, Lawrence Livermore National Laboratory
11 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
12 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
13 //
14 // File created by: ChangMo Yang, nichthierwohne@gmail.com, University of Illinois at Urbana-Champaign
15 //////////////////////////////////////////////////////////////////////////////////////
16 
17 
18 /**@file
19  */
20 #ifndef QMCPLUSPLUS_MULTIDIRACDETERMINANTCALCULATOR_H
21 #define QMCPLUSPLUS_MULTIDIRACDETERMINANTCALCULATOR_H
22 
23 #include "OhmmsPETE/OhmmsMatrix.h"
24 #include "Numerics/DeterminantOperators.h"
25 
26 namespace qmcplusplus
27 {
28 /** Function class calculates multi determinant ratio from matrix elements
29  * needs to be the size of your result matrix or larger
30  * includes manual expansions for smaller evaluations
31  */
32 template<typename T>
33 struct SmallMatrixDetCalculator
34 {
35  std::vector<T> M;
36  std::vector<int> Pivot;
37 
38  void resize(int n)
39  {
40  M.resize(n * n);
41  Pivot.resize(n);
42  }
43 
44  static T evaluate(T a11, T a12, T a21, T a22) { return a11 * a22 - a21 * a12; }
45 
46  static T evaluate(T a11, T a12, T a13, T a21, T a22, T a23, T a31, T a32, T a33)
47  {
48  return (a11 * (a22 * a33 - a32 * a23) - a21 * (a12 * a33 - a32 * a13) + a31 * (a12 * a23 - a22 * a13));
49  }
50 
51  static T evaluate(T a11,
52  T a12,
53  T a13,
54  T a14,
55  T a21,
56  T a22,
57  T a23,
58  T a24,
59  T a31,
60  T a32,
61  T a33,
62  T a34,
63  T a41,
64  T a42,
65  T a43,
66  T a44)
67  {
68  return (a11 * (a22 * (a33 * a44 - a43 * a34) - a32 * (a23 * a44 - a43 * a24) + a42 * (a23 * a34 - a33 * a24)) -
69  a21 * (a12 * (a33 * a44 - a43 * a34) - a32 * (a13 * a44 - a43 * a14) + a42 * (a13 * a34 - a33 * a14)) +
70  a31 * (a12 * (a23 * a44 - a43 * a24) - a22 * (a13 * a44 - a43 * a14) + a42 * (a13 * a24 - a23 * a14)) -
71  a41 * (a12 * (a23 * a34 - a33 * a24) - a22 * (a13 * a34 - a33 * a14) + a32 * (a13 * a24 - a23 * a14)));
72  }
73 
74  static T evaluate(T a11,
75  T a12,
76  T a13,
77  T a14,
78  T a15,
79  T a21,
80  T a22,
81  T a23,
82  T a24,
83  T a25,
84  T a31,
85  T a32,
86  T a33,
87  T a34,
88  T a35,
89  T a41,
90  T a42,
91  T a43,
92  T a44,
93  T a45,
94  T a51,
95  T a52,
96  T a53,
97  T a54,
98  T a55)
99  {
100  return (
101  a11 *
102  (a22 * (a33 * (a44 * a55 - a54 * a45) - a43 * (a34 * a55 - a54 * a35) + a53 * (a34 * a45 - a44 * a35)) -
103  a32 * (a23 * (a44 * a55 - a54 * a45) - a43 * (a24 * a55 - a54 * a25) + a53 * (a24 * a45 - a44 * a25)) +
104  a42 * (a23 * (a34 * a55 - a54 * a35) - a33 * (a24 * a55 - a54 * a25) + a53 * (a24 * a35 - a34 * a25)) -
105  a52 * (a23 * (a34 * a45 - a44 * a35) - a33 * (a24 * a45 - a44 * a25) + a43 * (a24 * a35 - a34 * a25))) -
106  a21 *
107  (a12 * (a33 * (a44 * a55 - a54 * a45) - a43 * (a34 * a55 - a54 * a35) + a53 * (a34 * a45 - a44 * a35)) -
108  a32 * (a13 * (a44 * a55 - a54 * a45) - a43 * (a14 * a55 - a54 * a15) + a53 * (a14 * a45 - a44 * a15)) +
109  a42 * (a13 * (a34 * a55 - a54 * a35) - a33 * (a14 * a55 - a54 * a15) + a53 * (a14 * a35 - a34 * a15)) -
110  a52 * (a13 * (a34 * a45 - a44 * a35) - a33 * (a14 * a45 - a44 * a15) + a43 * (a14 * a35 - a34 * a15))) +
111  a31 *
112  (a12 * (a23 * (a44 * a55 - a54 * a45) - a43 * (a24 * a55 - a54 * a25) + a53 * (a24 * a45 - a44 * a25)) -
113  a22 * (a13 * (a44 * a55 - a54 * a45) - a43 * (a14 * a55 - a54 * a15) + a53 * (a14 * a45 - a44 * a15)) +
114  a42 * (a13 * (a24 * a55 - a54 * a25) - a23 * (a14 * a55 - a54 * a15) + a53 * (a14 * a25 - a24 * a15)) -
115  a52 * (a13 * (a24 * a45 - a44 * a25) - a23 * (a14 * a45 - a44 * a15) + a43 * (a14 * a25 - a24 * a15))) -
116  a41 *
117  (a12 * (a23 * (a34 * a55 - a54 * a35) - a33 * (a24 * a55 - a54 * a25) + a53 * (a24 * a35 - a34 * a25)) -
118  a22 * (a13 * (a34 * a55 - a54 * a35) - a33 * (a14 * a55 - a54 * a15) + a53 * (a14 * a35 - a34 * a15)) +
119  a32 * (a13 * (a24 * a55 - a54 * a25) - a23 * (a14 * a55 - a54 * a15) + a53 * (a14 * a25 - a24 * a15)) -
120  a52 * (a13 * (a24 * a35 - a34 * a25) - a23 * (a14 * a35 - a34 * a15) + a33 * (a14 * a25 - a24 * a15))) +
121  a51 *
122  (a12 * (a23 * (a34 * a45 - a44 * a35) - a33 * (a24 * a45 - a44 * a25) + a43 * (a24 * a35 - a34 * a25)) -
123  a22 * (a13 * (a34 * a45 - a44 * a35) - a33 * (a14 * a45 - a44 * a15) + a43 * (a14 * a35 - a34 * a15)) +
124  a32 * (a13 * (a24 * a45 - a44 * a25) - a23 * (a14 * a45 - a44 * a15) + a43 * (a14 * a25 - a24 * a15)) -
125  a42 * (a13 * (a24 * a35 - a34 * a25) - a23 * (a14 * a35 - a34 * a15) + a33 * (a14 * a25 - a24 * a15))));
126  }
127 
128  /** default implementation of MultiDiracDeterminant::CustomizedMatrixDet
129  * If you don't have a perfect square below 25 of dots this is what you hit
130  * dots must be a 2n by 2n Matrix
131  * iter must be size n^2
132  * this object must have been resized to n
133  */
134  template<typename DT>
135  using OffloadMatrix = Matrix<DT, OffloadPinnedAllocator<DT>>;
136  template<typename ITER>
137  inline T evaluate(const OffloadMatrix<T>& dots, ITER it, int n)
138  {
139  typename std::vector<T>::iterator d = M.begin();
140  for (int i = 0; i < n; i++)
141  for (int j = 0; j < n; j++)
142  {
143  //performance through proper iterator indistiquishable from data pointer
144  *(d++) = dots(*(it + i), *(it + n + j));
145  }
146  return Determinant(M.data(), n, n, Pivot.data());
147  }
148 };
149 
150 template<unsigned NEXCITED>
151 class CustomizedMatrixDet;
152 
153 inline size_t flat_idx(const int i, const int a, const int n) { return a + i * n; }
154 
155 template<>
156 class CustomizedMatrixDet<1>
157 {
158 public:
159  template<typename VALUE>
160  static VALUE evaluate(const VALUE* table_matrix, const int* it, const size_t nb_cols)
161  {
162  const int i = *it;
163  const int a = *(it + 1);
164  const int n = nb_cols;
165  return table_matrix[flat_idx(i, a, nb_cols)];
166  }
167 };
168 
169 template<>
170 class CustomizedMatrixDet<2>
171 {
172 public:
173  template<typename VALUE>
174  static VALUE evaluate(const VALUE* table_matrix, const int* it, const size_t nb_cols)
175  {
176  const int i = *it;
177  const int j = *(it + 1);
178  const int a = *(it + 2);
179  const int b = *(it + 3);
180  const int n = nb_cols;
181  return SmallMatrixDetCalculator<VALUE>::evaluate(table_matrix[flat_idx(i, a, nb_cols)],
182  table_matrix[flat_idx(i, b, nb_cols)],
183  table_matrix[flat_idx(j, a, nb_cols)],
184  table_matrix[flat_idx(j, b, nb_cols)]);
185  }
186 };
187 
188 template<>
189 class CustomizedMatrixDet<3>
190 {
191 public:
192  template<typename VALUE>
193  static VALUE evaluate(const VALUE* table_matrix, const int* it, const size_t nb_cols)
194  {
195  const int i1 = *it;
196  const int i2 = *(it + 1);
197  const int i3 = *(it + 2);
198  const int a1 = *(it + 3);
199  const int a2 = *(it + 4);
200  const int a3 = *(it + 5);
201  return SmallMatrixDetCalculator<VALUE>::evaluate(table_matrix[flat_idx(i1, a1, nb_cols)],
202  table_matrix[flat_idx(i1, a2, nb_cols)],
203  table_matrix[flat_idx(i1, a3, nb_cols)],
204  table_matrix[flat_idx(i2, a1, nb_cols)],
205  table_matrix[flat_idx(i2, a2, nb_cols)],
206  table_matrix[flat_idx(i2, a3, nb_cols)],
207  table_matrix[flat_idx(i3, a1, nb_cols)],
208  table_matrix[flat_idx(i3, a2, nb_cols)],
209  table_matrix[flat_idx(i3, a3, nb_cols)]);
210  }
211 };
212 
213 template<>
214 class CustomizedMatrixDet<4>
215 {
216 public:
217  template<typename VALUE>
218  static VALUE evaluate(const VALUE* table_matrix, const int* it, const size_t nb_cols)
219  {
220  const int i1 = *it;
221  const int i2 = *(it + 1);
222  const int i3 = *(it + 2);
223  const int i4 = *(it + 3);
224  const int a1 = *(it + 4);
225  const int a2 = *(it + 5);
226  const int a3 = *(it + 6);
227  const int a4 = *(it + 7);
228  return SmallMatrixDetCalculator<VALUE>::evaluate(table_matrix[flat_idx(i1, a1, nb_cols)],
229  table_matrix[flat_idx(i1, a2, nb_cols)],
230  table_matrix[flat_idx(i1, a3, nb_cols)],
231  table_matrix[flat_idx(i1, a4, nb_cols)],
232  table_matrix[flat_idx(i2, a1, nb_cols)],
233  table_matrix[flat_idx(i2, a2, nb_cols)],
234  table_matrix[flat_idx(i2, a3, nb_cols)],
235  table_matrix[flat_idx(i2, a4, nb_cols)],
236  table_matrix[flat_idx(i3, a1, nb_cols)],
237  table_matrix[flat_idx(i3, a2, nb_cols)],
238  table_matrix[flat_idx(i3, a3, nb_cols)],
239  table_matrix[flat_idx(i3, a4, nb_cols)],
240  table_matrix[flat_idx(i4, a1, nb_cols)],
241  table_matrix[flat_idx(i4, a2, nb_cols)],
242  table_matrix[flat_idx(i4, a3, nb_cols)],
243  table_matrix[flat_idx(i4, a4, nb_cols)]);
244  }
245 };
246 
247 template<>
248 class CustomizedMatrixDet<5>
249 {
250 public:
251  template<typename VALUE>
252  static VALUE evaluate(const VALUE* table_matrix, const int* it, const size_t nb_cols)
253  {
254  const int i1 = *it;
255  const int i2 = *(it + 1);
256  const int i3 = *(it + 2);
257  const int i4 = *(it + 3);
258  const int i5 = *(it + 4);
259  const int a1 = *(it + 5);
260  const int a2 = *(it + 6);
261  const int a3 = *(it + 7);
262  const int a4 = *(it + 8);
263  const int a5 = *(it + 9);
264  return SmallMatrixDetCalculator<
265  VALUE>::evaluate(table_matrix[flat_idx(i1, a1, nb_cols)], table_matrix[flat_idx(i1, a2, nb_cols)],
266  table_matrix[flat_idx(i1, a3, nb_cols)], table_matrix[flat_idx(i1, a4, nb_cols)],
267  table_matrix[flat_idx(i1, a5, nb_cols)], table_matrix[flat_idx(i2, a1, nb_cols)],
268  table_matrix[flat_idx(i2, a2, nb_cols)], table_matrix[flat_idx(i2, a3, nb_cols)],
269  table_matrix[flat_idx(i2, a4, nb_cols)], table_matrix[flat_idx(i2, a5, nb_cols)],
270  table_matrix[flat_idx(i3, a1, nb_cols)], table_matrix[flat_idx(i3, a2, nb_cols)],
271  table_matrix[flat_idx(i3, a3, nb_cols)], table_matrix[flat_idx(i3, a4, nb_cols)],
272  table_matrix[flat_idx(i3, a5, nb_cols)], table_matrix[flat_idx(i4, a1, nb_cols)],
273  table_matrix[flat_idx(i4, a2, nb_cols)], table_matrix[flat_idx(i4, a3, nb_cols)],
274  table_matrix[flat_idx(i4, a4, nb_cols)], table_matrix[flat_idx(i4, a5, nb_cols)],
275  table_matrix[flat_idx(i5, a1, nb_cols)], table_matrix[flat_idx(i5, a2, nb_cols)],
276  table_matrix[flat_idx(i5, a3, nb_cols)], table_matrix[flat_idx(i5, a4, nb_cols)],
277  table_matrix[flat_idx(i5, a5, nb_cols)]);
278  }
279 };
280 
281 template<typename VALUE>
282 inline VALUE calcSmallDeterminant(size_t n, const VALUE* table_matrix, const int* it, const size_t nb_cols)
283 {
284  switch (n)
285  {
286  case 1:
287  return CustomizedMatrixDet<1>::evaluate(table_matrix, it, nb_cols);
288  case 2:
289  return CustomizedMatrixDet<2>::evaluate(table_matrix, it, nb_cols);
290  case 3:
291  return CustomizedMatrixDet<3>::evaluate(table_matrix, it, nb_cols);
292  case 4:
293  return CustomizedMatrixDet<4>::evaluate(table_matrix, it, nb_cols);
294  case 5:
295  return CustomizedMatrixDet<5>::evaluate(table_matrix, it, nb_cols);
296  default:
297  throw std::runtime_error("calculateSmallDeterminant only handles the number of excitations <= 5.");
298  }
299 }
300 
301 
302 } // namespace qmcplusplus
303 #endif
qmcplusplus::SmallMatrixDetCalculator::evaluate
static T evaluate(T a11, T a12, T a13, T a21, T a22, T a23, T a31, T a32, T a33)
Definition: SmallMatrixDetCalculator.h:46
qmcplusplus::SmallMatrixDetCalculator
Function class calculates multi determinant ratio from matrix elements needs to be the size of your r...
Definition: SmallMatrixDetCalculator.h:33
qmcplusplus
helper functions for EinsplineSetBuilder
Definition: Configuration.h:43
qmcplusplus::Determinant
T Determinant(T *restrict x, int n, int m, int *restrict pivot)
determinant of a matrix
Definition: DeterminantOperators.h:141
qmcplusplus::CustomizedMatrixDet< 2 >::evaluate
static VALUE evaluate(const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:174
qmcplusplus::calcSmallDeterminant
VALUE calcSmallDeterminant(size_t n, const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:282
qmcplusplus::SmallMatrixDetCalculator::evaluate
T evaluate(const OffloadMatrix< T > &dots, ITER it, int n)
Definition: SmallMatrixDetCalculator.h:137
qmcplusplus::CustomizedMatrixDet< 4 >::evaluate
static VALUE evaluate(const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:218
qmcplusplus::SmallMatrixDetCalculator::evaluate
static T evaluate(T a11, T a12, T a13, T a14, T a15, T a21, T a22, T a23, T a24, T a25, T a31, T a32, T a33, T a34, T a35, T a41, T a42, T a43, T a44, T a45, T a51, T a52, T a53, T a54, T a55)
Definition: SmallMatrixDetCalculator.h:74
qmcplusplus::SmallMatrixDetCalculator::evaluate
static T evaluate(T a11, T a12, T a21, T a22)
Definition: SmallMatrixDetCalculator.h:44
qmcplusplus::CustomizedMatrixDet< 3 >::evaluate
static VALUE evaluate(const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:193
DeterminantOperators.h
Define determinant operators.
qmcplusplus::SmallMatrixDetCalculator::evaluate
static T evaluate(T a11, T a12, T a13, T a14, T a21, T a22, T a23, T a24, T a31, T a32, T a33, T a34, T a41, T a42, T a43, T a44)
Definition: SmallMatrixDetCalculator.h:51
qmcplusplus::evaluate
void evaluate(Matrix< T, Alloc > &lhs, const Op &op, const Expression< RHS > &rhs)
Definition: OhmmsMatrix.h:514
qmcplusplus::CustomizedMatrixDet< 1 >::evaluate
static VALUE evaluate(const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:160
qmcplusplus::flat_idx
size_t flat_idx(const int i, const int a, const int n)
Definition: SmallMatrixDetCalculator.h:153
qmcplusplus::CustomizedMatrixDet< 5 >::evaluate
static VALUE evaluate(const VALUE *table_matrix, const int *it, const size_t nb_cols)
Definition: SmallMatrixDetCalculator.h:252