BsplineOneDimSolvers.h File Reference

Include dependency graph for BsplineOneDimSolvers.h:

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Classes
struct	SolvePeriodicInterp1D< T >
	dummy declaration to be specialized More...
struct	SolvePeriodicInterp1D< double >
struct	SolvePeriodicInterp1D< std::complex< double > >
struct	SolveFirstDerivInterp1D< T >
struct	SolveFirstDerivInterp1D< double >
struct	SolveFirstDerivInterp1D< std::complex< double > >
struct	SolveFirstDerivInterp1D< float >
struct	SolveFirstDerivInterp1D< std::complex< float > >

Functions
template<typename T >
void	FuncSolvePeriodicInterp1D (const std::vector< T > &data, std::vector< T > &p)

---

Class Documentation

SolvePeriodicInterp1D

struct SolvePeriodicInterp1D

template<typename T>
struct SolvePeriodicInterp1D< T >

dummy declaration to be specialized

Specializations for double and std::complex<double>

Definition at line 24 of file BsplineOneDimSolvers.h.

Collaboration diagram for SolvePeriodicInterp1D< T >:

SolveFirstDerivInterp1D

struct SolveFirstDerivInterp1D

template<typename T>
struct SolveFirstDerivInterp1D< T >

Definition at line 130 of file BsplineOneDimSolvers.h.

Collaboration diagram for SolveFirstDerivInterp1D< T >:

Function Documentation

FuncSolvePeriodicInterp1D()

void FuncSolvePeriodicInterp1D ( const std::vector< T > &  data, std::vector< T > &  p  )

inline

Definition at line 28 of file BsplineOneDimSolvers.h.

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

Referenced by SolvePeriodicInterp1D< double >::apply().

{
  const T ratio = 0.25;
  int N         = data.size();
  std::vector<T> d(N), gamma(N), mu(N);
  //d = 1.5*data;
  for (int i = 0; i < N; i++)
    d[i] = 1.5 * data[i];
  // First, eliminate leading coefficients
  gamma[0]     = ratio;
  mu[0]        = ratio;
  mu[N - 1]    = ratio;
  gamma[N - 1] = 1.0;
  for (int row = 1; row < (N - 1); row++)
  {
    double diag    = 1.0 - mu[row - 1] * ratio;
    double diagInv = 1.0 / diag;
    gamma[row]     = -ratio * gamma[row - 1] * diagInv;
    mu[row]        = diagInv * ratio;
    d[row]         = diagInv * (d[row] - ratio * d[row - 1]);
    // Last row
    d[N - 1] -= mu[N - 1] * d[row - 1];
    gamma[N - 1] -= mu[N - 1] * gamma[row - 1];
    mu[N - 1] = -mu[N - 1] * mu[row - 1];
  }
  // Last row:  gamma(N-1) hold diagonal element
  mu[N - 1] += ratio;
  gamma[N - 1] -= mu[N - 1] * (mu[N - 2] + gamma[N - 2]);
  d[N - 1] -= mu[N - 1] * d[N - 2];
  p[N - 1] = d[N - 1] / gamma[N - 1];
  // Now go back upward, back substituting
  for (int row = N - 2; row >= 0; row--)
    p[row] = d[row] - mu[row] * p[row + 1] - gamma[row] * p[N - 1];
}