Skip to content

trunc_svd duplicates its kept-singular-value sum across dtype branches #19

Description

@ultimatile

In include/tci/cytnx_typed_tensor_impl.h, the trunc_svd overload taking chi_min, chi_max, target_trunc_err and s_min writes its trunc_err out-parameter as (frobenius_sq - kept_s2) / frobenius_sq, where frobenius_sq is the squared Frobenius norm taken before truncation and kept_s2 is the sum of squares of the retained singular values. Only kept_s2 needs the singular-value array, and because ptr_as<T> requires a concrete element type, that one accumulation is written once per dtype:

{
  double kept_s2 = 0.0;
  bool computed = false;
  if (frobenius_sq > 0.0 && bond_dim > 0) {
    if (s_backend.dtype() == cytnx::Type.Double) {
      auto* s_data = s_backend.template ptr_as<double>();
      for (bond_dim_t<TenT> i = 0; i < bond_dim; ++i) {
        kept_s2 += static_cast<double>(s_data[i]) * static_cast<double>(s_data[i]);
      }
      computed = true;
    } else if (s_backend.dtype() == cytnx::Type.Float) {
      auto* s_data = s_backend.template ptr_as<float>();
      for (bond_dim_t<TenT> i = 0; i < bond_dim; ++i) {
        kept_s2 += static_cast<double>(s_data[i]) * static_cast<double>(s_data[i]);
      }
      computed = true;
    }
  }
  if (computed) {
    trunc_err = std::clamp((frobenius_sq - kept_s2) / frobenius_sq, 0.0, 1.0);
  } else {
    trunc_err = 0.0;
  }
}

The two arms are character-for-character identical apart from the pointer type. Keeping them in step means any change to the summation — its order, its accumulator type, its bounds — has to be made twice, or it silently applies to one dtype only.

Proposal

Nothing in the loop depends on the element type except through s_data, so a generic lambda can deduce it and each arm shrinks to a call:

auto sum_kept_squares = [&](const auto* s_data) {
  double acc = 0.0;
  for (bond_dim_t<TenT> i = 0; i < bond_dim; ++i) {
    acc += static_cast<double>(s_data[i]) * static_cast<double>(s_data[i]);
  }
  return acc;
};
if (frobenius_sq > 0.0 && bond_dim > 0) {
  if (s_backend.dtype() == cytnx::Type.Double) {
    kept_s2 = sum_kept_squares(s_backend.template ptr_as<double>());
    computed = true;
  } else if (s_backend.dtype() == cytnx::Type.Float) {
    kept_s2 = sum_kept_squares(s_backend.template ptr_as<float>());
    computed = true;
  }
}

Assigning rather than accumulating into kept_s2 is equivalent here because the arms are mutually exclusive and nothing writes kept_s2 between its zero-initialization and this block.

The computed flag stays as is: the dispatch covers only Type.Double and Type.Float, so a singular-value tensor of any other dtype leaves both arms unrun, and the flag is what routes that case to trunc_err = 0.0 instead of dividing by an unset kept_s2.

Summation order, accumulator type and bounds are all unchanged, so no behavioural difference is intended.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions