Skip to content

Conversation

@maki49
Copy link
Collaborator

@maki49 maki49 commented Nov 14, 2024

fix #5376

"precondition_funcs.h" is designed as a function library with following function types:

  1. single-value transformer: none / qe-pw hard-coded formula)
  2. vector-vector operator: div(prevec) / div(trans(prevec-eig)
  3. (In the future, at one's need) Matrix-vector operator

Two ways to use this library:

  • Construct your own lambda / funciton with the library functions (or not), and pass it directly into the eigensolver (e.g. in module_lr)
  • Construct a PreOP object to preallocate the device memory at its construction, then use .get() interface to get the function to be called by the eigensolver (e.g. in hsolver_pw)

Now this library is applied to "dav" and "dav_subspace". The precondition method for "cg" and "bpcg" seems more complicated, so I will not change them before further investigation.

@maki49 maki49 changed the title Refactor: extract preconditioner as functions Refactor: Preconditioner function library Nov 14, 2024
@maki49 maki49 requested review from Cstandardlib, a1henu, denghuilu and haozhihan and removed request for Cstandardlib November 14, 2024 20:20
@maki49 maki49 force-pushed the precondition branch 4 times, most recently from d88022c to 4cc41f5 Compare November 14, 2024 21:00
public:

DiagoDavid(const Real* precondition_in,
using PreFunc = fvec::Div<T>;
Copy link
Collaborator

@Cstandardlib Cstandardlib Nov 15, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is it possible to give PrecondFunc the interface like HPsi or SPsi?

// in precondition lib
using Div = std::function<void(T*, const size_t&, const size_t&)>;

Maybe we can give it the form

using HPsiFunc = std::function<void(T*, T*, const int, const int)>;
using SPsiFunc = std::function<void(T*, T*, const int, const int)>;
using PrecndFunc = std::function<void(T*, T*, const int, const int)>;
// maybe it's better to change int to size_t, I'll do it on current HPsiFunc and SPsiFunc

and pass the preconditioners to this universal function wrapper.

Sometimes the source and destination of preconditioning may be different. For example from a temporary block residual r to a block of total search space X[:, nband:nband+notconv]. The results will be directly stored to a new array in such case.

Copy link
Collaborator Author

@maki49 maki49 Nov 15, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good point, see the latest change.

And seems on HPsiFunc and SPsiFunc can "const" also be added like this:

using HPsiFunc = std::function<void(T*const, const T*const, const int, const int)>;

@kirk0830 kirk0830 removed the request for review from jieli-matrix November 15, 2024 10:53
@mohanchen mohanchen added the Feature Discussed The features will be discussed first but will not be implemented soon label Nov 17, 2024
std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
};

hsolver::PreOP<std::complex<double>, base_device::DEVICE_CPU, hsolver::fvec::DivTransMinusEigKernel<std::complex<double>, base_device::DEVICE_CPU>>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The name "DivTransMinusEigKernel" is a disaster.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry but do you have any better idea?

Copy link
Collaborator

@kirk0830 kirk0830 Nov 17, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The function should be named with following considerations: precisely describe what it does. In present case, you mix the mathematical operation with some rare numerical trick description like "kernel". So there will be quite few people can know what this function does from its name. You should either name the function purely with names mathematical operations, or use human natural language describe the action that what quantity is calculated, additionally from what quantity.

syncmem_op()(this->ctx, this->ctx, spsi_out, psi_in, static_cast<size_t>(nbands * nrow));
};

hsolver::PreOP<std::complex<double>> pre_op(precond_vec);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

其它地方命名都是下划线命名法,PreOP这个名字起得能不能有点默契!

Copy link
Collaborator Author

@maki49 maki49 Nov 17, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

https://github.com/MCresearch/abacus-user-guide/blob/master/develop-C++.md
根据这个文档,变量名和命名空间用下划线命名法,类型名用的是每个单词首字母大写的命名法。
image

@mohanchen
Copy link
Collaborator

如果编程风格不改到我满意,我拒绝所有PR。

@mohanchen mohanchen closed this Nov 17, 2024
@maki49 maki49 reopened this Nov 17, 2024
@mohanchen mohanchen closed this Nov 17, 2024
@kirk0830
Copy link
Collaborator

I will carefully review this PR tomorrow.

Copy link
Collaborator

@kirk0830 kirk0830 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have reviewed only parts of this PR, there are still many lines I have not enough time to review, I will continue later.

namespace fval
{
template <typename T> T none(const T& x) { return x; }
template <typename T> T qe_pw(const T& x) { return 0.5 * (1.0 + x + sqrt(1 + (x - 1.0) * (x - 1.0))); }
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am a little bit confused about the name of this function, also I am curious if it can be accelerated automatically with SIMD, can you answer my question?

/// type 1: directly divide each vector by the precondition vector
///---------------------------------------------------------------------------------------------
template <typename T, typename Device = base_device::DEVICE_CPU>
void div_prevec(T* const dst, const T* const src, const size_t& dim, const size_t& nvec,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is not a good idea to use a term "prevec". Although it is scarcely seen the abbreviation of it directly or easily, at least you can know the word "condition" has abbr. "cond", so I think "precond" is okay. Then I will question "what is precondition vector"? is it a formal name of some numerical/technical operation? I guess not. Then I will ask you what is the meaning of "vector" in your code? Is it about the data structure? the way to store data which indicates the memory is continuous, or a mathematical description? Only the last one is acceptable.

for (size_t m = 0; m < nvec; m++)
{
const size_t offset = m * dim;
vector_div_vector_op<T, Device>()(ctx, dim, dst + offset, src + offset, pre);
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

actually the name vector_div_vector_op is not a good name, so please do not get any idea about naming functions from it. But you are not expected to improve its name in this PR

std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
};

hsolver::PreOP<std::complex<double>, base_device::DEVICE_CPU, hsolver::fvec::DivTransMinusEigKernel<std::complex<double>, base_device::DEVICE_CPU>>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I cannot see any significance of the necessarity of instantiating the hsolver::PreOP. The name of PreOP is also confusing quite much. From the rest of your PR, I can understand it is about preconditioning, but unfortunately you did not name them with much care and enough considerations.
You should spend at least 60% of your time on coding to carefully consider the name of variable, function, class, namespace, and the design of interface.


hsolver::PreOP<std::complex<double>, base_device::DEVICE_CPU, hsolver::fvec::DivTransMinusEigKernel<std::complex<double>, base_device::DEVICE_CPU>>
pre_op(precond_vec, hsolver::fvec::div_trans_prevec_minus_eigen<std::complex<double>>, hsolver::fval::qe_pw<double>);
obj = std::make_unique<hsolver::Diago_DavSubspace<std::complex<double>, base_device::DEVICE_CPU>>(
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it is the first time that I find there is an instance even named obj. Informativeless.

std::copy(hpsi_ptr, hpsi_ptr + nvec * ld_psi, hpsi_out);
};

hsolver::PreOP<std::complex<double>, base_device::DEVICE_CPU, hsolver::fvec::DivTransMinusEigKernel<std::complex<double>, base_device::DEVICE_CPU>>
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

presently the lines 133, 134 are quite confusing and it is quite difficult for others to know what the operation is. It wastes other's time. You should always consider a clearer way to organize your code like:

// instantiate hsolver::PreOP
hsolver::PreOP<std::complex<double>, 
               base_device::DEVICE_CPU, 
               hsolver::fvec::DivTransMinusEigKernel<
                   std::complex<double>, 
                   base_device::DEVICE_CPU>
               > pre_op(precond_vec, 
                        hsolver::fvec::div_trans_prevec_minus_eigen<std::complex<double>>, 
                        hsolver::fval::qe_pw<double>);

however, do you think naming the instance of PreOP as pre_op is a good idea?


template <typename T, typename Device>
Diago_DavSubspace<T, Device>::Diago_DavSubspace(const std::vector<Real>& precondition_in,
Diago_DavSubspace<T, Device>::Diago_DavSubspace(PreFunc&& precondition_in,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

could you explain the necessarity of using rvalue instead of lvalue?

Comment on lines +579 to 582
DiagoDavid<T, Device> david(pre_op.get(),
nband,
dim,
PARAM.inp.pw_diag_ndim,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

besides the PreOP, now I realize the DiagoDavid class has the similar problem. It is confusing to instantiate the class DiagoDavid then name as david. A factory function will be better and more consistent, but will not be work of this PR:

std::unique_ptr<Diago> dav = make_diago("dav", /* other args */); // or use object directly without pointer

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Feature Discussed The features will be discussed first but will not be implemented soon

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Numerical instability of Davidson iterative eigensolver in solving LR-TDDFT excited states

4 participants