Feature: Use pesudo inverse to avoid singularly values and use `dsysv` to replace direct inverse in Broyden mixing

[WHUweiqingzhou]

Fix #4802.

List of Changes

1. replace replace dsytrf & dsytri & dgemv by dsysv in Broyden mixing. Plz note that The pass of PR #4842 CI tests show that dsysv is OK to replace dsytrf & dsytri & dgemv.
2. implement a pesudo-inverse in Broyden mixing.
3. set threshold of singularly value as 1e-12, waiting for more tests.
4. update the refdata in 201_NO_15_f_pseudopots.

update the refdata in 201_NO_15_f_pseudopots since this example triggers the pseudo-inverse at 7-th iteration, namely:

* * * * * *
 << Start SCF iteration.
 ITER       ETOT/eV          EDIFF/eV         DRHO     TIME/s
 GE1     -3.43887808e+03   0.00000000e+00   6.2153e-02   0.36
 GE2     -3.43901104e+03  -1.32966465e-01   2.4813e-02   0.06
 GE3     -3.43901115e+03  -1.06514536e-04   5.5111e-03   0.06
 GE4     -3.43900800e+03   3.14573854e-03   4.1309e-04   0.06
 GE5     -3.43900789e+03   1.09726259e-04   1.1831e-04   0.06
 GE6     -3.43900793e+03  -3.70139444e-05   1.2857e-05   0.06
 GE7     -3.43900793e+03  -1.40156736e-06   3.0691e-07   0.06
eigenvalues(sv): 
5.74409e-16 
 GE8     -3.43900793e+03   1.19930268e-08   1.2789e-08   0.06
 GE9     -3.43900793e+03  -4.62296497e-09   5.1121e-10   0.06
 >> Leave SCF iteration.
 * * * * * *

leading to a tiny difference in final DMR.

To prove the implementation of pesudo-inverse in this PR is correct

I choose a simple example at abacus-develop/examples/scf/lcao_Si2, and set INPUT as

INPUT_PARAMETERS
#Parameters     (General)
pseudo_dir              ../../../tests/PP_ORB
orbital_dir             ../../../tests/PP_ORB
#Parameters (Accuracy)
ecutwfc                 50
scf_nmax                100
scf_thr                 1e-6
basis_type              lcao
mixing_beta 0.4

During SCF, I output beta_tmp and inverse_beta_wosv at each iteration. In this calculation, I set thershold for singularly valuable as 1e-8. At 4-th iteration, the beta_matrix and inverse_beta_wosv is

beta_tmp: 
0.30991644980730138 0.45375775064234292 -0.0063860954924527128 0.0025808819908687709 
0.45375775064234292 0.66535882042578864 -0.0076176091637253866 0.0037556401927858148 
-0.0063860954924527128 -0.0076176091637253866 0.0031428520369101002 -9.3922204323843109e-05 
0.0025808819908687709 0.0037556401927858148 -9.3922204323843109e-05 2.2285714311579603e-05 

inverse_beta_wosv: 
1329781.8370982155 -907612.25488413975 538726.20938949555 1222909.9542547921 
-907612.25488413963 619596.45690384216 -368073.12804228161 -857566.09836292267 
538726.20938949543 -368073.12804228161 219744.12776514384 565398.73159492889 
1222909.9542547921 -857566.09836292267 565398.73159492889 5323003.029010226 

In this matrix, there is no singularly valuable. I use python to inverse the matrix by:

import numpy as np
mat = np.loadtxt('m4.dat')
# direct inverse (trf+tri)
inv1 = np.linalg.inv(mat)
print("fully inverse beta")
print(inv1)

and get the inverse matrix:

fully inverse beta
[[1329781.83747472 -907612.25514676  538726.20955929 1222909.95563534]
 [-907612.25514676  619596.45708704 -368073.12816078 -857566.09932965]
 [ 538726.2095593  -368073.12816078  219744.12784191  565398.73222901]
 [1222909.95563613 -857566.09933019  565398.73222932 5323003.03476567]]

You can see that the result is the same as inverse_beta_wosv in the calculation. So the implementation of pseudo-inverse can reproduce the full inverse if there is no singularly valuable in the matrix.

At 7-th iteration, the beta_matrix and inverse_beta_wosv is

