implements delayed update on CPU using BLAS More...

Collaboration diagram for DelayedUpdate< T, T_FP >:

Public Member Functions
	DelayedUpdate ()
	default constructor More...

void	resize (int norb, int delay)
	resize the internal storage More...

template<typename TREAL >
void	invert_transpose (const Matrix< T > &logdetT, Matrix< T > &Ainv, std::complex< TREAL > &log_value)
	compute the inverse of the transpose of matrix A More...

void	initializeInv (const Matrix< T > &Ainv)
	initialize internal objects when Ainv is refreshed More...

int	getDelayCount () const

template<typename VVT >
void	getInvRow (const Matrix< T > &Ainv, int rowchanged, VVT &invRow)
	compute the row of up-to-date Ainv More...

template<typename VVT , typename RATIOT >
void	acceptRow (Matrix< T > &Ainv, int rowchanged, const VVT &psiV, const RATIOT ratio_new)
	accept a move with the update delayed More...

void	updateInvMat (Matrix< T > &Ainv)
	update the full Ainv and reset delay_count More...

Private Attributes
Matrix< T >	U
	orbital values of delayed electrons More...

Matrix< T >	V
	rows of Ainv corresponding to delayed electrons More...

Matrix< T >	Binv
	Matrix inverse of B, at maximum KxK. More...

Matrix< T >	tempMat
	scratch space, used during inverse update More...

Vector< T >	temp
	temporal scratch space used by SM-1 More...

Vector< T >	p
	new column of B More...

std::vector< int >	delay_list
	list of delayed electrons More...

int	delay_count
	current number of delays, increase one for each acceptance, reset to 0 after updating Ainv More...

DiracMatrix< T_FP >	detEng
	matrix inversion engine More...

Detailed Description

template<typename T, typename T_FP>
class qmcplusplus::DelayedUpdate< T, T_FP >

implements delayed update on CPU using BLAS

Template Parameters

T	base precision for most computation
T_FP	high precision for matrix inversion, T_FP >= T

Definition at line 29 of file DelayedUpdate.h.

Constructor & Destructor Documentation

◆ DelayedUpdate()

DelayedUpdate ( )

inline

default constructor

Definition at line 52 of file DelayedUpdate.h.

52 : delay_count(0) {}

qmcplusplus::DelayedUpdate::delay_count

int delay_count

current number of delays, increase one for each acceptance, reset to 0 after updating Ainv ...

Definition: DelayedUpdate.h:46

Member Function Documentation

◆ acceptRow()

void acceptRow	(	Matrix< T > &	Ainv,
		int	rowchanged,
		const VVT &	psiV,
		const RATIOT	ratio_new
	)

inline

accept a move with the update delayed

Parameters

Ainv	inverse matrix
rowchanged	the row id corresponding to the proposed electron
psiV	new orbital values

Before delay_count reaches the maximum delay, only Binv is updated with a recursive algorithm

Definition at line 125 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::Binv, Matrix< T, Alloc >::cols(), BLAS::cone, qmcplusplus::syclBLAS::copy_n(), BLAS::czero, Matrix< T, Alloc >::data(), Vector< T, Alloc >::data(), DelayedUpdate< T, T_FP >::delay_count, DelayedUpdate< T, T_FP >::delay_list, BLAS::gemv(), BLAS::ger(), DelayedUpdate< T, T_FP >::p, Matrix< T, Alloc >::rows(), DelayedUpdate< T, T_FP >::U, DelayedUpdate< T, T_FP >::updateInvMat(), and DelayedUpdate< T, T_FP >::V.

Referenced by qmcplusplus::TEST_CASE().

   {
     constexpr T cone(1);
     constexpr T czero(0);
     const int norb     = Ainv.rows();
     const int lda_Binv = Binv.cols();
     std::copy_n(Ainv[rowchanged], norb, V[delay_count]);
     std::copy_n(psiV.data(), norb, U[delay_count]);
     delay_list[delay_count] = rowchanged;
     // the new Binv is [[X Y] [Z sigma]]
     BLAS::gemv('T', norb, delay_count + 1, -cone, V.data(), norb, psiV.data(), 1, czero, p.data(), 1);
     // sigma
     const T sigma                  = static_cast<T>(RATIOT(1) / ratio_new);
     Binv[delay_count][delay_count] = sigma;
     // Y
     BLAS::gemv('T', delay_count, delay_count, sigma, Binv.data(), lda_Binv, p.data(), 1, czero,
                Binv.data() + delay_count, lda_Binv);
     // X
     BLAS::ger(delay_count, delay_count, cone, Binv[delay_count], 1, Binv.data() + delay_count, lda_Binv, Binv.data(),
               lda_Binv);
     // Z
     for (int i = 0; i < delay_count; i++)
       Binv[delay_count][i] *= sigma;
     delay_count++;
     // update Ainv when maximal delay is reached
     if (delay_count == lda_Binv)
       updateInvMat(Ainv);
   }

