O2scl User's Guide

0.805

Download O2scl

O2scl is a C++ class library for object-oriented numerical programming. It includes

See licensing information at License Information.

Download O2scl      PDF documentation

Quick Reference to User's Guide

Installation

The rules for installation are generally the same as that for other GNU libraries. The file INSTALL has some details on this procedure. Generally, you should be able to run ./configure and then type make and make install. More information on the configure command can also be obtained from ./configure --help. O2scl requires the GSL libraries. If the configure script cannot find them, you may have to specify their location in the CPPFLAGS and LDFLAGS environment variables (./configure --help shows some information on this). The documentation is automatically installed by make install.

After make install, you may test the library with make o2scl-test.

This library requires GSL and is designed to work with GSL versions 1.9 or 1.10. Some classes may work with older versions of GSL, but this cannot be guaranteed.

Range-checking for vectors and matrices is performed similar to the GSL approach, and is turned on by default. You can disable range-checking by defining -DO2SCL_NO_RANGE_CHECK 

    CPPFLAGS="-DO2SCL_NO_RANGE_CHECK" ./configure

The separate libraries o2scl_eos and o2scl_part are installed by default. To disable the installation of these libraries and their associated documentation, run configure with the flags --disable-eoslib or --disable-partlib. Note that o2scl_eos depends on o2scl_part so using --disable-partlib without --disable-eoslib may not work.

There are several warning flags that are useful when configuring and compiling with g++. (See the GSL documentation for an excellent discussion.) For running configure, I use something like 

    CFLAGS="" CXXFLAGS="" CPPFLAGS="-ansi -pedantic -Wall -W           \
    -Wconversion -Wno-unused -Wshadow -Wpointer-arith -Wcast-align     \
    -Wwrite-strings -fshort-enums -ggdb -O3 -DGSL_RANGE_CHECK=1        \
    -DHAVE_INLINE -DO2SCL_ARRAY_ABORT                                  \
    -I/home/asteiner/install/include"                                  \
    LDFLAGS="-L/home/asteiner/install/lib" ./configure -C              \
    --prefix=/home/asteiner/install
In this example, specifying -I/home/asteiner/install/include and -L/home/asteiner/install/lib above ensures that the GSL libraries can be found (this is where they are installed on my machine). The --prefix=/home/asteiner/install argument to ./configure ensures that O2scl is installed there as well.

with the references to the directory /home/asteiner/install in order to install somewhere in my user directory, and because my copy of GSL is installed there, rather than in /usr/local.

The documentation is generated with Doxygen. In principle, the documentation can be regenerated by the end-user, but this is not supported and may require several external applications not included in the distribution.

Un-installation: While there is no explicit "uninstall" procedure, there are only a couple places to check. Installation creates directories named o2scl in the include, doc and shared files directory (which default to /usr/local/include, /usr/local/doc, and /usr/local/share) which can be removed. Finally, all of the libraries are named with the prefix libo2scl and are created by default in /usr/local/lib. As configured with the settings above, the files are in /home/asteiner/install/include/o2scl, /home/asteiner/install/lib, /home/asteiner/install/share/o2scl, and /home/asteiner/install/doc/o2scl.

General Usage

Namespaces

Most of the classes reside in the namespace o2scl (removed from the documentation). Numerical constants (many of them based on the GSL constants) are placed in separate namespaces (gsl_cgs, gsl_cgsm, gsl_mks, gsl_mksa, gsl_num, o2scl_const, and o2scl_fm). The CBLAS functions are in the o2scl_cblas namespace and vector and complex number arithmetic is in o2scl_arith. There are also two namespaces which hold integration coefficients, o2scl_inte_qag_coeffs and o2scl_inte_qng_coeffs.

Documentation conventions

In the following documentation, function parameters are denoted by parameter, except when used in mathematical formulas as in $ \mathrm{variable} $ .

Nomenclature

Classes directly derived from the GNU Scientific Library are preceeded by the prefix gsl_ and classes derived from CERNLIB are preceeded by the prefix cern_. Some of those classes derived from GSL and CERNLIB operate slightly differently from the original versions. The differences are detailed in the corresponding class documentation.

Error handling

Error handling is GSL-like with functions returning 0 for success and calling a GSL-like error handler. The error handler, err_hnd is a global pointer to an object of type err_class.

The default behaviour for all errors is simply store the error code and information. You can require the default error handler to print out any error (or even exit) using err_class::set_mode(). Also, if O2SCL_ARRAY_ABORT is defined, then exit() will be called any time a gsl_index error is called, e.g. an out-of-bounds array access.

Errors can be set through the macros set_err, set_err_ret, add_err, and add_err_ret, which are defined in the file err_hnd.h.

Functionality similar to assert() is provided with the macro err_assert, which exits if its argument is non-zero, and bool_assert which exits if its argument is false.

Note that the library functions do not typically reset the error handler. For this reason the user may need to reset the error handler before calling functions which do not return an integer error value so that they can easily test to see if an error occured, e.g. using err_hnd::get_errno().

The default error handler can be replaced by simply assigning the address of a descendant of err_class to err_hnd.

Library dependencies

All of the libraries need libo2scl_base, libgsl, and either libgslcblas or some externally-specified libcblas to operate. Other additional dependencies are listed below:

Compiling examples

A few example programs are in the examples directory. After installation, they can be compiled and executed by running make o2scl-examples in that directory. This will also test the output of the examples to make sure it is correct and summarize these tests in examples-summary.txt. The output for each example is placed in the corresponding file with a .scr extension.

Alternatively, you can make the executable for each example in the examples directory individually using, e.g. make ex_mroot.

See Example source code for the documented source code of the individual examples

Also, the testing code for each class is sometimes useful for providing examples of their usage. The testing source code for each source file is named with an _ts.cpp prefix in the same directory as the class source.

Related projects

Several noteworthy related projects:

Complex Numbers

Some rudimentary arithmetic operators for gsl_complex are defined in cx_arith.h, but no constructor has been written. The object gsl_complex is still a struct, not a class. For example, 
    gsl_complex a={{1,2}}, b={{3,4}};
    gsl_complex c=a+b;
    cout << GSL_REAL(c) << " " << GSL_IMAG(C) << endl;
In case the user needs to convert between gsl_complex and std::complex<double>, two conversion functions gsl_to_complex() and complex_to_gsl() are provided in ovector_cx_tlate.h.

Arrays, Vectors, Matrices and Tensors

Introduction

The O2scl library uses a standard nomenclature to distinguish a couple different concepts. The word "array" is always used to refer to C-style arrays, i.e. double[]. If there are two dimensions in the array, it is a "two-dimensional array", i.e. double[][] . The word "vector" is reserved generic objects with array-like semantics, i.e. any type of object or class which can be treated similar to an array in that it has an function operator[] which gives array-like indexing. Thus arrays are vectors, but not all vectors are arrays. There are a couple specific vector types defined in O2scl like ovector and uvector . The STL vector std::vector<double> is also a vector type in this language. The word "matrix" is reserved for the a generic object which has matrix-like semantics and can be accessed using either operator() or two successive applications of operator[] (or sometimes both). The O2scl objects omatrix and umatrix are matrix objects, as is a C-style two-dimensional array, double[][] . The header files are named in this spirit also: functions and classes which operate on generic vector types are in vector.h and functions and classes which work only with C-style arrays are in array.h . The word "tensor" is used for a generic object which has rank n and then has n associated indices. A vector is just a tensor of rank 1 and a matrix is just a tensor of rank 2.

Most of the classes in O2scl which use vectors and/or matrices are designed so that they can be used with any appropriately-defined vector or matrix types. This is a major part of the design goals for O2scl and most of the classes are compatible with matrix and vector objects from GSL, TNT, MV++, uBlas, and Blitz++.

The first index of a matrix type is defined always to be the index associated with "rows" and the second is associated with "columns". The O2scl matrix output functions respect this notation as well, so that all of the elements of the first row is sent to the screen, then all of the elements in the second row, and so on. With this in mind, one can make the distinction between "row-major" and "column-major" matrix storage. C-style two-dimensional arrays are "row-major" in that the elements of the first row occupy the first locations in memory as opposed "column-major" where the first column occupies the first locations in memory. It is important to note that the majority of the classes in O2scl do not care about the details of the underlying memory structure, so long as two successive applications of operator[] (or in some cases operator(,) ) works properly. The storage format used by omatrix and umatrix is row-major, and there are no column-major matrix classes in O2scl (yet).

Matrix indexing with [][] vs. (,)

While vector indexing with operator[] is very compatible with almost any vector type, matrix indexing is a bit more difficult. There are two options: assume matrix objects provide an operator[] method which can be applied twice, i.e. m[i][j], or assume that matrix elements should be referred to with m(i,j). Most of the O2scl classes use the former approach so that they are also compatible with two-dimensional arrays. However, there are sometimes good reasons to want to use operator() for matrix-intensive operations from linear algebra. For this reason, some of the functions given in the linalg directory are specified in two forms: the first default form which assumes [][], and the second form with the same name, but in a namespace which has a suffix _paren and assumes matrix types use (,).

O2scl matrix and vector types

The rest of this section in the User's guide is dedicated to describing the O2scl implementations of vector, matrix, and tensor types.

Vectors and matrices are designed using the templates ovector_tlate and omatrix_tlate, which are compatible with gsl_vector and gsl_matrix. Vectors and matrices with unit stride are provided in uvector_tlate and umatrix_tlate. The most commonly used double-precision versions of these template classes are ovector, omatrix, uvector and umatrix, and their associated "views" (analogous to GSL vector and matrix views) which are named with a _view suffix. The classes ovector_tlate and omatrix_tlate offer the syntactic simplicity of array-like indexing, which is easy to apply to vectors and matrices created with GSL.

The following sections primarily discuss the operation objects of type ovector. The generalizations to objects of type uvector, omatrix, and the complex vector and matrix objects ovector_cx, omatrix_cx, uvector_cx, and umatrix_cx are straightforward.

Views

Vector and matrix views are provided as parents of the vector and matrix classes which do not have methods for memory allocation. As in GSL, it is simple to "view" normal C-style arrays or pointer arrays (see ovector_array), parts of vectors (ovector_subvector), and rows and columns of matrices (omatrix_row and omatrix_col). Several operations are defined, including addition, subtraction, and dot products.

Typedefs

Several typedefs are used to give smaller names to often used templates. Vectors and matrices of double-precision numbers all begin with the prefixes ovector and omatrix. Integer versions begin with the prefixes ovector_int and omatrix_int. Complex versions have an additional "_cx", e.g. ovector_cx. See ovector_tlate.h, ovector_cx_tlate.h, omatrix_tlate.h, and omatrix_cx_tlate.h.

Unit-stride vectors

The uvector_tlate objects are naturally somewhat faster albeit less flexible than their finite-stride counterparts. Conversion to GSL vectors and matrices is not trivial for uvector_tlate objects, but demands copying the entire vector. Vector views, operators, and the nomenclature is similar to that of ovector.

Memory allocation

Memory for vectors can be allocated using ovector::allocate() and ovector::free(). Allocation can also be performed by the constructor, and the destructor automatically calls free() if necessary. In contrast to gsl_vector_alloc(), ovector::allocate() will call ovector::free(), if necessary, to free previously allocated space. Allocating memory does not clear the recently allocated memory to zero. You can use ovector::set_all() with an argument of 0.0 to clear a vector (and similarly for a matrix).

If the memory allocation fails, either in the constructor or in allocate(), then the error handler will be called, partially allocated memory will be freed, and the size will be reset to zero. You can either use the error handler or test to see if the allocation succeeded by examining the size argument, e.g. 

    const size_t n=10;
    ovector x(10);
    if (x.size()==0) cout << "Failed." << endl;

Get and set

