MRKISS Library
MR RK (Runge-Kutta) KISS (Keep It Super Simple)

MRKISS is a simple, tiny library with zero dependencies[*] that aims to make it easy to use and experiment with explicit Runge-Kutta methods.

From a feature standpoint this library doesn't do anything you can't find in more comprehensive packages like SUNDIALS or Julia's DifferentialEquations.jl; however, it has the rather charming feature of being super simple and self contained. This makes it ideal for sharing results and code with others without asking them to recreate my technical computing environment or learn about new software. It also makes it easy to get things up and and running quickly when moving to a new supercomputer, cloud, or HPC cluster.

Here is a complete example in 23 lines that solves Lorenz's system and generates a CSV file of the results.

program minimal

  use :: mrkiss_config,      only: rk, istats_size
  use :: mrkiss_solvers_nt,  only: fixed_t_steps
  use :: mrkiss_utils,       only: print_solution
  use :: mrkiss_erk_kutta_4, only: a, b, c

  implicit none (type, external)

  real(kind=rk), parameter :: y_iv(3)  = [1.0_rk, 0.0_rk, 0.0_rk]
  real(kind=rk), parameter :: param(3) = [10.0_rk, 28.0_rk, 8.0_rk/3.0_rk]
  real(kind=rk), parameter :: t_end    = 50.0_rk

  real(kind=rk)            :: solution(7, 10000)
  integer                  :: status, istats(istats_size)

  call fixed_t_steps(status, istats, solution, eq, y_iv, param, a, b, c, t_end_o=t_end)
  call print_solution(status, solution, filename_o="minimal.csv")

contains

  subroutine eq(status, dydt, y, param)
    integer,          intent(out) :: status
    real(kind=rk),    intent(out) :: dydt(:)
    real(kind=rk),    intent(in)  :: y(:), param(:)
    dydt = [ param(1)*(y(2)-y(1)), y(1)*(param(2)-y(3))-y(2), y(1)*y(2)-param(3)*y(3) ]
    status = 0
  end subroutine eq

end program minimal

You can check out a couple larger examples:

• Three Body Problem
• Lornez Attracter

If you just want to jump in, then take a look at the Quick Start section.

Before discussing the software, it is best if we define the problem we are solving and fix some definitions and notation – notation that is used throughout the source code and documentation.

Please note that the material below deviates from the standard presentation of Runge-Kutta methods and the Butcher tableau. In particular we dispense with separate \(\mathbf{b}\) vectors per method, the silly hats, and express most of the computations in terms of linear algebra.

Solving an ODE (Ordinary Differential Equation):

We wish to compute values for an unknown function \(\mathbf{y}:\mathbb{R}\rightarrow\mathbb{R}^n\) on some interval \([t_\mathrm{iv}, t_\mathrm{max}]\) given an initial value \(\mathbf{y}(t_\mathrm{iv})=\mathbf{y_\mathrm{iv}}\) and an implicit equation involving \(t\), \(\mathbf{y}\), and it's derivatives:

\[ \mathbf{y}' = \mathbf{f}(t, \mathbf{y}) \]

The expression above may be written in component form as:

\[ \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \mathbf{f}(t, \mathbf{y}) = \left[\begin{array}{c} \frac{\mathrm{d}y_1}{\mathrm{d}t} \\ \vdots \\ \frac{\mathrm{d}y_n}{\mathrm{d}t} \\ \end{array}\right] = \left[\begin{array}{c} f_1(t, \mathbf{y}) \\ \vdots \\ f_n(t, \mathbf{y}) \\ \end{array}\right] = \left[\begin{array}{c} f_1(t, [y_1, \cdots, y_n]^\mathrm{T}) \\ \vdots \\ f_n(t, [y_1, \cdots, y_n]^\mathrm{T}) \\ \end{array}\right] \]

If the function \(\mathbf{f}(t, \mathbf{y}) \) doesn't actually depend on \(t\), then the ODE is said to be Autonomous.

We define an \(s\)-stage embedded explicit Runge-Kutta method via a set of coefficients organized into a Butcher tableau:

\[ \begin{array}{l|llll} c_1 & a_{1,1} & a_{1,2} & \dots & a_{1,s} \\ c_2 & a_{2,1} & a_{2,2} & \dots & a_{2,s} \\ c_3 & a_{3,1} & a_{3,2} & \dots & a_{3,s} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s,1} & a_{s,2} & \dots & a_{s,s} \\ \hline \rule{0pt}{12pt} & b_{1,1} & b_{1,2} & \dots & b_{1,s} \\ & \vdots & \vdots & \ddots & \vdots \\ & b_{m,1} & b_{m,2} & \dots & b_{m,s} \\ \end{array} \]

