M. T. Homer Reid MIT Home Page
Physics Problems Research teaching My Music About Me Miscellany


C/C++ Libraries for Computational Math and Physics

Here's a compendium of C/C++ libraries that I've put together to carry out various computational tasks that seem to arise frequently in my research. Some of these are C/C++ wrappers around legacy Fortran codes; one of them is a hacked version of another open-source C/C++ code; and the rest are my own creations. (Which is which is indicated in the table below.)

Table Of Contents
1. One-Line Summaries of the Libraries
2. General Installation and Usage Notes
3. Detailed Overviews of The Libraries

1. One-Line Summaries of the Libraries

In this table,

  • “wrapper” denotes a Fortran code that I found online and repackaged as a C/C++ library, sometimes with slight modifications or extensions
  • “hack” denotes a C/C++ code that I found online and modified or extended
  • “original” denotes an original creation, i.e. a homemade code that I created from scratch.

Library Purpose Provenance
Libraries for fitting and interpolation
libPolyFit 1D polynomial fitting of data Original
libMDInterp Interpolation of functions in 1, 2, 3, or 4 dimensions Original
Libraries for special functions
libAmosBessel Bessel functions of complex arguments Wrapper
Libraries for numerical integration
libvvqag Adaptive quadrature of one-dimensional functions Hack
libTSQ Tanh-sinh quadrature of one-dimensional functions Original
libDCUHRE Adaptive cubature of multidimensional functions over hypercubic regions Wrapper
libDCUTRI Adaptive cubature of functions over triangles Wrapper
libTriInt Fixed-order cubature of functions over triangles Original
Libraries for Matrix and Vector Manipulations
libhmat Linear algebra operations on matrices and vectors Original
Libraries for Electromagnetic Scattering
libRWG A boundary-element field solver for 3D scattering problems Original
libTDRT A boundary-element field solver for 2D scattering problems Original
libIncField Calculation of electromagnetic fields due to various types of known sources Original
libMatProp Efficient database maintenance for frequency-dependent material properties Original
libSpherical Convenience routines for computations in spherical coordinates Original
libSphericalScattering Electromagnetic scattering from spheres and spherical shells Original
Libraries for Quantum Chemistry
libHQC An infrastructure library for quantum chemistry calculations with contracted Gaussian orbitals Original
libHHF A Hartree-Fock solver for quantum chemistry Original
Miscellaneous Utility Libraries
libhrutil Convenient utility functions for various purposes Original

2. Installation and usage notes

  • In general, you should be able to just download and unpack the source packages and do a simple make to compile each library. (There is no make install; each library here has just one .h and one .a file, and it is up to you whether you want to copy them somewhere else or leave them where they are.)

  • In addition to the code snippets below, all source packages come with a test/demonstration program named tlibName (where libName is the library in question; thus, tlibRWG, tlibPolyFit, etc.) You can build this program by doing a make test, and then (if this succeeds) review its command-line options by doing a tlibName --help.

  • Many of the libraries here are simple enough that you can understand what they offer just by looking at the code snippets, the example programs, and the header files. For the larger-scale libraries, I have provided API documentation.

  • Some of the libraries on this page depend on other libraries on this page (“internal” dependencies) or on open-source libraries available elsewhere (“external” dependencies). These dependencies are specified in the summaries below.

  • By default, the Makefile attempts to resolve internal dependencies by assuming that you have unpacked and built any other required library packages in the same directory as the package you are trying to build. Thus, if you are trying to compile libSpherical (which depends on libAmosBessel) and you have unpacked the libSpherical source package into the directory ~HRLibs/libSpherical, then the Makefile will attempt to reference the file ~HRlibs/libAmosBessel/libAmosBessel.h (and also ~HRlibs/libAmosBessel/libAmosBessel.a if you attempt to build the test program tlibSpherical.)

  • Send questions, bug reports, job offers, and (particularly) reports of scenarios in which you found any of these useful to: .

3. Detailed overviews of the libraries

libPolyFit: Polynomial Fitting of 1D Data

The Setting:
  • You have a collection of data points (xi,yi) that you want to fit to a one-dimensional polynomial y=f(x).
  • You choose the order of the polynomial, which may be greater than, equal to, or less than the number of data points.

