Commit 8f4d6a56 authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Use ParticleDescriptor in CLUSTER

parent c87351ce

src/cluster/cluster.cpp

+78 −78
Original line number Diff line number Diff line
@@ -200,7 +200,7 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      char virtual_line[256];

      // Create and initialise pristine cid for MPI proc 0 and thread 0

      ClusterIterationData *cid = new ClusterIterationData(gconf, sconf, mpidata, device_count);

      const int ndi = cid->c4->nsph * cid->c4->nlim;

      const int ndi = cid->c1->nsph * cid->c1->nlim;

      np_int ndit = 2 * ndi;

      logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");

      time_logger->log("INFO: Size of matrices to invert: " + to_string((int64_t)ndit) + " x " + to_string((int64_t)ndit) +".\n");
@@ -212,12 +212,12 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      p_output->append_line(virtual_line);

#ifdef USE_ILP64

      sprintf(virtual_line, " %5d%5d%5d%5ld%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#else

      sprintf(virtual_line, " %5d%5d%5d%5d%5d%5d%5d%5d%5d\n",

	      nsph, cid->c4->li, cid->c4->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      nsph, cid->c1->li, cid->c1->le, gconf->mxndm, gconf->in_pol, gconf->npnt,

	      gconf->npntts, gconf->iavm, gconf->iavm

	      );

#endif
@@ -266,8 +266,8 @@ void cluster(const string& config_file, const string& data_file, const string& o 
      p_output->append_line(virtual_line);

      sprintf(virtual_line, " \n");

      p_output->append_line(virtual_line);

      str(sconf, cid->c1, cid->c1ao, cid->c3, cid->c4, cid->c6);

      thdps(cid->c4->lm, cid->zpv);

      str(sconf, cid->c1, cid->c3, cid->c6);

      thdps(cid->c1->lm, cid->zpv);

      double exdc = sconf->exdc;

      double exri = sqrt(exdc);

      sprintf(virtual_line, "  REFR. INDEX OF EXTERNAL MEDIUM=%15.7lE\n", exri);
@@ -546,6 +546,8 @@ void cluster(const string& config_file, const string& data_file, const string& o 
    logger->log(message);

    logger->log("Finished: output written to " + output_path + "/c_OCLU\n");

    time_logger->log(message);

    fclose(timing_file);

    delete time_logger;

  } // end instructions block of MPI process 0

  

    //===============================
@@ -661,8 +663,6 @@ void cluster(const string& config_file, const string& data_file, const string& o 
    logger->log("INFO: Process " + to_string(mpidata->rank) + " finalizes MAGMA.\n");

    magma_finalize();

#endif

    // fclose(timing_file);

    // delete time_logger;

    delete logger;

#ifdef MPI_VERSION

  }
@@ -694,7 +694,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  int npntts = gconf->npntts;

  int isam = gconf->iavm;

  int jwtm = gconf->jwtm;

  np_int ndit = 2 * nsph * cid->c4->nlim;

  np_int ndit = 2 * nsph * cid->c1->nlim;

  int isq, ibf;

  int last_configuration;
@@ -719,7 +719,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
    sprintf(virtual_line, "  XI=%15.7lE\n", xi);

    output->append_line(virtual_line);

  }

  hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1, cid->c1ao, cid->c4);

  hjv(exri, vkarg, jer, lcalc, cid->arg, cid->c1);

  if (jer != 0) {

    sprintf(virtual_line, "  STOP IN HJV\n");

    output->append_line(virtual_line);
@@ -749,7 +749,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
      }

      if (nsh % 2 == 0) cid->c2->dc0[ici] = exdc;

      dme(

	  cid->c4->li, i132, npnt, npntts, vkarg, exdc, exri,

	  cid->c1->li, i132, npnt, npntts, vkarg, exdc, exri,

	  cid->c1, cid->c2, jer, lcalc, cid->arg, last_configuration

	  );

      if (jer != 0) {
@@ -783,7 +783,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  outam0->write_to_disk(outam0_name);

  delete outam0;

#endif

  cms(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6);

  cms(cid->am, cid->c1, cid->c6);

#ifdef DEBUG_AM

  VirtualAsciiFile *outam1 = new VirtualAsciiFile();

  string outam1_name = output_path + "/c_AM1_JXI" + to_string(jxi488) + ".txt";
@@ -835,7 +835,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
#ifdef USE_NVTX

  nvtxRangePush("Average calculation");

#endif

  ztm(cid->am, cid->c1, cid->c1ao, cid->c4, cid->c6, cid->c9);

  ztm(cid->am, cid->c1, cid->c6, cid->c9);

#ifdef DEBUG_AM

  VirtualAsciiFile *outam3 = new VirtualAsciiFile();

  string outam3_name = output_path + "/c_AM3_JXI" + to_string(jxi488) + ".txt";
@@ -849,11 +849,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
#endif

  if (idfc >= 0) {

    if (jxi488 == jwtm) {

      int nlemt = 2 * cid->c4->nlem;

      int nlemt = 2 * cid->c1->nlem;

      string ttms_name = output_path + "/c_TTMS.hd5";

      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m, "HDF5");

      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1->am0m, "HDF5");

      ttms_name = output_path + "/c_TTMS";

      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1ao->am0m);

      TransitionMatrix::write_binary(ttms_name, nlemt, lm, cid->vk, exri, cid->c1->am0m);

    }

  }

  // label 156: continue from here
