diff options
Diffstat (limited to 'src/hhtop.c')
| -rw-r--r-- | src/hhtop.c | 597 |
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, ¶ms, 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(¶ms_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); +} + |