Download: libPolyFit.tar.gz
Dependencies: libhrutil
Provenance: Original.
Code snippet:

          int NumPts=6;
          double X[NumPts], Y[NumPts]; 
          PolyFit *PF4, *PF8;

          /* fill in data arrays */
          GetExperimentalData(X, Y, NumPts);

          /* fit the data to an 4th order polynomial */
          PF4=new PolyFit(X,Y,NumPts,4);

          /* fit the data to an 8th order polynomial */
          PF8=new PolyFit(X,Y,NumPts,8);

          /* plot the fits with the data points */
          PF4->PlotFit(X, Y, NumPts);
          PF8->PlotFit(X, Y, NumPts);

          /* use the fits to extrapolate to X=0 */
          printf("Extrapolated to X=0 (4th order): %e\n",PF4->f(0.0));
          printf("Extrapolated to X=0 (8th order): %e\n",PF8->f(0.0));


libMDInterp: Interpolation of functions in 1, 2, 3, or 4 dimensions

The Setting:
  • You have a collection of functions of 1, 2, 3, or 4 variables that are expensive to compute. (A common example is a definite integral, which depends on a number of external parameters, and which cannot be evaluated in closed form but which can be evaluated numerically for arbitrary values of the parameters.)
  • You need to evaluate your function(s) at a large number of points within some 1, 2, 3, or 4 dimensional hypercubic region (which may be infinite), and you don't want to pay the cost of recomputing the function at each separate point.

  • Instead, you want to precompute function values at the points of a 1, 2, 3, or 4-dimensional grid spanning your hypercube, store these data in a big table, and then obtain approximate values of your function at arbitrary points inside the hypercube by interpolating among the known function values at the grid points.

  • Often, you will want to do the precomputation step ahead of time and save the grid data in a big data file, which you can then subsequently reread into memory whenever you have a code that needs to compute values of your function(s).

Download: libMDInterp.tar.gz
Dependencies: libhmat, libhrutil
Provenance: Original.
Code snippet:

         /* 3D function to be interpolated ******************************/
         void My3DFunction(double X1, double X2, double X3, double *PhiVD)
             PhiVD[0] = exp(-X1*X1-X2*X2-X3*X3);               /* value */
             PhiVD[1] = -2.0*X1*exp(-X1*X1-X2*X2-X3*X3);       /* d/dX */
             PhiVD[7] = -8.0*X1*X2*X3*exp(-X1*X1-X2*X2-X3*X3); /* d^3/dXdYdZ*/

          int main()
              int N1=4;
              double X1Points[4]={0.0, 0.1, 0.2, 0.3};
              int N1=5;
              double X2Points[5]={0.2, 0.4, 0.6, 0.8, 1.0};
              int N1=7;
              double X3Points[7]={-3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0};
              /* initialize the interpolator */
              Interp3D *I3D = new Interp3D(X1Points, N1, 
                                           X2Points, N2, 
                                           X3Points, N3,
                                           1, nThread, My3DFunction);

              /* get interpolated function values */
              double F;
              I3D->Evaluate(0.123, 0.234, 0.345, &F);




libAmosBessel: Bessel functions of general complex arguments

The Setting:
  • You need to evaluate Bessel-type functions (cylindrical Bessel functions, spherical Bessel functions, or Airy functions) at general complex arguments.
  • You need Bessel functions of arbitrary real (not just integer or half-integer) order.

Download: libAmosBessel.tar.bz2
Dependencies: None
Provenance: C++ wrapper around a Fortran code by D. E. Amos. This is the code used by Matlab for evaluating Bessel and Airy functions.
Code snippet:

          cdouble z;
          double MinOrder = 0.25;
          int NumOrders=3;
          cdouble JValues[NumOrders], iValues[NumOrders], AValue;
          z = 1.2 + 3.4fi;   // arbitrary complex number

          AmosBessel('J', z, MinOrder, NumOrders, 0, JValues);

          // at this point i have 
          //  JValues[0] = regular cylindrical bessel func J_{0.25}(z)
          //  JValues[1] = regular cylindrical bessel func J_{1.25}(z)
          //  JValues[2] = regular cylindrical bessel func J_{2.25}(z)

          AmosBessel('i', z, MinOrder, NumOrders, 0, iValues);

          // at this point i have 
          //  iValues[0] = regular modified spherical bessel func i_{0.25}(z)
          //  iValues[1] = regular modified spherical bessel func i_{1.25}(z)
          //  iValues[2] = regular modified spherical bessel func i_{2.25}(z)

          AmosAiry('A', z, 0, &AValue);
          // at this point i have AValue = airy function Ai(z)


