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

Merge branch 'split_in_out_orders' into 'master' 

Split inner and outer orders

See merge request giacomo.mulas/np_tmcode!59
parents a474a1a3 5f071869

.gitignore

+0 −1
Original line number Diff line number Diff line 
build/aclocal.m4

build/autom4te.cache

build/config.*

build/cluster/*

build/aclocal.m4

0 → 100644
+10298 −0

File added.

Preview size limit exceeded, changes collapsed.

src/cluster/cluster.cpp

+5 −1
Original line number Diff line number Diff line
@@ -536,6 +536,10 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      

    delete sconf;

    delete gconf;

#ifdef USE_MAGMA

    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");

    magma_finalize();

#endif

    chrono::time_point<chrono::high_resolution_clock> t_end = chrono::high_resolution_clock::now();

    elapsed = t_end - t_start;

    string message = "INFO: Calculation lasted " + to_string(elapsed.count()) + "s.\n";
@@ -682,7 +686,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  int lm = 0;

  if (le > lm) lm = le;

  if (li > lm) lm = li;

  int nsph = gconf->number_of_spheres;

  int nsph = sconf->number_of_spheres;

  np_int mxndm = gconf->mxndm;

  int iavm = gconf->iavm;

  int inpol = gconf->in_pol;

src/include/Commons.h

+40 −35
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,8 +74,8 @@ protected: 
public:

  //! \brief Number of spheres.

  int nsph;

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

  int nlmmt;

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

  int nlemt;

  //! \brief Number of configurations

  int configurations;

  //! \brief QUESTION: definition?
@@ -84,21 +88,21 @@ public: 
  dcomplex *fsas;

  //! \brief QUESTION: definition?

  dcomplex *vint;

  //! \brief QUESTION: definition?

  //! \brief Components of the scattered field intensities.

  dcomplex **vints;

  //! \brief QUESTION: definition?

  //! \brief Single sphere scattering cross-sections.

  double *sscs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere extinction cross-sections.

  double *sexs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere absorption cross-sections.

  double *sabs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere scattering efficiencies.

  double *sqscs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere extinction efficiencies.

  double *sqexs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere absorption efficiencies.

  double *sqabs;

  //! \brief QUESTION: definition?

  //! \brief Single sphere geometric cross-section vector.

  double *gcsv;

  //! \brief Vector of sphere X coordinates.

  double *rxx;
@@ -165,7 +169,7 @@ public: 
  dcomplex *ris;

  //! \brief QUESTION: definition?

  dcomplex *dlri;

  //! \brief QUESTION: definition?

  //! \brief Vector of dielectric constants.

  dcomplex *dc0;

  //! \brief QUESTION: definition?

  dcomplex *vkt;
@@ -212,13 +216,13 @@ public: 
  dcomplex tfsas;

  //! \brief QUESTION: definition?

  dcomplex **tsas;

  //! \brief QUESTION: definition?

  //! \brief Total geometric cross-section.

  double gcs;

  //! \brief QUESTION: definition?

  //! \brief Total scattering cross-section.

  double scs;

  //! \brief QUESTION: definition?

  //! \brief Total extinction cross-section.

  double ecs;

  //! \brief QUESTION: definition?

  //! \brief Total absorption cross-section.

  double acs;



  /*! \brief C3 instance constructor.
@@ -265,17 +269,17 @@ public: 
  int lmtpos;

  //! \brief Internal field expansion order.

  int li;

  //! \brief QUESTION: definition?

  //! \brief NLIM = LI * (LI + 2)

  int nlim;

  //! \brief External field expansion order.

  int le;

  //! \brief QUESTION: definition?

  //! \brief NLEM = LE * (LE + 2)

  int nlem;

  //! \brief Maximum field expansion order.

  int lm;

  //! \brief Number of spheres.

  int nsph;

  //! \brief QUESTION: definition?

  //! \brief NV3J = (LM * (LM + 1) * (2 * LM + 7)) / 6

  int nv3j;



  /*! \brief C4 instance constructor.
@@ -317,21 +321,21 @@ class C1_AddOns { 
protected:

  //! \brief Number of spheres.

  int nsph;

  //! \brief QUESTION: definition?

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

  int nlemt;

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

  //! \brief LMPO = LM + 1

  int lmpo;

  //! \brief QUESTION: definition?

  //! \brief LITPO = 2 * LI + 1

  int litpo;

  //! \brief QUESTION: definition?

  //! \brief LMTPO = 2 * LM + 1

  int lmtpo;

  //! \brief QUESTION: definition?

  //! \brief LITPOS = LITPO * LITPO

  int litpos;

  //! \brief QUESTION: definition?

  //! \brief LMTPOS = LMTPO * LMTPO

  int lmtpos;

  //! \brief QUESTION: definition?

  //! \brief NV3J = (LM * (LM + 1) * (2 * LM + 7)) / 6

  int nv3j;

  //! \brief QUESTION: definition?

  //! \brief Maximum field expansion order.

  int lm;



public:
@@ -347,7 +351,7 @@ public: 
  dcomplex *vyj0;

  //! \brief QUESTION: definition?

  dcomplex **am0m;

  //! \brief QUESTION: definition?

  //! \brief Vectorized AM0M matrix.

  dcomplex *am0v;

  //! \brief QUESTION: definition?

  dcomplex *vintm;
@@ -465,9 +469,9 @@ protected: 
  int sam_size_0;



public:

  //! \brief QUESTION: definition?

  //! \brief NLEM = LE * (LE + 2)

  int nlem;

  //! \brief QUESTION: definition?

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

  int nlemt;

  //! \brief QUESTION: definition?

  dcomplex **gis;
@@ -478,10 +482,10 @@ public: 


  /*! \brief C9 instance constructor.

   *

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

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

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

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

   * \param ndi: `int` NDI = NSPH * LI * (LI + 2)

   * \param nlem: `int` NLEM = LE * (LE + 2)

   * \param ndit: `int` NDIT = 2 * NSPH * LI * (LI + 2)

   * \param nlemt: `int` NLEMT = 2 * LE * (LE + 2)

   */

  C9(int ndi, int nlem, int ndit, int nlemt);