beta_tmp: 
0.30991644980730138 0.45375775064234292 -0.0063860954924527128 0.0025808819908687709 0.0017322960687262837 -4.8537701873346579e-05 -1.8086621143195072e-05 
0.45375775064234292 0.66535882042578864 -0.0076176091637253866 0.0037556401927858148 0.0025223399442106708 -7.0396554149326749e-05 -2.6318512057081064e-05 
-0.0063860954924527128 -0.0076176091637253866 0.0031428520369101002 -9.3922204323843109e-05 -6.191416912745716e-05 2.2149545751773874e-06 6.753640391042234e-07 
0.0025808819908687709 0.0037556401927858148 -9.3922204323843109e-05 2.2285714311579603e-05 1.5632503314831404e-05 -4.3451231126606036e-07 -1.5795089413856005e-07 
0.0017322960687262837 0.0025223399442106708 -6.191416912745716e-05 1.5632503314831404e-05 1.2915074979472969e-05 -3.301103803883504e-07 -1.1471873560213929e-07 
-4.8537701873346579e-05 -7.0396554149326749e-05 2.2149545751773874e-06 -4.3451231126606036e-07 -3.301103803883504e-07 9.0275512101134368e-09 3.2312200552144175e-09 
-1.8086621143195072e-05 -2.6318512057081064e-05 6.753640391042234e-07 -1.5795089413856005e-07 -1.1471873560213929e-07 3.2312200552144175e-09 1.1881888757367936e-09 
eigenvalues(sv): 
5.9739721575589226e-15 3.0988131226578748e-12 1.6513760321049919e-10 
inverse_beta_wosv: 
1308970.4144005312 -888307.38733344967 514581.87891458248 57982.548945154551 313776.90349846962 18417.665282821468 13250.160174952882 
-888307.38733344967 602841.8572745705 -349230.27878850576 -39764.622848515246 -214427.28045297775 -12479.048973192026 -8989.0780916783679 
514581.87891458254 -349230.27878850576 202687.97516340908 24028.235162993187 127789.90231289201 7182.9336764697991 5201.151515010587 
57982.548945154551 -39764.622848515246 24028.235162993187 24941.633835385637 93992.558019673539 -246.06824438769928 431.8048211896371 
313776.90349846962 -214427.28045297775 127789.90231289201 93992.558019673554 361942.13980911387 613.53845310958309 2620.9416704590385 
18417.665282821468 -12479.048973192026 7182.9336764697991 -246.06824438769928 613.53845310958366 309.54833017925165 193.80012808272735 
13250.160174952882 -8989.0780916783679 5201.151515010587 431.80482118963704 2620.9416704590385 193.80012808272738 135.20347640814734 

there are 3 singularly valuables in the matrix. I use python to inverse the matrix by:

import numpy as np
thr = 1e-8
mat = np.loadtxt('m7.dat')
# pseudo-inverse
pinv2 = np.linalg.pinv(mat, rcond=thr)
print("pseudo inverse beta")
print(pinv2)

and get the inverse matrix:

pseudo inverse beta
[[ 1.30897041e+06 -8.88307387e+05  5.14581879e+05  5.79825489e+04
   3.13776904e+05  1.84176650e+04  1.32501602e+04]
 [-8.88307387e+05  6.02841857e+05 -3.49230279e+05 -3.97646228e+04
  -2.14427280e+05 -1.24790488e+04 -8.98907809e+03]
 [ 5.14581879e+05 -3.49230279e+05  2.02687975e+05  2.40282352e+04
   1.27789902e+05  7.18293356e+03  5.20115151e+03]
 [ 5.79825489e+04 -3.97646228e+04  2.40282352e+04  2.49416338e+04
   9.39925580e+04 -2.46068261e+02  4.31804821e+02]
 [ 3.13776904e+05 -2.14427280e+05  1.27789902e+05  9.39925580e+04
   3.61942140e+05  6.13538369e+02  2.62094167e+03]
 [ 1.84176650e+04 -1.24790488e+04  7.18293356e+03 -2.46068261e+02
   6.13538369e+02  3.09548322e+02  1.93800125e+02]
 [ 1.32501602e+04 -8.98907809e+03  5.20115151e+03  4.31804821e+02
   3.62094167e+03  1.93800125e+02  1.35203476e+02]]

You can see that the result is the same as inverse_beta_wosv in the calculation.