◆ getDelayCount()

int getDelayCount ( ) const

inline

Definition at line 90 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::delay_count.

90 { return delay_count; }

qmcplusplus::DelayedUpdate::delay_count

int delay_count

current number of delays, increase one for each acceptance, reset to 0 after updating Ainv ...

Definition: DelayedUpdate.h:46

◆ getInvRow()

void getInvRow	(	const Matrix< T > &	Ainv,
		int	rowchanged,
		VVT &	invRow
	)

inline

compute the row of up-to-date Ainv

Parameters

Ainv	inverse matrix
rowchanged	the row id corresponding to the proposed electron

Definition at line 97 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::Binv, Matrix< T, Alloc >::cols(), BLAS::cone, qmcplusplus::syclBLAS::copy_n(), BLAS::czero, Matrix< T, Alloc >::data(), Vector< T, Alloc >::data(), DelayedUpdate< T, T_FP >::delay_count, BLAS::gemv(), DelayedUpdate< T, T_FP >::p, Matrix< T, Alloc >::rows(), Matrix< T, Alloc >::size(), DelayedUpdate< T, T_FP >::U, and DelayedUpdate< T, T_FP >::V.

Referenced by qmcplusplus::TEST_CASE().

   {
     if (delay_count == 0)
     {
       // Ainv is fresh, directly access Ainv
       std::copy_n(Ainv[rowchanged], invRow.size(), invRow.data());
       return;
     }
     constexpr T cone(1);
     constexpr T czero(0);
     const int norb     = Ainv.rows();
     const int lda_Binv = Binv.cols();
     // save Ainv[rowchanged] to invRow
     std::copy_n(Ainv[rowchanged], norb, invRow.data());
     // multiply V (NxK) Binv(KxK) U(KxN) invRow right to the left
     BLAS::gemv('T', norb, delay_count, cone, U.data(), norb, invRow.data(), 1, czero, p.data(), 1);
     BLAS::gemv('N', delay_count, delay_count, -cone, Binv.data(), lda_Binv, p.data(), 1, czero, Binv[delay_count], 1);
     BLAS::gemv('N', norb, delay_count, cone, V.data(), norb, Binv[delay_count], 1, cone, invRow.data(), 1);
   }

◆ initializeInv()

void initializeInv ( const Matrix< T > & Ainv )

inline

initialize internal objects when Ainv is refreshed

Parameters

Ainv	inverse matrix

Definition at line 84 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::delay_count.

   {
     // safe mechanism
     delay_count = 0;
   }

◆ invert_transpose()

void invert_transpose	(	const Matrix< T > &	logdetT,
		Matrix< T > &	Ainv,
		std::complex< TREAL > &	log_value
	)

inline

compute the inverse of the transpose of matrix A

Parameters

logdetT	orbital value matrix
Ainv	inverse matrix

Definition at line 74 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::delay_count, and DelayedUpdate< T, T_FP >::detEng.

   {
     detEng.invert_transpose(logdetT, Ainv, log_value);
     // safe mechanism
     delay_count = 0;
   }

◆ resize()

void resize	(	int	norb,
		int	delay
	)

inline

resize the internal storage

Parameters

norb	number of electrons/orbitals
delay,maximum	delay 0<delay<=norb

Definition at line 58 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::Binv, DelayedUpdate< T, T_FP >::delay_list, DelayedUpdate< T, T_FP >::p, Matrix< T, Alloc >::resize(), Vector< T, Alloc >::resize(), DelayedUpdate< T, T_FP >::temp, DelayedUpdate< T, T_FP >::tempMat, DelayedUpdate< T, T_FP >::U, and DelayedUpdate< T, T_FP >::V.

Referenced by qmcplusplus::TEST_CASE().

   {
     V.resize(delay, norb);
     U.resize(delay, norb);
     p.resize(delay);
     temp.resize(norb);
     tempMat.resize(norb, delay);
     Binv.resize(delay, delay);
     delay_list.resize(delay);
   }