@@ -560,7 +564,7 @@ public: 
  C6 *c6;

  //! \brief Pointer to a C9 structure.

  C9 *c9;

  //! \brief Pointer to a formatted output file.

  //! \brief Vector of geometric asymmetry factors.

  double *gaps;

  double **tqse;

  dcomplex **tqspe;
@@ -595,6 +599,7 @@ public: 
  double *unsmp;

  double *upmp;

  double *upsmp;

  //! \brief Scattering angle.

  double scan;

  double cfmp;

  double sfmp;
@@ -605,7 +610,7 @@ public: 
  dcomplex *am_vector;

  dcomplex **am;

  dcomplex arg;

  //! \brief Wave vector.

  //! \brief Vacuum magnitude of wave vector.

  double vk;

  //! \brief Wave number.

  double wn;

src/libnptm/Commons.cpp

+53 −43
Original line number Diff line number Diff line
@@ -36,25 +36,28 @@ 


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

  lm = gconf->l_max;

  int li = gconf->li;

  int le = gconf->le;

  li = gconf->li;

  le = gconf->le;

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

  if (lm == 0) {

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

  } else {

    if (li == 0) li = lm;

    if (le == 0) nlemt = 2 * (lm * (lm + 2));

  }

  nsph = gconf->number_of_spheres;

  nlmmt = 2 * (lm * (lm + 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]();
@@ -107,29 +110,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 */

  }
@@ -197,24 +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 */

  }
@@ -278,12 +285,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);
@@ -755,7 +764,8 @@ C4::C4(GeometryConfiguration *gconf) { 
  // The following is needed to initialize C1_AddOns

  litpo = li + li + 1;

  litpos = litpo * litpo;

  lmtpo = li + le + 1;

  //lmtpo = li + le + 1;

  lmtpo = lm + lm + 1;

  lmtpos = lmtpo * lmtpo;

  nlim = li * (li + 2);

  nlem = le * (le + 2);