Vectors and matrices can be modified using ovector::get() and ovector::set() methods analogous to gsl_vector_get() and gsl_vector_set(), or they can be modified through ovector::operator[] (or ovector::operator() ), e.g. 

    ovector a(4);
    a.set(0,2.0);
    a.set(1,3.0);
    a[2]=4.0;
    a[3]=2.0*a[1];

If you want to set all of the values in an ovector or an omatrix at the same time, then use ovector::set_all() or omatrix::set_all().

Range checking

Range checking is performed depending on whether or not O2SCL_NO_RANGE_CHECK is defined. It can be defined in the arguments to ./configure upon installation to turn off range checking. Note that this is completely separate from the GSL range checking mechanism, so range checking may be on in O2scl even if it has been turned off in GSL. Range checking is used primarily in the vector, matrix, and tensor get() and set() methods.

To see if range checking was turned on during installation (without calling the error handler), use lib_settings_class::range_check().

Note that range checking in O2scl code is most often present in header files, rather than in source code. This means that range checking can be turned on or off in user-defined functions separately from whether or not it was used in the library classes and functions.

Shallow and deep copy

Copying O2scl vectors using constructors or the = operator is performed according to what kind of object is on the left-hand side (LHS) of the equals sign. If the LHS is a view, then a shallow copy is performed, and if the LHS is a ovector, then a deep copy is performed. If an attempt is made to perform a deep copy onto a vector that has already been allocated, then that previously allocated memory is automatically freed. The user must be careful to ensure that information is not lost this way, even though no memory leak will occur.

For generic deep vector and matrix copying, you can use the template functions vector_copy(), matrix_copy(). These would allow you, for example, to copy an ovector to a std::vector<double> object (assuming the memory allocation has already been taken care of).

Vector and matrix arithmetic

Several operators are available as member functions of the corresponding template:

Vector_view unary operators:

Matrix_view unary operators:

Binary operators like addition, subtraction, and matrix multiplication are also defined for ovector, uvector, and related objects in the o2scl_arith namespace. The generic template for a binary operator, e.g. 

    template<class vec_t> vec_t &operator+(vec_t &v1, vec_t &v2);
is difficult because the compiler has no way of distinguishing vector and non-vector classes. At the moment, this is solved by creating a define macro for the binary operators. In addition to the predefined operators for native classes, the user may also define binary operators for other classes using the same macros. For example, 
    O2SCL_OP_VEC_VEC_ADD(o2scl::ovector,std::vector<double>,
    std::vector<double>)
would provide an addition operator for ovector and vectors from the Standard Template Library. The macros are detailed in vec_arith.h.

The GSL BLAS routines can also be used directly with ovector and omatrix objects.

Note that some of these arithmetic operations succeed even with non-matching vector and matrix sizes. For example, adding a 3x3 matrix to a 4x4 matrix will result in a 3x3 matrix and the 7 outer elements of the 4x4 matrix are ignored.

Converting to and from GSL forms

Because of the way ovector is constructed, you may use type conversion to convert to and from objects of type gsl_vector. 

    ovector a(2);
    a[0]=1.0;
    a[1]=2.0;
    gsl_vector *g=(gsl_vector *)(&a);
    cout << gsl_vector_get(g,0) << " " << gsl_vector_get(g,1) << endl;
Or, 
    gsl_vector *g=gsl_vector_alloc(2);
    gsl_vector_set(0,1.0);
    gsl_vector_set(1,2.0);
    ovector &a=(ovector &)(*g);
    cout << a[0] << " " << a[1] << endl;
This sort of type-casting is discouraged among unrelated classes, but is permissible here because ovector_tlate is a descendant of gsl_vector. In particular, this will not generate "type-punning" warnings in later gcc versions. If this bothers your sensibilities, however, then you can use the following approach: 
    ovector a(2);
    gsl_vector *g=a.get_gsl_vector();
The ease of converting between these two kind of objects makes it easy to use gsl functions on objects of type ovector, i.e. 
    ovector a(2);
    a[0]=2.0;
    a[1]=1.0;
    gsl_vector_sort((gsl_vector *)(&a));
    cout << a[0] << " " << a[1] << endl;

Converting from STL form

To "view" a std::vector<double>, you can use ovector_array 

    std::vector<double> d;
    d.push_back(1.0);
    d.push_back(3.0);
    ovector_array aa(d.size,&(d[0]));
    cout << aa[0] << " " << aa[1] << endl;

However, you should note that if the memory for the std::vector is reallocated (for example because of a call to push_back()), then a previously created ovector_view will be incorrect.

push_back() and pop() methods

These two functions give a behavior similar to the corresponding methods for std::vector<>. This will work in O2scl classes, but may not be compatible with all of the GSL functions. This will break if the address of a ovector_tlate is given to a GSL function which accesses the block->size parameter instead of the size parameter of a gsl_vector. Please contact the author of O2scl if you find a GSL function with this behavior.

Views

Views are slightly different than in GSL in that they are now implemented as parent classes. The code 

    double x[2]={1.0,2.0};
    gsl_vector_view_array v(2,x);
    gsl_vector *g=&(v.vector);
    gsl_vector_set(g,0,3.0);
    cout << gsl_vector_get(g,0) << " " << gsl_vector_get(g,1) << endl;
can be replaced by 
    double x[2]={1.0,2.0};
    ovector_array a(2,x);
    a[0]=3.0;
    cout << a[0] << " " << a[1] << endl;

Passing ovector parameters

It is often best to pass an ovector as a const reference to an ovector_view, i.e. 

    void function(const ovector_view &a);
If the function may change the values in the ovector, then just leave out const 
    void function(ovector_view &a);
This way, you ensure that the function is not allowed to modify the memory for the vector argument.

If you intend for a function (rather than the user) to handle the memory allocation, then some care is necessary. The following code 

