Commit cb809b82 authored by Mulas, Giacomo's avatar Mulas, Giacomo
Browse files

- add constructors to the "Common" objects allowing them to be copied 

- turn C4 from a structure to a class (to do the above), move its initialisation to its constructor
- eliminate "allocate_vectors()", redundant and including a (wrong) double allocation
parent f47b2156

src/cluster/cluster.cpp

+1 −14
Original line number Diff line number Diff line
@@ -107,20 +107,7 @@ void cluster(string config_file, string data_file, string output_path) { 
      c1->rzz[c1i] = gconf->sph_z[c1i];

      c1->ros[c1i] = sconf->radii_of_spheres[c1i];

    }

    C4 *c4 = new C4;

    c4->li = gconf->li;

    c4->le = gconf->le;

    c4->lm = (c4->li > c4-> le) ? c4->li : c4->le;

    c4->nv3j = (c4->lm * (c4->lm  + 1) * (2 * c4->lm + 7)) / 6;

    c4->nsph = nsph;

    // The following is needed to initialize C1_AddOns

    c4->litpo = c4->li + c4->li + 1;

    c4->litpos = c4->litpo * c4->litpo;

    c4->lmtpo = c4->li + c4->le + 1;

    c4->lmtpos = c4->lmtpo * c4->lmtpo;

    c4->nlim = c4->li * (c4->li + 2);

    c4->nlem = c4->le * (c4->le + 2);

    c4->lmpo = c4->lm + 1;

    C4 *c4 = new C4(gconf->li, gconf->le, nsph);

    C1_AddOns *c1ao = new C1_AddOns(c4);

    // End of add-ons initialization

    C6 *c6 = new C6(c4->lmtpo);

src/include/Commons.h

+290 −229
Original line number Diff line number Diff line
@@ -42,6 +42,8 @@ protected: 
  int lm;

  //! \brief NLMMT. QUESTION: definition?

  int nlmmt;

  //! \brief Number of configurations

  int configurations;

public:

  //! \brief QUESTION: definition?

  dcomplex **rmi;
