Changing the Linear Algebra Library

Kriging, which is at the root of many geostatistical algorithms requires basic linear algebra facilities, essentially matrix inversion. Most of the computing time in kriging is actually spent building and solving the kriging system, hence the importance of using an efficient linear algebra library.

The default library used by G_sTL is a slightly modified version of TNT, the Template Numerical Toolkit. It is a public domain library written by Roldan Pozo, Mathematical and Computational Sciences Division, National Institute of Standards and Technology. This library was chosen because it is relatively efficient, it is easy to use and modify, and it is public domain.

However, it is easy to change the linear algebra library used by the G_sTL algorithms without having to modify existing code. The procedure is twofold. The first step is to wrap all the needed functionalities of the library into a single structure, call it new_matrix_lib. This structure will contain nested types (a matrix type, a vector type, ...) and static functions (for example static void inverse(matrix& A)), which must honor the requirements detailed in section 2.4.1.

The following is an extract from the TNT wrapper:

template<class T>
struct TNT_lib{

  typedef TNT::Matrix<T> tnt_Matrix;
  typedef TNT::Vector<T> tnt_Vector;

  // Cholesky factorization.
  static inline int cholesky(TNT::Matrix<T>& A, TNT::Matrix<T>& B){
    return Cholesky_upper_factorization(A,B);
  }

  // LU factorization.
  static inline int LU_factor(TNT::Matrix<T>& A, TNT::Vector<int>& index){
    return TNT::LU_factor(A,index);
  }


It has two nested types: Matrix, and Vector which are defined by TNT, and two functions: a Cholesky factorization, a LU factorization, which simply call existing TNT functions.

The second step is to specialize the matrix library trait class for the new library wrapper new_matrix_lib. The trait class is defined in <matrix_lib_traits.h>. What has been done in the case of TNT can serve as a model.