Commit 318ba1da authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Add inclu_subs to Makefile and C++ porting of CNF()

parent 0b28e145

build/Makefile.am

+1 −1
Original line number Diff line number Diff line 
LDADD=libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}

lib_LTLIBRARIES=libnptm/libnptm.la

libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

libnptm_libnptm_la_SOURCES=../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/inclu_subs.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

if BUILDFORTRAN

PROGS=cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu inclusion/np_inclusion sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_ParticleDescriptor testing/test_TEDF testing/test_TTMS

bin_PROGRAMS=$(PROGS)

build/Makefile.in

+12 −6
Original line number Diff line number Diff line
@@ -160,11 +160,11 @@ am__dirstamp = $(am__leading_dot)dirstamp 
am_libnptm_libnptm_la_OBJECTS = ../src/libnptm/algebraic.lo \

	../src/libnptm/clu_subs.lo ../src/libnptm/Commons.lo \

	../src/libnptm/Configuration.lo ../src/libnptm/file_io.lo \

	../src/libnptm/lapack_calls.lo ../src/libnptm/logging.lo \

	../src/libnptm/magma_calls.lo ../src/libnptm/Parsers.lo \

	../src/libnptm/sph_subs.lo ../src/libnptm/utils.lo \

	../src/libnptm/tfrfme.lo ../src/libnptm/TransitionMatrix.lo \

	../src/libnptm/tra_subs.lo

	../src/libnptm/inclu_subs.lo ../src/libnptm/lapack_calls.lo \

	../src/libnptm/logging.lo ../src/libnptm/magma_calls.lo \

	../src/libnptm/Parsers.lo ../src/libnptm/sph_subs.lo \

	../src/libnptm/utils.lo ../src/libnptm/tfrfme.lo \

	../src/libnptm/TransitionMatrix.lo ../src/libnptm/tra_subs.lo

libnptm_libnptm_la_OBJECTS = $(am_libnptm_libnptm_la_OBJECTS)

AM_V_lt = $(am__v_lt_@AM_V@)

am__v_lt_ = $(am__v_lt_@AM_DEFAULT_V@)
@@ -348,6 +348,7 @@ am__depfiles_remade = ../src/cluster/$(DEPDIR)/cluster.Po \ 
	../src/libnptm/$(DEPDIR)/algebraic.Plo \

	../src/libnptm/$(DEPDIR)/clu_subs.Plo \

	../src/libnptm/$(DEPDIR)/file_io.Plo \

	../src/libnptm/$(DEPDIR)/inclu_subs.Plo \

	../src/libnptm/$(DEPDIR)/lapack_calls.Plo \

	../src/libnptm/$(DEPDIR)/logging.Plo \

	../src/libnptm/$(DEPDIR)/magma_calls.Plo \
@@ -637,7 +638,7 @@ top_builddir = @top_builddir@ 
top_srcdir = @top_srcdir@

LDADD = libnptm/libnptm.la -L/usr/lib64 ${USER_LDFLAGS} ${HDF5_LDFLAGS} ${LAPACKLDFLAGS} ${MAGMALDFLAGS}

lib_LTLIBRARIES = libnptm/libnptm.la

libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

libnptm_libnptm_la_SOURCES = ../src/libnptm/algebraic.cpp ../src/libnptm/clu_subs.cpp ../src/libnptm/Commons.cpp ../src/libnptm/Configuration.cpp ../src/libnptm/file_io.cpp ../src/libnptm/inclu_subs.cpp ../src/libnptm/lapack_calls.cpp ../src/libnptm/logging.cpp ../src/libnptm/magma_calls.cpp ../src/libnptm/Parsers.cpp ../src/libnptm/sph_subs.cpp ../src/libnptm/utils.cpp ../src/libnptm/tfrfme.cpp ../src/libnptm/TransitionMatrix.cpp ../src/libnptm/tra_subs.cpp

@BUILDFORTRAN_FALSE@PROGS = cluster/np_cluster inclusion/np_inclusion sphere/np_sphere trapping/np_trapping testing/test_ParticleDescriptor testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_TRUE@PROGS = cluster/edfb_clu cluster/clu cluster/np_cluster inclusion/edfb_inclu inclusion/inclu inclusion/np_inclusion sphere/edfb_sph sphere/sph sphere/np_sphere trapping/frfme trapping/lffft trapping/np_trapping testing/test_ParticleDescriptor testing/test_TEDF testing/test_TTMS

@BUILDFORTRAN_TRUE@EDFBCLUSOURCES = ../src/cluster/edfb_clu.f
@@ -807,6 +808,8 @@ clean-libLTLIBRARIES: 
	../src/libnptm/$(DEPDIR)/$(am__dirstamp)

