Code Quality: improve numerical stability of `base_mixing` class

[WHUweiqingzhou]

Describe the Code Quality Issue

Thanks to @jinzx10 mentined in #4058

module_base/module_mixing/broyden_mixing.cpp (line number 140-174):

double* work = new double[ndim_cal_dF];
int* iwork = new int[ndim_cal_dF];
char uu = 'U';
int info;
dsytrf_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &ndim_cal_dF, &info);
if (info != 0)
    ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when factorizing beta.");
dsytri_(&uu, &ndim_cal_dF, beta_tmp.c, &ndim_cal_dF, iwork, work, &info);
if (info != 0)
    ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYTRI beta.");
for (int i = 0; i < ndim_cal_dF; ++i)
{
    for (int j = i + 1; j < ndim_cal_dF; ++j)
    {
        beta_tmp(i, j) = beta_tmp(j, i);
    }
}
for (int i = 0; i < ndim_cal_dF; ++i)
{
    FPTYPE* dFi = FP_dF + i * length;
    work[i] = inner_product(dFi, FP_F);
}
// gamma[i] = \sum_j beta_tmp(i,j) * work[j]
std::vector<double> gamma(ndim_cal_dF);
container::BlasConnector::gemv('N',
                               ndim_cal_dF,
                               ndim_cal_dF,
                               1.0,
                               beta_tmp.c,
                               ndim_cal_dF,
                               work,
                               1,
                               0.0,
                               gamma.data(),
                               1);

It seems to me that the above code is solving a linear equation Ax=b by computing inv(A) followed by a multiplication. While this procedure should behave as expected when A is a well-conditioned matrix, it could lead to severe numerical instability when A is poorly conditioned (i.e., error vectors have near linear dependency); dsysv should be preferred in this case (https://netlib.org/lapack/explore-html-3.6.1/d6/d0e/group__double_s_ysolve_gac0d0ad0edaa9d2014264c78874055db1.html).

Additional Context

No response

Task list for Issue attackers (only for developers)

• Identify the specific code file or section with the code quality issue.
• Investigate the issue and determine the root cause.
• Research best practices and potential solutions for the identified issue.
• Refactor the code to improve code quality, following the suggested solution.
• Ensure the refactored code adheres to the project's coding standards.
• Test the refactored code to ensure it functions as expected.
• Update any relevant documentation, if necessary.
• Submit a pull request with the refactored code and a description of the changes made.

[jinzx10]

Just to add some information:

While dsysv should in general performs better than dsytrf+dsytri+dgemm, it still assumes that the problem is not singular. In case the matrix is very close to singular, extra steps might be necessary, and I think there's no unique solution. For example, pyscf would do an eigen-decomposition and filters out tiny eigenvalues (more or less the pseudo-inverse); some quantum chemistry softwares chooses to shrink the mixing subspace until the problem is well-conditioned; using least square solvers might be another option.