We reference the blocks of the tableau as (Note \(\mathbf{b}\) is transposed compared to the tableau):

\[ \mathbf{a} = \begin{bmatrix} a_{1,1} & a_{1,2} & \dots & a_{1,s} \\ a_{2,1} & a_{2,2} & \dots & a_{2,s} \\ a_{3,1} & a_{3,2} & \dots & a_{3,s} \\ \vdots & \vdots & \ddots & \vdots \\ a_{s,1} & a_{s,2} & \dots & a_{s,s} \end{bmatrix} \,\,\,\,\,\,\,\, \mathbf{c} = \begin{bmatrix} c_{1} \\ c_{2} \\ c_{3} \\ \vdots \\ c_{s} \end{bmatrix} \,\,\,\,\,\,\,\, \mathbf{b} = \begin{bmatrix} b_{1,1} & \dots & b_{1,m} \\ b_{2,1} & \dots & b_{2,m} \\ b_{3,1} & \dots & b_{3,m} \\ \vdots & \ddots & \vdots \\ b_{s,1} & \dots & b_{s,m} \end{bmatrix} = \begin{bmatrix} \vert & & \vert \\ \mathbf{b}_1 & \dots & \mathbf{b}_m \\ \vert & & \vert \\ \end{bmatrix} \]

Note we use \(\mathbf{b}_i\) with an index to indicate the \(i\)-th column of the matrix, and \(b_{i,j}\) to indicate the \(i,j\)-th element of the matrix.

The word embedded is used to indicate that the tableau actually defines \(m\) Runge-Kutta methods – one for each column of \(\mathbf{b}\). The most common cases are \(m=1\) and \(m=2\) with the first being most commonly used for simple, fixed step size solvers and second most commonly used for adaptive step size solvers.

We say the overall number of stages for the entire embedded method is \(s\), but the number of stages required for an individual method associated with a column of \(\mathbf{b}\) might be lower if that column has trailing zeros. If a column has \(k\) trailing zeros the Runge-Kutta method associated with that column has\(s-k\) stages.

Our goal is to compute values of \(\mathbf{y}(t)\) for \(t\in[t_\mathrm{iv}, t_\mathrm{max}]\). We are given the value of \(\mathbf{y}(t_\mathrm{iv})\). So we define a small, positive value \(\Delta{t}\in\mathbb{R}\), and attempt to find \(\mathbf{\Delta{y}}\) such that \(\mathbf{y}(t_\mathrm{iv})+\mathbf{\Delta{y}}\approx\mathbf{y}(t_\mathrm{iv}+\Delta{t})\).

We begin by defining an \(s\times s\) matrix \(\mathbf{k}\):

\[ \mathbf{a} = \begin{bmatrix} k_{1,1} & k_{1,2} & \dots & k_{1,s} \\ k_{2,1} & k_{2,2} & \dots & k_{2,s} \\ k_{3,1} & k_{3,2} & \dots & k_{3,s} \\ \vdots & \vdots & \ddots & \vdots \\ k_{s,1} & k_{s,2} & \dots & k_{s,s} \end{bmatrix} = \begin{bmatrix} \vert & & \vert \\ \mathbf{k}_1 & \dots & \mathbf{k}_s \\ \vert & & \vert \\ \end{bmatrix} \]

Using

\[ \mathbf{k}_i = \mathbf{f}\left(t + c_i \Delta{t},\, \mathbf{y} + \Delta{t} \mathbf{k}_{1..s,1..i} \mathbf{a}_{1..i,1..i} \right) \]

Explicit methods, which are the focus of MRKISS, have \(c_1=0\) and \(a_{ij}=0\) for \(i\le j\). Therefore \(\mathbf{k}_1 = \mathbf{f}(t,\, \mathbf{y})\), and each subsequent \(\mathbf{k}\) value may be computed in sequence.

Each column of \(\mathbf{k}\mathbf{b}\) forms a \(\mathbf{\Delta{y}}\) we can use to approximate \(\mathbf{y}(t_\mathrm{iv}+\Delta{t})\):

\[ \mathbf{k}\mathbf{b} = \begin{bmatrix} \vert & & \vert \\ \mathbf{\Delta{y}}_1 & \dots & \mathbf{\Delta{y}}_m \\ \vert & & \vert \\ \end{bmatrix} \]