@@ -868,11 +868,11 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  double cs0 = 0.25 * cid->vk * cid->vk * cid->vk / acos(0.0);

  double csch = 0.0, qschu = 0.0, pschu = 0.0, s0mag = 0.0;

  dcomplex s0 = 0.0 + 0.0 * I;

  scr0(cid->vk, exri, cid->c1, cid->c1ao, cid->c3, cid->c4);

  scr0(cid->vk, exri, cid->c1, cid->c3);

  double sqk = cid->vk * cid->vk * exdc;

  aps(cid->zpv, cid->c4->li, nsph, cid->c1, sqk, cid->gaps);

  rabas(inpol, cid->c4->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);

  if (cid->c4->li != cid->c4->le) {

  aps(cid->zpv, cid->c1->li, nsph, cid->c1, sqk, cid->gaps);

  rabas(inpol, cid->c1->li, nsph, cid->c1, cid->tqse, cid->tqspe, cid->tqss, cid->tqsps);

  if (cid->c1->li != cid->c1->le) {

    sprintf(virtual_line, "     SPHERES; LMX=LI\n");

    output->append_line(virtual_line);

  }
@@ -936,8 +936,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
  output->append_line(virtual_line);

  // tppoan.write(reinterpret_cast<char *>(&(cid->vk)), sizeof(double));

  vtppoanp->append_line(VirtualBinaryLine(cid->vk));

  pcrsm0(cid->vk, exri, inpol, cid->c1, cid->c1ao, cid->c4);

  apcra(cid->zpv, cid->c4->le, cid->c1ao->am0m, inpol, sqk, cid->gapm, cid->gappm);

  pcrsm0(cid->vk, exri, inpol, cid->c1);

  apcra(cid->zpv, cid->c1->le, cid->c1->am0m, inpol, sqk, cid->gapm, cid->gappm);

#ifdef USE_NVTX

  nvtxRangePop();

#endif
@@ -959,19 +959,19 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	upvmp(th, ph, 0, cost, sint, cosp, sinp, cid->u, cid->upmp, cid->unmp);

	if (isam >= 0) {

	  wmamp(

		0, cost, sint, cosp, sinp, inpol, cid->c4->le, 0,

		0, cost, sint, cosp, sinp, inpol, cid->c1->le, 0,

		nsph, cid->argi, cid->u, cid->upmp, cid->unmp, cid->c1

		);

	  // label 182

	  apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	  raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	  apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	  raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	  jw = 1;

	}

      } else { // label 180, NK == 1 AND JXI488 == 1

	if (isam >= 0) {

	  // label 182

	  apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	  raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	  apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	  raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	  jw = 1;

	}

      }