libhmat: Manipulation of Matrices and Vectors

The Setting:
  • Your C/C++ programs need to work with vectors and matrices, which may be ordinary, symmetric, sparse, or otherwise, and you can't be bothered to worry about efficient storage and data retrieval implementations for all the various cases.
  • Your C/C++ programs need to use LAPACK routines for linear algebra, and you would prefer to write C++ syntax that looks like this:

    instead of this:
                dgetrf_(&NR, &NC, DM, &NR, ipiv, &info); 
                dgetrs_("N", &NR, &iOne, DM, &NR, ipiv, X->DV, &NR, &info); 
Download: libhmat.tar.bz2
Dependencies: libC2ML
Provenance: Original.
Code snippet:

          // create a matrix and fill in its entries by hand 
          int m, n;
          HMatrix *M1=new HMatrix(DIM, DIM, LHM_SYMMETRIC);
          for(m=0; m<DIM; m++)
           for(n=m; n<DIM; n++)
            M1->SetEntry(m,n, MyMatrixEntryFunction(m,n) );

          // read in a matrix from a matlab-style ascii data file 
          HMatrix *M2=new HMatrix("MyMatrixDataFile.mat");

          // create and initialize a vector (with random entries in this case)
          HVector *V=new HVector(DIM);
          for(m=0; m<DIM; m++)
           V->SetEntry(m, lrand48() );

          // now factorize the matrix and solve the linear system 
          if ( Success )
              // if M is positive-definite, then i can solve the linear 
              // system using cholesky factorization
           {  // otherwise, i fall back on solving the system 
              // by LU factorization

libMatProp: Database Maintenance for Frequency-Dependent Material Properties

The Setting:
  • You want to do computational electromagnetism (possibly with one or more of my software packages) involving materials with frequency-dependent material properties (permittivity ε and permeability μ).
  • For each of the various materials you will be working with, you have either (a) tabulated data for ε and μ at a range of frequencies, or (b) functional forms (such as a sum of lorentzians) for ε(ω) and μ(ω), which you would like to store somewhere in a universal database file which can be used over and over again by various software tools.

Download: libMatProp.tar.bz2
Dependencies: libhrutil, muParser
Provenance: Original.
Code snippet:

        // material properties specified by functional forms 
        // tabulated in a global database file 
        MatProp *MP1 = new MatProp("Gold");
        MatProp *MP2 = new MatProp("Bromobenzene");

        // permittivity of gold at omega = 1e14 radians / second
        cdouble Eps = MP1->GetEps(1.0e14, REAL_FREQ);

        // permeability of bromobenzene at \omega = 1e14 on imaginary axis
        double Mu = MP2->GetMu(1.0e14, IMAG_FREQ);

        // material properties obtained by interpolating among
        // values in a data file
        MatProp *MP3 = new MatProp("FILE_MyInterpolatedDataFile.dat");

        // often it's more convenient to work in terms of 
        // units in which c==1 and the default length unit is 1 um
        MatProp::SetFreqUnit( SPEED_OF_LIGHT / 1.0e-6 );

        MP3->GetEpsMu(0.12, &Eps, &Mu);


libSpherical: Utilities for Computations in 3D Spherical Polar Coordinates

The Setting:
  • You need to convert coordinates, vectors, finite-difference expressions for div, grad, and curl, and other entities back and forth between cartesian and spherical polar coordinates.
  • You need to do computations involving solutions of the scalar or vector Helmholtz equation, with either real or imaginary values of the k parameter, i.e. spherical Bessel functions and spherical harmonics.

  • You want to do electromagnetic scattering from spherical objects using my libSphericalScattering code.

Download: libSpherical.tar.bz2
Dependencies: libC2ML
Provenance: Original.
Code snippet:

         * given the cartesian coordinates of a point in space,
         * compute the polar coordinates of the same point
         * note: C2S = 'cartesian to spherical'
         *       S2C = 'spherical to cartesian'
        double xC[3], r, Theta, Phi;
        xC[0] = 0.1; // x coord
        xC[1] = 0.2; // y coord
        xC[2] = 0.3; // z coord
        CoordinateC2S(xC, &r, &Theta, &Phi);
         * given the x, y, z components of a vector at some 
         * point in space, compute the r, theta, phi components
         * of the same vector 
         * note: on entry to VectorC2S i have 
         *  vC[0,1,2] = x, y, z components of vector 
         * on return i have 
         *  vS[0,1,2] = r, theta, phi components of vector 
        double vC[3], vS[3];
        vC[0] = 0.1; // x component
        vC[1] = 0.2; // y component
        vC[2] = 0.3; // z componet
        VectorC2S(Theta, Phi, vC, vS);
         * get the values of the 'M' and 'N' vector-valued
         * Helmholtz solutions for a range of l,m values 
        int lMax=3;
        int NAlpha=(lMax+1)*(lMax+1); // total # of functions with l<= lMax
        cdouble MInt[3*NAlpha];       // interior 'M' functions
        cdouble NInt[3*NAlpha];       // interior 'N' functions
        cdouble MExt[3*NAlpha];       // exterior 'M' functions
        cdouble NExt[3*NAlpha];       // exterior 'N' functions
        GetMNlmArray(lMax, k, 1, r, Theta, Phi, INTERIOR, MInt, NInt);
        GetMNlmArray(lMax, k, 1, r, Theta, Phi, EXTERIOR, MExt, NExt);

libSphericalScattering: Electromagnetic Scattering from Spherical Geometries

The Setting:
  • You need to solve problems involving the scattering of electromagnetic radiation from a spherical object (either a solid sphere or a spherical shell).

  • You set up the problem by specifying

    • the incident field (which may be a plane wave, the field of a point source, or any other arbitrary incident field you like)

    • the frequency (which may be real or imaginary)

    • the radius (for a sphere) or the inner and outer radii (for a spherical shell)

    • the material properties (permittivity and permeability) of all regions: the exterior region, the scatterer, and the innermost region in the case of a spherical shell. Lossy materials are allowed; permittivities may have negative real parts and nonzero imaginary parts.

  • Then, libSphericalScattering solves the problem for you, after which point you can call a GetFields() routine to compute the scattered E and H fields at arbitrary points in space (both inside and outside the scatterer).

Download: libSphericalScattering.tar.bz2
API Documentation: libSphericalScattering.pdf
Dependencies: libSpherical, libAmosBessel, libhmat, libC2ML, libIncField, libSGJV, libhrutil
Provenance: Original.
Code snippet:

          // create an instance of the SSGeometry class
          SSGeometry *SSG = new SSGeometry();

          // set frequency, geometry, material properties, etc.
          SSG->SetParameters("--freq 0.1 --R 0.2 --rEpsI 13.0 --MuI 0.4");  

          // prepare to solve scattering problems using the given parameters 

          // solve a scattering problem involving a particular incident field
          SSG->Solve(MyIncidentFieldRoutine, MyUserData, EXTERIOR);

          // now compute scattered E and H fields at arbitrary points 
          SSG->GetFields(X1, EH1);
          SSG->GetFields(X2, EH2);

libhrutil: Miscellaneous Utility Functions

The Setting:
  • You want to streamline your C++ programming process by automating a number of simple rote tasks that arise over and over again.
  • More likely, you just need to compute and link one or more of my original codes, which pretty much all use functions contained in this library.

Download: libhrutil.tar.bz2
Dependencies: None.
Provenance: Original.
Code snippet:


          int main(int argc, char *argv[])
              // Vararg versions of common functions 
              char *FB=GetFileBase(argv[0]);
              FILE *f=vfopen("%s.out","w",FB);
              int err=vsystem("grep -iq %s.out 'MySymbol'",FB);

              // Matlab-style timing
              double SecondsElapsed=Toc();
              // ... and many more similar utility functions


Homer Reid's C/C++ Libraries, by Homer Reid
Last Modified: 11/16/16