Numerical library

MADNESS provides a high-level environment for the solution of integral and differential equations in many dimensions using adaptive, fast methods with guaranteed precision based on multi-resolution analysis and novel separated representations.

Useful examples can be found in the src/examples directory, where the use of the numerical library can be practiced by example.

• sininteg.cc: Create a function and integrate

• hatom_energy.cc: Compute the energy of the hydrogen atom

• heat.cc: apply the Greens function to the heat equation

• 3dharmonic.cc: solve the 3D harmonic oscillator

• h2.cc: solve the Hartree-Fock equations for the H2 molecule

• nonlinschro.cc: use the KAIN solver for accelerating the solution of a system of non-linear equations

• hedft.cc: use density functional theory

Some (dated) documentation for the examples and the API is here.

Example for MADNESS as an external library

To use MADNESS as an external library in a code we recommend cmake. Build and install MADNESS to an install directory. Set

export MADNESS_DIR=/path/to/madness/install/directory/

In file CMakeLists.txt:

cmake_minimum_required(VERSION 3.22)
project(yourbinary)
set(CMAKE_CXX_STANDARD 17)
find_package(MADNESS CONFIG REQUIRED)
add_executable(yourbinary main.cc)
target_link_libraries(yourbinary madness)

In file main.cc:

#include <madness.h>
using namespace madness;
int main(int argc, char* argv[]) {
    World& world=initialize(argc,argv);
    startup(world,argc,argv,true);
    FunctionDefaults<1>::set_cubic_cell(-10,10);
    FunctionDefaults<1>::set_k(8);
    FunctionDefaults<1>::set_thresh(1.e-6);
    try {
        auto gaussian = [](const Vector<double,1>& r){return std::exp(-r[0]*r[0]);};
        Function<double,1> f = FunctionFactory<double,1>(world).f(gaussian);
        double I = f.trace();
        double value = f(0.7);
        if (world.rank() == 0) {
          print("trace(exp(-r^2)",I,"error",I-std::sqrt(M_PI));
          print("f(0.7)",value,"error",value-exp(-0.7*0.7));
        }
    } catch (...) {
         std::cout << "caught an error " << std::endl;
    } 
    finalize();
    return 0;
}
}

Going through the code line by line:

#include <madness.h>

Includes the MADNESS library.

using namespace madness;
int main(int argc, char* argv[]) {
    World& world=initialize(argc,argv)

World contains the MPI communicator and must always be set, even if the code will run only on one node.

    startup(world,argc,argv,true)

Reads all necessary numerical data (e.g. the wavelet twoscale coefficients) and prints out compiler flags.

    FunctionDefaults<1>::set_cubic_cell(-10,10)

Defines the computation cell/intervall, here in 1 dimension. MADNESS is templated with respect to the dimensions.

    FunctionDefaults<1>::set_k(8);

Sets the wavelet order to 8. Anything between 2 and 30 will work, the optimal choice depends on the specific problem, 8 is usually a good guess.

    FunctionDefaults<1>::set_thresh(1.e-7);

Sets the precision threshold to \(\epsilon=10^{-6}\) in the \(L^2\) norm. Unless the function is singular or has other pathological features the threshold will be met. Up to the limits of numerical precision, it is possible to tighten the threshold to secure more digits. There are also different truncation modes. The default is mode=0 which is designed to yield accurate function values and norms, whereas mode=1 aims to also yield accurate derivatives.

    try {
        auto functor=[](const Vector<double,1>& r){return std::exp(-r[0]*r[0]);}

Defines the function to be represented in MRA. Vector is a MADNESS class defining coordinates.

        Function<double,1> f=FunctionFactory<double,1>(world).f(functor);

Projects the function \(e^{-r^2}\) onto the MRA representation. A FunctionFactory is used to define various properties of the Function, it requires world as input, a function (or functor) to compute the value. Additional optional arguments can be provided using the named-parameter idiom.

        double I=f.trace();

Integrates the function \(\int_{-10}^{-10} e^{-r^2}\mathrm dx\). It is a collective operation.

        double value = f(0.7);

evalues the numerical function at point `x=0.7. This is executed by a single process but might involve remote communication. If you wish to evaluate at many points (e.g., a line or a cube, there are more efficient interfaces).

        if (world.rank() == 0) {
          print("trace(exp(-r^2)",I,"error",I-std::sqrt(M_PI));
          print("f(0.7)",value,"error",value-exp(-0.7*0.7));
        }

Prints out the values using the MADNESS Python-like print function. Be sure to add the if block to avoid verbose output if run on many MPI ranks. Also make sure that the trace or any other collective operation is not called inside this if block. A collective operation called by only one rank will cause the program to hang. This is a common error.

    } catch (...) {
         std::cout << "caught an error " << std::endl;
    } 
    finalize();

Finalizes the communicator. It is important that all MRA objects (e.g. Function<double,1>) are destructed before finalize() is called, otherwise segmentation faults might occur since the destructor for these objects will erroneously be called at the very end of the program after the runtime has been dismantled. Thus, for simple programs enclose all MRA code after initialize in a sub-scope (e.g., using braces or inside a try/catch block), or after obtaining World pass it (by reference) into another procedure that contains your MRA code.

    return 0;
}