If we use \(\mathbf{\Delta{y}}_1\) to approximate \(\mathbf{y}(t_\mathrm{iv}+\Delta{t})\), then we can use \(\mathbf{\Delta{y}}_2\) to approximate the error like this:

\[\vert\mathbf{\Delta{y}}_1 - \mathbf{\Delta{y}}_2 \vert\]

We say that a Runge-Kutta method is of order \(p\in\mathbb{Z_+}=\{i\in\mathbb{Z}\,\,\vert\,\,i>0\}\) if the error at each step is on the order of \(\mathcal{O}(\Delta{t}^{p+1})\) and the accumulated error of multiple steps is of order \(\mathcal{O}(\Delta{t}^{p})\). This abuse of big-O terminology is traditionally also extended to an abuse of big-O notation by saying the Runge-Kutta method is \(\mathcal{O}(p)\).

The IVPs I work with are generally pretty well behaved

• Non-stiff
• Time forward (\(\Delta{t} \ge 0\))
• Defined by a small (<50) set of equations expressible in closed form.

Typical examples are strange attractors and systems related to chaotic science models from celestial/classical mechanics, population dynamics, oscillating chemical reactions, and electronic circuits. My motivation for solving IVPs generally revolve around generative art and visualization. You will actually see this in the code and feature set of the library.

Things I care about

• Simple to use for simple problems.
• Graceful response to evaluation failure in derivative functions
• A good selection of predefined RK methods
• Standard, adaptive step solver with programmable step processing:
• Easy deployment & sharing
  • Easy to compile and tune for a new hardware architectures.
  • Zero external dependencies[*] except a Fortran compiler.
  • 100% standard Fortran that works with various compilers.
  • Simple text output that can be compressed and sent back home or shared with others.
• Runge-Kutta Research
  • Try out new RK methods by simply feeding the solvers a Butcher tableau.
  • Directly accessible one step routines for assembling custom solvers.
  • Simple code flow to facilitate instrumentation and deep runtime analysis and reporting.
  • Individual access to each method in an embedded tableau, and control over how each is used.
  • Machine readable butcher tableau data and Maple worksheets to process that data.
  • A few of the RK methods included are of historical or research interest – not necessarily very practical.
