Computational Physics Library (cpl)

A numerical analysis library written in C. Interface inspired by the GNU Scientific Library.

Monte Carlo simulation of a heavy ion collision in the Glauber model.

Include the header file cpl.h to get started.

#include "path/to/code/cpl.h"

Vector and matrices

This library introduces new data types cpl_vector and cpl_matrix which are heap allocated (and therefore must be memory managed manually). A new 3-dimensional vector can be initialized by calling cpl_vector_alloc, which returns a pointer to the allocated memory block.

cpl_vector *v = cpl_vector_alloc(3);

However, the vector has not been initialized and accessing its components would give junk values. Components can be assigned values all together

cpl_vector_build(v, 5, 6, 7);

or one by one with

cpl_set(v, 1, 5); // sets the first component to 5
cpl_set(v, 2, 6); // sets the second component to 6
cpl_set(v, 3, 7); // sets the third component to 7

which both change v in-place. Indexing starts at 1.

There are analogous functions for defining matrices.

cpl_matrix *A = cpl_matrix_alloc(3, 3);
cpl_matrix_build(A,  5,  6,  7,
                     8,  9, 10,
                    11, 12, 13);

Components of matrices can also be set one by one.

cpl_set(A, 1, 1, 5); // sets the (1, 1) component to 5
cpl_set(A, 1, 2, 6); // sets the (1, 2) component to 6
cpl_set(A, 1, 3, 7); // sets the (1, 3) component to 7

Components of vectors and matrices can be accessed by cpl_get.

scalar v1 = cpl_get(v, 1);

scalar A11 = cpl_get(A, 1, 1);

Scalars. The keyword scalar is set in cpl_defines.h to double but it can be changed to float or complex if needed.

The library also includes the cpl_print function which prints well-formatted vectors and matrices.

cpl_print(v);

outputs

3-dimensional vector (double)
 5.000
 6.000
 7.000

and similarly

cpl_print(A);

outputs

3×3 matrix (double)
 5.000     6.000    7.000
 8.000     9.000    10.000
 11.000    12.000   13.000

Note. cpl_set, cpl_get and cpl_print are generics which get dispatched according to the type of the first argument.

Utility functions

In the following code v and w are (pointers to) cpl_vectors and A and B are (pointers to) cpl_matrixs.

int d = cpl_vector_dim(v);   // returns the dimension on v
int m = cpl_matrix_rows(A);  // returns the number of rows of A
int n = cpl_matrix_cols(A);  // returns the number of columns of A

cpl_vector_isequal(v, w);    // returns 1 if v equals w component-wise, 0 otherwise
cpl_matrix_isequal(A ,B);    // returns 1 if A equals B component-wise, 0 otherwise  

cpl_vector_overwrite(v, w);  // overwrites v with w
cpl_matrix_overwrite(A, B);  // overwrites A with B

cpl_matrix *id = cpl_matrix_id(d);  // returns a d×d identity matrix

Inner products and norms.

cpl_vector_l2norm(v);     // returns the Euclidean norm of v
cpl_vector_inner(v, w);   // returns the Euclidean inner product of v and w
cpl_vector_l2dist(v, w);  // returns the Euclidean norm of v - w

cpl_vector_normalize(v);  // normalizes v so that its norm is 1

Vector and matrix algebra

cpl_scale(v, 5.0);  // multiplies each entry of v by 5.0
cpl_scale(A, 0.2);  // multiplies each entry of v by 0.2

cpl_add(v, w);      // overwrites w with the sum v + w
cpl_add(A, B);      // overwrites B with the sum A + B

cpl_mult(A, v);     // overwrites v with the product A*v
cpl_mult(A, B);     // overwrites B with the product A*B

cpl_vector *Av = cpl_mult_alloc(A, v); // returns a newly allocated vector for the product A*v
cpl_matrix *AB = cpl_mult_alloc(A, B); // returns a newly allocated matrix for the product A*B

Freeing allocated memory

Heap allocations should be freed when no longer in use to avoid memory leakage.

cpl_free(v);
cpl_free(w);
cpl_free(A);
cpl_free(B);

Solving linear systems

Suppose that A is an n × n square matrix and B is a n × 1 matrix, then the solution to the linear system A*X = B---where X is an unknown n × 1 matrix---can be computed using Gauss--Jordan elimination as follows.

cpl_linalg_gaussjordan(A, B);  // overwrites B with the solution

Decompositions

For a symmetric square matrix A, its Cholesky decomposition A = L*Lᵀ can be computed as follows.

cpl_matrix *L = cpl_linalg_cholesky(A);   // allocates new memory

For a square matrix A, its LU decomposition A = L*U can be computed as follows.

cpl_matrix *LU = cpl_linalg_ludecomp(A);  // allocates new memory

LU is a square matrix such that the upper triangular matrix U[i, j] = LU[i, j] for i <= j and the lower triangular matrix L[i, j] = LU[i, j] for i > j and L[i, i] = 1. There is a utility function to get L and U from LU separately.

cpl_linalg_luseparate(LU, L, U);          // overwrites L and U

There is a solver that uses LU decomposition and forward/back substitutions for linear systems. Let A be an n × n matrix and b be an n-dimensional vector, then the solution to the linear system A*x = b can be computed as follows.

cpl_linalg_lusolve(A, b);                 // overwrites b with the solution

The LU decomposition can also be used to compute the determinant of a matrix.

scalar det = cpl_linalg_ludet(A);

Iterative solvers

The linear system A*X = B---where A is an n × n matrix, B is an n × m matrix, and X is the unknown n × m matrix---can be solved by Jacobi iteration as follows. Iteration starts with the passed matrix X.

cpl_linalg_jacobi(A, B, X, res);  // overwrites X with the solution

The final argument res is a 0-dimensional cpl_vector that keeps track of residues as iteration progresses. NULL can be passed instead of a cpl_vector if this information is not relevant.

The inverse of A can be computed using the above method if B is set to the identity matrix.

int dim = cpl_matrix_rows(A);

cpl_matrix *id = cpl_matrix_id(dim);
cpl_matrix *Ainv = cpl_matrix_calloc(dim, dim);

cpl_linalg_jacobi(A, id, Ainv, NULL);

The system can also be solved by the Gauss--Seidel method.

cpl_linalg_seidel(A, id, Ainv, NULL);

And the conjugate gradient method.

cpl_linalg_conjgrad(A, id, Ainv);

Tolerance and maximum iterations. The tolerance TOL and the maximum number of iterations MAX_ITERS for iterative methods can be set in cpl_defines.h. By default they are set to 1e-6 and 20000 respectively.