class my_class {
  int afunction(ovector &a) {
    a.allocate(1);
    // do something with a
    return 0;
  }
};
is confusing because the user may have already allocated memory for a. To avoid this, you may want to ensure that the user sends an empty vector. For example, 
class my_class {
  int afunction(ovector &a) {
    if (a.get_size()>0 && a.is_owner()==true) {
      set_err("Unallocated vector not sent.",1);
      return 1;
    } else {
      a.allocate(1);
      // do something with a
      return 0;
    }
  }
};
In lieu of this, it is often preferable to use a local variable for the storage and offer a get() function, 
class my_class {
protected:
  ovector a;
public:
  int afunction() {
    a.allocate(1);
    // do something with a
    return 0;
  }
  int get_result(const ovector_view &av) { av=a; return 0; }
};
The O2scl classes run into this situation quite frequently, but the vector type is specified through a template 
template<class vec_t> class my_class {
protected:
  vec_t a;
public:
  int afunction(vec_t &a) {
    a.allocate(1);
    // do something with a
    return 0;
  }
};

Vectors and operator=()

An "operator=(value)" method for setting all vector elements to the same value is not included because it leads to confusion between, ovector_tlate::operator=(const data_t &val) and ovector_tlate::ovector_tlate(size_t val) For example, after implementing operator=() and executing the following 

    ovector_int o1=2;
    ovector_int o2;
    o2=2;
o1 will be a vector of size two, and o2 will be an empty vector!

To set all of the vector elements to the same value, use ovector_tlate::set_all(). Because of the existence of constructors like ovector_tlate::ovector_tlate(size_t val), the following code 

    ovector_int o1=2;
still compiles, and is equivalent to 
    ovector_int o1(2);
while the code 
    ovector_int o1;
    o1=2;
will not compile. As a matter of style, ovector_int o1(2); is preferable to ovector_int o1=2; to avoid confusion.

Matrix structure

The matrices from omatrix_tlate are structured in exactly the same way as in GSL. For a matrix with 2 rows, 4 colums, and a "tda" or "trailing dimension" of 7, the memory for the matrix is structured in the following way:

    00 01 02 03 XX XX XX
    10 11 12 13 XX XX XX

where XX indicates portions of memory that are unused. The tda can be accessed through, for example, the method omatrix_view_tlate::tda(). The get(size_t,size_t) methods always take the row index as the first argument and the column index as the second argument. The matrices from umatrix_tlate have a trailing dimension which is always equal to the number of columns.

Reversing the order of vectors

You can get a reversed vector view from ovector_reverse_tlate, or uvector_reverse_tlate. For these classes, operator[] and related methods are redefined to perform the reversal. If you want to make many calls to these indexing methods for a reversed vector, then simply copying the vector to a reversed version may be faster.

Const-correctness with vectors

There are several classes named with "_const" to provide different kinds of const views of const vectors. The keyword const still ought to be included to ensure that the object is treated properly. For example, 

    ovector o(2);
    o[0]=3.0;
    o[1]=-1.0;
    const ovector_const_subvector ocs(o,1,1);

At present, const-correctness in O2scl can be improperly removed, if the const keyword is not properly included. For example, the following code will compile, violated the const-correctness of the ocs variable. 

    ovector o(2);
    o[0]=3.0;
    o[1]=-1.0;
    ovector_const_subvector ocs(o,1,1);
    ovector_view ov(ocs);
    ov[0]=2.0;

Tensors

Some preliminary support is provided for tensors of arbitrary rank and size in the class tensor. Classes tensor1, tensor2, tensor3, and tensor4 are rank-specific versions for 1-, 2-, 3- and 4-rank tensors. For n-dimsional data defined on a grid, tensor_grid provides a space to define a hyper-cubic grid in addition to the the tensor data. This class tensor_grid also provides n-dimensional interpolation of the data defined on the specified grid.

Permutations

Permutations are implemented through the permutation class. This class is fully compatible with gsl_permutation objects since it is inherited from gsl_permutation_struct. The class also contains no new data members, so upcasting and downcasting can always be performed. It is perfectly permissible to call GSL permutation functions from permutation objects by simply passing the address of the permutation, i.e. 
    permutation p(4);
    p.init();
    gsl_permutation_swap(&p,2,3);

The functions which apply a permutation to a user-specified vector are member template functions in the permutation class (see permutation::apply() ).

Memory allocation/deallocation between the class and the gsl_struct is compatible in many cases, but mixing these forms is strongly discouraged, i.e. avoid using gsl_permutation_alloc() on a permutation object, but rather use permutation::allocate() instead. The use of permutation::free() is encouraged, but any remaining memory is deallocated in the object destructor.

Linear algebra

There is a small set of linear algebra routines. These are not intended to be a replacement for LAPACK, but are designed for use by O2scl routines so that they work for generic matrix and vector types. For vector and matrix types using operator[], the BLAS and linear algebra routines are inside the o2scl_cblas and o2scl_linalg namespaces. For vector and matrix types using operator(), the BLAS and linear algebra routines routines are inside the o2scl_cblas_paren and o2scl_linalg_paren namespaces.

The linear algebra classes and functions include:

Interpolation

The classes o2scl_interp and o2scl_interp_vec allow basic interpolation, lookup, differentiation, and integration of data given in two ovectors or ovector views. In contrast to the GSL routines, data which is presented with a decreasing independent variable is handled automatically. For interpolation with arrays rather than ovectors, use array_interp or array_interp_vec.

For fast interpolation of arrays where the independent variable is strictly increasing and without error-checking, you can directly use the children of base_interp.

The two interpolation interfaces

The difference between the two classes, o2scl_interp and o2scl_interp_vec, analogous to the difference between using gsl_interp_eval() and gsl_spline_eval() in GSL. You can create a o2scl_interp object and use it to interpolate among any pair of chosen vectors, i.e. 

    ovector x(20), y(20);
    // fill x and y with data
    o2scl_interp oi;
    double y_o2sclf=oi.interp(0.5,20,x,y);
Alternatively, you can create a o2scl_interp_vec object which can be optimized for a pair of vectors that you specify in advance 
    ovector x(20), y(20);
    // fill x and y with data
    o2scl_interp_vec oi(20,x,y);
    double y_o2sclf=oi.interp(0.5);

Lookup and binary search

