• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

adc-connect / adcc / 5362537629

pending completion
5362537629

push

github

mfherbst
Bump version: 0.15.16 → 0.15.17

1487 of 2318 branches covered (64.15%)

Branch coverage included in aggregate %.

1 of 1 new or added line in 1 file covered. (100.0%)

7130 of 9470 relevant lines covered (75.29%)

191956.13 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

88.52
/libadcc/TensorImpl/ExpressionTree.cc
1
//
2
// Copyright (C) 2020 by the adcc authors
3
//
4
// This file is part of adcc.
5
//
6
// adcc is free software: you can redistribute it and/or modify
7
// it under the terms of the GNU General Public License as published
8
// by the Free Software Foundation, either version 3 of the License, or
9
// (at your option) any later version.
10
//
11
// adcc is distributed in the hope that it will be useful,
12
// but WITHOUT ANY WARRANTY; without even the implied warranty of
13
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
14
// GNU General Public License for more details.
15
//
16
// You should have received a copy of the GNU General Public License
17
// along with adcc. If not, see <http://www.gnu.org/licenses/>.
18
//
19

20
// This needs to be the first include ... libtensor reasons
21
#pragma GCC visibility push(default)
22
#include <libtensor/expr/eval/eval_exception.h>
23
#pragma GCC visibility pop
24
//
25
#include "../exceptions.hh"
26
#include "ExpressionTree.hh"
27

28
// Change visibility of libtensor singletons to public
29
#pragma GCC visibility push(default)
30
#include <libtensor/block_tensor/btod_add.h>
31
#include <libtensor/block_tensor/btod_set.h>
32
#include <libtensor/expr/btensor/eval_btensor_double.h>
33
#include <libtensor/expr/btensor/impl/eval_tree_builder_btensor.h>
34
#include <libtensor/expr/btensor/impl/tensor_from_node.h>
35
#include <libtensor/expr/dag/node_add.h>
36
#include <libtensor/expr/dag/node_assign.h>
37
#include <libtensor/expr/dag/node_ident.h>
38
#include <libtensor/expr/dag/node_transform.h>
39
#include <libtensor/expr/iface/node_ident_any_tensor.h>
40
#include <libtensor/expr/opt/opt_add_before_transf.h>
41
#include <libtensor/expr/opt/opt_merge_adjacent_add.h>
42
#include <libtensor/expr/opt/opt_merge_adjacent_transf.h>
43
#include <libtensor/expr/opt/opt_merge_equiv_ident.h>
44
#include <libtensor/linalg/BlasSequential.h>
45
#include <libtensor/symmetry/so_copy.h>
46
#pragma GCC visibility pop
47

48
namespace libadcc {
49
namespace lt = libtensor;
50

51
template <size_t N>
52
lt::expr::expr_rhs<N, scalar_type> ExpressionTree::attach_letters(
2,076,184✔
53
      const std::vector<std::shared_ptr<const lt::letter>>& letters) const {
54

55
  if (N != permutation.size()) {
2,076,184✔
56
    throw runtime_error("Internal error: Mismatch between permutation.size() == " +
×
57
                        std::to_string(permutation.size()) +
58
                        " and expr_rhs dimensionality " + std::to_string(N) + ".");
59
  }
60
  if (letters.size() != permutation.size()) {
2,076,184✔
61
    throw runtime_error("Internal error: Mismatch between permutation.size() == " +
×
62
                        std::to_string(permutation.size()) + " and letters size " +
63
                        std::to_string(letters.size()) + ".");
64
  }
65

66
  // The permutation is to be applied *after* the expression is evaluated,
67
  // so these letters need to pick up the inverse permutation.
68
  std::vector<const lt::letter*> let_ptrs;
2,076,184✔
69
  for (size_t i = 0; i < N; ++i) {
8,294,608✔
70
    auto it = std::find(permutation.begin(), permutation.end(), i);
6,218,424✔
71
    if (it == permutation.end()) {
6,218,424✔
72
      throw runtime_error("Internal error: Could not build inverse permutation");
×
73
    }
74
    const size_t inv_i = static_cast<size_t>(it - permutation.begin());
6,218,424✔
75

76
    let_ptrs.push_back(letters[inv_i].get());
6,218,424✔
77
  }
78

79
  // This copies both the expression as well as the letters to internal storage.
80
  return lt::expr::expr_rhs<N, scalar_type>(*tree_ptr, let_ptrs);
4,152,368✔
81
}
2,076,184✔
82

83
namespace {
84
lt::expr::expr_tree optimise_tree(const lt::expr::expr_tree& in) {
200,845✔
85
  lt::expr::expr_tree ret(in);
200,845✔
86
  lt::expr::opt_merge_equiv_ident(ret);
200,845✔
87
  lt::expr::opt_merge_adjacent_transf(ret);
200,845✔
88
  lt::expr::opt_add_before_transf(ret);
200,845✔
89
  lt::expr::opt_merge_adjacent_transf(ret);
200,845✔
90
  lt::expr::opt_merge_adjacent_add(ret);
200,845✔
91
  return ret;
200,845✔
92
}
×
93

94
bool is_linear_combination(const lt::expr::expr_tree& assigntree) {
214,493✔
95
  using namespace libtensor::expr;
96

97
  auto is_valid_summand = [&assigntree](size_t id) {
1,296,941✔
98
    if (assigntree.get_vertex(id).check_type<node_ident>()) return true;
578,298✔
99

100
    if (assigntree.get_vertex(id).check_type<node_transform_base>()) {
236,315✔
101
      const node_transform<scalar_type>& nt =
102
            assigntree.get_vertex(id).template recast_as<node_transform<scalar_type>>();
165,210✔
103
      for (size_t i = 0; i < nt.get_perm().size(); ++i) {
491,732✔
104
        if (i != nt.get_perm()[i]) return false;  // Need identity permutation
386,026✔
105
      }
106

107
      // Only one subnode for transformation.
108
      if (assigntree.get_edges_out(id).size() != 1) return false;
105,706✔
109
      const size_t tid = assigntree.get_edges_out(id)[0];
105,706✔
110

111
      if (assigntree.get_vertex(tid).check_type<node_ident>()) {
105,706✔
112
        return true;
83,742✔
113
      }
114
    }
115
    return false;
93,069✔
116
  };
214,493✔
117

118
  bool found_add = false;
214,493✔
119
  for (auto it = assigntree.begin(); it != assigntree.end(); ++it) {
820,211✔
120
    if (assigntree.get_vertex(it).check_type<node_add>()) {
758,291✔
121
      if (found_add) return false;  // There should be exactly one add
133,024✔
122
      auto add_edges = assigntree.get_edges_out(it);
80,160✔
123
      for (auto aed = add_edges.begin(); aed != add_edges.end(); ++aed) {
141,956✔
124
        if (!is_valid_summand(*aed)) return false;
114,660✔
125
      }
126
      found_add = true;
27,296✔
127
    } else if (assigntree.get_vertex(it).check_type<node_assign>()) {
758,291✔
128
      if (assigntree.get_root() != assigntree.get_id(it)) return false;
214,493✔
129
    } else if (!is_valid_summand(assigntree.get_id(it))) {
463,638✔
130
      return false;
99,709✔
131
    }
132
  }
133

134
  return found_add;
61,920✔
135
}
136

137
template <size_t N>
138
void evaluate_linear_combination(const lt::expr::expr_tree& assigntree,
13,648✔
139
                                 lt::btensor<N, scalar_type>& result, bool add) {
140
  using namespace libtensor::expr;
141
  if (!is_linear_combination(assigntree)) {
13,648✔
142
    throw invalid_argument(
×
143
          "evaluate_linear_combination got a tree which is not a linear combination "
144
          "tree");
145
  }
146

147
  lt::btensor<N, scalar_type>* first_btensor = nullptr;
13,648✔
148
  std::unique_ptr<lt::btod_add<N>> operator_ptr;
13,648✔
149
  auto add_summand = [&assigntree, &operator_ptr, &first_btensor](size_t id) {
41,870✔
150
    eval_btensor_double::btensor_from_node<N, scalar_type> extract(assigntree, id);
28,222✔
151
    auto& btensor = dynamic_cast<lt::btensor<N, scalar_type>&>(extract.get_btensor());
28,222✔
152
    // Notice: No permutation allowed here
153
    const lt::permutation<N>& perm = extract.get_transf().get_perm();
28,222✔
154
    if (!perm.is_identity()) {
28,222✔
155
      throw runtime_error("Internal error: Caught non-identity permutation.");
×
156
    }
157

158
    const double c = extract.get_transf().get_scalar_tr().get_coeff();
28,222✔
159
    if (!operator_ptr) {
28,222✔
160
      // Operator not yet set up -> do so with first summand
161
      operator_ptr.reset(new lt::btod_add<N>(btensor, c));
13,648✔
162
      first_btensor = &btensor;
13,648✔
163
    } else {
164
      // Operator is set up -> just add summand
165
      operator_ptr->add_op(btensor, c);
14,574✔
166
    }
167
  };
168

169
  auto root_edges = assigntree.get_edges_out(assigntree.get_root());
13,648✔
170
  for (auto it = root_edges.begin(); it != root_edges.end(); ++it) {
27,296✔
171
    if (assigntree.get_vertex(*it).check_type<node_add>()) {
27,296✔
172
      auto add_edges = assigntree.get_edges_out(*it);
13,648✔
173
      for (auto aed = add_edges.begin(); aed != add_edges.end(); ++aed) {
41,870✔
174
        add_summand(*aed);
28,222✔
175
      }
176
      break;  // Just one add tree
13,648✔
177
    }
13,648✔
178
  }
179

180
  // Copy symmmetry over
181
  if (first_btensor == nullptr) {
13,648✔
182
    throw runtime_error("Internal error: Got nullptr where set pointer was expected.");
×
183
  }
184
  lt::block_tensor_ctrl<N, scalar_type> ctrl_to(result);
13,648✔
185
  lt::block_tensor_ctrl<N, scalar_type> ctrl_from(*first_btensor);
13,648✔
186
  lt::so_copy<N, scalar_type>(ctrl_from.req_const_symmetry())
13,648✔
187
        .perform(ctrl_to.req_symmetry());
13,648✔
188

189
  if (!add) lt::btod_set<N>(0.0).perform(result);
13,648✔
190
  operator_ptr->perform(result, 1.0);
13,648✔
191
}
13,648✔
192

193
/** Return the expression tree resulting from assigning tree to
194
 *  the passed result. Optionally the assignment adds elements instead of setting them.
195
 */
196
template <size_t N>
197
libtensor::expr::expr_tree assignment_tree(const ExpressionTree& tree,
200,845✔
198
                                           libtensor::btensor<N, scalar_type>& result,
199
                                           bool add) {
200
  using namespace libtensor::expr;
201
  if (add) throw not_implemented_error("add = true not tested so far");
200,845✔
202

203
  // Form the RHS of the assignment
204
  std::vector<std::shared_ptr<const lt::letter>> label_safe;
200,845✔
205
  for (size_t i = 0; i < N; ++i) {
753,875✔
206
    label_safe.push_back(std::make_shared<lt::letter>());
553,030✔
207
  }
208
  expr_rhs<N, scalar_type> rhs = tree.attach_letters<N>(label_safe);
200,845✔
209

210
  // Form the labels for libtensor
211
  std::vector<const lt::letter*> label_unsafe;
200,845✔
212
  for (size_t i = 0; i < N; ++i) {
753,875✔
213
    label_unsafe.push_back(label_safe[i].get());
553,030✔
214
  }
215
  label<N> label(label_unsafe);
200,845✔
216

217
  // Build the assignment tree, starting from the assignment node
218
  node_assign nassign(N, add);
200,845✔
219
  expr_tree assigntree(nassign);
200,845✔
220
  expr_tree::node_id_t id = assigntree.get_root();
200,845✔
221
  node_ident_any_tensor<N, scalar_type> nresult(result);  // LHS of assignment
200,845✔
222
  assigntree.add(id, nresult);
200,845✔
223

224
  // Check if permutations are needed
225
  lt::permutation<N> px = label.permutation_of(rhs.get_label());
200,845✔
226
  if (!px.is_identity()) {
200,845✔
227
    std::vector<size_t> perm(N);
55,522✔
228
    for (size_t i = 0; i < N; i++) perm[i] = px[i];
207,630✔
229

230
    node_transform<scalar_type> n3(perm, lt::scalar_transf<scalar_type>());
55,522✔
231
    id = assigntree.add(id, n3);
55,522✔
232
  }
55,522✔
233

234
  // Add permuted or unpermuted RHS
235
  assigntree.add(id, rhs.get_expr());
200,845✔
236

237
  return assigntree;
401,690✔
238
}
200,845✔
239
}  // namespace
240

241
/** Return an optimised form of the internal expression tree */
242
lt::expr::expr_tree ExpressionTree::optimised_tree() { return optimise_tree(*tree_ptr); }
×
243

244
template <size_t N>
245
libtensor::expr::expr_tree ExpressionTree::evaluation_tree(
×
246
      libtensor::btensor<N, scalar_type>& result, bool add) const {
247
  // This does exactly the same stuff as
248
  // lt::expr::eval_btensor<double>().evaluate(assigntree) except doing the evaluation.
249
  lt::expr::expr_tree assigntree = assignment_tree(*this, result, add);
×
250
  lt::expr::eval_tree_builder_btensor bld(assigntree);
×
251
  bld.build();
×
252
  return bld.get_tree();
×
253
}
×
254

255
template <size_t N>
256
void ExpressionTree::evaluate_to(lt::btensor<N, scalar_type>& result, bool add) const {
200,845✔
257
  lt::expr::expr_tree assigntree = optimise_tree(assignment_tree(*this, result, add));
200,845✔
258

259
  lt::BlasSequential seq;  // Switch to sequential BLAS for evaluations
200,845✔
260

261
  if (is_linear_combination(assigntree)) {
200,845✔
262
    evaluate_linear_combination(assigntree, result, add);
13,648✔
263
    return;
13,648✔
264
  }
265

266
  // Normal evaluation via libtensor
267
  lt::expr::eval_btensor<double>().evaluate(assigntree);
187,197✔
268
}
214,493✔
269

270
//
271
// Explicit instantiation
272
//
273

274
#define INSTANTIATE(DIM)                                               \
275
  template void ExpressionTree::evaluate_to(                           \
276
        libtensor::btensor<DIM, scalar_type>& result, bool add) const; \
277
                                                                       \
278
  template libtensor::expr::expr_tree ExpressionTree::evaluation_tree( \
279
        libtensor::btensor<DIM, scalar_type>& result, bool add) const;
280

281
INSTANTIATE(1)
282
INSTANTIATE(2)
283
INSTANTIATE(3)
284
INSTANTIATE(4)
285

286
#undef INSTANTIATE
287

288
}  // namespace libadcc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc