unstable bug of "Assertion `psi_norm > 0.0' failed" in Davidson method

[pxlxingliang]

Describe the Bug

In this LiSi system, when using davidson method, ABACUS sometimes throw the error of "Assertion `psi_norm > 0.0' failed". This bug randomly occur when using 16 threads on 16 physical cores' machine.

abacus: /root/abacus-develop/source/module_hsolver/diago_david.cpp:544: void hsolver::DiagoDavid::SchmitOrth(hamilt::Hamilt *, const int &, int, int, const ModuleBase::ComplexMatrix &, std::complex *, std::complex *): Assertion `psi_norm > 0.0' failed.

dav.zip

[pxlxingliang]

This bug is occurred in the ABACUS complied by intel, and ABACUS complied by GNU has no problem.
A detail checking, I find that the function Parallel_Reduce::reduce_complex_double_pool in DiagoDavid::cal_elem sometimes "is bad" for hc matrix. This function should plus the values in each core and than assign it back to each core, while sometimes some elements in hc matrix in some cores are not same.

Interestingly, I try to running this function by each element, and the bug is not occurred. Even I don't know why, but it works, at least not occurred in 100 repeat running on below Al example.
Modify
Parallel_Reduce::reduce_complex_double_pool(hc.c + offset_h, notconv * hc.nr);
to be
for(int i=0;i<notconv * hc.nr;i++) { Parallel_Reduce::reduce_complex_double_pool(hc.c + offset_h + i, 1); }

al.zip

[caic99]

#1221 is a temporary fix; it will cause significant communication overhead. Fixing MKL linking issue should be the best solution.

My profiling result of P000_si16_pw with 32 processes:

• 22.15% hsolver::DiagoDavid::diag_zhegvx
• 21.52% Parallel_Common::bcast_complex_double
• 0.60% zhegvx_

[pxlxingliang]

The davidson method need do cal_elem and then do diag_zhegvx.
The cal_elem is to construct the HC/SC matrix, where the element is calculated by parallel and followed the MPI_REDUCE.
The diag_zhegvx is to dragonize the HC/SC matrix, and get the eigenvalues and eigenvectors.

The old version (before #1221 ) do the ALL TO ALL MPI_REDUCE in cal_elem, and each process do the diag_zhegvx. But the ALL TO ALL MPI_REDUCE has some problems when complied with intel.

#1221 is a different solution that doing the MPI_Reduce to rank0, and then do diag_zhegvx in rank0, and then do bcast_complex_double for eigenvalues and eigenvectors.

So in #1221 , there have an extra cost in hsolver::DiagoDavid::diag_zhegvx and Parallel_Common::bcast_complex_double, but a smaller cost in cal_elem.