The class search_vec contains a searching functions for objects of type ovector which are monotonic. Note that if you want to find the index of an ovector where a particular value is located without any assumptions with regard to the ordering, you can use ovector::lookup() which performs an exhaustive search.

"Smart" interpolation

The classes smart_interp and smart_interp_vec allow interpolation, lookup, differentiation, and integration of data which is non-monotonic or multiply-valued outside the region of interest. As with o2scl_interp above, the corresponding array versions are given in sma_interp and sma_interp_vec.

Two and higher dimensional interpolation

Preliminary support for two-dimensional interpolation is given in twod_intp, and n-dimensional interpolation in tensor_grid.

Physical constants

The constants from GSL are reworked with the type const double and placed in namespaces called gsl_cgs, gsl_cgsm, gsl_mks, gsl_mksa, and gsl_num. Some additional constants are given in the namespace o2scl_const. Some of the numerical values have been updated from recently released data from NIST.

Function Objects

Functions are passed to numerical routines using template-based function classes. There are several basic kinds of function objects:

For each of these classes (except funct), there is a version named _vfunct instead of _funct which is designed to be used with C-style arrays instead.

The class name suffixes denote children of a generic function type which are created using different kinds of inputs:

There is a small overhead associated with the indirection: a "user class" accesses the function class which then calls function which was specified in the constructor of the function class. In many problems, the overhead associated with the indirection is small. Some of this overhead can always be avoided by inheriting directly from the function class and thus the user class will make a direct virtual function call. To eliminate the overhead entirely, one can specifying a new type for the template parameter in the user class.

Note that virtual functions can be specified through this mechanism as well. For example, if cern_mroot is used to solve a set of equations specified as 

class my_type_t {
  virtual member_func();
};
my_type_t my_instance;
class my_derived_type_t : public my_type_t {
  virtual member_func();
};
my_derived_type_t my_inst2;
mm_funct_mfptr<my_type_t> func(&my_inst2,&my_instance::member_func);
Then the solver will solve the member function in the derived type, not the parent type.

Data tables

The class table is a container to hold and perform operations on related columns of data. It supports column operations, interpolation, column reference by either name or index, binary searching (in the case of ordered columns), sorting, and fitting two columns to a user-specified function.

String manipulation

There are a couple classes and functions to help manipulate strings of text. Conversion routines for std::string objects are given in string_conv.h and include

See also size_of_exponent(), double_to_latex(), double_to_html, and double_to_ieee_string().

There is a class called columnify, which converts a set of strings into nicely formatted columns by padding with the necessary amount of spaces. This class operates on string objects of type std::string, and also works will for formatting columns of floating-point numbers. This class is used to provide output for matrices in the functions matrix_out(), matrix_out_paren(), and matrix_cx_out_paren(). For output of vectors, see vector_out() in array.h.

A related function, screenify(), reformats a column of strings into many columns stored row-by-row in a new string array. It operates very similar to the way the classic Unix command ls organizes files and directories in multiple columns in order to save screen space.

The function count_words() counts the number of "words" in a string, which are delimited by whitespace.

Differentiation

Differentiation is performed by descendants of deriv and the classes are provided. These allow one to calculate either first, second, and third derivatives. The GSL approach is used in gsl_deriv, and the CERNLIB routine is used in cern_deriv. Both of these compute derivatives for a function specified using a descendant of funct. For functions which are tabulated over equally-spaced abscissas, the class eqi_deriv is provided which applies the formulas from Abramowitz and Stegun at a specified order.

Warning: For gsl_deriv and cern_deriv, the second and third derivatives are calculated by naive repeated application of the code for the first derivative and can be particularly troublesome if the function is not sufficiently smooth. Error estimation is also incorrect for second and third derivatives.

Integration

Integration is performed by descendants of inte and is provided in the library o2scl_inte.

There are several routines for one-dimensional integration.

For the GSL-based integration routines, the variables inte::tolx and inte::tolf have the same role as the quantities usually denoted in the GSL integration routines by epsabs and epsrel. In particular, the integration classes attempt to ensure that

\[ |\mathrm{result}-I| \leq \mathrm{Max}(\mathrm{tolx}, \mathrm{tolf}|I|) \]

and returns an error to attempt to ensure that

\[ |\mathrm{result}-I| \leq \mathrm{abserr} \leq \mathrm{Max}(\mathrm{tolx},\mathrm{tolf}|I|) \]

where I is the integral to be evaluated. Even when the corresponding descendant of inte::integ() returns success, these inequalities may fail for sufficiently difficult functions. All of the GSL integration routines except for gsl_inte_qng use a workspace given in gsl_inte_table which holds the results of the various subdivisions of the original interval. The size of this workspace (and thus then number of subdivisions) can be controlled with gsl_inte_table::set_wkspace().

The GSL routines were originally based on QUADPACK, which is available at http://www.netlib.org/quadpack .

Multi-dimensional hypercubic integration is performed by composite_inte, the sole descendant of multi_inte. composite_inte allows you to specify a set of one-dimensional integration routines (objects of type inte) and apply them to a multi-dimensional problem.

General multi-dimensional integration is performed by comp_gen_inte, the sole descendant of gen_inte. The user is allowed to specify a upper and lower limits which are functions of the variables for integrations which have not yet been performed, i.e. the n-dimensional integral

\[ \int_{x_0=a_0}^{x_0=b_0} f(x_0) \int_{x_1=a_1(x_0)}^{x_1=b_1(x_0)} f(x_0, x_1) ... \int_{x_{\mathrm{n}-1}=a_{\mathrm{n}-1}(x_0,x_1,..,x_{\mathrm{n}-2})}^ {x_{\mathrm{n}-1}=b_{\mathrm{n}-1}(x_0,x_1,..,x_{\mathrm{n}-2})} f(x_0,x_1,...,x_{\mathrm{n-1}})~d x_{\mathrm{n}-1}~...~d x_1~d x_0 \]

Again, one specifies a set of inte objects to apply to each variable to be integrated over.

Monte Carlo integration is also provided (see Monte Carlo Integration).

Roots of Polynomials

Classes are provided for solving quadratic, cubic, and quartic equations as well as general polynomials. There is a standard nomenclature: classes which handle polynomials with real coefficients and real roots end with the suffix "_real" (quadratic_real, cubic_real and quartic_real), classes which handle real coefficients and complex roots end with the suffix "_real_coeff" (quadratic_real_coeff, cubic_real_coeff, quartic_real_coeff, and poly_real_coeff), and classes which handle complex polynomials with complex coefficients (quadratic_complex, cubic_complex, quartic_complex, and poly_complex). As a reminder, complex roots may not occur in conjugate pairs if the coefficients are not real. Most of these routines do not separately handle cases where the leading coefficient is zero.

In the public interfaces to the polynomial solvers, the complex type std::complex<double> is used. These can be converted to and from the GSL complex type using the complex_to_gsl() and gsl_to_complex() functions.

At present, the polynomial routines work with complex numbers as objects of type std::complex<double> and are located in library o2scl_other.

For quadratics, gsl_quadratic_real_coeff is the best if the coefficients are real, while if the coefficients are complex, use quadratic_std_complex. For cubics with real coefficients, cern_cubic_real_coeff is the best, while if the coefficients are complex, use cubic_std_complex.

For a quartic polynomial with real coefficients, cern_quartic_real_coeff is the best, unless the coefficients of odd powers happen to be small, in which case, gsl_quartic_real2 tends to work better. For quartics, generic polynomial solvers such as gsl_poly_real_coeff can provide more accurate (but slower) results. If the coefficients are complex, then you can use simple_quartic_complex.

Equation Solving

Equation solving classes are stored in the library o2scl_root.

One-dimensional solvers

Solution of one equation in one variable is accomplished by children of the class root. This base class provides the structure for three different solving methods:

For one-dimensional solving, use cern_root or gsl_root_brent if you have the root bracketed, or gsl_root_stef if you have the derivative available. If you have neither a bracket or a derivative, you can use cern_mroot_root.

If not all of these three functions are overloaded, then the source code in root is designed to try to automatically provide the solution using the remaining functions. Most of the one-dimensional solving routines, in their original form, are written in the second or third form above. For example, gsl_root_brent is originally a bracketing routine of the form root::solve_bkt(), but calls to either root::solve() or root::solve_de() will attempt to automatically bracket the function given the initial guess that is provided. Also, gsl_root_stef is a "root-polishing" routine given derivatives of the form root::solve_de(). If either root::solve() or root::solve_bkt() are called, then root::solve_de() will be called with finite-differencing used to estimate the derivative. Of course, it is frequently most efficient to use the solver in the way it was intended.

Todo:
Double check this documentation above
Multi-dimensional solvers

Solution of more than one equation is accomplished by descendants of the class mroot. There are two basic functions

For multi-dimensional solving, you can use either cern_mroot or gsl_mroot_hybrids. While cern_mroot does not use user-supplied derivatives, gsl_mroot_hybrids can use user-supplied derivative information (as in the GSL hybridsj method).

Minimization

One-dimensional minimization is performed by descendants of minimize and provided in the library o2scl_minimize. There are two one-dimensional minimization algorithms, cern_minimize and gsl_min_brent, and they are both bracketing algorithms type where an interval and an initial guess must be provided. If only an initial guess and no bracket is given, these two classes will attempt to find a suitable bracket from the initial guess. While the minimize base class is designed to allow future descendants to optionally use derivative information, this is not yet supported for any one-dimensional minimizers.

Multi-dimensional minimization is performed by descendants of multi_min: gsl_mmin_simp, gsl_mmin_conp, gsl_mmin_conf, and gsl_mmin_bfgs2. The class multi_min_fix is a convenient way to perform a minimization while fixing some of the original parameters. The class gsl_mmin_simp does not require or use any derivative information, but the other minimization classes are intended for use when derivatives are available.

Simulated annealing methods are also provided (see Simulated Annealing).

It is important to note that not all of the minimization routines test the second derivative to ensure that it doesn't vanish to ensure that we have found a true minimum.

A naive way of implementing constraints is to add a function to the original which increases the value outside of the allowed region. This can be done with the functions constraint() and lower_bound. There are two analogous functions, cont_constraint() and cont_lower_bound(), which continuous and differentiable versions. Where possible, it is better to use the constrained minimization routines described below.

Constrained Minimization

O2scl reimplements the Open Optimization Library (OOL) available at http://ool.sourceforge.net. The associated classes allow constrained minimization when the constraint can be expressed as a hyper-cubic constraint on all of the independent variables. The routines have been rewritten and reformatted for C++ in order to facilitate the use of member functions and user-defined vector types as arguments. The base class is ool_constr_mmin and there are two different constrained minimzation algorithms implemented in ool_mmin_pgrad, ool_mmin_spg. (The ool_mmin_gencan minimizer is not yet finished). The O2scl implementation should be essentially identical to the most recently released version of OOL.

The constrained minimization classes operate in a similar way to the other multi-dimensional minimization classes (which are derived from multi_min). The constraints are specified with the function 

    ool_constr_mmin::set_constraints(size_t nc, vec_t &lower, 
    vec_t &upper);
and the minimization can be performed by calling either multi_min::mmin() or multi_min::mmin_de() (if the gradient is provided by the user). The "GENCAN" method requires a Hessian vector product and the user can specify this product for the minimization by using ool_constr_mmin::mmin_hess(). The Hessian product function can be specified as an object of type ool_hfunct or ool_hvfunct in a similar way to the other function objects in O2scl .

There are five error codes defined in ool_constr_mmin which are specific to the OOL classes.

Monte Carlo Integration

Monte Carlo integration is performed by descendants of mcarlo_inte in the library o2scl_mcarlo (gsl_monte, gsl_miser, and gsl_vegas). These routines are generally superior to the direct methods for integrals over regions with large numbers of spatial dimensions.

Simulated Annealing

