Skip to content

Fix: memory access out of bound #3205

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

Merged
dyzheng merged 1 commit into from
Nov 15, 2023
Merged

Fix: memory access out of bound #3205

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
43 changes: 22 additions & 21 deletions source/module_base/math_sphbes.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -369,30 +369,31 @@ void Sphbes::Spherical_Bessel_Roots
Sphbes::Spherical_Bessel(msh, r, 1, l, jl);

int n=0;
for (int i=0; i<msh && n<num; i++)
for (int i=0; i<msh-1 && n<num; i++)
{
if (jl[i]*jl[i+1] < 0.0)
{
double y_1 = jl[i];
double y_2 = jl[i+1];
double x_1 = r[i];
double x_2 = r[i+1];
double acc = std::fabs( y_2 - y_1 );
double acc = std::fabs(y_2 - y_1);

const int grid=100;
double *rad = new double[grid];
double *jl_new = new double[grid];
while (acc > epsilon)
{
double *rad = new double[100];
double *jl_new = new double[100];

// if not enough accurate, divide again.
const double delta2 = (x_2 - x_1)/99.0;
for (int j=0;j<100;j++)
const double delta2 = (x_2 - x_1)/(grid-1);
for (int j=0;j<grid;j++)
{
rad[j] = x_1 + j*delta2;
}
Sphbes::Spherical_Bessel(100,rad,1,l,jl_new);
Sphbes::Spherical_Bessel(grid,rad,1,l,jl_new);

int j=0;
for (;j<100;j++)
for (;j<grid-1;j++)
Comment thread
dyzheng marked this conversation as resolved.
{
if (jl_new[j]*jl_new[j+1]<0)break;
}
Expand All @@ -402,9 +403,9 @@ void Sphbes::Spherical_Bessel_Roots
y_1 = jl_new[j];
y_2 = jl_new[j+1];
acc = std::fabs( y_2 - y_1 );
delete[] rad;
delete[] jl_new;
}
delete[] rad;
delete[] jl_new;
eigenvalue[n]=(x_2 + x_1)*0.5/rcut;
n++;
}
Expand Down Expand Up @@ -638,13 +639,13 @@ void Sphbes::dSpherical_Bessel_dx
std::cout << "We temporarily only calculate derivative of l >= 0." << std::endl;
exit(0);
}

double djl0 = 0;
if(l == 1)
{
djl0 = 1.0/3.0;
}

if(l == 0 )
{
for (int ir = 0;ir < msh; ir++)
Expand All @@ -658,7 +659,7 @@ void Sphbes::dSpherical_Bessel_dx
{
djl[ir] = (x1 * std::cos(x1) - std::sin(x1)) / (x1*x1);
}

}
}
else
Expand All @@ -669,7 +670,7 @@ void Sphbes::dSpherical_Bessel_dx
for (int ir = 0;ir < msh; ir++)
{
double x1 = q * r[ir];
if(x1 < 1e-8)
if(x1 < 1e-8)
{
djl[ir] = djl0;
}
Expand All @@ -687,7 +688,7 @@ void Sphbes::dSpherical_Bessel_dx
double Sphbes::_sphbesj_ascending_recurrence(int l, double x) {

// should be used when x > l && l > 0

double invx = 1.0 / x;
double j0 = std::sin(x) * invx;
double j1 = ( j0 - std::cos(x) ) * invx;
Expand Down Expand Up @@ -753,21 +754,21 @@ double Sphbes::sphbesj(const int l, const double x)
case 0:
return std::sin(x) * invx;
case 1:
return ( std::sin(x) * invx - std::cos(x) ) * invx;
return ( std::sin(x) * invx - std::cos(x) ) * invx;
// NOTE: the following explicit expressions are not necessarily faster than ascending recurrence,
// but we keep them just in case we need them in the future.
//case 2:
// return ( (3.0 * invx - x) * std::sin(x) - 3.0 * std::cos(x) ) * (invx * invx);
//case 3:
// return ( std::sin(x) * (15.0 * invx - 6.0 * x) + std::cos(x) * (x * x - 15.0) ) * std::pow(invx, 3);
// return ( std::sin(x) * (15.0 * invx - 6.0 * x) + std::cos(x) * (x * x - 15.0) ) * std::pow(invx, 3);
//case 4:
// return ( std::sin(x) * (std::pow(x,3) - 45.0 * x + 105.0 * invx)
// return ( std::sin(x) * (std::pow(x,3) - 45.0 * x + 105.0 * invx)
// + std::cos(x) * (10.0 * x * x - 105.0) ) * std::pow(invx, 4);
//case 5:
// return ( std::sin(x) * (15.0 * std::pow(x,3) - 420.0 * x + 945.0 * invx)
// return ( std::sin(x) * (15.0 * std::pow(x,3) - 420.0 * x + 945.0 * invx)
// + std::cos(x) * (-std::pow(x, 4) + 105.0 * x * x - 945.0) ) * std::pow(invx, 5);
//case 6:
// return ( std::sin(x) * (-std::pow(x, 5) + 210.0 * std::pow(x, 3) - 4725.0 * x + 10395.0 * invx)
// return ( std::sin(x) * (-std::pow(x, 5) + 210.0 * std::pow(x, 3) - 4725.0 * x + 10395.0 * invx)
// + std::cos(x) * (-21.0 * std::pow(x, 4) + 1260.0 * x * x - 10395.0) ) * std::pow(invx, 6);
default:
return _sphbesj_ascending_recurrence(l, x);
Expand Down