Method Smooth

linalg/sparsecholesky.cpp:2073–2106

  template <class TM, class TV_ROW, class TV_COL>
  void SparseCholesky<TM, TV_ROW, TV_COL> :: 
  Smooth (BaseVector & u, const BaseVector & f, BaseVector & y) const
  {
    static Timer t("SparseCholesky::Smooth");
    RegionTimer reg(t);

    if(dynamic_pointer_cast<const SparseMatrixSymmetric<TM,TV>>(this->matrix.lock()))
      {
        // use the original one ...
        SparseFactorization::Smooth(u,f,y);
        return;
      }

    const FlatVector<TVX> fu = u.FV<TVX> ();
    FlatVector<TVX> fy = y.FV<TVX> ();
    
    Vector<TVX> hy(this->nused);
    auto spmat = dynamic_pointer_cast<const SparseMatrix<TM,TV,TV>> (this->matrix.lock());
    if(!spmat)
      throw Exception("A matrix not available any more, needed for Smooth!");
    auto & hmat = *spmat;
    
    ParallelFor (this->nused, [&] (int i)
                 {
                   hy(i) = fy(inv_order[i]) - hmat.RowTimesVector(inv_order[i], fu);
                 });
    
    SolveReordered(hy);

    ParallelFor (this->nused, [&] (int i)
                 {
                   fu(inv_order[i]) += hy(i);
                 });
  }