Minimization by simulated annealing is performed by descendants of sim_anneal (see gsl_anneal). The annealing schedule is given by a descendant of tptr_schedule (see tptr_geoseries).

Non-linear Least-Squares Fitting

Fitting is performed by descendants of fit_base and fitting functions can be specifed using fit_funct. The GSL fitting routines (scaled and unscaled) are implemented in gsl_fit. A generic fitting routine using a minimizer object specified as a child of multi_min is implemented in min_fit. When the multi_min object is (for example) a sim_anneal object, min_fit can avoid local minima which can occur when fitting noisy data.

Solution of Ordinary Differential Equations

Classes for non-adaptive integration are provided as descendents of odestep and classes for adaptive integration are descendants of adapt_step. To specify a set of functions to these classes, use a child of ode_funct for a generic vector type or a child of ode_vfunct when using arrays.

Solution of simple initial value problems is performed by ode_iv_solve.

Preliminary support for boundary value problems is given in children of ode_bv_solve.

Random Number Generation

Random number generators are descendants of rnga and are provided in the library o2scl_rnga. While the base object rnga is created to allow user-defined random number generators, the only random number generator presently included are from GSL. The GSL random number generator code is reimplemented in the class gsl_rnga to avoid an additional performance penalty. This may not be a truly "object-oriented" interface in that it does not use virtual functions, but it avoids any possible performance penalty. Random number generators are implemented as templates in sim_anneal and mcarlo_inte. In these classes, the random number generator is a template type, rather than a member data pointer, in order to ensure fast execution.

Two-dimensional Interpolation

Successive use of smart_interp is implemented in twod_intp. Also, see planar_intp and quad_intp and the computation of contour lines in contour. These latter three classes are somewhat experiemental at present.

Other Routines

(These are all experimental)

Fourier transforms - see gsl_fft

Series acceleration - see gsl_series

Chebyshev approximations - see gsl_chebapp

Timing execution - see timer_gettod and timer_clock

Polylogarithms - see polylog

Library settings

There are a couple library settings which are handled by a global object lib_settings of type lib_settings_class.

There are several data files that are used by various classes in the library. The installation procedure should ensure that these files are automatically found. However, if these data files are moved after installation, then a call to lib_settings_class::set_data_dir() can adjust the library to use the new directory. It is assumed that the directory structure within the data directory has not changed.

Object I/O

The I/O portion of the library is still experimental.

Collections of objects can be stored in a collection class, and these collections can be written to or read from text or binary files. User-defined classes may be added to the collections and may be read and written to files as long as a descendant of io_base is provided.