../src/libnptm/file_io.lo: ../src/libnptm/$(am__dirstamp) \

	../src/libnptm/$(DEPDIR)/$(am__dirstamp)

../src/libnptm/inclu_subs.lo: ../src/libnptm/$(am__dirstamp) \

	../src/libnptm/$(DEPDIR)/$(am__dirstamp)

../src/libnptm/lapack_calls.lo: ../src/libnptm/$(am__dirstamp) \

	../src/libnptm/$(DEPDIR)/$(am__dirstamp)

../src/libnptm/logging.lo: ../src/libnptm/$(am__dirstamp) \
@@ -1006,6 +1009,7 @@ distclean-compile: 
@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/algebraic.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/clu_subs.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/file_io.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/inclu_subs.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/lapack_calls.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/logging.Plo@am__quote@ # am--include-marker

@AMDEP_TRUE@@am__include@ @am__quote@../src/libnptm/$(DEPDIR)/magma_calls.Plo@am__quote@ # am--include-marker
@@ -1382,6 +1386,7 @@ distclean: distclean-am 
	-rm -f ../src/libnptm/$(DEPDIR)/algebraic.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/clu_subs.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/file_io.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/inclu_subs.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/lapack_calls.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/logging.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/magma_calls.Plo
@@ -1455,6 +1460,7 @@ maintainer-clean: maintainer-clean-am 
	-rm -f ../src/libnptm/$(DEPDIR)/algebraic.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/clu_subs.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/file_io.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/inclu_subs.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/lapack_calls.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/logging.Plo

	-rm -f ../src/libnptm/$(DEPDIR)/magma_calls.Plo

src/include/inclu_subs.h

+12 −1
Original line number Diff line number Diff line
@@ -29,6 +29,16 @@ 
#ifndef INCLUDE_INCLU_SUBS_H_

#define INCLUDE_INCLU_SUBS_H_



/*! \brief C++ porting of CNF.

 *

 * \param n: `int` Bessel y function order.

 * \param z: `dcomplex` Argument of Bessel y function.

 * \param nm: `int` Maximum computed order.

 * \param csj: `dcomplex *` TBD.

 * \param csy: `dcomplex *` Complex spherical Bessel functions up to desired order.

 */

void cnf(int n, dcomplex z, int nm, dcomplex *csj, dcomplex *csy);



/*! \brief C++ porting of EXMA.

 *

 * \param am: `dcomplex **` Field transition coefficients matrix.
@@ -43,8 +53,9 @@ void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1); 
 * \param at: `dcomplex **` TBD.

 * \param enti: `double` TBD.

 * \param c1: `ParticleDescriptor *` Pointer to a ParticleDescriptor instance.

 * \param c6: `C6 *` Pointer to a C6 instance.

 */

void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1);

void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1, C6 *c6);



/*! \brief C++ porting of INDME.

 *

src/libnptm/inclu_subs.cpp

+70 −26
Original line number Diff line number Diff line
@@ -31,30 +31,71 @@ 
#include "../include/Commons.h"

#endif



#ifndef INCLUDE_SPH_SUBS_H_

#include "../include/sph_subs.h"

#endif



#ifndef INCLUDE_CLU_SUBS_H_

#include "../include/clu_subs.h"

#endif



using namespace std;



void cnf(int n, dcomplex z, int nm, dcomplex *csj, dcomplex *csy) {

  /*

  FROM CSPHJY OF LIBRARY specfun 

    

  ==========================================================

    Purpose: Compute spherical Bessel functions y

    Input :  Z --- Complex argument of y

             N --- Order of y ( N = 0,1,2,... )

             CSJ(N+1) --- j

             NM --- Highest order computed

    Output:  CSY(N+1) --- y

  ==========================================================

  */

  double a0 = cabs(z);

  if (a0 < 1.0e-60) {

    for (int k = 0; k <= n; k++) csy[k] = -1.0e300;

  } else {

    csy[0] = -ccos(z) / z;

    if (n > 0) {

      csy[1] = (csy[0] - csin(z)) / z;

      if (n > 1) {

	for (int k = 2; k <= nm; k++) {

	  double absjk = cabs(csj[k - 1]);

	  double absjkmo = cabs(csj[k - 2]);

	  csy[k] = (absjk > absjkmo) ?

	    (csj[k] * csy[k - 1] - 1.0 / (z * z)) / csj[k - 1] :

	    (csj[k] * csy[k - 2] - (2.0 * k - 1.0) / (z * z * z)) / csj[k - 2];

	}

      }

    }

  }

}



