Commit efa73552 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Enable independent internal and external orders in C1

parent e50f1ece

src/include/Commons.h

+6 −0
Original line number Diff line number Diff line
@@ -58,6 +58,10 @@ class C1 { 
protected:

  //! \brief Maximum order of field expansion.

  int lm;

  //! \brief Maximum order of internal field expansion.

  int li;

  //! \brief Maximum order of external field expansion.

  int le;

  //! \brief Contiguous RMI vector.

  dcomplex *vec_rmi;

  //! \brief Contiguous REI vector.
@@ -70,6 +74,8 @@ protected: 
public:

  //! \brief Number of spheres.

  int nsph;

  //! \brief NLEMT = 2 * LE * (LE + 2).

  int nlemt;

  //! \brief NLMMT = 2 * LM * (LM + 2).

  int nlmmt;

  //! \brief Number of configurations

src/libnptm/Commons.cpp

+48 −38
Original line number Diff line number Diff line
@@ -36,25 +36,26 @@ 


C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) {

  lm = gconf->l_max;

  int li = gconf->li;

  int le = gconf->le;

  li = gconf->li;

  le = gconf->le;

  if (lm == 0) {

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

  }

  nsph = gconf->number_of_spheres;

  nlmmt = 2 * (lm * (lm + 2));

  nlemt = 2 * (le * (le + 2));



  vec_rmi = new dcomplex[nsph * lm]();

  vec_rei = new dcomplex[nsph * lm]();

  rmi = new dcomplex*[lm];

  rei = new dcomplex*[lm];

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

  vec_rmi = new dcomplex[nsph * li]();

  vec_rei = new dcomplex[nsph * li]();

  rmi = new dcomplex*[li];

  rei = new dcomplex*[li];

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

    rmi[ri] = &(vec_rmi[nsph * ri]);

    rei[ri] = &(vec_rei[nsph * ri]);

  }

  vec_w = new dcomplex[4 * nlmmt]();

  w = new dcomplex*[nlmmt];

  for (int wi = 0; wi < nlmmt; wi++) w[wi] = &(vec_w[4 * wi]);

  vec_w = new dcomplex[4 * nlemt]();

  w = new dcomplex*[nlemt];

  for (int wi = 0; wi < nlemt; wi++) w[wi] = &(vec_w[4 * wi]);

  configurations = sconf->configurations;

  vint = new dcomplex[16]();

  vec_vints = new dcomplex[nsph * 16]();
@@ -108,28 +109,31 @@ C1::C1(GeometryConfiguration *gconf, ScattererConfiguration *sconf) { 


C1::C1(const C1& rhs) {

  nlmmt = rhs.nlmmt;

  nlemt = rhs.nlemt;

  nsph = rhs.nsph;

  lm = rhs.lm;

  li = rhs.li;

  le = rhs.le;



  vec_rmi = new dcomplex[nsph * lm]();

  vec_rei = new dcomplex[nsph * lm]();

  for (long li=0; li<(nsph*lm); li++) {

    vec_rmi[li] = rhs.vec_rmi[li];

    vec_rei[li] = rhs.vec_rei[li];

  vec_rmi = new dcomplex[nsph * li]();

  vec_rei = new dcomplex[nsph * li]();

  for (long ili=0; ili < (nsph * li); ili++) {

    vec_rmi[ili] = rhs.vec_rmi[ili];

    vec_rei[ili] = rhs.vec_rei[ili];

  }

  rmi = new dcomplex*[lm];

  rei = new dcomplex*[lm];

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

  rmi = new dcomplex*[li];

  rei = new dcomplex*[li];

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

    rmi[ri] = &(vec_rmi[nsph * ri]);

    rei[ri] = &(vec_rei[nsph * ri]);

    /*! The contents were already copied via vec_rmi and vec_rei */

  }

  vec_w = new dcomplex[4 * nlmmt]();

  for (long li=0; li<(4*nlmmt); li++) {

    vec_w[li] = rhs.vec_w[li];

  vec_w = new dcomplex[4 * nlemt]();

  for (long ili=0; ili < (4*nlemt); ili++) {

    vec_w[ili] = rhs.vec_w[ili];

  }    

  w = new dcomplex*[nlmmt];

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

  w = new dcomplex*[nlemt];

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

    w[wi] = &(vec_w[4 * wi]);

    /*! The contents were already copied via vec_w */

  }
@@ -198,23 +202,26 @@ C1::C1(const C1& rhs) { 
#ifdef MPI_VERSION

C1::C1(const mixMPI *mpidata) {

  MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);

  vec_rmi = new dcomplex[nsph * lm]();

  vec_rei = new dcomplex[nsph * lm]();

  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  rmi = new dcomplex*[lm];

  rei = new dcomplex*[lm];

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

  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);

  vec_rmi = new dcomplex[nsph * li]();

  vec_rei = new dcomplex[nsph * li]();

  MPI_Bcast(vec_rmi, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rei, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  rmi = new dcomplex*[li];

  rei = new dcomplex*[li];

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

    rmi[ri] = &(vec_rmi[nsph * ri]);

    rei[ri] = &(vec_rei[nsph * ri]);

    /*! The contents were copied already via vec_rmi and vec_rei */

  }

  vec_w = new dcomplex[4 * nlmmt]();

  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  w = new dcomplex*[nlmmt];

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

  vec_w = new dcomplex[4 * nlemt]();

  MPI_Bcast(vec_w, 4*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  w = new dcomplex*[nlemt];

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

    w[wi] = &(vec_w[4 * wi]);

    /*! The contents were copied already via vec_w */

  }
@@ -279,11 +286,14 @@ C1::C1(const mixMPI *mpidata) { 


void C1::mpibcast(const mixMPI *mpidata) {

  MPI_Bcast(&nlmmt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&nlemt, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&nsph, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&lm, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rmi, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rei, nsph*lm, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_w, 4*nlmmt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&li, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(&le, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rmi, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_rei, nsph*li, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_w, 4*nlemt, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(&configurations, 1, MPI_INT, 0, MPI_COMM_WORLD);

  MPI_Bcast(vint, 16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);

  MPI_Bcast(vec_vints, nsph*16, MPI_C_DOUBLE_COMPLEX, 0, MPI_COMM_WORLD);