Skip to content
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.

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

Should solve integer overflows in SAPT code #619

Merged
merged 2 commits into from Feb 17, 2017
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion psi4/src/psi4/libsapt_solver/amplitudes.cc
Expand Up @@ -620,7 +620,7 @@ void SAPT2::ijkl_to_ikjl(double *tARAR, int ilength, int jlength, int klength,

for(int i=0; i<ilength; i++) {
for(int l=0; l<llength; l++) {
long int ijkl = i*jlength*klength*llength + l;
long int ijkl = (long int)i*jlength*klength*llength + l;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I might remember wrong, but don't you need to also cast jlength, klength and llength, since the order of evaluation isn't guaranteed...? What this is doing is typecasting i to long int, but jlengthklengthllength is still int and might overflow.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nevermind, the compiler should promote all of the operands to long int before the multiplication takes place.

C_DCOPY(jlength*klength,&(tARAR[ijkl]),llength,X,1);
for(int j=0; j<jlength; j++) {
for(int k=0; k<klength; k++) {
Expand Down
34 changes: 17 additions & 17 deletions psi4/src/psi4/libsapt_solver/disp2ccd.cc
Expand Up @@ -240,7 +240,7 @@ void SAPT2p::r_ccd_prep(const char *TARBS, const char *ARBS, const char *CA_RBS,

double **tARBS = block_matrix(occA*virA,occB*virB); //!

C_DCOPY(occA*virA*occB*virB,&(vARBS[0][0]),1,&(tARBS[0][0]),1);
C_DCOPY((size_t)occA*virA*occB*virB,&(vARBS[0][0]),1,&(tARBS[0][0]),1);
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same thing here.


double **thARAR = block_matrix(occA*virA,occA*virA); //!

Expand Down Expand Up @@ -394,7 +394,7 @@ double SAPT2p::r_ccd_energy(const char *TARBS, const char *ARBS, int occA, int v
psio_->read_entry(PSIF_SAPT_CCD,TARBS,(char *) &(tARBS[0][0]),
occA*virA*occB*virB*(ULI) sizeof(double));

double energy = C_DDOT(occA*virA*occB*virB,&(vARBS[0][0]),1,
double energy = C_DDOT((size_t)occA*virA*occB*virB,&(vARBS[0][0]),1,
&(tARBS[0][0]),1);

free_block(vARBS); //!
Expand All @@ -419,7 +419,7 @@ double SAPT2p::r_ccd_iterate(const char *TARBS, const char *TARBSerr, const char
int iter = 1;
double E_old=0.0, E_new=0.0, RMS=0.0;

SAPTDIIS diis(PSIF_SAPT_CCD,TARBS,TARBSerr,occA*virA*occB*virB,
SAPTDIIS diis(PSIF_SAPT_CCD,TARBS,TARBSerr,(size_t)occA*virA*occB*virB,
max_ccd_vecs_,psio_);

do {
Expand Down Expand Up @@ -569,10 +569,10 @@ double SAPT2p::r_ccd_amplitudes(const char *TARBS, const char *TARBSerr, const c
psio_->write_entry(PSIF_SAPT_CCD,TARBS,(char *) &(t2ARBS[0][0]),
occA*virA*occB*virB*(ULI) sizeof(double));

C_DAXPY(occA*virA*occB*virB,-1.0,tARBS[0],1,t2ARBS[0],1);
C_DAXPY((size_t)occA*virA*occB*virB,-1.0,tARBS[0],1,t2ARBS[0],1);

double RMS = C_DDOT(occA*virA*occB*virB,t2ARBS[0],1,t2ARBS[0],1);
RMS /= (double) (occA*virA*occB*virB);
double RMS = C_DDOT((size_t)occA*virA*occB*virB,t2ARBS[0],1,t2ARBS[0],1);
RMS /= (double) ((size_t)occA*virA*occB*virB);

psio_->write_entry(PSIF_SAPT_CCD,TARBSerr,(char *) &(t2ARBS[0][0]),
occA*virA*occB*virB*(ULI) sizeof(double));
Expand Down Expand Up @@ -712,7 +712,7 @@ double SAPT2p::s_ccd_iterate(const char *SARAR, const char *SARARerr, const char
int iter = 1;
double E_old=0.0, E_new=0.0, RMS=0.0;

SAPTDIIS diis(PSIF_SAPT_CCD,SARAR,SARARerr,occA*virA*occA*virA,
SAPTDIIS diis(PSIF_SAPT_CCD,SARAR,SARARerr,(size_t)occA*virA*occA*virA,
max_ccd_vecs_,psio_);

do {
Expand Down Expand Up @@ -1077,9 +1077,9 @@ double SAPT2p::s_ccd_amplitudes(const char *SARAR, const char *SARARerr, const c
psio_->write_entry(PSIF_SAPT_CCD,SARAR,(char *) &(s2ARAR[0][0]),
occA*virA*occA*virA*(ULI) sizeof(double));

C_DAXPY(occA*virA*occA*virA,-1.0,sARAR[0],1,s2ARAR[0],1);
double RMS = C_DDOT(occA*virA*occA*virA,s2ARAR[0],1,s2ARAR[0],1);
RMS /= (double) (occA*virA*occA*virA);
C_DAXPY((size_t)occA*virA*occA*virA,-1.0,sARAR[0],1,s2ARAR[0],1);
double RMS = C_DDOT((size_t)occA*virA*occA*virA,s2ARAR[0],1,s2ARAR[0],1);
RMS /= (double) ((size_t)occA*virA*occA*virA);

psio_->write_entry(PSIF_SAPT_CCD,SARARerr,(char *) &(s2ARAR[0][0]),
occA*virA*occA*virA*(ULI) sizeof(double));
Expand Down Expand Up @@ -1661,7 +1661,7 @@ void SAPT2p::ccd_iterate(const char *TARAR, const char *TARARerr, const char *Th
int iter = 1;
double E_old=0.0, E_new=0.0, RMS=0.0;

SAPTDIIS diis(PSIF_SAPT_CCD,TARAR,TARARerr,occA*virA*occA*virA,
SAPTDIIS diis(PSIF_SAPT_CCD,TARAR,TARARerr,(size_t)occA*virA*occA*virA,
max_ccd_vecs_,psio_);

do {
Expand Down Expand Up @@ -1710,7 +1710,7 @@ double SAPT2p::ccd_energy(const char *TARAR, const char *GARAR, int occA, int vi
psio_->read_entry(PSIF_SAPT_CCD,TARAR,(char *) &(tARAR[0][0]),
occA*virA*occA*virA*(ULI) sizeof(double));

double energy = C_DDOT(occA*virA*occA*virA,&(gARAR[0][0]),1,
double energy = C_DDOT((size_t)occA*virA*occA*virA,&(gARAR[0][0]),1,
&(tARAR[0][0]),1);

free_block(gARAR); //!
Expand All @@ -1730,7 +1730,7 @@ double SAPT2p::ccd_amplitudes(const char *TARAR, const char *TARARerr, const cha
psio_->read_entry(PSIF_SAPT_CCD,ARAR,(char *) &(t2ARAR[0][0]),
occA*virA*occA*virA*(ULI) sizeof(double));

C_DSCAL(occA*virA*occA*virA,0.5,&(t2ARAR[0][0]),1);
C_DSCAL((size_t)occA*virA*occA*virA,0.5,&(t2ARAR[0][0]),1);

double **gARRA = block_matrix(occA*virA,occA*virA); //!
double **tARAR = block_matrix(occA*virA,occA*virA); //!
Expand Down Expand Up @@ -1993,10 +1993,10 @@ double SAPT2p::ccd_amplitudes(const char *TARAR, const char *TARARerr, const cha

free_block(th2ARAR); //!

C_DAXPY(occA*virA*occA*virA,-1.0,t2ARAR[0],1,tARAR[0],1);
C_DAXPY((size_t)occA*virA*occA*virA,-1.0,t2ARAR[0],1,tARAR[0],1);

double RMS = C_DDOT(occA*virA*occA*virA,tARAR[0],1,tARAR[0],1);
RMS /= (double) (occA*virA*occA*virA);
double RMS = C_DDOT((size_t)occA*virA*occA*virA,tARAR[0],1,tARAR[0],1);
RMS /= (double) ((size_t)occA*virA*occA*virA);

psio_->write_entry(PSIF_SAPT_CCD,TARARerr,(char *) &(tARAR[0][0]),
occA*virA*occA*virA*(ULI) sizeof(double));
Expand Down Expand Up @@ -2029,7 +2029,7 @@ void SAPT2p::write_IJKL(double **A, int filenum, const char *label, int length_I
free_block(A);
}

SAPTDIIS::SAPTDIIS(int ampfile, const char *amplabel, const char *errlabel, int length,
SAPTDIIS::SAPTDIIS(int ampfile, const char *amplabel, const char *errlabel, size_t length,
int maxvec, std::shared_ptr<PSIO> psio) : psio_(psio), vec_label_(amplabel), err_label_(errlabel)
{
diis_file_ = 56;
Expand Down
30 changes: 15 additions & 15 deletions psi4/src/psi4/libsapt_solver/sapt.h
Expand Up @@ -64,19 +64,19 @@ class SAPT : public Wavefunction {
std::shared_ptr<BasisSet> elstbasis_;
std::shared_ptr<BasisSet> zero_;

int nsoA_;
int nmoA_;
int nsoB_;
int nmoB_;
int ndf_;
int noccA_;
int foccA_;
int aoccA_;
int noccB_;
int foccB_;
int aoccB_;
int nvirA_;
int nvirB_;
size_t nsoA_;
size_t nmoA_;
size_t nsoB_;
size_t nmoB_;
size_t ndf_;
size_t noccA_;
size_t foccA_;
size_t aoccA_;
size_t noccB_;
size_t foccB_;
size_t aoccB_;
size_t nvirA_;
size_t nvirB_;
int NA_;
int NB_;
int natomsA_;
Expand Down Expand Up @@ -110,7 +110,7 @@ class SAPT : public Wavefunction {

std::shared_ptr<SAPTDenominator> denom_;

int nvec_;
size_t nvec_;

double **dAR_;
double **dBS_;
Expand All @@ -130,7 +130,7 @@ class CPHFDIIS {

private:
int max_diis_vecs_;
int vec_length_;
size_t vec_length_;

int curr_vec_;
int num_vecs_;
Expand Down
14 changes: 7 additions & 7 deletions psi4/src/psi4/libsapt_solver/sapt0.h
Expand Up @@ -179,11 +179,11 @@ struct SAPTDFInts {
bool dress_disk_;
bool active_;

int i_length_;
int j_length_;
int ij_length_;
int i_start_;
int j_start_;
size_t i_length_;
size_t j_length_;
size_t ij_length_;
size_t i_start_;
size_t j_start_;

double **B_p_;
double **B_d_;
Expand All @@ -206,10 +206,10 @@ struct SAPTDFInts {

struct Iterator {

int num_blocks;
size_t num_blocks;
int* block_size;

int curr_block;
size_t curr_block;
long int curr_size;

~Iterator() { free(block_size); };
Expand Down
6 changes: 4 additions & 2 deletions psi4/src/psi4/libsapt_solver/sapt2.h
Expand Up @@ -38,15 +38,17 @@ class SAPT2 : public SAPT {
virtual void print_results();

protected:
int no_nvirA_;
int no_nvirB_;
size_t no_nvirA_;
size_t no_nvirB_;

double *no_evalsA_;
double *no_evalsB_;

double **no_CA_;
double **no_CB_;

// These ints should not overflow below 46000 basis functions
// They contain indices up to nso_ * (nso_ + 1) / 2
int *ioff_;
int *index2i_;
int *index2j_;
Expand Down
4 changes: 2 additions & 2 deletions psi4/src/psi4/libsapt_solver/sapt2p.h
Expand Up @@ -171,7 +171,7 @@ class SAPTDIIS {
int max_diis_vecs_;

int diis_file_;
int vec_length_;
size_t vec_length_;

int curr_vec_;
int num_vecs_;
Expand All @@ -183,7 +183,7 @@ class SAPTDIIS {
std::shared_ptr<PSIO> psio_;

public:
SAPTDIIS(int, const char *, const char *, int, int, std::shared_ptr<PSIO>);
SAPTDIIS(int, const char *, const char *, size_t, int, std::shared_ptr<PSIO>);
~SAPTDIIS();

void store_vectors();
Expand Down