void exma(dcomplex **am, dcomplex **at, ParticleDescriptor *c1) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  const ndit = c1->nsph * c1->nlim * 2;

  const int ndit = c1->nsph * c1->nlim * 2;

  const int ndm = ndit + c1->nlemt;

  for (int j20 = 1; j20 <= c1->nlemt; j20++) {

    int j0 = j20 + ndit;

    for (int i20 = 1; i20 <= c1->nlemt; i20++) {

      dcomplex sum = cc0;

      for (int k = 0; k < ndm; k++) sum += at[i20 - 1][k] * am[k][j0 - 1];

    }

      c1->am0m[i20 - 1][j20 - 1] = sum;

    }

  } // j20 loop

}



void incms(dcomplex **am, dcomplex **at, double enti) {

void incms(dcomplex **am, dcomplex **at, double enti, ParticleDescriptor *c1, C6 *c6) {

  const dcomplex cc0 = 0.0 + I * 0.0;

  int nbl;

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

  const int ndit = ndi + ndi;

  const int ndm = ndit + c1->nlemt;

  nbl = 0;

  for (int n1 = 1, n1 < c1->nsph; n1++) {

  for (int n1 = 1; n1 < c1->nsph; n1++) {

    int in1 = (n1 - 1) * c1->nlim;

    int n1po = n1 + 1;

    for (int n2 = n1po; n2 <= c1->nsph; n2++) {
@@ -86,8 +127,8 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
	      int i2e = in2 + ilm2e;

	      int j2 = in1 + ilm2;

	      int j2e = in1 + ilm2e;

	      dcomplex cgh = ghit(0, 0, nbl, l1, m1, l2, m2);

	      dcomplex cgk = ghit(0, 1, nbl, l1, m1, l2, m2);

	      dcomplex cgh = ghit(0, 0, nbl, l1, m1, l2, m2, c1, c6);

	      dcomplex cgk = ghit(0, 1, nbl, l1, m1, l2, m2, c1, c6);

	      am[i1 - 1][i2 - 1] = cgh;

	      am[i1 - 1][i2e - 1] = cgk;

	      am[i1e - 1][i2 - 1] = cgk;
@@ -113,10 +154,11 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
      for (int im1 = 1; im1 <= l1tpo; im1++) {

	int m1 = im1 - l1po;

	int ilm1 = il1 + m1;

	int i1 = in1 + ilm1;

	int i1e = i1 + ndi;

	for (int ilm2 = 1; ilm2 <= c1->nlim; ilm2++) {

	  int i2 = in1 + ilm2;

	  i2e = i2 + ndi;

	  int i2e = i2 + ndi;

	  am[i1 - 1][i2 - 1] = cc0;

	  am[i1 - 1][i2e - 1] = cc0;

	  am[i1e - 1][i2 - 1] = cc0;
@@ -189,8 +231,8 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
	      int j1 = in1 + jlm1;

	      int j1e = j1 + ndi;

	      int isil = ((m2 + m1) % 2 == 0) ? 1 : -1;

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2);

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2, c1, c6);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2, c1, c6);

	      am[i1 - 1][i3 - 1] = cgi;

	      am[i1 - 1][i3e - 1] = cgl;

	      am[i1e - 1][i3 - 1] = cgl;
@@ -224,18 +266,18 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
	i2++;

	int i2e = i2 + c1->nlem;

	int i3 = i2 + ndit;

	int i3e = i3 + nlem;

	int i3e = i3 + c1->nlem;

	int i1 = 0;

	for (int n1 = 1; n1 <= c1->nsph; n1++) {

	  for (int l1 = 1; l1 <= c1->li; l1++) {

	    l1tpo = l1 + l1 + 1;

	    int l1tpo = l1 + l1 + 1;

	    int m1 = -l1 - 1;

	    for (int im1 = 1; im1 <= l1tpo; im1++) {

	      m1++;

	      i1++;

	      int i1e = i1 + ndi;

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2);

	      dcomplex cgi = ghit(2, 0, n1, l1, m1, l2, m2, c1, c6);

	      dcomplex cgl = ghit(2, 1, n1, l1, m1, l2, m2, c1, c6);

	      am[i1 - 1][i3 - 1] = cgi;

	      am[i1 - 1][i3e - 1] = cgl;

	      am[i1e - 1][i3 - 1] = cgl;
@@ -249,7 +291,7 @@ void incms(dcomplex **am, dcomplex **at, double enti) { 
	      at[i2 - 1][i1 - 1] = cgi * ftm;

	      at[i2 - 1][i1e - 1] = cgl * ftm;

	      at[i2e - 1][i1 - 1] = cgl * fte;

	      at[i2e - 1][i1e - 1] cgi * fte;

	      at[i2e - 1][i1e - 1] = cgi * fte;

	    } // im1 loop 60

	  } // l1 loop 60

	} // n1 loop 60
@@ -279,8 +321,8 @@ void indme( 
  double sz = vk * c1->ros[i - 1];

  c2->vsz[i - 1] = sz;

  double vkr1 = vk * c1->rc[i - 1][0];

  int nsh = c1->nshl[i . 1];

  c2->vkt[i - 1] = sqrt(c2->dc0[0]);

  int nsh = c1->nshl[i - 1];

  c2->vkt[i - 1] = csqrt(c2->dc0[0]);

  arg = vkr1 * c2->vkt[i - 1];

  dcomplex arin = arg;

  if (imag(arg) != 0.0) {
@@ -374,6 +416,7 @@ void indme( 
  dcomplex *drmf = new dcomplex[c1->li]();

  dcomplex *ref = new dcomplex[c1->li]();

  dcomplex *dref = new dcomplex[c1->li]();

  int ic = 0;

  if (nsh < 2) { // nsh == 1

    dcomplex cri = c2->dc0[0] / ent;

    for (int l160 = 1; l160 <= c1->li; l160++) {
@@ -387,8 +430,8 @@ void indme( 
      dcomplex ccnb = fn[lpo - 1] * dfbi;

      dcomplex ccnc = fbi[lpo - 1] * dfb;

      dcomplex ccnd = fb[lpo - 1] * dfbi;

      rmi[l60 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);

      rei[l60 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);

      c1->rmi[l160 - 1][i - 1] = 1.0 + uim * (ccna - ccnb) / (ccnc - ccnd);

      c1->rei[l160 - 1][i - 1] = 1.0 + uim * (cri * ccna - ccnb) / (cri * ccnc - ccnd);

    } // l160 loop

  } else { // label 165: nsh > 1

    for (int l180 = 1; l180 <= c1->li; l180++) {
@@ -400,18 +443,19 @@ void indme( 
      dcomplex dy1 = (l180 * fbi[l180 - 1] - lpo * fbi[lpt - 1]) * c2->vkt[i - 1] / dltpo;

      dcomplex y2 = y1;

      dcomplex dy2 = dy1;

      int ic = 1;

      ic = 0;

      for (int ns = 2; ns <= nsh; ns++) {

	int nsmo = ns - 1;

	double vkr = vk * c1->rc[i - 1][nsmo - 1];

	if {ns % 2 != 0) {

	if (ns % 2 != 0) {

	  // ic is incremented before being read in this loop.

	  int step = vk * (c1->rc[i - 1][ns - 1] - c1->rc[i - 1][nsmo - 1]) / nstp;

	  arg = c2->dc0[ic++];

	  arg = c2->dc0[++ic];

	  rkc(nstp, step, arg, vkr, lpo, y1, y2, dy1, dy2);

	} else { // label 170

	  diel(nstpts, nsmo, i, ic, vk);

	  diel(nstpts, nsmo, i, ic, vk, c1, c2);

	  double stepts = vk * (c1->rc[i - 1][ns - 1] - c1->rc[i - 1][nsmo - 1]) / nstpts;

	  rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2);

	  rkt(nstpts, stepts, vkr, lpo, y1, y2, dy1, dy2, c2);

	}

      } // ns loop 176

      rmf[l180 - 1] = y1 * sz;
@@ -465,10 +509,10 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6) { 
    }

  } // i18 loop

  int i = 0;

  for (int l1po28 = 1; l1po28 <= c1->lmpo; l1po28++) {

  for (int l1po = 1; l1po <= c1->lmpo; l1po++) {

    int l1 = l1po - 1;

    for (int l2 = 1; l2 <= c1->lm; l2++) {

      r3j000(l1, l2);

      r3j000(l1, l2, c6);

      c1->ind3j[l1po - 1][l2 - 1] = i;

      int lmnpo = 1 + ((l2 - l1 > 0) ? l2 - l1 : l1 - l2);

      int lmxpo = l2 + l1 + 1;
@@ -479,7 +523,7 @@ void instr(ScattererConfiguration *sconf, ParticleDescriptor *c1, C6 *c6) { 
	lpo += 2;

      }

    } // l2 loop

  } // l1po28 loop

  } // l1po loop 28

  int nsphmo = c1->nsph - 1;

  int lit = c1->li + c1->li;

  ivy = 0;