Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 4 commits into from
Closed

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

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5ce8a2e
add a pre-condition to filter bad condition numbers during mixing
WHUweiqingzhou Jul 31, 2024
dcd9c08
use dsysv if non-sv is found, while use pseudo-inverse if sv is found
WHUweiqingzhou Aug 1, 2024
26af7d9
update the refdate in 201_NO_15_f_pseudopots
WHUweiqingzhou Aug 1, 2024
46eaadc
add a input mixing_eps for test
WHUweiqingzhou Aug 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions source/module_base/lapack_connector.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,6 +174,10 @@ extern "C"

// solve a tridiagonal linear system
void dgtsv_(int* N, int* NRHS, double* DL, double* D, double* DU, double* B, int* LDB, int* INFO);

// solve Ax = b
void dsysv_(const char* uplo, const int* n, const int* m, double * a, const int* lda,
int *ipiv, double * b, const int* ldb, double *work, const int* lwork ,int *info);
}

#ifdef GATHER_INFO
Expand Down
150 changes: 118 additions & 32 deletions source/module_base/module_mixing/broyden_mixing.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,7 @@ void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function<double
FPTYPE* FP_F = static_cast<FPTYPE*>(F);
if (ndim_cal_dF > 0)
{
ModuleBase::matrix beta_tmp(ndim_cal_dF, ndim_cal_dF);
ModuleBase::matrix beta_tmp(ndim_cal_dF, ndim_cal_dF); // store <dF|dF>
// beta(i, j) = <dF_i, dF_j>
for (int i = 0; i < ndim_cal_dF; ++i)
{
Expand All @@ -137,50 +137,120 @@ void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function<double
}
}
}
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.");

// diagonalization for Singularly Valuable Decomposition (SVD)
ModuleBase::matrix beta_diag(ndim_cal_dF, ndim_cal_dF); // for diagonalization
beta_diag = beta_tmp;
double* val = new double[ndim_cal_dF]; // store eigenvalues

diag(ndim_cal_dF, beta_diag.c, val); // diagonalize beta_tmp

double eps = 1e-12 * val[ndim_cal_dF - 1]; // threshold for filtering eigenvalues
int sv_num = 0; // number of singularly valuables
for (int i = 0; i < ndim_cal_dF; ++i)
{
for (int j = i + 1; j < ndim_cal_dF; ++j)
if (val[i] < eps)
{
beta_tmp(i, j) = beta_tmp(j, i);
sv_num++; // find a sv
}
}
for (int i = 0; i < ndim_cal_dF; ++i)

std::vector<double> gamma(ndim_cal_dF); // store gamma coeficients for Broyden mixing

