[Why we need a subroutine library]

Due to the nature of the math in our rheological update, we need to incorporate new libraries into DynEarthSol capable of solving the newly introduced program. Our project uses the effective medium theory, a series of formulae to estimate elastic parameters of some material based on the abundance and orientation of cracks within that material. Our math routine begins with the construction of the stiffness matrix for an intact rock using Lame’s first parameter (\lambda) and the shear modulus of the rock; from there, we manipulate the variables such that we can estimate the new stiffness matrix for a cracked (damaged) rock and eventually the incremental stress acting on the material. However, prior to this implementation, DynEarthSol did not require math routines capable of manipulating rank-4 tensors such as the elastic stiffness tensor. Although we temporarily incorporated the math directly into the program, we have since worked to find a more permanent subroutine library specifically for these high-order structures.

A dedicated subroutine library for mathematical computations is highly important. Years worth of resources have already been invested in developing well-known algorithms (specifically for numerical software such as DynEarthSol) and they have been rigorously tested for accuracy and efficiency. This makes them much less likely to contain bugs than if someone naively re-implements the routines directly. Furthermore, incorporating a subroutine library will generally result in a more readable script.

[Subroutines Considered]

We considered three different subroutine libraries specializing in linear algebra operations: ALGLIB, BLAS/LAPACK, and Eigen. The oldest of these is the Basic Linear Algebra Subprogram (BLAS). Initially released in 1979, this library contains low-level linear algebra computations, subdivided into three levels. Generally, Level 1 contains vector operations; Level 2 matrix-vector operations; and Level 3 matrix-matrix operations. BLAS has since been incorporated into newer algebra packages such as LINPACK, MATLAB, and NumPy. Relevant to our work, BLAS is incorporated into both the Linear Algebra PACKage (LAPACK) and Eigen. LAPACK was initially released in 1992 and written in the Fortran coding language. LAPACK is a successor to both LINPACK and BLAS, building upon both to create efficient routines for updated computer software systems. A few years following LAPACK’s debut, ALGLIB was released in 1999. Although ALGLIB is capable of linear algebra, it also contains functions for solving fast fourier transforms, ordinary differential equations, data analysis, and more. It is available as a free, public library or as an optimized commercial version. Unlike the other libraries, ALGLIB can boast that it is available in multiple programming languages, including C++, C#, and Python. Despite being touted for its portability and ease of use, we were dissatisfied with ALGLIB because the standard routine required importing many header and C preprocessor files into the project folder for compilation. Supposedly, this was to ensure that ALGLIB’s original dependent files were not corrupted, but it created clutter in the project tree which we did not appreciate. The third and final library was Eigen, initially released in 2006. Eigen is an open-source library written in C++ which also provides all of the BLAS library and some routines from LAPACK. Like LAPACK, Eigen is known for its linear algebraic and vector-matrix operative capabilities.

[The Evaluation/ The Chosen One]

We evaluated the libraries by providing the same task to each and timing how long it took each to perform 100,000 iterations. They were to construct the same 6x6 matrix composed of real positive numbers less than 1, overwrite the diagonal elements with random numbers within the range of (0,100), and solve for the inverse of the new matrix. Repeat. Without any optimizations, ALGLIB performed the worst taking 2.836 seconds, followed by LAPACK at 0.174 seconds, leaving Eigen as best-performing with 0.085 seconds. ALGLIB’s disappointing performance combined with the hassle involved with running it led us to stop considering it at this point. Then we tested LAPACK and Eigen’s performances using the -O2 optimizer. We concluded that LAPACK as it is installed on our machine is already optimized as the difference in run-time was minimal. On the other hand, Eigen ran an order faster, clocking in at 0.0089 seconds, hence we decided to use Eigen for DynEarthSol’s rheological update.