@@ -92,7 +94,11 @@ public: 
   * \param iog: `int *` Vector of spherical units ID numbers.

   */

  C1(int ns, int l_max, int *nshl, int *iog);



  /*! \brief C1 instance constructor copying all contents from a preexisting template

   *

   * \param rhs: `C1` preexisting template.

   */

  C1(const C1& rhs);

  //! \brief C1 instance destroyer.

  ~C1();

};
@@ -105,6 +111,8 @@ class C2 { 
  int nsph;

  //! \brief Number of required orders.

  int nhspo;

  //! \brief QUESTION: what is nl?

  int nl;

public:

  //! \brief QUESTION: definition?

  dcomplex *ris;
@@ -125,6 +133,11 @@ public: 
   * \param npntts: `int`

   */

  C2(int ns, int nl, int npnt, int npntts);

  /*! \brief C2 instance constructor copying its contents from preexisting instance.

   *

   * \param rhs: `C2` object to copy contents from

   */

  C2(const C2& rhs);



  //! \brief C2 instance destroyer.

  ~C2();
@@ -150,6 +163,9 @@ public: 
  /*! \brief C3 instance constructor.

   */

  C3();

  /*! \brief C3 instance constructor copying its contents from a preexisting object.

   */

  C3(const C3& rhs);



  /*! \brief C3 instance destroyer.

   */
@@ -158,7 +174,8 @@ public: 


/*! \brief Representation of the FORTRAN C4 blocks.

 */

struct C4 {

class C4 {

public:

  //! \brief QUESTION: definition?

  int litpo;

  //! \brief QUESTION: definition?
@@ -183,6 +200,17 @@ struct C4 { 
  int nsph;

  //! \brief QUESTION: definition?

  int nv3j;



  /*! \brief C3 instance constructor.

   */

  C4(int li, int le, int nsph);

  /*! \brief C3 instance constructor copying its contents from a preexisting object.

   */

  C4(const C4& rhs);



  /*! \brief C3 instance destroyer.

   */

  ~C4();

};



/*! \brief Vectors and matrices that are specific to cluster C1 blocks.
@@ -196,6 +224,18 @@ protected: 
  int nlemt;

  //! \brief Maximum expansion order plus one. QUESTION: correct?

  int lmpo;

  //! \brief QUESTION: definition?

  int litpo;

  //! \brief QUESTION: definition?

  int lmtpo;

  //! \brief QUESTION: definition?

  int litpos;

  //! \brief QUESTION: definition?

  int lmtpos;

  //! \brief QUESTION: definition?

  int nv3j;

  //! \brief QUESTION: definition?

  int lm;



  /*! \brief Allocate the necessary common vectors depending on configuration.

   *
@@ -209,7 +249,7 @@ protected: 
   *

   * \param c4: `C4 *` Pointer to a C4 structure.

   */

	void allocate_vectors(C4 *c4);

  //void allocate_vectors(C4 *c4);

public:

  //! \brief QUESTION: definition?

  dcomplex *vh;
@@ -265,6 +305,11 @@ public: 
   * \param c4: `C4 *` Pointer to a C4 structure.

   */

  C1_AddOns(C4 *c4);

  /*! \brief C1_AddOns instance constructor copying contents from a preexisting object.

   *

   * \param rhs: `C1_AddOns` preexisting object to copy from.

   */

  C1_AddOns(const C1_AddOns& rhs);



  //! \brief C1_AddOns instance destroyer.

  ~C1_AddOns();
@@ -274,6 +319,8 @@ public: 
 */

class C6 {

public:

  //! \brief QUESTION: definition?

  int lmtpo;

  //! \brief QUESTION: definition?

  double *rac3j;
@@ -282,6 +329,11 @@ public: 
   * \param lmtpo: `int` QUESTION: definition?

   */

  C6(int lmtpo);

  /*! \brief C6 instance constructor copying contents from preexisting object.

   *

   * \param lmtpo: `int` QUESTION: definition?

   */

  C6(const C6& rhs);



  /*! \brief C6 instance destroyer.

   */
@@ -296,6 +348,10 @@ protected: 
  int gis_size_0;

  //! \brief Number of rows in the SAM matrix

  int sam_size_0;

  //! \brief QUESTION: definition?

  int nlem;

  //! \brief QUESTION: definition?

  int nlemt;

public:

  //! \brief QUESTION: definition?

  dcomplex **gis;
@@ -312,6 +368,11 @@ public: 
   * \param nlemt: `int` QUESTION: definition?

   */

  C9(int ndi, int nlem, int ndit, int nlemt);

  /*! \brief C9 instance constructor copying contents from preexisting object.

   *

   * \param rhs: `C9` preexisting object to copy from

   */

  C9(const C9& rhs);



  /*! \brief C9 instance destroyer.

   */

src/libnptm/Commons.cpp

+293 −24
Original line number Diff line number Diff line
@@ -72,6 +72,78 @@ C1::C1(int ns, int l_max, int *_nshl, int *_iog) { 
  }

}



C1::C1(const C1& rhs) {

  nlmmt = rhs.nlmmt;

  nsph = rhs.nsph;

  lm = rhs.lm;



  rmi = new dcomplex*[lm];

  rei = new dcomplex*[lm];

  for (int ri = 0; ri < lm; ri++) {

    rmi[ri] = new dcomplex[nsph]();

    rei[ri] = new dcomplex[nsph]();

    /*! Copy the contents from the template */

    for (int rj=0; rj<nsph; rj++) {

      rmi[ri][rj] = rhs.rmi[ri][rj];

      rei[ri][rj] = rhs.rei[ri][rj];

    }

  }

  w = new dcomplex*[nlmmt];

  for (int wi = 0; wi < nlmmt; wi++) {

    w[wi] = new dcomplex[4]();

    for (int wj=0; wj<4; wj++) w[wi][wj] = rhs.w[wi][wj];

  }

  configurations = rhs.configurations;

  vints = new dcomplex*[nsph];

  rc = new double*[configurations];

  nshl = new int[configurations]();

  for (int ni=0; ni<configurations; ni++) nshl[ni] = rhs.nshl[ni];

  iog = new int[nsph]();

  for (int ni=0; ni<nsph; ni++) iog[ni] = rhs.iog[ni];

  int conf_index = 0;

  for (int vi = 0; vi < nsph; vi++) {

    vints[vi] = new dcomplex[16]();

    for (int vj=0; vj<16; vj++) vints[vi][vj] = rhs.vints[vi][vj];

  }

  for (int ri=0; ri<configurations; ri++) {

    rc[ri] = new double[nshl[ri]]();

    for (int rj=0; rj<nshl[ri]; rj++) rc[ri][rj] = rhs.rc[ri][rj];

  }

  fsas = new dcomplex[nsph]();

  sscs = new double[nsph]();

  sexs = new double[nsph]();

  sabs = new double[nsph]();

  sqscs = new double[nsph]();

  sqexs = new double[nsph]();

  sqabs = new double[nsph]();

  gcsv = new double[nsph]();

  rxx = new double[nsph]();

  ryy = new double[nsph]();

  rzz = new double[nsph]();

  ros = new double[nsph]();



  sas = new dcomplex**[nsph];

  for (int si = 0; si < nsph; si++) {

    sas[si] = new dcomplex*[2];

    for (int sj=0; sj<2; sj++) {

      sas[si][sj] = new dcomplex[2]();

      for (int sk=0; sk<2; sk++) sas[si][sj][sk] = rhs.sas[si][sj][sk];

    }

    fsas[si] = rhs.fsas[si];

    sscs[si] = rhs.sscs[si];

    sexs[si] = rhs.sexs[si];

    sabs[si] = rhs.sabs[si];

    sqscs[si] = rhs.sqscs[si];

    sqexs[si] = rhs.sqexs[si];

    sqabs[si] = rhs.sqabs[si];

    gcsv[si] = rhs.gcsv[si];

    rxx[si] = rhs.rxx[si];

    ryy[si] = rhs.ryy[si];

    rzz[si] = rhs.rzz[si];

    ros[si] = rhs.ros[si];    

  }

}



C1::~C1() {

  delete[] rmi;

  delete[] rei;
@@ -114,11 +186,17 @@ C1_AddOns::C1_AddOns(C4 *c4) { 
  nsph = c4->nsph;

  lmpo = c4->lmpo;

  nlemt = 2 * c4->nlem;

  vh = new dcomplex[(nsph * nsph - 1) * c4->litpo]();

  vj0 = new dcomplex[nsph * c4->lmtpo]();

  litpo = c4->litpo;

  lmtpo = c4->lmtpo;

  litpos = c4->litpos;

  lmtpos = c4->lmtpos;

  nv3j = c4->nv3j;

  lm = c4->lm;

  vh = new dcomplex[(nsph * nsph - 1) * litpo]();

  vj0 = new dcomplex[nsph * lmtpo]();

  vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?

  vyhj = new dcomplex[(nsph * nsph - 1) * c4->litpos]();

  vyj0 = new dcomplex[nsph * c4->lmtpos]();

  vyhj = new dcomplex[(nsph * nsph - 1) * litpos]();

  vyj0 = new dcomplex[nsph * lmtpos]();

  am0m = new dcomplex*[nlemt];

  for (int ai = 0; ai < nlemt; ai++) {

    am0m[ai] = new dcomplex[nlemt]();
@@ -140,12 +218,103 @@ C1_AddOns::C1_AddOns(C4 *c4) { 
  ecscp = new dcomplex[2]();

  scscpm = new dcomplex[2]();

  ecscpm = new dcomplex[2]();

  allocate_vectors(c4);

  v3j0 = new double[nv3j]();

  ind3j = new int*[lmpo];

  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[lm]();

  //allocate_vectors(c4);

  sscs = new double[nsph]();

  ecscm = new double[2]();

  scscm = new double[2]();

  scsc = new double[2]();

  ecsc = new double[2]();

}



C1_AddOns::C1_AddOns(const C1_AddOns& rhs) {

  nsph = rhs.nsph;

  lmpo = rhs.lmpo;

  nlemt = rhs.nlemt;

  litpo = rhs.litpo;

  lmtpo = rhs.lmtpo;

  litpos = rhs.litpos;

  lmtpos = rhs.lmtpos;

  nv3j = rhs.nv3j;

  lm = rhs.lm;

  int vhsize=(nsph * nsph - 1) * litpo;

  vh = new dcomplex[vhsize]();

  for (int hi=0; hi<vhsize; hi++) vh[hi] = rhs.vh[hi];

  int vj0size=nsph * lmtpo;

  vj0 = new dcomplex[vj0size]();

  for (int hi=0; hi<vj0size; hi++) vj0[hi] = rhs.vj0[hi];

  vj = new dcomplex[1](); // QUESTION: is 1 really enough for a general case?

  vj[0] = rhs.vj[0];

  int vyhjsize=(nsph * nsph - 1) * litpos;

  vyhj = new dcomplex[vyhjsize]();

  for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi];

  int vyj0size = nsph * lmtpos;

  vyj0 = new dcomplex[vyj0size]();

  for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];

  am0m = new dcomplex*[nlemt];

  for (int ai = 0; ai < nlemt; ai++) {

    am0m[ai] = new dcomplex[nlemt]();

    for (int aj = 0; aj < nlemt; aj++) am0m[ai][aj] = rhs.am0m[ai][aj];

  }

  vint = new dcomplex[16](); // QUESTION: is dimension 16 generally fixed?

  vintm = new dcomplex[16]();

  vintt = new dcomplex[16]();

  for (int hi=0; hi<16; hi++) {

    vint[hi] = rhs.vint[hi];

    vintm[hi] = rhs.vintm[hi];

    vintt[hi] = rhs.vintt[hi];

  }

  vints = new dcomplex*[nsph];

  for (int vi = 0; vi < nsph; vi++) {

    vints[vi] = new dcomplex[16]();

    for (int hj=0; hj<16; hj++) vints[vi][hj] = rhs.vints[vi][hj];

  }

  fsac = new dcomplex*[2];

  sac = new dcomplex*[2];

  fsacm = new dcomplex*[2];

  for (int fi = 0; fi < 2; fi++) {

    fsac[fi] = new dcomplex[2]();

    sac[fi] = new dcomplex[2]();

    fsacm[fi] = new dcomplex[2]();

    for (int fj = 0; fj < 2; fj++) {

      fsac[fi][fj] = rhs.fsac[fi][fj];

      sac[fi][fj] = rhs.sac[fi][fj];

      fsacm[fi][fj] = rhs.fsacm[fi][fj];

    }

  }

  scscp = new dcomplex[2]();

  ecscp = new dcomplex[2]();

  scscpm = new dcomplex[2]();

  ecscpm = new dcomplex[2]();

  v3j0 = new double[nv3j]();

  for (int hi=0; hi<nv3j; hi++) v3j0[hi] = rhs.v3j0[hi];

  ind3j = new int*[lmpo];

  for (int ii = 0; ii < lmpo; ii++) {

    ind3j[ii] = new int[lm]();

    for (int ij=0; ij<lm; ij++) ind3j[ii][ij] = rhs.ind3j[ii][ij];

  }

  // vyhj = new dcomplex[int vyhjsize = (nsph * nsph - 1) * litpos]();

  // for (int hi=0; hi<vyhjsize; hi++) vyhj[hi] = rhs.vyhj[hi]

  // vyj0 = new dcomplex[int vyj0size=lmtpos * nsph]();

  // for (int hi=0; hi<vyj0size; hi++) vyj0[hi] = rhs.vyj0[hi];

  sscs = new double[nsph]();

  for (int hi=0; hi<nsph; hi++) sscs[hi] = rhs.sscs[hi];

  ecscm = new double[2]();

  scscm = new double[2]();

  scsc = new double[2]();

  ecsc = new double[2]();

  for (int hi=0; hi<2; hi++) {

    scscp[hi] = rhs.scscp[hi];

    ecscp[hi] = rhs.ecscp[hi];

    scscpm[hi] = rhs.scscpm[hi];

    ecscpm[hi] = rhs.ecscpm[hi];

    ecscm[hi] = rhs.ecscm[hi];

    scscm[hi] = rhs.scscm[hi];

    scsc[hi] = rhs.scsc[hi];

    ecsc[hi] = rhs.ecsc[hi];

  }

}



C1_AddOns::~C1_AddOns() {
@@ -184,25 +353,26 @@ C1_AddOns::~C1_AddOns() { 
  delete[] ecsc;

}



void C1_AddOns::allocate_vectors(C4 *c4) {

  // This calculates the size of v3j0

  int lm = (c4->li > c4->le) ? c4->li : c4->le;

  const int nv3j = c4->nv3j;

  v3j0 = new double[nv3j]();

  ind3j = new int*[lmpo];

  for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();

  // This calculates the size of vyhj

  int ivy = (nsph * nsph - 1) * c4->litpos;

  vyhj = new dcomplex[ivy]();

  // This calculates the size of vyj0

  ivy = c4->lmtpos * c4->nsph;

  vyj0 = new dcomplex[ivy]();

}

// void C1_AddOns::allocate_vectors(C4 *c4) {

//   // This calculates the size of v3j0

//   //int lm = (c4->li > c4->le) ? c4->li : c4->le;

//   const int nv3j = c4->nv3j;

//   v3j0 = new double[nv3j]();

//   ind3j = new int*[lmpo];

//   for (int ii = 0; ii < lmpo; ii++) ind3j[ii] = new int[c4->lm]();

//   // This calculates the size of vyhj

//   // int ivy = (nsph * nsph - 1) * c4->litpos;

//   // vyhj = new dcomplex[ivy]();

//   // This calculates the size of vyj0

//   // ivy = c4->lmtpos * c4->nsph;

//   // vyj0 = new dcomplex[ivy]();

// }



C2::C2(int ns, int nl, int npnt, int npntts) {

C2::C2(int ns, int _nl, int npnt, int npntts) {

  nsph = ns;

  int max_n = (npnt > npntts) ? npnt : npntts;

  nhspo = 2 * max_n - 1;

  nl = _nl;

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

  vkt = new dcomplex[nsph]();
@@ -210,6 +380,26 @@ C2::C2(int ns, int nl, int npnt, int npntts) { 
  vsz = new double[nsph]();

}



C2::C2(const C2& rhs) {

  nsph = rhs.nsph;

  nhspo = rhs.nhspo;

  nl = rhs.nl;

  ris = new dcomplex[nhspo]();

  dlri = new dcomplex[nhspo]();

  for (int ind=0; ind<nhspo; ind++) {

    ris[ind] = rhs.ris[ind];

    dlri[ind] = rhs.dlri[ind];

  }

  vkt = new dcomplex[nsph]();

  vsz = new double[nsph]();

  for (int ind=0; ind<nsph; ind++) {

    vkt[ind] = rhs.vkt[ind];

    vsz[ind] = rhs.vsz[ind];

  }

  dc0 = new dcomplex[nl]();

  for (int ind=0; ind<nl; ind++) dc0[ind] = rhs.dc0[ind];

}



C2::~C2() {

  delete[] ris;

  delete[] dlri;
@@ -229,21 +419,80 @@ C3::C3() { 
  acs = 0.0;

}



C3::C3(const C3& rhs) {

  tsas = new dcomplex*[2];

  for (int ti=0; ti<2; ti++) {

    tsas[ti] = new dcomplex[2];

    for (int tj=0; tj<2; tj++) tsas[ti][tj] = rhs.tsas[ti][tj];

  }

  tfsas = rhs.tfsas;

  gcs = rhs.gcs;

  scs = rhs.scs;

  ecs = rhs.ecs;

  acs = rhs.acs;

}



C3::~C3() {

  delete[] tsas[1];

  delete[] tsas[0];

  delete[] tsas;

}



C6::C6(int lmtpo) {

C4::C4(int _li, int _le, int _nsph) {

    li = _li;

    le = _le;

    lm = (li > le) ? li : le;

    nv3j = (lm * (lm  + 1) * (2 * lm + 7)) / 6;

    nsph = _nsph;

    // The following is needed to initialize C1_AddOns

    litpo = li + li + 1;

    litpos = litpo * litpo;

    lmtpo = li + le + 1;

    lmtpos = lmtpo * lmtpo;

    nlim = li * (li + 2);

    nlem = le * (le + 2);

    lmpo = lm + 1;

}



C4::C4(const C4& rhs) {

    li = rhs.li;

    le = rhs.le;

    lm = rhs.lm;

    nv3j = rhs.nv3j;

    nsph = rhs.nsph;

    // The following is needed to initialize C1_AddOns

    litpo = rhs.litpo;

    litpos = rhs.litpos;

    lmtpo = rhs.lmtpo;

    lmtpos = rhs.lmtpos;

    nlim = rhs.nlim;

    nlem = rhs.nlem;

    lmpo = rhs.lmpo;

}



C4::~C4() {

}



C6::C6(int _lmtpo) {

  lmtpo = _lmtpo;

  rac3j = new double[lmtpo]();

}



C6::C6(const C6& rhs) {

  lmtpo = rhs.lmtpo;

  rac3j = new double[lmtpo]();

  for (int ri=0; ri<lmtpo; ri++) rac3j[ri] = rhs.rac3j[ri];

}



C6::~C6() {

  delete[] rac3j;

}



C9::C9(int ndi, int nlem, int ndit, int nlemt) {

C9::C9(int ndi, int _nlem, int ndit, int _nlemt) {

  gis_size_0 = ndi;

  sam_size_0 = ndit;

  nlem = _nlem;

  nlemt = _nlemt;

  gis = new dcomplex*[ndi];

  gls = new dcomplex*[ndi];

  for (int gi = 0; gi < ndi; gi++) {
@@ -252,8 +501,28 @@ C9::C9(int ndi, int nlem, int ndit, int nlemt) { 
  }

  sam = new dcomplex*[ndit];

  for (int si = 0; si < ndit; si++) sam[si] = new dcomplex[nlemt]();

  gis_size_0 = ndi;

  sam_size_0 = ndit;

}



C9::C9(const C9& rhs) {

  gis_size_0 = rhs.gis_size_0;

  sam_size_0 = rhs.sam_size_0;

  nlem = rhs.nlem;

  nlemt = rhs.nlemt;

  gis = new dcomplex*[gis_size_0];

  gls = new dcomplex*[gis_size_0];

  for (int gi = 0; gi < gis_size_0; gi++) {

    gis[gi] = new dcomplex[nlem]();

    gls[gi] = new dcomplex[nlem]();

    for (int gj=0; gj<nlem; gj++) {

      gis[gi][gj] = rhs.gis[gi][gj];

      gls[gi][gj] = rhs.gls[gi][gj];

    }

  }

  sam = new dcomplex*[sam_size_0];

  for (int si = 0; si < sam_size_0; si++) {

    sam[si] = new dcomplex[nlemt]();

    for (int sj=0; sj<nlemt; sj++) sam[si][sj] = rhs.sam[si][sj];

  }

}



C9::~C9() {