Initial commitHEADv1.0master

1 files changed, 597 insertions, 0 deletions

diff --git a/src/hhtop.c b/src/hhtop.c
new file mode 100644
index 0000000..0889b76
--- /dev/null
+++ b/src/hhtop.c
@@ -0,0 +1,597 @@
+/*
+Copyright 2016 Ian Jauslin
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+
+    http://www.apache.org/licenses/LICENSE-2.0
+
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+*/
+
+#include <stdio.h>
+#include <stdarg.h>
+// define MPFR_USE_VA_LIST to enable the use of mpfr_inits and mpfr_clears
+#define MPFR_USE_VA_LIST
+#include <mpfr.h>
+#include <math.h>
+#include <libinum.h>
+
+#include "types.h"
+#include "hh_integral.h"
+#include "hh_root.h"
+#include "hh_integral_double.h"
+#include "hh_root_double.h"
+#include "ss_integral_double.h"
+#include "zz_integral.h"
+#include "zz_integral_double.h"
+#include "parser.h"
+#include "definitions.h"
+
+// usage message
+int print_usage();
+// read arguments
+int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double);
+
+// compute the first loop correction to the topological phase diagram
+int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order);
+// using doubles instead of mpfr
+int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order);
+
+// compute the sunrise diagram
+int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads);
+// using doubles instead of mpfr
+int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads);
+
+// codes for possible computations
+#define COMPUTATION_PHASE 1
+#define COMPUTATION_ZZ 2
+#define COMPUTATION_ZZZZ 3
+
+int main(int argc, char* argv[]){
+  hh_params params;
+  mpfr_t tolerance;
+  unsigned int maxiter;
+  unsigned int order;
+  int ret;
+  unsigned int computation_nr;
+  unsigned int threads;
+  unsigned int use_double;
+
+  // default computation: phase diagram
+  computation_nr=COMPUTATION_PHASE;
+
+  mpfr_inits(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
+  // read command line arguments
+  ret=read_args(argc, argv, &params, tolerance, &maxiter, &order, &computation_nr, &threads, &use_double);
+  if(ret<0){
+    mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
+    return(-1);
+  }
+  if(ret>0){
+    mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
+    return(0);
+  }
+
+  // phase diagram
+  if(computation_nr==COMPUTATION_PHASE){
+    // compute the first-loop correction to the phase diagram
+    if(use_double==0){
+      compute_hh(params, tolerance, maxiter, order);
+    }
+    else{
+      compute_hh_double(params, tolerance, maxiter, order);
+    }
+  }
+  else if(computation_nr==COMPUTATION_ZZ){
+    // compute second-order correction to z1-z2
+    if(use_double==0){
+      compute_ss(&zz_I, params, tolerance, maxiter, order, threads);
+    }
+    else{
+      compute_ss_double(&zz_I_double, params, tolerance, maxiter, order, threads);
+    }
+  }
+  else if(computation_nr==COMPUTATION_ZZZZ){
+    // compute second-order correction to z1+z2
+    if(use_double==0){
+      compute_ss(&ZZ_I, params, tolerance, maxiter, order, threads);
+    }
+    else{
+      compute_ss_double(&ZZ_I_double, params, tolerance, maxiter, order, threads);
+    }
+  }
+
+  mpfr_clears(params.W, params.sinphi, params.phi, params.t1, params.t2, params.lambda, tolerance, NULL);
+
+  return(0);
+}
+
+// usage message
+int print_usage(){
+  fprintf(stderr, "usage:\n       hhtop phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax]\n       hhtop z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n       hhtop z1+z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-P precision] [-E emax] [-C threads]\n\n       hhtop -D phase [-p params] [-v] [-O order] [-t tolerance] [-N maxiter]\n       hhtop -D z1-z2 [-p params] [-v] [-O order] [-t tolerance] [-N maxiter] [-C threads]\n\n       hhtop -V [-v]\n\n");
+  return(0);
+}
+
+// read command line arguments
+#define CP_FLAG_PARAMS 1
+#define CP_FLAG_ORDER 2
+#define CP_FLAG_TOLERANCE 3
+#define CP_FLAG_MAXITER 4
+#define CP_FLAG_MPFR_PREC 5
+#define CP_FLAG_MPFR_EXP 6
+#define CP_FLAG_THREADS 7
+int read_args(int argc, char* argv[], hh_params* params, mpfr_t tolerance, unsigned int* maxiter, unsigned int* order, unsigned int* computation_nr, unsigned int* threads, unsigned int* use_double){
+  int i;
+  int ret;
+  // temporary long int
+  long int tmp_lint;
+  // temporary unsigned int
+  unsigned int tmp_uint;
+  // pointers
+  char* ptr;
+  // flag that indicates what argument is being read
+  int flag=0;
+  // pointer to various arguments
+  char* tolerance_str=NULL;
+  char* params_str=NULL;
+  // whether to print the variables after they are read
+  int print_vars=0;
+  // keep track of which flags were used (to check for incompatibilities)
+  unsigned char pflag=0;
+  unsigned char Oflag=0;
+  unsigned char tflag=0;
+  unsigned char Nflag=0;
+  unsigned char Pflag=0;
+  unsigned char Eflag=0;
+  unsigned char vflag=0;
+  unsigned char Cflag=0;
+  unsigned char Dflag=0;
+  unsigned char Vflag=0;
+
+  // defaults
+  mpfr_set_d(params->t1, 1., MPFR_RNDN);
+  mpfr_set_d(params->t2, .1, MPFR_RNDN);
+  mpfr_set_d(params->lambda, .01, MPFR_RNDN);
+  mpfr_set_d(params->sinphi, 1., MPFR_RNDN);
+  mpfr_set_d(tolerance, 1e-11, MPFR_RNDN);
+  *maxiter=1000000;
+  *order=10;
+  params->omega=+1;
+  *threads=1;
+  *use_double=0;
+
+  mpfr_const_pi(params->phi, MPFR_RNDN);
+  mpfr_div_ui(params->phi, params->phi, 2, MPFR_RNDN);
+
+  // default W=-3*sqrt(3)*t2*sin(phi)
+  mpfr_sqrt_ui(params->W, 3, MPFR_RNDN);
+  mpfr_mul_si(params->W, params->W, -3, MPFR_RNDN);
+  mpfr_mul(params->W, params->W, params->sinphi, MPFR_RNDN);
+  mpfr_mul(params->W, params->W, params->t2, MPFR_RNDN);
+
+  // loop over arguments
+  for(i=1;i<argc;i++){
+    // flag
+    if(argv[i][0]=='-'){
+      for(ptr=((char*)argv[i])+1;*ptr!='\0';ptr++){
+	switch(*ptr){
+	// parameters
+	case 'p':
+	  flag=CP_FLAG_PARAMS;
+	  pflag=1;
+	  break;
+	// order of the integration
+	case 'O':
+	  flag=CP_FLAG_ORDER;
+	  Oflag=1;
+	  break;
+	// tolerance
+	case 't':
+	  flag=CP_FLAG_TOLERANCE;
+	  tflag=1;
+	  break;
+	// maximal number of Newton steps
+	case 'N':
+	  flag=CP_FLAG_MAXITER;
+	  Nflag=1;
+	  break;
+	// mpfr precision
+	case 'P':
+	  flag=CP_FLAG_MPFR_PREC;
+	  Pflag=1;
+	  break;
+	// mpfr emax
+	case 'E':
+	  flag=CP_FLAG_MPFR_EXP;
+	  Eflag=1;
+	  break;
+	// print value of variables
+	case 'v':
+	  print_vars=1;
+	  vflag=1;
+	  break;
+	// number of threads
+	case 'C':
+	  flag=CP_FLAG_THREADS;
+	  Cflag=1;
+	  break;
+	// use doubles instead of mpfr
+	case 'D':
+	  *use_double=1;
+	  Dflag=1;
+	  // set prec to that of long doubles (for consistency when reading cli arguments with many digits)
+	  mpfr_set_default_prec(64);
+	  break;
+	// print version
+	case 'V':
+	  Vflag=1;
+	  break;
+	default:
+	  fprintf(stderr, "unrecognized option '-%c'\n", *ptr);
+	  print_usage();
+	  return(-1);
+	  break;
+	}
+      }
+    }
+    // parameters
+    else if(flag==CP_FLAG_PARAMS){
+      // read str later (after having set the MPFR precision and emax)
+      params_str=argv[i];
+      flag=0;
+    }
+    // order of the integration
+    else if(flag==CP_FLAG_ORDER){
+      ret=sscanf(argv[i],"%u",&tmp_uint);
+      if(ret!=1){
+	fprintf(stderr, "error: '-O' should be followed by an unsigned int\n       got '%s'\n",argv[i]);
+	return(-1);
+      }
+      *order=tmp_uint;
+      flag=0;
+    }
+    // tolerance
+    else if(flag==CP_FLAG_TOLERANCE){
+      // read str later (after having set the MPFR precision and emax)
+      tolerance_str=argv[i];
+      flag=0;
+    }
+    // maximal number of Newton steps
+    else if(flag==CP_FLAG_MAXITER){
+      ret=sscanf(argv[i],"%u",maxiter);
+      if(ret!=1){
+	fprintf(stderr, "error: '-N' should be followed by a positive int\n       got '%s'\n",argv[i]);
+	return(-1);
+      }
+      flag=0;
+    }
+    // mpfr precision
+    else if(flag==CP_FLAG_MPFR_PREC){
+      ret=sscanf(argv[i],"%ld",&tmp_lint);
+      if(ret!=1){
+	fprintf(stderr, "error: '-P' should be followed by a long int\n       got '%s'\n",argv[i]);
+	return(-1);
+      }
+      mpfr_set_default_prec(tmp_lint);
+      flag=0;
+    }
+    // mpfr emax
+    else if(flag==CP_FLAG_MPFR_EXP){
+      ret=sscanf(argv[i],"%ld",&tmp_lint);
+      if(ret!=1){
+	fprintf(stderr, "error: '-E' should be followed by a long int\n       got '%s'\n",argv[i]);
+	return(-1);
+      }
+      mpfr_set_emax(tmp_lint);
+      flag=0;
+    }
+    // number of threads to use for the computation
+    else if(flag==CP_FLAG_THREADS){
+      ret=sscanf(argv[i],"%u",threads);
+      if(ret!=1){
+	fprintf(stderr, "error: '-C' should be followed by a positive int\n       got '%s'\n",argv[i]);
+	return(-1);
+      }
+      flag=0;
+    }
+    // computation to run
+    else{
+      if(str_cmp(argv[i], "phase")==1){
+	*computation_nr=COMPUTATION_PHASE;
+      }
+      else if(str_cmp(argv[i], "z1-z2")==1){
+	*computation_nr=COMPUTATION_ZZ;
+      }
+      else if(str_cmp(argv[i], "z1+z2")==1){
+	*computation_nr=COMPUTATION_ZZZZ;
+      }
+      else{
+	fprintf(stderr, "error: unrecognized computation: '%s'\n",argv[i]);
+	print_usage();
+	return(-1);
+      }
+      flag=0;
+    }
+  }
+  if(tolerance_str!=NULL){
+    ret=mpfr_set_str(tolerance, tolerance_str, 10, MPFR_RNDN);
+    if(ret<0){
+      fprintf(stderr, "error: '-t' should be followed by an MPFR floating point number\n       got '%s'\n", tolerance_str);
+      return(-1);
+    }
+  }
+  if(params_str!=NULL){
+    ret=read_params(params, params_str);
+    if(ret<0){
+      return(ret);
+    }
+  }
+
+  // check for incompatible flags
+  if((Vflag==1 && (pflag!=0 || Oflag!=0 || tflag!=0 || Nflag!=0 || Pflag!=0 || Eflag!=0 || Cflag!=0 || Dflag!=0)) || \
+     (*computation_nr!=COMPUTATION_ZZ && *computation_nr!=COMPUTATION_ZZZZ && Cflag==1) || \
+     (Dflag==1 && (Pflag==1 || Eflag==1)) \
+     ){
+    print_usage();
+    return(-1);
+  }
+
+  // print version and exit
+  if(Vflag==1){
+    printf("hhtop " VERSION "\n");
+    printf("libinum " LIBINUM_VERSION "\n");
+    if(vflag==1){
+      // print datatype information
+      printf("\n\n");
+      print_datatype_info(stdout);
+    }
+    return(1);
+  }
+
+  // print variables
+  if(print_vars==1){
+    fprintf(stderr, "t1=");
+    fprint_mpfr(stderr,params->t1);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "t2=");
+    fprint_mpfr(stderr,params->t2);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "lambda=");
+    fprint_mpfr(stderr,params->lambda);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "phi=");
+    fprint_mpfr(stderr,params->phi);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "sinphi=");
+    fprint_mpfr(stderr,params->sinphi);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "W=");
+    fprint_mpfr(stderr,params->W);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "omega=%d\n",params->omega);
+
+    fprintf(stderr, "\ntolerance=");
+    fprint_mpfr(stderr,tolerance);
+    fprintf(stderr, "\n");
+
+    fprintf(stderr, "order of integration=%d\n", *order);
+
+    fprintf(stderr, "\nMPFR precision=%ld\n", mpfr_get_default_prec());
+    fprintf(stderr, "MPFR emax=%ld\n", mpfr_get_emax());
+
+    fprintf(stderr, "\n");
+  }
+
+  return(0);
+}
+
+// compute the first loop correcion to the topological phase diagram
+int compute_hh(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){
+  mpfr_t val;
+  int ret;
+  args_integration args_int;
+
+  // compute weights
+  ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &(args_int.abcissa), &(args_int.weights));
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the precision or the tolerance\n");
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
+    return(ret);
+  }
+
+  // set args
+  args_int.params=params;
+  mpfr_init(args_int.params.W);
+
+  mpfr_init(val);
+  // initial value
+  mpfr_sqrt_ui(val, 3, MPFR_RNDN);
+  mpfr_mul_ui(val, val, 3, MPFR_RNDN);
+  mpfr_mul(val, val, params.sinphi, MPFR_RNDN);
+  mpfr_mul(val, val, params.t2, MPFR_RNDN);
+  if(params.omega==1){
+    mpfr_neg(val, val, MPFR_RNDN);
+  }
+
+  // compute root
+  ret=root_newton_inplace_mpfr(&val, &integration_wrapper, &d_integration_wrapper, tolerance, maxiter, &args_int);
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n       try increasing the precision or the tolerance\n");
+    mpfr_clear(val);
+    mpfr_clear(args_int.params.W);
+    array_mpfr_free(args_int.abcissa);
+    array_mpfr_free(args_int.weights);
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n");
+    mpfr_clear(val);
+    mpfr_clear(args_int.params.W);
+    array_mpfr_free(args_int.abcissa);
+    array_mpfr_free(args_int.weights);
+    return(ret);
+  }
+
+  fprint_mpfr(stdout, val);
+  printf("\n");
+
+  mpfr_clear(val);
+  mpfr_clear(args_int.params.W);
+  array_mpfr_free(args_int.abcissa);
+  array_mpfr_free(args_int.weights);
+
+  return(0);
+}
+// using double instead of mpfr
+int compute_hh_double(hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order){
+  long double tolerance_d;
+  hh_args_integration_double args_int;
+  long double val;
+  int ret;
+
+  // convert mpfr to double
+  tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN);
+  hh_params_todouble(&(args_int.params), params);
+
+  // compute weights
+  ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &(args_int.abcissa), &(args_int.weights));
+
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n");
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
+    return(ret);
+  }
+
+
+  // initial value
+  val=-args_int.params.omega*3*sqrtl(3)*args_int.params.t2*args_int.params.sinphi;
+
+  ret=root_newton_inplace_ldouble(&val, &hh_integration_wrapper_double, &hh_d_integration_wrapper_double, tolerance_d, maxiter, &args_int);
+
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the solution to m=0\n");
+    array_ldouble_free(args_int.abcissa);
+    array_ldouble_free(args_int.weights);
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered: either the integrand is singular or the derivative of the integral vanishes at a point of the Newton iteration\n");
+    array_ldouble_free(args_int.abcissa);
+    array_ldouble_free(args_int.weights);
+    return(ret);
+  }
+
+  printf("% .19Le\n",val);
+
+  array_ldouble_free(args_int.abcissa);
+  array_ldouble_free(args_int.weights);
+
+  return(0);
+}
+
+// compute the sunrise diagram
+int compute_ss(TYPE_I, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){
+  mpfr_t val;
+  int ret;
+  array_mpfr abcissa;
+  array_mpfr weights;
+
+  // compute weights
+  ret=gauss_legendre_weights_mpfr(order, tolerance, maxiter, &abcissa, &weights);
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the precision, the tolerance, or the maximal number of Newton steps\n");
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
+    return(ret);
+  }
+
+  mpfr_init(val);
+
+  // compute integral
+  ret=ss_integrate(&val, I, params, abcissa, weights, threads);
+  // return codes
+  if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered: the integrand is singular\n");
+    mpfr_clear(val);
+    array_mpfr_free(abcissa);
+    array_mpfr_free(weights);
+    return(ret);
+  }
+
+  fprint_mpfr(stdout, val);
+  printf("\n");
+
+  mpfr_clear(val);
+  array_mpfr_free(abcissa);
+  array_mpfr_free(weights);
+
+  return(0);
+}
+// using doubles instead of mpfr
+int compute_ss_double(TYPE_I_DOUBLE, hh_params params, mpfr_t tolerance, unsigned int maxiter, unsigned int order, unsigned int threads){
+  long double tolerance_d;
+  hh_params_double params_d;
+  long double val;
+  int ret;
+  array_ldouble abcissa;
+  array_ldouble weights;
+
+  // convert mpfr to double
+  tolerance_d=mpfr_get_ld(tolerance, MPFR_RNDN);
+  hh_params_todouble(&params_d, params);
+
+  // compute weights
+  ret=gauss_legendre_weights_ldouble(order, tolerance_d, maxiter, &abcissa, &weights);
+  // return codes
+  if(ret==LIBINUM_ERROR_MAXITER){
+    fprintf(stderr, "error: maximum number of iterations reached when computing the integration abcissa\n       try increasing the tolerance, or the maximal number of Newton steps\n");
+    return(ret);
+  }
+  else if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered when computing the integration abcissa\n");
+    return(ret);
+  }
+
+  // compute integral
+  ret=ss_integrate_double(&val, I, params_d, abcissa, weights, threads);
+  // return codes
+  if(ret==LIBINUM_ERROR_NAN){
+    fprintf(stderr, "error: infinity encountered: the integrand is singular\n");
+    array_ldouble_free(abcissa);
+    array_ldouble_free(weights);
+    return(ret);
+  }
+
+  printf("% .19Le\n",val);
+
+  array_ldouble_free(abcissa);
+  array_ldouble_free(weights);
+
+  return(0);
+}
+