Merge pull request #44 from EXP-code/SignConvention

Code to enforce sign convention for basis creation
EXP-code · Sep 26, 2023 · 9a12d9b · 9a12d9b
2 parents 2bb4fa2 + 83a4ed2
commit 9a12d9b
Show file tree

Hide file tree

Showing 6 changed files with 96 additions and 9 deletions.
diff --git a/exputil/EmpCylSL.cc b/exputil/EmpCylSL.cc
@@ -3336,6 +3336,20 @@ void EmpCylSL::eigen_problem(int request_id, int M, Timer& timer)
       efE = esE.eigenvectors().real();
       efO = esO.eigenvectors().real();
 
+      // Choose sign conventions for the ef table
+      //
+      int nfid = std::min<int>(4, efE.rows()) - 1;
+
+      // Even case
+      for (int j=0; j<efE.cols(); j++) {
+	if (efE(nfid, j) < 0.0) efE.col(j) *= -1;
+      }
+
+      // Odd case
+      for (int j=0; j<efO.cols(); j++) {
+	if (efO(nfid, j) < 0.0) efO.col(j) *= -1;
+      }
+
       if (VFLAG & 32) {
 
 	std::ostringstream sout;
@@ -3381,6 +3395,13 @@ void EmpCylSL::eigen_problem(int request_id, int M, Timer& timer)
       ev = es.eigenvalues().real();
       ef = es.eigenvectors().real();
 
+      // Sign convention
+      //
+      int nfid = std::min<int>(4, efE.rows()) - 1;
+      for (int j=0; j<ef.cols(); j++) {
+	if (ef(nfid, j) < 0.0) ef.col(j) *= -1;
+      }
+
       if (VFLAG & 32) {
 
 	std::ostringstream sout;
@@ -3554,6 +3575,23 @@ void EmpCylSL::eigen_problem(int request_id, int M, Timer& timer)
       efE = esE.eigenvectors().real();
       efO = esO.eigenvectors().real();
 
+
+      // Sign convention
+      //
+      int nfid = std::min<int>(4, efE.rows()) - 1;
+
+      // Even case
+      //
+      for (int j=0; j<efE.cols(); j++) {
+	if (efE(nfid, j) < 0.0) efE.col(j) *= -1;
+      }
+
+      // Odd case
+      //
+      for (int j=0; j<efO.cols(); j++) {
+	if (efO(nfid, j) < 0.0) efO.col(j) *= -1;
+      }
+
       if (VFLAG & 32) {
 
 	std::ostringstream sout;
@@ -3586,6 +3624,13 @@ void EmpCylSL::eigen_problem(int request_id, int M, Timer& timer)
       ev = es.eigenvalues().real();
       ef = es.eigenvectors().real();
 
+      // Sign convention
+      //
+      int nfid = std::min<int>(4, efE.rows()) - 1;
+      for (int j=0; j<ef.cols(); j++) {
+	if (ef(nfid, j) < 0.0) ef.col(j) *= -1;
+      }
+
       if (VFLAG & 32) {
 
 	std::ostringstream sout;
@@ -4424,6 +4469,14 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp)
 
 	(*pb)[mm]->evalJK = es.eigenvalues().real();
 	(*pb)[mm]->evecJK = es.eigenvectors().real();
+
+	// Sign convention
+	//
+	int nfid = std::min<int>(4, efE.rows()) - 1;
+	for (int j=0; j<(*pb)[mm]->evecJK.cols(); j++) {
+	  if ((*pb)[mm]->evecJK(nfid, j) < 0.0) (*pb)[mm]->evecJK.col(j) *= -1;
+	}
+
       }
 
       // Transformation output
@@ -4479,6 +4532,14 @@ void EmpCylSL::pca_hall(bool compute, bool subsamp)
 	  evalVar = es.eigenvalues().real();
 	  evecVar = es.eigenvectors().real();
 
+	  // Sign convention for eigenvectors
+	  //
+	  int nfid = std::min<int>(4, evecVar.rows()) - 1;
+	  for (int j=0; j<evecVar.cols(); j++) {
+	    if (evecVar(nfid, j) < 0.0) evecVar.col(j) *= -1;
+	  }
+
+
 	  mout << "# EOF eigenvalues" << std::endl;
 	  double total = 0.0;
 	  for (int nn=0; nn<rank3; nn++) {

diff --git a/exputil/SLGridMP2.cc b/exputil/SLGridMP2.cc
@@ -2679,9 +2679,17 @@ void SLGridSph::compute_table(struct TableSph* table, int l)
 
   table->ef.resize(N, numr);
 
+  // Choose sign conventions for the ef table
+  //
+  int nfid = std::min<int>(nevsign, N) - 1;
+  Eigen::VectorXi sgn = Eigen::VectorXi::Ones(N);
+  for (int j=0; j<N; j++) {
+    if (ef[j*NUM+nfid]<0.0) sgn(j) = -1;
+  }
+
   for (int i=0; i<numr; i++) {
     for(int j=0; j<N; j++) 
-      table->ef(j, i) = ef[j*NUM+i];
+      table->ef(j, i) = ef[j*NUM+i] * sgn(j);
   }
 
   table->l = l;
@@ -2862,15 +2870,24 @@ void SLGridSph::compute_table_worker(void)
 	}
       }
     }
