#include <scconfig.h>
#include <util/misc/exenv.h>

#undef SCF_CHECK_BOUNDS

#ifdef SCF_CHECK_BOUNDS
#define CHECK(ival,pval,ij,kl,con) check(ival,pval,ij,kl,con)
#else
#define CHECK(ival,pval,ij,kl,con)
#endif

namespace sc {

class LocalCLHFContribution {
  private:
    double * const restrictxx gmat;
    double * const restrictxx pmat;

    double ibound_;
    double pbound_;

  public:
    LocalCLHFContribution(double *g, double *p) : gmat(g), pmat(p) {}
    ~LocalCLHFContribution() {}

    void set_bound(double i, double p) { ibound_ = i; pbound_ = p; }
    void check(double ival, double pval, int ij, int kl, const char *contrib)
        {
          int bad = 0;
          if ( 1.000001 * ibound_ < (ival > 0 ? ival : -ival)) {
              ExEnv::errn() << "BAD INTEGRAL BOUND" << std::endl;
              ExEnv::errn() << " bound = " << ibound_ << std::endl;
              ExEnv::errn() << " value = " << ival << std::endl;
              bad = 1;
            }
          if ( 1.000001 * pbound_ < (pval > 0 ? pval : -pval)) {
              ExEnv::errn() << "BAD DENSITY BOUND" << std::endl;
              ExEnv::errn() << " bound = " << pbound_ << std::endl;
              ExEnv::errn() << " value = " << pval << std::endl;
              bad = 1;
            }
          if (bad) {
              ExEnv::errn() << " ij    = " << ij << std::endl;
              ExEnv::errn() << " kl    = " << kl << std::endl;
              ExEnv::errn() << " cont  = " << contrib << std::endl;
              abort();
            }
        }

    inline void cont1(int ij, int kl, double val) {
      gmat[ij] += val*pmat[kl]; CHECK(val,pmat[kl],ij,kl,"cont1a");
      gmat[kl] += val*pmat[ij]; CHECK(val,pmat[ij],ij,kl,"cont1b");
    }
    
    inline void cont2(int ij, int kl, double val) {
      val *= -0.25;
      gmat[ij] += val*pmat[kl]; CHECK(4*val,0.25*pmat[kl],ij,kl,"cont2a");
      gmat[kl] += val*pmat[ij]; CHECK(4*val,0.25*pmat[ij],ij,kl,"cont2b");
    }
    
    inline void cont3(int ij, int kl, double val) {
      val *= -0.5;
      gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont3a");
      gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont3b");
    }
    
    inline void cont4(int ij, int kl, double val) {
      val *= 0.75;
      gmat[ij] += val*pmat[kl]; CHECK(4./3.*val,0.75*pmat[kl],ij,kl,"cont4a");
      gmat[kl] += val*pmat[ij]; CHECK(4./3.*val,0.75*pmat[ij],ij,kl,"cont4b");
    }
    
    inline void cont5(int ij, int kl, double val) {
      val *= 0.5;
      gmat[ij] += val*pmat[kl]; CHECK(2*val,0.5*pmat[kl],ij,kl,"cont5a");
      gmat[kl] += val*pmat[ij]; CHECK(2*val,0.5*pmat[ij],ij,kl,"cont5b");
    }
};

class LocalCLHFEnergyContribution {
  private:
    double * const pmat;

  public:
    double ec;
    double ex;

    void set_bound(double,double) {}
    
    LocalCLHFEnergyContribution(double *p) : pmat(p) {
      ec=ex=0;
    }
    ~LocalCLHFEnergyContribution() {}

    inline void cont1(int ij, int kl, double val) {
      ec += val*pmat[ij]*pmat[kl];
    }
    
    inline void cont2(int ij, int kl, double val) {
      ex -= 0.25*val*pmat[ij]*pmat[kl];
    }
    
    inline void cont3(int ij, int kl, double val) {
      ex -= 0.5*val*pmat[ij]*pmat[kl];
    }
    
    inline void cont4(int ij, int kl, double val) {
      ec += val*pmat[ij]*pmat[kl];
      ex -= 0.25*val*pmat[ij]*pmat[kl];
    }
    
    inline void cont5(int ij, int kl, double val) {
      ec += val*pmat[ij]*pmat[kl];
      ex -= 0.5*val*pmat[ij]*pmat[kl];
    }
};

class LocalCLHFGradContribution {
  private:
    double * const pmat;

  public:
    LocalCLHFGradContribution(double *p) : pmat(p) {}
    ~LocalCLHFGradContribution() {}

    inline double cont1(int ij, int kl) {
      return pmat[ij]*pmat[kl];
    }

    inline double cont2(int ij, int kl) {
      return pmat[ij]*pmat[kl];
    }
};

}