@@ -1003,7 +1003,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    upvmp(ths, phs, icspnv, costs, sints, cosps, sinps, cid->us, cid->upsmp, cid->unsmp);

	    if (isam >= 0)

	      wmamp(

		    2, costs, sints, cosps, sinps, inpol, cid->c4->le,

		    2, costs, sints, cosps, sinps, inpol, cid->c1->le,

		    0, nsph, cid->args, cid->us, cid->upsmp, cid->unsmp, cid->c1

		    );

	  }
@@ -1016,7 +1016,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    if (isam < 0) {

	      wmasp(

		    cost, sint, cosp, sinp, costs, sints, cosps, sinps,

		    cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c4->le,

		    cid->u, cid->up, cid->un, cid->us, cid->ups, cid->uns, isq, ibf, inpol, cid->c1->le,

		    0, nsph, cid->argi, cid->args, cid->c1

		    );

	    } else { // label 192
@@ -1029,10 +1029,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    }

	  }

	  // label 194

	  if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c1ao, cid->c4, cid->c6);

	  if (iavm == 1) crsm1(cid->vk, exri, cid->c1, cid->c6);

	  if (isam < 0) {

	    apc(cid->zpv, cid->c4->le, cid->c1ao->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	    raba(cid->c4->le, cid->c1ao->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	    apc(cid->zpv, cid->c1->le, cid->c1->am0m, cid->c1->w, sqk, cid->gap, cid->gapp);

	    raba(cid->c1->le, cid->c1->am0m, cid->c1->w, cid->tqce, cid->tqcpe, cid->tqcs, cid->tqcps);

	    jw = 1;

	  }

	  // label 196
@@ -1048,7 +1048,7 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	  vtppoanp->append_line(VirtualBinaryLine(cid->scan));

	  if (jaw != 0) {

	    jaw = 0;

	    mextc(cid->vk, exri, cid->c1ao->fsacm, cid->cextlr, cid->cext);

	    mextc(cid->vk, exri, cid->c1->fsacm, cid->cextlr, cid->cext);

	    // We now have some implicit loops writing to binary

	    for (int i = 0; i < 4; i++) {

	      for (int j = 0; j < 4; j++) {
@@ -1058,22 +1058,22 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      }

	    }

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

	      double value = cid->c1ao->scscm[i];

	      double value = cid->c1->scscm[i];

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = real(cid->c1ao->scscpm[i]);

	      value = real(cid->c1->scscpm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = imag(cid->c1ao->scscpm[i]);

	      value = imag(cid->c1->scscpm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = cid->c1ao->ecscm[i];

	      value = cid->c1->ecscm[i];

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = real(cid->c1ao->ecscpm[i]);

	      value = real(cid->c1->ecscpm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = imag(cid->c1ao->ecscpm[i]);

	      value = imag(cid->c1->ecscpm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	    }
@@ -1096,10 +1096,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    for (int ilr210 = 1; ilr210 <= 2; ilr210++) {

	      int ipol = (ilr210 % 2 == 0) ? 1 : -1;

	      if (ilr210 == 2) jlr = 1;

	      double extsm = cid->c1ao->ecscm[ilr210 - 1];

	      double extsm = cid->c1->ecscm[ilr210 - 1];

	      double qextm = extsm * cid->sqsfi / cid->c3->gcs;

	      double extrm = extsm / cid->c3->ecs;

	      double scasm = cid->c1ao->scscm[ilr210 - 1];

	      double scasm = cid->c1->scscm[ilr210 - 1];

	      double albdm = scasm / extsm;

	      double qscam = scasm * cid->sqsfi / cid->c3->gcs;

	      double scarm = scasm / cid->c3->scs;
@@ -1108,12 +1108,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      double absrm = abssm / cid->c3->acs;

	      double acsecs = cid->c3->acs / cid->c3->ecs;

	      if (acsecs >= -1.0e-6 && acsecs <= 1.0e-6) absrm = 1.0;

	      dcomplex s0m = cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1] * exri;

	      dcomplex s0m = cid->c1->fsacm[ilr210 - 1][ilr210 - 1] * exri;

	      double qschum = imag(s0m) * csch;

	      double pschum = real(s0m) * csch;

	      double s0magm = cabs(s0m) * cs0;

	      double rfinrm = real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas);

	      double extcrm = imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas);

	      double rfinrm = real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]) / real(cid->c3->tfsas);

	      double extcrm = imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]) / imag(cid->c3->tfsas);

	      if (inpol == 0) {

		sprintf(virtual_line, "   LIN %2d\n", ipol);

		output->append_line(virtual_line);
@@ -1136,9 +1136,9 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      output->append_line(virtual_line);

	      sprintf(

		      virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",

		      ilr210, ilr210, real(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]),

		      imag(cid->c1ao->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,

		      real(cid->c1ao->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1ao->fsacm[jlr - 1][ilr210 - 1])

		      ilr210, ilr210, real(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]),

		      imag(cid->c1->fsacm[ilr210 - 1][ilr210 - 1]), jlr, ilr210,

		      real(cid->c1->fsacm[jlr - 1][ilr210 - 1]), imag(cid->c1->fsacm[jlr - 1][ilr210 - 1])

		      );

	      output->append_line(virtual_line);

	      sprintf(
@@ -1148,8 +1148,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      output->append_line(virtual_line);

	      sprintf(virtual_line, "  QSCHU=%15.7lE, PSCHU=%15.7lE, S0MAG=%15.7lE\n", qschum, pschum, s0magm);

	      output->append_line(virtual_line);

	      double rapr = cid->c1ao->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];

	      double cosav = cid->gapm[2][ilr210 - 1] / cid->c1ao->scscm[ilr210 - 1];

	      double rapr = cid->c1->ecscm[ilr210 - 1] - cid->gapm[2][ilr210 - 1];

	      double cosav = cid->gapm[2][ilr210 - 1] / cid->c1->scscm[ilr210 - 1];

	      double fz = rapr;

	      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);

	      output->append_line(virtual_line);