-				// Load table
 
+    // Load table
+    //
     table.ev.resize(N);
     for (int i=0; i<N; i++) table.ev[i] = ev[i];
 
+    // Choose sign conventions for the ef table
+    //
+    int nfid = std::min<int>(nevsign, N) - 1;
+    Eigen::VectorXi sgn = Eigen::VectorXi::Ones(N);
+    for (int j=0; j<N; j++) {
+      if (ef[j*NUM+nfid]<0.0) sgn(j) = -1;
+    }
+
     table.ef.resize(N, numr);
     for (int i=0; i<numr; i++) {
       for (int j=0; j<N; j++) 
-	table.ef(j, i) = ef[j*NUM+i];
+	table.ef(j, i) = ef[j*NUM+i] * sgn(j);
     }
 
     table.l = L;

diff --git a/exputil/libvars.cc b/exputil/libvars.cc
@@ -30,6 +30,10 @@ namespace __EXP__
 
   //! Source/line info exception strings
   bool             sourceline       = false;
+
+  //! Sign convention element rank for eigenfunctions and eigenvectors
+  int              nevsign          = 4;
+
 };
 
 
diff --git a/include/SLGridMP2.H b/include/SLGridMP2.H
@@ -15,6 +15,9 @@
 #include <sltableMP2.H>
 #include <yaml-cpp/yaml.h>
 
+#include <libvars.H>
+using namespace __EXP__;
+
 #if HAVE_LIBCUDA==1
 #include <cudaUtil.cuH>
 #include <cudaMappingConstants.cuH>

diff --git a/include/libvars.H b/include/libvars.H
@@ -36,6 +36,9 @@ namespace __EXP__
   //! Control flag for source/line debug info
   extern bool sourceline;
 
+  //! Sign convention element rank for eigenfunctions and eigenvectors
+  extern int nevsign;
+
 };
 
 #endif	// END _LIBVARS_H
diff --git a/src/expand.cc b/src/expand.cc
@@ -636,12 +636,9 @@ main(int argc, char** argv)
     // Initialize slurm api
     //=====================
 #ifdef HAVE_LIBSLURM
-    {
-      int rc = slurm_init(0);
-      if (rc != SLURM_SUCCESS)
-	std::cerr << "EXP [" << myid << "]: error initializing Slurm API: "
-		  << slurm_strerror(rc) << std::endl;
-    }
+#if SLURM_VERSION_NUMBER > SLURM_VERSION_NUM(21,10,0)
+    slurm_init(0);
+#endif
 #endif
 
     //==============================================
@@ -795,7 +792,9 @@ main(int argc, char** argv)
   //======================
 
 #ifdef HAVE_LIBSLURM
+#if SLURM_VERSION_NUMBER > SLURM_VERSION_NUM(21,10,0)
   slurm_fini();
+#endif
 #endif
 
   //=================