Function CreateSparseMatrixInverse

linalg/sparsematrix.cpp:2346–2401

shared_ptr<BaseMatrix> CreateSparseMatrixInverse(shared_ptr<const BaseSparseMatrix> baseA,
                                                 shared_ptr<BitArray> subset,
                                                 shared_ptr<const Array<int>> clusters)
{
  if (baseA->invcreator)
    return baseA->invcreator(const_pointer_cast<BaseSparseMatrix>(baseA), subset, clusters);

  auto index = BaseMatrix::invcreators.CheckIndex(baseA->GetInverseType());
  if(index != -1)
    return BaseMatrix::invcreators[index](const_pointer_cast<BaseSparseMatrix>(baseA), subset, clusters);
  
  auto [dynH, dynW] = baseA->EntrySizes();

  shared_ptr<BaseMatrix> inv = nullptr;

  if ( dynH == dynW )
  {
    Switch<MAX_SYS_DIM>(dynH - 1, [&](auto HM)
    {
      constexpr int H = HM + 1;

      if constexpr(H == 1)
      {
        if ( auto A = dynamic_pointer_cast<const SparseMatrix<double>>(baseA) )
        {
          inv = InverseSparseMatrixTM<double>(A, subset, clusters);
        }
        else
        {
          inv = InverseSparseMatrixTM<Complex>(dynamic_pointer_cast<const SparseMatrix<Complex>>(baseA), subset, clusters);
        }
      }
      else
      {
        if (auto A = dynamic_pointer_cast<const SparseMatrix<Mat<H, H, double>>>(baseA))
        {
          inv = InverseSparseMatrixTM(A, subset, clusters);
        }
        else
        {
          inv = InverseSparseMatrixTM(dynamic_pointer_cast<const SparseMatrix<Mat<H, H, Complex>>>(baseA), subset, clusters);
        }
      }
    });

    if ( inv == nullptr )
    {
      throw Exception("SparseMatrix inverse not available for block-size " + std::to_string(dynH) + " matrix");
    }
  }
  else
  {
    throw Exception("SparseMatrix inverse not available for non-square " + std::to_string(dynH) + "x" + std::to_string(dynW) + " matrix");
  }
  return inv;
}