A potential bug for BlockedTensor batched implementation

[lcyyork]

Here is the problem: suppose I want to compute C2["ijrs"] += batched("r", alpha * B["gar"] * B["gbs"] * T2["ijab"]); where r and s are composite indices. Let's say they are either block c or v. The above expression works fine when C2 contains all the blocks (i.e., x-x-c-c, x-x-c-v, x-x-v-c, x-x-v-v; x used for indices ij). If part of the C2 blocks are missing (e.g., C2 contains only x-x-core-core), the expression generate some strange results that do not seem to be random.

Here is the code that shows the problem.

#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>
#include "ambit/tensor.h"
#include "ambit/blocked_tensor.h"

using namespace std;
using namespace ambit;

ambit::TensorType tensor_type = CoreTensor;

void init_random(ambit::Tensor A) {
    std::vector<double>& data = A.data();
    for (size_t i = 0, size = data.size(); i < size; ++i) {
        data[i] = double(std::rand()) / double(RAND_MAX);
    }
}

std::vector<std::string> od_two_labels_hhpp() {
    std::vector<std::string> labels;
    for (const std::string& p : {"c", "a"}) {
        for (const std::string& q : {"c", "a"}) {
            for (const std::string& r : {"a", "v"}) {
                for (const std::string& s : {"a", "v"}) {
                    if (p == "a" && q == "a" && r == "a" && s == "a") {
                        continue;
                    }
                    labels.push_back(p + q + r + s);
                }
            }
        }
    }
    return labels;
}

int main() {
    ambit::initialize();

    ambit::BlockedTensor::reset_mo_spaces();
    ambit::BlockedTensor::set_expert_mode(true);

    std::vector<size_t> core_mos{0, 1};
    std::vector<size_t> actv_mos{2, 3};
    std::vector<size_t> virt_mos{4, 5, 6, 7};

    std::vector<size_t> aux_mos = std::vector<size_t>(10);
    std::iota(aux_mos.begin(), aux_mos.end(), 0);

    ambit::BlockedTensor::add_mo_space("c", "m,n", core_mos, NoSpin);
    ambit::BlockedTensor::add_mo_space("a", "u,v,w,x,y,z", actv_mos, NoSpin);
    ambit::BlockedTensor::add_mo_space("v", "e,f", virt_mos, NoSpin);

    ambit::BlockedTensor::add_composite_mo_space("h", "i,j,k,l", {"c", "a"});
    ambit::BlockedTensor::add_composite_mo_space("p", "a,b,c,d", {"a", "v"});
    ambit::BlockedTensor::add_composite_mo_space("g", "p,q,r,s,t,o", {"c", "a", "v"});

    ambit::BlockedTensor::add_mo_space("L", "g", aux_mos, NoSpin);

    auto B = ambit::BlockedTensor::build(tensor_type, "B", {"Lgg"});
    auto T2 = ambit::BlockedTensor::build(tensor_type, "T2", {"hhpp"});
    auto C2 = ambit::BlockedTensor::build(tensor_type, "C2", od_two_labels_hhpp());
    auto X2 = ambit::BlockedTensor::build(tensor_type, "X2", od_two_labels_hhpp());

    for (const std::string& block : B.block_labels()) {
        init_random(B.block(block));
    }
    for (const std::string& block : T2.block_labels()) {
        init_random(T2.block(block));
    }
    for (const std::string& block : C2.block_labels()) {
        init_random(C2.block(block));
    }
    X2["pqrs"] = C2["pqrs"];

    C2["ijes"] += batched("e", B["gae"] * B["gbs"] * T2["ijab"]);
    C2["ijus"] += batched("u", B["gau"] * B["gbs"] * T2["ijab"]);
    C2["ijms"] += batched("m", B["gam"] * B["gbs"] * T2["ijab"]);

    X2["ijrs"] += batched("r", B["gar"] * B["gbs"] * T2["ijab"]);

    std::cout << std::setprecision(15) << C2.norm() << std::endl;
    std::cout << std::setprecision(15) << X2.norm() << std::endl;

    return 0;
}