added help/usage messages with -h ; dsyevd falls back to dsyev if con…

…vergence fails ; use smart pointers in diag.h

c++/diag.h

@@ -43,9 +43,9 @@ size_t diagonalise_dsyev(Matrix &m, Eigen &d, char jobz = 'V') {
   LAPACK_dsyev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK0, &LWORK0, &INFO);
   my_assert(INFO == 0);
   int LWORK    = int(WORK0[0]);
-  auto *WORK = new double[LWORK];
+  auto WORK = std::make_unique<double[]>(LWORK);
   // Step 2: perform the diagonalisation
-  LAPACK_dsyev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK, &LWORK, &INFO);
+  LAPACK_dsyev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK.get(), &LWORK, &INFO);
   if (INFO != 0) my_error("dsyev failed. INFO=%i", INFO);
   d = Eigen(dim, dim);
   copy_values(eigenvalues, d.value, dim);
@@ -54,7 +54,6 @@ size_t diagonalise_dsyev(Matrix &m, Eigen &d, char jobz = 'V') {
       for (size_t j = 0; j < dim; j++) d.vektor(r, j) = ham[dim * r + j];
     d.perform_checks();
   }
-  delete[] WORK;
   return dim;
 }
 #endif
@@ -78,20 +77,23 @@ size_t diagonalise_dsyevd(Matrix &m, Eigen &d, char jobz = 'V')
   my_assert(INFO == 0);
   LWORK = int(WORK0[0]);
   LIWORK = IWORK0[0];
-  auto *WORK = new double[LWORK];
-  auto *IWORK = new int[LIWORK];
-  LAPACK_dsyevd(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK, &LWORK,
-                IWORK, &LIWORK, &INFO);
-  if (INFO != 0) my_error("dsyevd failed. INFO=%i", INFO);
+  auto WORK = std::make_unique<double[]>(LWORK);
+  auto IWORK = std::make_unique<int[]>(LIWORK);
+  LAPACK_dsyevd(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK.get(), &LWORK,
+                IWORK.get(), &LIWORK, &INFO);
+  if (INFO != 0) {
+    // dsyevd sometimes fails to converge (INFO>0). In such cases we do not trigger
+    // an error but return 0, to permit error recovery.
+    if (INFO > 0) return 0;
+    my_error("dsyevd failed. INFO=%i", INFO);
+  }
   d = Eigen(dim, dim);
   copy_values(eigenvalues, d.value, dim);
   if (jobz == 'V') {
     for (size_t r = 0; r < dim; r++)
       for (size_t j = 0; j < dim; j++) d.vektor(r, j) = ham[dim * r + j];
     d.perform_checks();
   }
-  delete[] WORK;
-  delete[] IWORK;
   return dim;
 }
 #endif  
@@ -141,11 +143,11 @@ size_t diagonalise_dsyevr(Matrix &m, Eigen &d, char jobz = 'V',
   my_assert(INFO == 0);
   int LWORK    = int(WORK0[0]);
   int LIWORK   = IWORK0[0];
-  auto *WORK = new double[LWORK];
-  int *IWORK   = new int[LIWORK];
+  auto WORK = std::make_unique<double[]>(LWORK);
+  auto IWORK  = std::make_unique<int[]>(LIWORK);
   // Step 2: perform the diagonalisation
-  LAPACK_dsyevr(&jobz, &RANGE, &UPLO, &NN, ham, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &MM, (double *)eigenvalues, &Z[0], &LDZ, ISUPPZ, WORK, &LWORK,
-                IWORK, &LIWORK, &INFO);
+  LAPACK_dsyevr(&jobz, &RANGE, &UPLO, &NN, ham, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &MM, (double *)eigenvalues, &Z[0], &LDZ, ISUPPZ, WORK.get(), &LWORK,
+                IWORK.get(), &LIWORK, &INFO);
   if (INFO != 0) my_error("dsyevr failed. INFO=%i", INFO);
   if (MM != int(M)) {
     cout << "dsyevr computed " << MM << "/" << M << endl;
@@ -159,8 +161,6 @@ size_t diagonalise_dsyevr(Matrix &m, Eigen &d, char jobz = 'V',
       for (size_t j = 0; j < dim; j++) d.vektor(r, j) = Z[dim * r + j];
     d.perform_checks();
   }
-  delete[] WORK;
-  delete[] IWORK;
   return M;
 }
 #endif
@@ -182,9 +182,9 @@ size_t diagonalise_zheev(Matrix &m, Eigen &d, char jobz = 'V') {
   LAPACK_zheev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK0, &LWORK0, RWORK, &INFO);
   my_assert(INFO == 0);
   int LWORK  = int(WORK0[0].real);
-  lapack_complex_double *WORK = new lapack_complex_double[LWORK];
+  auto WORK = std::make_unique<lapack_complex_double[]>(LWORK);
   // Step 2: perform the diagonalisation
-  LAPACK_zheev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK, &LWORK, RWORK, &INFO);
+  LAPACK_zheev(&jobz, &UPLO, &NN, ham, &LDA, (double *)eigenvalues, WORK.get(), &LWORK, RWORK, &INFO);
   if (INFO != 0) my_error("zheev failed. INFO=%i", INFO);
   d = Eigen(dim, dim);
   copy_values(eigenvalues, d.value, dim);
