lup Struct Reference

Collaboration diagram for lup:

Static Public Member Functions
template<class Matrix , class Permutation >
static auto	decompose (Matrix &&A, Permutation &&P, double tol=std::numeric_limits< double >::epsilon())
template<class Matrix , class Permutation , class VectorSol >
static auto	solve (Matrix const &LU, Permutation const &P, VectorSol &&x) -> VectorSol &&

Static Private Member Functions
template<class Matrix , class Permutation , class Index >
static auto	permute_max_diagonal (Matrix &&LU, Permutation &&P, Index i)
template<class Permutation , class Vector >
static auto	permute (Permutation const &p, Vector &&data) -> Vector &&
template<class LUMatrix , class Vector >
static auto	lower_solve (LUMatrix const &LU, Vector &&x) -> Vector &&
template<class LUMatrix , class Vector >
static auto	upper_solve (LUMatrix const &LU, Vector &&x) -> Vector &&

Detailed Description

Definition at line 32 of file MatrixOps.h.

Member Function Documentation

decompose()

static auto decompose ( Matrix &&  A, Permutation &&  P, double  tol = std::numeric_limits<double>::epsilon()  )

inlinestatic

Definition at line 36 of file MatrixOps.h.

References qmcplusplus::Units::distance::A, omptarget::min(), qmcplusplus::Units::force::N, and permute_max_diagonal().

                                                                                                   {
  throw;
  std::iota(begin(P), end(P), typename std::decay_t<Permutation>::value_type{0});
  auto const N = std::min(size(A), size(~A));
  assert( P.size() >= N );

  auto&& ret = A({0, N}, {0, N});
  for(auto i : extension(ret)){
    if(lup::permute_max_diagonal(A, P, i) < tol) return A({0, i}, {0, i});

    for(auto&& row : A({i + 1, N})){
      auto&& urow = row({i + 1, N});
      std::transform(
        cbegin(urow), cend(urow), cbegin(A[i]({i + 1, N})), begin(urow), 
        [f = row[i] /= A[i][i]](auto const& a, auto const& b){return a - f*b;}
      );
    }
  }
  return A({0, N}, {0, N});
}

lower_solve()

static auto lower_solve ( LUMatrix const &  LU, Vector &&  x  ) -> Vector&&

inlinestaticprivate

Definition at line 92 of file MatrixOps.h.

References qmcplusplus::Units::force::N.

Referenced by solve().

                                                                   {
  assert(size(LU) <= size(x));
  auto const N = size(LU);
  for(typename LUMatrix::size_type i = 0; i != N; ++i){
    auto const& Lrowi = LU[i]({0, i});
    x[i] -= std::inner_product(begin(Lrowi), end(Lrowi), cbegin(x), 0.);
  }
  return std::forward<Vector>(x);
}

permute()

static auto permute ( Permutation const &  p, Vector &&  data  ) -> Vector&&

inlinestaticprivate

Definition at line 73 of file MatrixOps.h.

Referenced by solve().

                                                                    {
  assert(size(p) <= size(data));
  using index = typename Permutation::size_type;
  for(index i = 0; i != size(p); ++i){
    index k = p[i];
    for( ; k > i; k = p[k]){}
    index pk = p[k];
    if(k >=i and pk != i){
      auto const t = data[i];
      for( ; pk != i; k = pk, pk = p[k]){
        data[k] = data[pk];
      };
      data[k] = t;
    }
  }
  return std::forward<Vector>(data);
}

permute_max_diagonal()

static auto permute_max_diagonal ( Matrix &&  LU, Permutation &&  P, Index  i  )

inlinestaticprivate

Definition at line 65 of file MatrixOps.h.

References qmcplusplus::abs().

Referenced by decompose().

                                                                       {
  auto mi = std::max_element(begin(LU) + i, end(LU), [i](auto const& a, auto const& b){return std::abs(a[i]) < std::abs(b[i]);}) - begin(LU);
       swap(LU[i], LU[mi]); 
  std::swap(P [i], P [mi]);
  return std::abs(LU[i][i]);
}

solve()

static auto solve ( Matrix const &  LU, Permutation const &  P, VectorSol &&  x  ) -> VectorSol&&

inlinestatic

Definition at line 58 of file MatrixOps.h.

References lower_solve(), permute(), and upper_solve().

                                                                                       {
  return upper_solve(LU, lower_solve(LU, permute(P, x)));
}

upper_solve()

static auto upper_solve ( LUMatrix const &  LU, Vector &&  x  ) -> Vector&&

inlinestaticprivate

Definition at line 103 of file MatrixOps.h.

References qmcplusplus::Units::force::N.

Referenced by solve().

                                                                   {
  assert(size(LU) <= size(x));
  auto const N = size(LU);
  for(typename LUMatrix::size_type i = N - 1; i >= 0; --i){
    auto const& Urowi = LU[i]({i + 1, N});
    (x[i] -= std::inner_product(begin(Urowi), end(Urowi), cbegin(x) + i + 1, 0.)) /= LU[i][i];
  }
  return std::forward<Vector>(x);
}

---

The documentation for this struct was generated from the following file:

• /home/pk7/projects/qmc/for_cron_doxygen/qmcpack/src/QMCTools/ppconvert/src/common/MatrixOps.h