◆ updateInvMat()

void updateInvMat ( Matrix< T > & Ainv )

inline

update the full Ainv and reset delay_count

Parameters

Ainv	inverse matrix

Definition at line 157 of file DelayedUpdate.h.

References DelayedUpdate< T, T_FP >::Binv, Matrix< T, Alloc >::cols(), BLAS::cone, BLAS::czero, Matrix< T, Alloc >::data(), Vector< T, Alloc >::data(), DelayedUpdate< T, T_FP >::delay_count, DelayedUpdate< T, T_FP >::delay_list, BLAS::gemm(), BLAS::gemv(), BLAS::ger(), getNextLevelNumThreads(), omptarget::min(), BlasThreadingEnv::NestedThreadingSupported(), Matrix< T, Alloc >::rows(), DelayedUpdate< T, T_FP >::temp, DelayedUpdate< T, T_FP >::tempMat, DelayedUpdate< T, T_FP >::U, and DelayedUpdate< T, T_FP >::V.

Referenced by DelayedUpdate< T, T_FP >::acceptRow().

   {
     if (delay_count == 0)
       return;
     // update the inverse matrix
     constexpr T cone(1);
     constexpr T czero(0);
     const int norb = Ainv.rows();
     if (delay_count == 1)
     {
       // this is a special case invoking the Fahy's variant of Sherman-Morrison update.
       // Only use the first norb elements of tempMat as a temporal array
       BLAS::gemv('T', norb, norb, cone, Ainv.data(), norb, U[0], 1, czero, temp.data(), 1);
       temp[delay_list[0]] -= cone;
       BLAS::ger(norb, norb, -Binv[0][0], V[0], 1, temp.data(), 1, Ainv.data(), norb);
     }
     else
     {
       const int lda_Binv = Binv.cols();
       // number of threads at the next level, forced to 1 if the problem is small.
       const int num_threads = (norb < 256 ? 1 : getNextLevelNumThreads());
       if (num_threads == 1 || BlasThreadingEnv::NestedThreadingSupported())
       {
         // threading depends on BLAS
         BlasThreadingEnv knob(num_threads);
         BLAS::gemm('T', 'N', delay_count, norb, norb, cone, U.data(), norb, Ainv.data(), norb, czero, tempMat.data(),
                    lda_Binv);
         for (int i = 0; i < delay_count; i++)
           tempMat(delay_list[i], i) -= cone;
         BLAS::gemm('N', 'N', norb, delay_count, delay_count, cone, V.data(), norb, Binv.data(), lda_Binv, czero,
                    U.data(), norb);
         BLAS::gemm('N', 'N', norb, norb, delay_count, -cone, U.data(), norb, tempMat.data(), lda_Binv, cone,
                    Ainv.data(), norb);
       }
       else
       {
         // manually threaded version of the above GEMM calls
 #pragma omp parallel
         {
           const int block_size = getAlignedSize<T>((norb + num_threads - 1) / num_threads);
           int num_block        = (norb + block_size - 1) / block_size;
 #pragma omp for
           for (int ix = 0; ix < num_block; ix++)
           {
             int x_offset = ix * block_size;
             BLAS::gemm('T', 'N', delay_count, std::min(norb - x_offset, block_size), norb, cone, U.data(), norb,
                        Ainv[x_offset], norb, czero, tempMat[x_offset], lda_Binv);
           }
 #pragma omp master
           for (int i = 0; i < delay_count; i++)
             tempMat(delay_list[i], i) -= cone;
 #pragma omp for
           for (int iy = 0; iy < num_block; iy++)
           {
             int y_offset = iy * block_size;
             BLAS::gemm('N', 'N', std::min(norb - y_offset, block_size), delay_count, delay_count, cone,
                        V.data() + y_offset, norb, Binv.data(), lda_Binv, czero, U.data() + y_offset, norb);
           }
 #pragma omp for collapse(2) nowait
           for (int iy = 0; iy < num_block; iy++)
             for (int ix = 0; ix < num_block; ix++)
             {
               int x_offset = ix * block_size;
               int y_offset = iy * block_size;
               BLAS::gemm('N', 'N', std::min(norb - y_offset, block_size), std::min(norb - x_offset, block_size),
                          delay_count, -cone, U.data() + y_offset, norb, tempMat[x_offset], lda_Binv, cone,
                          Ainv[x_offset] + y_offset, norb);
             }
         }
       }
     }
     delay_count = 0;
   }