@@ -196,7 +196,6 @@ size_t diagonalise_zheev(Matrix &m, Eigen &d, char jobz = 'V') {
       }
     d.perform_checks();
   }
-  delete[] WORK;
   return dim;
 }
 #endif
@@ -245,14 +244,14 @@ size_t diagonalise_zheevr(Matrix &m, Eigen &d, char jobz = 'V', double ratio = 1
                 RWORK0, &LRWORK0, IWORK0, &LIWORK0, &INFO);
   my_assert(INFO == 0);
   int LWORK   = int(WORK0[0].real);
-  lapack_complex_double *WORK  = new lapack_complex_double[LWORK];
+  auto WORK  = std::make_unique<lapack_complex_double[]>(LWORK);
   int LRWORK  = int(RWORK0[0]);
-  double *RWORK = new double[LRWORK];
+  auto RWORK = std::make_unique<double[]>(LRWORK);
   int LIWORK  = IWORK0[0];
-  int *IWORK  = new int[LIWORK];
+  auto IWORK  = std::make_unique<int[]>(LIWORK);
   // Step 2: perform the diagonalisation
-  LAPACK_zheevr(&jobz, &RANGE, &UPLO, &NN, ham, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &MM, (double *)eigenvalues, &Z[0], &LDZ, ISUPPZ, WORK, &LWORK,
-                RWORK, &LRWORK, IWORK, &LIWORK, &INFO);
+  LAPACK_zheevr(&jobz, &RANGE, &UPLO, &NN, ham, &LDA, &VL, &VU, &IL, &IU, &ABSTOL, &MM, (double *)eigenvalues, &Z[0], &LDZ, ISUPPZ, WORK.get(), &LWORK,
+                RWORK.get(), &LRWORK, IWORK.get(), &LIWORK, &INFO);
   if (INFO != 0) my_error("zheevr failed. INFO=%i", INFO);
   if (MM != int(M)) {
     cout << "zheevr computed " << MM << "/" << M << endl;
@@ -269,9 +268,6 @@ size_t diagonalise_zheevr(Matrix &m, Eigen &d, char jobz = 'V', double ratio = 1
       }
     d.perform_checks();
   }
-  delete[] WORK;
-  delete[] RWORK;
-  delete[] IWORK;
   return M;
 }
 #endif
@@ -293,7 +289,13 @@ Eigen diagonalise(Matrix &m) {
   dr_value dr = sP.diagroutine;
 #ifdef NRG_REAL
   if (dr == diagdsyev) M = diagonalise_dsyev(m, d, 'V');
-  if (dr == diagdsyevd) M = diagonalise_dsyevd(m, d, 'V');
+  if (dr == diagdsyevd) {
+    M = diagonalise_dsyevd(m, d, 'V');
+    if (M == 0) {
+      std::cout << "dsyevd failed, falling back to dsyev" << std::endl;
+      M = diagonalise_dsyev(m, d, 'V');
+    }
+  }
   if (dr == diagdsyevr) M = diagonalise_dsyevr(m, d, 'V', sP.diagratio);
 #endif
 #ifdef NRG_COMPLEX

c++/nrg-lib.cc

@@ -2176,7 +2176,8 @@ void nrg_perform_measurements(DiagInfo &diag) {
  use the quantum number differences in array In[] (obtained by a call to
  function input_subspaces), make a list of all possible subspaces and
  remove duplicates. */
-void nrg_make_subspaces_list(list<Invar> &subspaces) {
+auto nrg_make_subspaces_list() {
+  list<Invar> subspaces;
   for (const auto &ii : diagprev)
     if (NRSTATES(ii)) {
       const Invar I = INVAR(ii);
@@ -2189,6 +2190,7 @@ void nrg_make_subspaces_list(list<Invar> &subspaces) {
     }
   subspaces.sort();
   subspaces.unique();
+  return subspaces;
 }
 
 /* Define recalculation strategy
@@ -2475,8 +2477,7 @@ void nrg_diagonalisations() {
 void nrg_determine_tasks(map<Invar, InvarVec> &ancestors) {
   nrglog('@', "@ nrg_determine_tasks()");
   // Make a list of all subspaces to consider.
-  list<Invar> subspaces;
-  nrg_make_subspaces_list(subspaces);
+  auto subspaces = nrg_make_subspaces_list();
   // Auxiliary information: ancestor subspaces and their dimensions.
   qsrmax.clear();
   ancestors.clear();

c++/nrg.cc

@@ -8,6 +8,15 @@ mpi::environment *mpienv;
 mpi::communicator *mpiw;
 #endif
 
+inline void help(int argc, char **argv, std::string help_message)
+{
+  std::vector<std::string> args(argv+1, argv+argc); // NOLINT
+  if (args.size() >= 1 && args[0] == "-h") {
+    std::cout << help_message << std::endl;
+    exit(EXIT_SUCCESS);
+  }
+}
+
 int main(int argc, char **argv) {
 #ifdef NRG_MPI
   mpi::environment env(argc, argv);
@@ -21,6 +30,7 @@ int main(int argc, char **argv) {
 #endif
     print_about_message();
     report_openMP();
+    help(argc, argv, "Usage: nrg [-h] [-w workdir]");
     set_workdir(argc, argv);
     run_nrg_master();
     time_mem::memory_report();

nrginit/nrginit

@@ -3,6 +3,12 @@
 # initialisation script
 # Rok Zitko, [email protected], 2009-2018
 
+if [ "$1" = "-h" ]
+then
+  echo Usage: nrginit
+  exit 0
+fi
+
 NRGDIR="$( cd "$( dirname "${BASH_SOURCE[0]}" )" >/dev/null && pwd )"
 NRG_INITSCRIPT=${NRGDIR}/nrginit.m
 ${1:-math} -batchinput -batchoutput <<EOF

tools/adapt/adapt.cc

@@ -338,7 +338,18 @@ double rhs_F(double x, double y) {
 
 void about(ostream &F = cout) {
   F << "# Discretization ODE solver" << endl;
-  F << "# Rok Zitko, [email protected], 2008-2019" << endl;
+  F << "# Rok Zitko, [email protected], 2008-2020" << endl;
+}
+
+const std::string usage{"Usage: adapt [-h] [P|N] [param_filename]"};
+
+inline void help(int argc, char **argv, std::string help_message)
+{
+  std::vector<std::string> args(argv+1, argv+argc); // NOLINT
+  if (args.size() >= 1 && args[0] == "-h") {
+    std::cout << help_message << std::endl;
+    exit(EXIT_SUCCESS);
+  }
 }
 
 void cmd_line(int argc, char *argv[]) {
@@ -347,7 +358,7 @@ void cmd_line(int argc, char *argv[]) {
     switch (first) {
       case 'P': sign = POS; break;
       case 'N': sign = NEG; break;
-      default: cerr << "Usage: adapt [P|N] [param_filename]" << endl; exit(1);
+      default: cerr << usage << endl; exit(1);
     }
   }
 
@@ -370,7 +381,7 @@ void load_init_rho() {
   string rhofn = Pstr("dos", "Delta.dat");
   vecrho       = load_rho(rhofn, sign);
   add_zero_point(vecrho);
-      
+
   rescalevecxy(vecrho, 1.0 / bandrescale, bandrescale);
   minmaxvec(vecrho, "rho");
   rho = LinInt(vecrho);
@@ -484,6 +495,7 @@ int main(int argc, char *argv[]) {
   clock_t start_clock = clock();
 
   about();
+  help(argc, argv, usage);
   cmd_line(argc, argv);
   parser(param_fn);
   set_parameters();

tools/binavg/binavg.cc

@@ -47,6 +47,7 @@ void usage(ostream &F = cout) {
   F << "Usage: binavg <name> <Nz>" << endl;
   F << endl;
   F << "Optional parameters:" << endl;
+  F << " -h -- show help" << endl;
   F << " -v -- verbose" << endl;
   F << " -o -- one .dat file" << endl;
   F << " -2 -- use the 2nd column for weight values (complex spectra)" << endl;
@@ -56,8 +57,12 @@ void usage(ostream &F = cout) {
 void cmd_line(int argc, char *argv[]) {
   char c;
 
-  while ((c = getopt(argc, argv, "vo23")) != -1) {
+  while ((c = getopt(argc, argv, "hvo23")) != -1) {
     switch (c) {
+      case 'h':
+       usage();
+       exit(EXIT_SUCCESS);
+
       case 'v': verbose = true; break;
 
       case 'o': one = true; break;

tools/broaden/broaden.cc

@@ -76,6 +76,7 @@ void usage(ostream &F = cout) {
   F << "Usage: broaden <name> <Nz> <alpha> <T> [omega0_ratio]" << endl;
   F << endl;
   F << "Optional parameters:" << endl;
+  F << " -h -- show help (when used as sole cmd line switch)" << endl;
   F << " -v -- verbose" << endl;
   F << " -m <min> -- minimal mesh frequency" << endl;
   F << " -M <max> -- maximal mesh frequency" << endl;
@@ -99,8 +100,11 @@ void usage(ostream &F = cout) {
 }
 
 void cmd_line(int argc, char *argv[]) {
+  if (argc == 2 && string(argv[1]) == "-h") {
+    usage();
+    exit(EXIT_SUCCESS);
+  }
   char c;
-
   while ((c = getopt(argc, argv, "vm:M:r:o23nscgf:x:a:l:h:PNAB")) != -1) {
     switch (c) {
       case 'c': cumulative = true; break;

tools/bw/bw.cc

@@ -68,13 +68,19 @@ vec mesh; // Frequency mesh
 vec a;                // Spectral function [current approximation]
 vec inta, intb, intc; // Integrated spectral functions
 
-void usage(ostream &F = cout) { F << "Usage: bw <name> <b0> <Nz>\n"; }
+void usage(ostream &F = cout) { 
+  F << "Usage: bw [-h] [-vVpsSo] [-m min] [-M max] [-r ratio] [-n nr_iter] [-t trim] [-q q] <name> <b0> <Nz>\n"; 
+}
 
 void cmd_line(int argc, char *argv[]) {
   char c;
 
-  while ((c = getopt(argc, argv, "vVm:M:r:n:t:pq:sSx:d:o")) != -1) {
+  while ((c = getopt(argc, argv, "hvVm:M:r:n:t:pq:sSx:d:o")) != -1) {
     switch (c) {
+      case 'h':
+        usage();
+        exit(EXIT_SUCCESS);
+
       case 'v': verbose = true; break;
 
       case 'V': veryverbose = true; break;

tools/diag/diag.cc

@@ -42,8 +42,6 @@ void about(ostream &OUT = cout) {
   }
 }
 
-void usage(ostream &OUT = cout) { OUT << "Usage: diag <input file>" << endl; }
-
 inline int MAX(int a, int b) { return (a > b ? a : b); }
 
 #include "lapack.h"
@@ -178,11 +176,19 @@ void diag_stream(istream &F) {
   diagonalize(N, data);
 }
 
+void usage(ostream &OUT = cout) { 
+  OUT << "Usage: diag [-h] [-t | -T] [-b | -B] [-vVq] [-o fn_val] [-O fn_vec] [-s scale] <input file>" << endl; 
+}
+
 void parse_param(int argc, char *argv[]) {
   char c;
 
-  while ((c = getopt(argc, argv, "tTbBvVqo:O:s:")) != -1) {
+  while ((c = getopt(argc, argv, "htTbBvVqo:O:s:")) != -1) {
     switch (c) {
+      case 'h':
+        usage();
+        exit(EXIT_SUCCESS);
+
       case 't': output_text = false; break;
 
       case 'T': output_text = true; break;

tools/hilb/hilb.cc

@@ -127,7 +127,7 @@ double calcim() {
     double lim2up   = -1;
 
     bool inside;
-    
+
     // x within the band
     if (W1 < 0 && W2 > 0) {
       lim1down = ln1016;
@@ -235,7 +235,7 @@ double calcre() {
     double lim2up   = -1;
 
     bool inside;
-    
+
     // x within the band
     if (W1 < 0 && W2 > 0) {
       lim1down = ln1016;
@@ -335,6 +335,7 @@ void usage() {
   cout << "Usage (3): hilb [options] <resigma.dat> <imsigma.dat> <reaw.dat> <imaw.dat>" << endl;
   cout << endl;
   cout << "Options:" << endl;
+  cout << "-h        show help" << endl;
   cout << "-d <dos>  Load the density of state data from file 'dos'" << endl;
   cout << "          If this option is not used, the Bethe lattice DOS is assumed." << endl;
   cout << "-v        Increase verbosity" << endl;
@@ -497,8 +498,12 @@ void parse_param_run(int argc, char *argv[]) {
   ostream *OUT = &cout;
   ofstream OUTFILE;
 
-  while ((c = getopt(argc, argv, "Gd:vVs:B:o:")) != -1) {
+  while ((c = getopt(argc, argv, "hGd:vVs:B:o:")) != -1) {
     switch (c) {
+      case 'h':
+        usage();
+        exit(EXIT_SUCCESS);
+
       case 'G': G = true; break;
 
       case 'd': load_dos(optarg); break;