• Generative art and visualization functionality
  • Solutions include derivative values for visualization tools.
  • Programmable step processing tuned for visualization problems.
    • Stopping integration. Examples:
      • If the solution curve is too long in \(\mathbf{y}\)-space
      • If the step delta, or some components of it, are too long in \(\mathbf{y}\)-space
      • If the solution has returned to the IV
      • If the solution intersects itself
    • Provide an alternate \(\Delta{t}\) and redo the step based on some condition.
    • Trigger a bisection search for a \(\Delta{t}\) fitting some condition based on t-space and/or \(\mathbf{y}\)-space. Examples:
      • Find \(\Delta{t}\) so that \(\Delta{\mathbf{y}}\), or some components of it, are the perfect length.
      • Find where a step crosses over a boundary in space (ex: root finding)
      • Find where a step approaches closest to a point (ex: like the problem's IV)
  • Easy to use, hardwired methods for fixed step size visualization use cases:
    • Fixed \(t\) step size solvers
    • Fixed \(\mathbf{y}\) step size solver with hard bounds on length error
    • Fixed \(\mathbf{y}\) step size solver with soft bounds on length error
  • Interpolate entire solutions to new time points (Hermite & linear).

Things that may weird you out

Derivative function computation dominates the compute time for many applications. This observation has had tremendous influence on the development of IVP numerical techniques and related software. For the IVPs we are interested in the derivative functions are quite simple and the Runge-Kutta computations themselves generally dominate the computation time. This changes things quite dramatically. For example, techniques like embedded polynomial schemes for intrastep interpolation designed to minimize function evaluation may actually perform worse than simply making more Runge-Kutta steps for our problems! As another example, performing iterative algorithms, like bisection, directly on the Runge-Kutta method itself is quite practical. Keeping this in mind will reduce the number of "WTH?!?!" moments you have while reading the code, and will keep you from applying some features, like fixed_y_steps(), to problems that do have difficult derivatives to compute.

While MRKISS provides comprehensive support for detecting and forwarding runtime errors, it doesn't do much programmer error checking. For example, the code doesn't double check that the programmer has provided appropriate arguments for routines.

In MRKISS explicit Runge-Kutta methods are specified by directly providing the Butcher tableau via arguments to subroutines.

MRKISS provides several predefined methods in modules found in the "MRKISS/lib/" directory. Each module defines a single tableau via parameters with names mirroring the Butcher Tableau arguments documented in the previous section.

The modules follow a simple naming conventions:

• They have one of two prefixes. mrkiss_eerk_ for embedded and mrkiss_erk_ for non-embedded.
• The names end with underscore separated integers indicating the orders of the first two methods.

In addition to the parameters, the comments in these files normally include at least the following three sections:

IMO

Personal commentary about the method in question. Please note this material is simply my personal opinion.

Known Aliases

These include names used in the literature as well as names in some common ODE software.

References

I try to include the original reference if I have it. I also frequently include discussions found in other texts.

To make all this concrete, here is what one of these modules looks like (mrkiss_erk_kutta_4.f90):

!> Butcher tableau for the classic 4 stage Runge-Kutta method of O(4)
!!
!! IMO: Useful for low accuracy applications; however, I find I rarely use it.
!!
!! Known Aliases: 'RK4' (OrdinaryDiffEq.jl), 'RK41' (Butcher), & 'The Runge-Kutta Method'.
!!
!! @image html erk_kutta_4-stab.png
!!
!! @par Stability Image Links
!! <a href="erk_kutta_4-stab.png">  <img src="erk_kutta_4-stab.png"  width="256px"> </a>
!! <a href="erk_kutta_4-astab.png"> <img src="erk_kutta_4-astab.png" width="256px"> </a>
!! <a href="erk_kutta_4-star1.png"> <img src="erk_kutta_4-star1.png" width="256px"> </a>
!! <a href="erk_kutta_4-star2.png"> <img src="erk_kutta_4-star2.png" width="256px"> </a>
!!
!! @par References:
!!  - Kutta (1901); Beitrag Zur N\"herungsweisen Integration Totaler Differentialgleichungen; Z. Math. Phys. 46; p435-53
!!  - Hairer, Norsett & Wanner (2009). Solving Ordinary Differential Equations. I: Nonstiff Problems. p138;
!!    zotero://select/items/0_VLZWN2CT
!!  - Butcher (2016); Numerical Methods for Ordinary Differential Equations. 3ed; p102
!!
module mrkiss_erk_kutta_4
  use mrkiss_config, only: rk
  implicit none
  public
  !> The order of the overall method
  integer,          parameter :: s      = 4
  !> Number of methods
  integer,          parameter :: m      = 1
  !> The @f$\mathbf{a}@f$ matrix for the Butcher Tableau. @hideinitializer @showinlinesource
  real(kind=rk),    parameter :: a(s,s) = reshape([ 0.0_rk, 0.0_rk, 0.0_rk, 0.0_rk,  &
       &                                            1.0_rk, 0.0_rk, 0.0_rk, 0.0_rk,  &
       &                                            0.0_rk, 1.0_rk, 0.0_rk, 0.0_rk,  &
       &                                            0.0_rk, 0.0_rk, 2.0_rk, 0.0_rk], [s, s]) / 2.0_rk
  !> The @f$\mathbf{b}@f$ matrix for the Butcher Tableau. @hideinitializer @showinlinesource
  real(kind=rk),    parameter :: b(s,m) = reshape([ 1.0_rk, 2.0_rk, 2.0_rk, 1.0_rk], [s, m]) / 6.0_rk
  !> The @f$\mathbf{c}@f$ matrix for the Butcher Tableau. @hideinitializer @showinlinesource
  real(kind=rk),    parameter :: c(s)   = [         0.0_rk, 1.0_rk, 1.0_rk, 2.0_rk]          / 2.0_rk
  !> The method orders
  integer,          parameter :: p(m)   = [4]
  !> Number of stages for each method
  integer,          parameter :: se(m)  = [4]
end module mrkiss_erk_kutta_4

Also note all the zeros. KISS! Seriously, it takes up a tiny bit of extra space and simplifies the code considerably…

Each embedded method defines at least two Runge-Kutta methods. Normally the first two methods are used in conjunction to simultaneously estimate the solution and the error. In this library, first method is recommended for approximating the solution while the second method should be used to estimate error. This is a recommendation, and is in no way enforced by the library. When the higher order method is used for the solution, we say we are using local extrapolation. Note that each of the methods in an embedded Butcher tableau may be used individually as a non-embedded method.

In addition to the module files, each tableau may be found in the "MRKISS/rk_methods_maple/" directory along with several Maple worksheets to manipulate them. The filenames mirror the names of the modules.

I use BOO to mean Best Of Order. When used with two integers it is for an embedded method. I use RC to indicate a method is a Release Candidate – i.e. it may well become a future BOO. Not the contents of the Status columns are simply my personal opinion based on my experience. YMMV

The equation to be solved is implemented in a user provided subroutine with one of the following two signatures:

For Non-autonomous (with t) problems (found in the module mrkiss_solvers_wt):

subroutine deq_iface(status, dydt, t, y, param) ! Non-Autonomous Case (with t)
  integer,          intent(out) :: status
  real(kind=rk),    intent(out) :: dydt(:)
  real(kind=rk),    intent(in)  :: t
  real(kind=rk),    intent(in)  :: y(:)
  real(kind=rk),    intent(in)  :: param(:)

For Autonomous (no t) problems (found in the module mrkiss_solvers_nt)::

subroutine deq_iface(status, dydt, y, param)    ! Autonomous Case (no t)
  integer,          intent(out) :: status
  real(kind=rk),    intent(out) :: dydt(:)
  real(kind=rk),    intent(in)  :: y(:)
  real(kind=rk),    intent(in)  :: param(:)

The arguments are as follows:

status ........ A status code. A positive value indicates failure.
                Do not return a value larger than 255!
dydt .......... The value of for f(t, y) is returned in this argument
t ............. The time (only for deq_iface_wt)
y ............. Values for the dependent variables
param ......... Constant parameters

This function should return the value for \( \mathbf{f}(t, \mathbf{y}) \) in dydt. The value of status should be non-positive, \((-\infty, 0]\), if everything worked, and a value between 1 and 255 inclusive, \([1, 255]\), if something went wrong. This value will be passed back via the status argument of higher level routines to indicate an error condition.

adaptive_steps() uses traditional adaptive step size

• This solver is similar to solvers found in other ODE packages.
• It provides somewhat finer control over step size adjustment than is typical.
• Programmable step processing is used to steer the algorithm instead of directly providing solution points.
• A programmable bisection option to solve for interesting \(\Delta{t}\) values.
• End precisely on a \(t\) value or when a maximum \(t\) is exceeded.

fixed_t_steps() uses fixed \(\Delta{t}\) steps

• With most modern ODE packages, this would be done with a "dense output" mode using embedded interpolation.
• This is a good place to start when writing a custom solver.
• Common use cases:
  • Curve evolution animations that naturally display velocity.

fixed_y_steps() uses fixed \(\Delta{\mathbf{y}}\)-space steps

• Produce solution points separated by fixed distance in \(\Delta{\mathbf{y}}\)-space.
• The components of \(\Delta{\mathbf{y}}\) used for length computations may be specified.
• Overall curve length is computed allowing exit options based on total solution length.
• This is a good place to start when writing a custom solver with a bisection step.
• With most modern ODE packages, this would be done with interpolation. Usually in combination with a step processing routine.
• This could be achieved in MRKISS with the programmable step processing and bisection features of adaptive_steps().
• Common use cases:
  • Parametric line/dot plots
  • Tube plots and sphere sweeps.

sloppy_fixed_y_steps() adjusts \(\Delta{t}\) to approximate \(\Delta{\mathbf{y}}\)-space steps

• Much like fixed_y_steps() in practice but faster and less precise.
• Computes a test step and then adjusts \(\Delta{t}\) in proportion to the ratio of the desired length vs test length.
• This method has no guarantee for correctness, but generally works well in practice.
• Can operate on every step or just steps that are too long. This allows for two handy behaviors:
  • A cheap kind of steps size error control
  • A way to bound the overall length of \(\Delta{\mathbf{y}}\)
• Common use cases:
  • fixed_t_steps() results have a few steps that are "too long"
  • Additional error control for systems with one or two "hot" coordinates.

fixed_t_steps_between() takes a set of \(t\) values at which to find solutions

• Uses several fixed Runge-Kutta steps to go from each \(t\) value to the next.
• Common use cases:
  • Produce a higher accuracy solution based on a previous solution.
  • Physical problems requiring solutions at specific time points.

interpolate_solution interpolates an existing solution to a new solution

• Not a Runge-Kutta method; however, it provides \(\mathcal{O}(3)\) accuracy when using Hermite interpolation.
• Derivative values are freshly computed for the interpolated points
• Common use cases:
  • Produce smooch plots from distant solution points.
  • To "line up" two solutions to the same set of \(t\) values.

Behind all of the above high level solvers are single step routines to carry out the step calculations. These are handy for creating DIY solvers.

• Generic RK methods:
  • step_one() – Use the first RK method in a tableau, returning \(\Delta{\mathbf{y}}\) as a vector.
  • step_all() – Use every RK method in a tableau, returning \(\Delta{\mathbf{y}}\) values as columns in a matrix .
• Hardwired RK methods (used mostly for unit tests):
  • step_one_rk4(): Corresponds to mrkiss_erk_kutta_4. Returns a single \(\Delta{\mathbf{y}}\) as a vector.
  • step_one_rkf45(): Corresponds to mrkiss_eerk_fehlberg_4_5, Returns a two \(\Delta{\mathbf{y}}\) values as columns of a matrix.
  • step_one_dp54(): Corresponds to mrkiss_eerk_dormand_prince_5_4. Returns a two \(\Delta{\mathbf{y}}\) values as columns of a matrix.

MRKISS provides a few utilities:

• Output
  • print_solution() Print a solution to a file or STDOUT.
  • print_istats() Print an istats array to a file or STDOUT.
• Miscellaneous
  • analyze_solution() Compute statistics related to the solution and print the results to STDOUT.
  • seq() Compute a fixed delta sequence of values in the same way steps are computed in fixed_t_steps()
• Status Codes
  • status_to_origin() Returns the subroutine or interface name assigned the given status code
  • status_to_message() Returns the error message for the given status code

If you are interested playing around with MRKISS as quickly as possible, then this section is for you.

All of the code is in the module source files with no external dependencies at all. So you just need to call the modules from your code, and then compile/link everything together.

You can do that by just listing all the source files on the command line with most Fortran compilers. For example, we could compile the lorenz.f90 example in the MRKISS/examples/ directly like this:

cd examples
gfortran.exe lorenz.f90 ../src/*.f90

That said, most people will probably want to use a build system. If GNU Make is your thing, then the files in the /MRKISS/make_include/ directory may be of help. In particular the makefile fragment include.mk provides useful targets and variables. The makefile in the MRKISS/examples/ directory is a good guide on how to use include.mk. In essence you do the following in your makefile:

1. Set MRKISS_PATH in your makefile to the path of the MRKISS source directory – that's the one with the include.mk file.
2. Set FC, FFLAGS, & AR if necessary – most of the time you can use the defaults.
3. Include the "include.mk" file in the MRKISS source directory.
4. Add a build rule for your program.

Your makefile will look something like this:

MRKISS_PATH = ../MRKISS

# Set FC, FFLAGS, & AR here.  The include below has the settings I use on my system.
include $(MRKISS_PATH)/tools_gfortran.mk

include $(MRKISS_PATH)/include.mk

your_program : your_program.f90 $(MRKISS_OBJ_FILES)
    $(FC) $(FFLAGS) $^ -o $@

Note the rule for your_program in the makefile above takes the lazy approach of adding every MRKISS module as a dependency regardless of if your program actually needs them all. This is how most people use the modules because it's simple. The cost might be a couple seconds of extra compile time. You can explicitly list out the modules in the makefile if you wish. Such a rule might look like the following:

your_program : your_program.f90 mrkiss_config$(OBJ_SUFFIX) mrkiss_solvers_wt(OBJ_SUFFIX) mrkiss_utils$(OBJ_SUFFIX)
    $(FC) $(FFLAGS) $^ -o $@

This section is about how I test MRKISS.

The MRKISS/tests/ directory contains code I primary use for testing MRKISS while the MRKISS/examples/ directory contains code I primarily use to demonstrate how to use MRKISS. The difference between a "test" and an "example" in MRKISS is a little bit slippery. Some of the tests, like tc1_* and tc2_*, could be considered demonstrations. In addition, I use all of the code in MRKISS/examples/ for tests.

The tests can be run by changing into the appropriate directory (MRKISS/tests/ or MRKISS/examples/), and building the make target tests. For example:

cd tests
make -j 16 tests

Note the -j 16 argument to make. When running all the tests, especially in the MRKISS/tests/ directory, I strongly recommend running in parallel.

In addition, the make files in MRKISS/tests/ and MRKISS/examples/ have numerous additional targets to run various individual tests or subsets of the test suite. These targets are documented in the subsections below.