if (sv_num > 0) // pesudo inverse is sv is found
{
FPTYPE* dFi = FP_dF + i * length;
work[i] = inner_product(dFi, FP_F);
// output eigenvalues, will be removed later
std::cout << "eigenvalues(sv): " << std::endl;
for (int i = 0; i < sv_num; ++i)
{
std::cout << val[i] << " ";
}
std::cout << std::endl;
// pesudo inverse eigenvalue matrix without sv
ModuleBase::matrix inverse_eigenval_wosv(ndim_cal_dF-sv_num, ndim_cal_dF-sv_num);
inverse_eigenval_wosv.zero_out();
for (int i = 0; i < ndim_cal_dF - sv_num; ++i)
{
inverse_eigenval_wosv(i, i) = 1.0 / val[i + sv_num];
}
//
ModuleBase::matrix inverse_tmp_wosv(ndim_cal_dF, ndim_cal_dF - sv_num); // store tmp matrix U * D^-1
ModuleBase::matrix inverse_beta_wosv(ndim_cal_dF, ndim_cal_dF); // store final matrix U * D^-1 * U^T
inverse_tmp_wosv.zero_out();
inverse_beta_wosv.zero_out();
// 1: inverse_beta = U * D^-1
const double alpha = 1.0;
const double beta = 0.0;
double* a_ptr = beta_diag.c + sv_num * ndim_cal_dF; // U
double* b_ptr = inverse_eigenval_wosv.c; // D^-1
double* c_ptr = inverse_tmp_wosv.c; // U * D^-1
const int ndim_cal_dF_wosv = ndim_cal_dF - sv_num;
dgemm_("N", "N", &ndim_cal_dF, &ndim_cal_dF_wosv, &ndim_cal_dF_wosv, &alpha,
a_ptr, &ndim_cal_dF,
b_ptr, &ndim_cal_dF_wosv, &beta,
c_ptr, &ndim_cal_dF);
// 2: (U * D^-1) * U^T
a_ptr = inverse_tmp_wosv.c; // U * D^-1
b_ptr = beta_diag.c + sv_num * ndim_cal_dF; // U
c_ptr = inverse_beta_wosv.c; // U * D^-1 * U^T
dgemm_("N", "T", &ndim_cal_dF, &ndim_cal_dF, &ndim_cal_dF_wosv, &alpha,
a_ptr, &ndim_cal_dF,
b_ptr, &ndim_cal_dF, &beta,
c_ptr, &ndim_cal_dF);

// get <dF|F>
double* dFF_vector = new double[ndim_cal_dF]; // <dF|F>
for (int i = 0; i < ndim_cal_dF; ++i)
{
FPTYPE* dFi = FP_dF + i * length;
dFF_vector[i] = inner_product(dFi, FP_F);
}
// gamma[i] = \sum_j beta_tmp(i,j) * dFF[j]
container::BlasConnector::gemv('N',
ndim_cal_dF,
ndim_cal_dF,
1.0,
inverse_beta_wosv.c,
ndim_cal_dF,
dFF_vector,
1,
0.0,
gamma.data(),
1);
}
else if (sv_num == 0)
{
double* work = new double[ndim_cal_dF]; // workspace
int* iwork = new int[ndim_cal_dF]; // ipiv
char uu = 'U';
int info;
int m = 1;
// gamma means the coeficients for mixing
// but now gamma store <dFi|Fm> since the gamma is input/output in DSYSV
for (int i = 0; i < ndim_cal_dF; ++i)
{
FPTYPE* dFi = FP_dF + i * length;
gamma[i] = inner_product(dFi, FP_F);
}
// solve aG = c
dsysv_(&uu, &ndim_cal_dF, &m, beta_tmp.c, &ndim_cal_dF, iwork, gamma.data(), &ndim_cal_dF, work, &ndim_cal_dF, &info);
if (info != 0)
ModuleBase::WARNING_QUIT("Charge_Mixing", "Error when DSYSV.");
// after solving, gamma store the coeficients for mixing
delete[] work;
delete[] iwork;
}
// 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);
else
{
ModuleBase::WARNING_QUIT("Charge_Mixing", "sv_num < 0");
}
// get the coeficients for mixing from gamma
coef[mdata.start] = 1 + gamma[dFindex_move(0)];
for (int i = 1; i < ndim_cal_dF; ++i)
{
coef[mdata.index_move(-i)] = gamma[dFindex_move(-i)] - gamma[dFindex_move(-i + 1)];
}
coef[mdata.index_move(-ndim_cal_dF)] = -gamma[dFindex_move(-ndim_cal_dF + 1)];

delete[] work;
delete[] iwork;
}
else
{
Expand All @@ -197,4 +267,20 @@ void Broyden_Mixing::tem_cal_coef(const Mixing_Data& mdata, std::function<double
}
ModuleBase::timer::tick("Charge", "Broyden_mixing");
};

void Broyden_Mixing::diag(int n, double* A, double* val) {
// diagonalize a real symmetric matrix A
// A: input matrix (n x n); overwritten with eigenvectors
// val: eigenvalues
// work space query
int lwork = -1;
int info = 0;
std::vector<double> work(1);
dsyev_("V", "U", &n, A, &n, val, work.data(), &lwork, &info);
lwork = work[0];
work.resize(lwork);
// actual computation
dsyev_("V", "U", &n, A, &n, val, work.data(), &lwork, &info);
}

} // namespace Base_Mixing
2 changes: 2 additions & 0 deletions source/module_base/module_mixing/broyden_mixing.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,6 +129,8 @@ class Broyden_Mixing : public Mixing
*/
template <class FPTYPE>
void tem_cal_coef(const Mixing_Data& mdata, std::function<double(FPTYPE*, FPTYPE*)> inner_product);
// to get eigenvalue and eigenvector
void diag(int n, double* A, double* val);

private:
// F = data_out - data_in
Expand Down
Loading
Loading