Every type has an associated I/O type which is a descendant of io_base. In order to perform any sort of input/output on any type, an object of the corresponding I/O type must be instantiated by the user. This is not done automatically by the library. (Since it doesn't know which objects are going to be used ahead of time, the library would have to instantiate all of the I/O objects, which is needlessly slow.) This makes the I/O slightly less user-friendly, but much more efficient. For convenience, each subsection of the library has a class (named with an _ioc suffix) which will automatically allocate all I/O types for that subsection.

Level 1 functions: Functions that input/output data from library-defined objects and internal types from files and combine these objects in collections. These are primarily member functions of the class collection.

Level 2 functions: Functions which are designed to allow the user to input or output data for user-generated objects. These are primarily member functions of classes cinput and coutput.

Level 3 functions: Functions which allow low-level modifications on how input and output is performed. Usage of level 3 functions is not immediately recommended for the casual user.

Level 1 usage:

For adding an object to a collection when you have a pointer to the I/O object for the associated type: 

    int collection::add(std::string name, io_base *tio, void *vec, 
    int sz=0, int sz2=0, bool overwrt=true, bool owner=false);
For adding an object to a collection otherwise: 
    int collection::add(std::string name, std::string stype,
    void *vec, int sz=0, int sz2=0,
    bool overwrt=true, bool owner=false);

To retrieve an object as a 

void *
from a collection use one of: 
    int get(std::string tname, void *&vec);
    int get(std::string tname, void *&vec, int &sz);
    int get(std::string tname, void *&vec, int &sz, int &sz2);
    int get(std::string tname, std::string &stype, void *&vec);
    int get(std::string tname, std::string &stype, void *&vec, int &sz);
    int get(std::string tname, std::string &stype, void *&vec, int &sz,
    int &sz2);
When retrieving a scalar object without error- and type-checking you can use the shorthand version: 
    void *get(std::string name);

To output one object to a file: 

    int collection::out_one(out_file_format *outs, std::string stype, 
    std::string name, void *vp, int sz=0, int sz2=0);

To input one object from a file with a given type and name: 

    int collection::in_one_name(in_file_format *ins, std::string stype, 
    std::string name, void *&vp, int &sz, int &sz2);
To input the first object of a given type from a file: 
    int collection::in_one(in_file_format *ins, std::string stype, 
    std::string &name, void *&vp, int &sz, int &sz2);

Level 2 usage (string-based):

If you don't have a pointer to the io_base child object corresponding to the type of subobject that you are manipulating, then you can use the following functions, which take the type name as a string.

To input a sub-object in an io_base template for which memory has already been allocated use one of: 

int collection::object_in(std::string type, in_file_format *ins, void *vp,
                          std::string &name);
int collection::object_in(std::string type, in_file_format *ins, void *vp,
                          int sz, std::string &name);
int collection::object_in(std::string type, in_file_format *ins, void *vp,
                          int sz, int sz2, std::string &name);

To automatically allocate memory and input a sub-object of a io_base template use one of: 

int collection::object_in_mem(std::string type, in_file_format *ins, 
                              void *vp, std::string &name);
int collection::object_in_mem(std::string type, in_file_format *ins, 
                              void *vp, int sz, std::string &name);
int collection::object_in_mem(std::string type, in_file_format *ins, 
                              void *vp, int sz, int sz2, std::string &name);

To output a subobject in an io_base template use: 

    int collection::object_out(std::string type, out_file_format *outs,
    void *op, int sz=0, int sz2=0, std::string name="");

Level 2 usage (with io_base pointer):

To input a sub-object in an io_base template for which memory has already been allocated use one of: 

virtual int object_in(cinput *cin, in_file_format *ins, object *op,
                      std::string &name);
virtual int object_in(cinput *cin, in_file_format *ins, object *op,
                      int sz, std::string &name);
virtual int object_in(cinput *cin, in_file_format *ins, object **op,
                      int sz, int sz2, std::string &name);
template<size_t N>
int object_in(cinput *co, in_file_format *ins,
              object op[][N], int sz, std::string &name);

To automatically allocate memory and input a sub-object of a io_base template use one of: 

virtual int object_in_mem(cinput *cin, in_file_format *ins,
                          object *&op, std::string &name);
virtual int object_in_mem(cinput *cin, in_file_format *ins, object *&op,
                          int &sz, std::string &name);
virtual int object_in_mem(cinput *cin, in_file_format *ins,
                          object **&op, int &sz, int &sz2,
                          std::string &name);
template<size_t N>
int object_in_mem(cinput *co, in_file_format *ins,
                  object op[][N], int &sz, std::string &name);

To output a subobject in an io_base template use: 

virtual int object_out(coutput *cout, out_file_format *outs, object *op,
                       int sz=0, std::string name="");
virtual int object_out(coutput *cout, out_file_format *outs, object **op,
                       int sz, int sz2, std::string name="");
template<size_t N>
int object_out(coutput *cout, out_file_format *outs,
               object op[][N], int sz, std::string name="");

To automatically allocate/deallocate memory for an object, use: 

    virtual int mem_alloc(object *&op);
    virtual int mem_alloc_arr(object *&op, int sz);
    virtual int mem_alloc_2darr(object **&op, int sz, int sz2);
    virtual int mem_free(object *op);
    virtual int mem_free_arr(object *op);
    virtual int mem_free_2darr(object **op, int sz);

Usage of io_tlate

The functions io_tlate::input() and io_tlate::output() need to be implemented for every class has information for I/O. For subobjects of the class, cinput::object_in() and cinput::object_out() can be called to input or output the information associated with the subobject. For input, cinput::object_in_name(), cinput::object_in_mem(), and cinput::object_in_mem_name() allow the freedom to input an object with a name or with memory allocation. The function coutput::object_out_name() allows one to output an object with a name. If the class contains a pointer to the subobject, then io_base::pointer_in() or io_base::pointer_out() can be used.

Example source code

Example list

Function and solver example

/* Example: ex_fptr.cpp
   -------------------------------------------------------------------
   This gives an example of the how member functions and external
   parameters are supplied to numerical routines. In this case, a
   member function with two parameters is passed to the gsl_root_brent
   class, which solves the equation. One of the parameters is member
   data, and the other is specified using the extra parameter argument
   to the function.
*/

#include <o2scl/funct.h>
#include <o2scl/gsl_root_brent.h>
#include <o2scl/test_mgr.h>

using namespace std;
using namespace o2scl;

class my_class {

private:

  double parameter;

public:
  
  void set_parameter() { parameter=0.01; }
  
  // A function demonstrating the different ways of implementing
  // function parameters
  double function_to_solve(double x, double &p) {
    return atan((x-parameter)*4)*(1.0+sin((x-parameter)*50.0)/p);
  }
  
};

// A simple function to make the plot 
int plot(double sol);

int main(void) {
  
  cout.setf(ios::scientific);
  
  test_mgr t;
  // Only print something out if one of the tests fails
  t.set_output_level(1);

  // The solver, specifying the type of the parameter (double)
  // and the function type (funct<double>)
  gsl_root_brent<double,funct<double> > solver;

  my_class c;
  c.set_parameter();
  
  // This is the "magic" that allows specification of class member
  // functions as functions to solve. This object-oriented approach
  // avoids the use of static variables and functions and multiple
  // inheritance at the expense of a little overhead. We need to
  // provide the address of an instantiated object and the address of
  // the member function.
  funct_mfptr_noerr<my_class,double> function(&c,&my_class::function_to_solve);
  
  double x1=-1;
  double x2=2;
  double p=1.1;
  
  // The value verbose=1 prints out iteration information
  // and verbose=2 requires a keypress between iterations.
  // The parameter p=0.1 is used.
  solver.verbose=1;
  solver.solve_bkt(x1,x2,p,function);

  // This is actually a somewhat difficult function to solve because
  // of the sinusoidal behavior.
  cout << "Solution: " << x1 
       << " Function value: " << c.function_to_solve(x1,p) << endl;

  // Write the function being solved to a file (see source code 
  // in examples directory for details)
  plot(x1);

  t.report();
  return 0;
}
// End of example

The image below shows how the solver progresses to the solution of the example function.

ex_fptr.png

ex_fptr.png

Multidimensional solver example

/* Example: ex_mroot.cpp
   -------------------------------------------------------------------
   Several ways to use an O2scl solver to solve a simple function
*/

#include <cmath>
#include <o2scl/test_mgr.h>
#include <o2scl/mm_funct.h>
#include <o2scl/gsl_mroot_hybrids.h>
#include <o2scl/cern_mroot.h>

using namespace std;
using namespace o2scl;

int gfn(size_t nv, const ovector_view &x, 
        ovector_view &y, void *&pa) {
  y[0]=sin(x[1]-0.2);
  y[1]=sin(x[0]-0.25);
  return 0;
}

class cl {

public:

  // Store the number of function and derivative evaluations
  int nf, nd;

  int mfn(size_t nv, const ovector_view &x, ovector_view &y, 
          void *&pa) {
    y[0]=sin(x[1]-0.2);
    y[1]=sin(x[0]-0.25);
    nf++;
    return 0;
  }

  int operator()(size_t nv, const ovector_view &x, 
                 ovector_view &y, void *&pa) {
    y[0]=sin(x[1]-0.2);
    y[1]=sin(x[0]-0.25);
    nf++;
    return 0;
  }

  int mfnd(size_t nv, ovector_view &x, ovector_view &y,
           omatrix_view &j, void *&pa) {
    j[0]