@@ -1161,8 +1161,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
		output->append_line(virtual_line);

	      }

	    } // ilr210 loop

	    double rmbrif = (real(cid->c1ao->fsacm[0][0]) - real(cid->c1ao->fsacm[1][1])) / real(cid->c1ao->fsacm[0][0]);

	    double rmdchr = (imag(cid->c1ao->fsacm[0][0]) - imag(cid->c1ao->fsacm[1][1])) / imag(cid->c1ao->fsacm[0][0]);

	    double rmbrif = (real(cid->c1->fsacm[0][0]) - real(cid->c1->fsacm[1][1])) / real(cid->c1->fsacm[0][0]);

	    double rmdchr = (imag(cid->c1->fsacm[0][0]) - imag(cid->c1->fsacm[1][1])) / imag(cid->c1->fsacm[0][0]);

	    sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rmbrif);

	    output->append_line(virtual_line);

	    sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rmdchr);
@@ -1197,8 +1197,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    output->append_line(virtual_line);

	  }

	  // label 224

	  scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c1ao, cid->c3, cid->c4);

	  if (cid->c4->li != cid->c4->le) {

	  scr2(cid->vk, vkarg, exri, cid->duk, cid->c1, cid->c3);

	  if (cid->c1->li != cid->c1->le) {

	    sprintf(virtual_line, "     SPHERES; LMX=MIN0(LI,LE)\n");

	    output->append_line(virtual_line);

	  }
@@ -1256,8 +1256,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	  output->append_line(virtual_line);

	  sprintf(virtual_line, "     CLUSTER\n");

	  output->append_line(virtual_line);

	  pcros(cid->vk, exri, cid->c1, cid->c1ao, cid->c4);

	  mextc(cid->vk, exri, cid->c1ao->fsac, cid->cextlr, cid->cext);

	  pcros(cid->vk, exri, cid->c1);

	  mextc(cid->vk, exri, cid->c1->fsac, cid->cextlr, cid->cext);

	  mmulc(cid->c1->vint, cid->cmullr, cid->cmul);

	  if (jw != 0) {

	    jw = 0;
@@ -1270,22 +1270,22 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      }

	    }

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

	      double value = cid->c1ao->scsc[i];

	      double value = cid->c1->scsc[i];

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = real(cid->c1ao->scscp[i]);

	      value = real(cid->c1->scscp[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = imag(cid->c1ao->scscp[i]);

	      value = imag(cid->c1->scscp[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = cid->c1ao->ecsc[i];

	      value = cid->c1->ecsc[i];

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = real(cid->c1ao->ecscp[i]);

	      value = real(cid->c1->ecscp[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = imag(cid->c1ao->ecscp[i]);

	      value = imag(cid->c1->ecscp[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	    }
@@ -1360,10 +1360,10 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	  for (int ilr290 = 1; ilr290 <= 2; ilr290++) {

	    int ipol = (ilr290 % 2 == 0) ? 1 : -1;

	    if (ilr290 == 2) jlr = 1;

	    double extsec = cid->c1ao->ecsc[ilr290 - 1];

	    double extsec = cid->c1->ecsc[ilr290 - 1];

	    double qext = extsec * cid->sqsfi / cid->c3->gcs;

	    double extrat = extsec / cid->c3->ecs;

	    double scasec = cid->c1ao->scsc[ilr290 - 1];

	    double scasec = cid->c1->scsc[ilr290 - 1];

	    double albedc = scasec / extsec;

	    double qsca = scasec * cid->sqsfi / cid->c3->gcs;

	    double scarat = scasec / cid->c3->scs;
@@ -1372,12 +1372,12 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    double absrat = 1.0;

	    double ratio = cid->c3->acs / cid->c3->ecs;

	    if (ratio < -1.0e-6 || ratio > 1.0e-6) absrat = abssec / cid->c3->acs;

	    s0 = cid->c1ao->fsac[ilr290 - 1][ilr290 - 1] * exri;

	    s0 = cid->c1->fsac[ilr290 - 1][ilr290 - 1] * exri;

	    double qschu = imag(s0) * csch;

	    double pschu = real(s0) * csch;

	    s0mag = cabs(s0) * cs0;

	    double refinr = real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);

	    double extcor = imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);

	    double refinr = real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / real(cid->c3->tfsas);

	    double extcor = imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]) / imag(cid->c3->tfsas);

	    if (inpol == 0) {

	      sprintf(virtual_line, "   LIN %2d\n", ipol);

	      output->append_line(virtual_line);
@@ -1409,14 +1409,14 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    output->append_line(virtual_line);

	    sprintf(

		    virtual_line, "  FSAC(%1d,%1d)=%15.7lE%15.7lE   FSAC(%1d,%1d)=%15.7lE%15.7lE\n",

		    ilr290, ilr290, real(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->fsac[ilr290 - 1][ilr290 - 1]),

		    jlr, ilr290, real(cid->c1ao->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->fsac[jlr - 1][ilr290 - 1])

		    ilr290, ilr290, real(cid->c1->fsac[ilr290 - 1][ilr290 - 1]), imag(cid->c1->fsac[ilr290 - 1][ilr290 - 1]),

		    jlr, ilr290, real(cid->c1->fsac[jlr - 1][ilr290 - 1]), imag(cid->c1->fsac[jlr - 1][ilr290 - 1])

		    );

	    output->append_line(virtual_line);

	    sprintf(

		    virtual_line, "   SAC(%1d,%1d)=%15.7lE%15.7lE    SAC(%1d,%1d)=%15.7lE%15.7lE\n",

		    ilr290, ilr290, real(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1ao->sac[ilr290 - 1][ilr290 - 1]),

		    jlr, ilr290, real(cid->c1ao->sac[jlr - 1][ilr290 - 1]), imag(cid->c1ao->sac[jlr - 1][ilr290 - 1])

		    ilr290, ilr290, real(cid->c1->sac[ilr290 - 1][ilr290 - 1]), imag(cid->c1->sac[ilr290 - 1][ilr290 - 1]),

		    jlr, ilr290, real(cid->c1->sac[jlr - 1][ilr290 - 1]), imag(cid->c1->sac[jlr - 1][ilr290 - 1])

		    );

	    output->append_line(virtual_line);

	    sprintf(
@@ -1441,8 +1441,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      cid->gapv[0] = cid->gap[0][ilr290 - 1];

	      cid->gapv[1] = cid->gap[1][ilr290 - 1];

	      cid->gapv[2] = cid->gap[2][ilr290 - 1];

	      double extins = cid->c1ao->ecsc[ilr290 - 1];

	      double scatts = cid->c1ao->scsc[ilr290 - 1];

	      double extins = cid->c1->ecsc[ilr290 - 1];

	      double scatts = cid->c1->scsc[ilr290 - 1];

	      double rapr, cosav, fp, fn, fk, fx, fy, fz;

	      rftr(cid->u, cid->up, cid->un, cid->gapv, extins, scatts, rapr, cosav, fp, fn, fk, fx, fy, fz);

	      sprintf(virtual_line, "  COSAV=%15.7lE, RAPRS=%15.7lE\n", cosav, rapr);
@@ -1475,8 +1475,8 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	      output->append_line(virtual_line);

	    }

	  } //ilr290 loop

	  double rbirif = (real(cid->c1ao->fsac[0][0]) - real(cid->c1ao->fsac[1][1])) / real(cid->c1ao->fsac[0][0]);

	  double rdichr = (imag(cid->c1ao->fsac[0][0]) - imag(cid->c1ao->fsac[1][1])) / imag(cid->c1ao->fsac[0][0]);

	  double rbirif = (real(cid->c1->fsac[0][0]) - real(cid->c1->fsac[1][1])) / real(cid->c1->fsac[0][0]);

	  double rdichr = (imag(cid->c1->fsac[0][0]) - imag(cid->c1->fsac[1][1])) / imag(cid->c1->fsac[0][0]);

	  sprintf(virtual_line, "  (RE(FSAC(1,1))-RE(FSAC(2,2)))/RE(FSAC(1,1))=%15.7lE\n", rbirif);

	  output->append_line(virtual_line);

	  sprintf(virtual_line, "  (IM(FSAC(1,1))-IM(FSAC(2,2)))/IM(FSAC(1,1))=%15.7lE\n", rdichr);
@@ -1500,14 +1500,14 @@ int cluster_jxi488_cycle(int jxi488, ScattererConfiguration *sconf, GeometryConf 
	    output->append_line(virtual_line);

	  }

	  if (iavm != 0) {

	    mmulc(cid->c1ao->vintm, cid->cmullr, cid->cmul);

	    mmulc(cid->c1->vintm, cid->cmullr, cid->cmul);

	    // Some implicit loops writing to binary.

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

	      double value;

	      value = real(cid->c1ao->vintm[i]);

	      value = real(cid->c1->vintm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	      value = imag(cid->c1ao->vintm[i]);

	      value = imag(cid->c1->vintm[i]);

	      // tppoan.write(reinterpret_cast<char *>(&value), sizeof(double));

	      vtppoanp->append_line(VirtualBinaryLine(value));

	    }

src/include/Commons.h

+86 −5
Original line number Diff line number Diff line
@@ -37,6 +37,7 @@ 
#include <mpi.h>

#endif



class ParticleDescriptor;

class mixMPI;



/*! \brief Representation of the FORTRAN C1 common blocks.
@@ -551,15 +552,11 @@ public: 
class ClusterIterationData {

public:

  //! \brief Pointer to a C1 structure.

  C1 *c1;

  //! \brief Pointer to a C1_AddOns structure.

  C1_AddOns *c1ao;

  ParticleDescriptor *c1;

  //! \brief Pointer to a C2 structure.

  C2 *c2;

  //! \brief Pointer to a C3 structure.

  C3 *c3;

  //! brief Pointer to a C4 structure.

  C4 *c4;

  //! \brief Pointer to a C6 structure.

  C6 *c6;

  //! \brief Pointer to a C9 structure.
@@ -698,6 +695,8 @@ protected: 
  int _litpo;

  //! \brief LITPOS = LITPO * LITPO

  int _litpos;

  //! \brief LMPO = LM + 1

  int _lmpo;

  //! \brief LMTPO = LI + LE + 1

  int _lmtpo;

  //! \brief LMTPOS = LMTPO * LMTPO
@@ -805,6 +804,8 @@ public: 
  const int& litpo = _litpo;

  //! \brief Read-only view of LITPOS.

  const int& litpos = _litpos;

  //! \brief Read-only view of LMPO.

  const int& lmpo = _lmpo;

  //! \brief Read-only view of LMTPO.

  const int& lmtpo = _lmtpo;

  //! \brief Read-only view of LMTPOS.
@@ -880,6 +881,26 @@ public: 
   */

  ParticleDescriptor(GeometryConfiguration *gconf, ScattererConfiguration *sconf);



  /*! \brief ParticleDescriptor copy constructor.

   *

   * \param rhs: `const ParticleDescriptor &` Reference to ParticleDescriptor object to be copied.

   */

  ParticleDescriptor(const ParticleDescriptor &rhs);



#ifdef MPI_VERSION

  /*! \brief ParticleDescriptor MPI constructor.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  ParticleDescriptor(const mixMPI *mpidata);

  

  /*! \brief Broadcast ParticleDescriptor instance over MPI.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  void mpibcast(const mixMPI *mpidata);

#endif // MPI_VERSION

  

  /*! \brief ParticleDescriptor instance destroyer.

   */

  ~ParticleDescriptor();
@@ -906,6 +927,26 @@ public: 
   */

  ParticleDescriptorCluster(GeometryConfiguration *gconf, ScattererConfiguration *sconf);



  /*! \brief ParticleDescriptorCluster copy constructor.

   *

   * \param rhs: `const ParticleDescriptorCluster &` Reference to ParticleDescriptorCluster object to be copied.

   */

  ParticleDescriptorCluster(const ParticleDescriptorCluster &rhs);



#ifdef MPI_VERSION

  /*! \brief ParticleDescriptorCluster MPI constructor.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  ParticleDescriptorCluster(const mixMPI *mpidata);

  

  /*! \brief Broadcast ParticleDescriptorCluster instance over MPI.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  void mpibcast(const mixMPI *mpidata);

#endif // MPI_VERSION

  

  /*! \brief Interface function to return the descriptor type as string.

   *

   * \return descriptor_type: `string` The descriptor type name.
@@ -928,6 +969,26 @@ public: 
   */

  ParticleDescriptorInclusion(GeometryConfiguration *gconf, ScattererConfiguration *sconf);



  /*! \brief ParticleDescriptorInclusion copy constructor.

   *

   * \param rhs: `const ParticleDescriptorInclusion &` Reference to ParticleDescriptorInclusion object to be copied.

   */

  ParticleDescriptorInclusion(const ParticleDescriptorInclusion &rhs);



#ifdef MPI_VERSION

  /*! \brief ParticleDescriptorInclusion MPI constructor.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  ParticleDescriptorInclusion(const mixMPI *mpidata);

  

  /*! \brief Broadcast ParticleDescriptorInclusion instance over MPI.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  void mpibcast(const mixMPI *mpidata);

#endif // MPI_VERSION

  

  /*! \brief Interface function to return the descriptor type as string.

   *

   * \return descriptor_type: `string` The descriptor type name.
@@ -950,6 +1011,26 @@ public: 
   */

  ParticleDescriptorSphere(GeometryConfiguration *gconf, ScattererConfiguration *sconf);



  /*! \brief ParticleDescriptorSphere copy constructor.

   *

   * \param rhs: `const ParticleDescriptorSphere &` Reference to ParticleDescriptorSphere object to be copied.

   */

  ParticleDescriptorSphere(const ParticleDescriptorSphere &rhs);



#ifdef MPI_VERSION

  /*! \brief ParticleDescriptorSphere MPI constructor.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  ParticleDescriptorSphere(const mixMPI *mpidata);

  

  /*! \brief Broadcast ParticleDescriptorSphere instance over MPI.

   *

   * \param mpidata: `const mixMPI *` Pointer to mixMPI instance.

   */

  void mpibcast(const mixMPI *mpidata);

#endif // MPI_VERSION

  

  /*! \brief Interface function to return the descriptor type as string.

   *

   * \return descriptor_type: `string` The descriptor type name.

src/include/clu_subs.h

+24 −46

File changed.

Preview size limit exceeded, changes collapsed.

src/libnptm/Commons.cpp

+768 −97

File changed.

Preview size limit exceeded, changes collapsed.

src/libnptm/clu_subs.cpp

+166 −166

File changed.

Preview size limit exceeded, changes collapsed.