"mini linear algebra"
A standalone Fortran implementation of common linear algebraic operations for dense matrices and vectors.
implemented (otherwise it has not yet been implemented)
(f) - function
(s) - subroutine
Simple linear algebra operations not included as Fortran intrinsics. Building blocks for other modules.
disp_mat(A)(s): displays matrix A
get_diag(A)(f): returns diagonal elements of matrix A as vector
set_diag(diag_vec)(f): returns square matrix of zeros with diagonal elements specified by diag_vec
identity(n)(f): returns the identity matrix of dimension n
rop(u1,u2)(f): rank one product of two vectors, produces matrix equal to u1 u2^T
diag_sub(D,b)(f): solve diagonal system Dx=b for x=b/get_diag(D)
bck_sub(T,b)(f): solve lower triangular system Tx=b for x using back substitution
fwd_sub(T,b)(f): solve upper triangular system Tx=b for x using forward substitution
strang(T,b)(f): solve triangular system Tx=b after detecting the triangularity type of T. Essentially wrapper around the *sub functions with O(n^2) checks and extra safeguards; if performance is critical use the *_sub functions directly. What a serindiptously great name for the function.
lupp(A,L,U,p)(s): LU decomposition of square matrix with partial pivoting, results satisfy PA=A(p,:)=LU where P is the matrix corresponding to pivot vector p
ldl(A,L,D)(s): LDL decomposition for symmetric matrix as A=L D L^T
qr(A,Q,R)(s): QR decomposition of matrix as A=QR
eig(A,V,d)(s): eigenvalue/spectral decomposition of matrix as A = V diag(d) V^-1
svd(A,U,d,V)(s): singular value decomposition of a matrix as A = U diag(d) V^T
linsol(A,B)(f): solves linear system AX=B via LU decomposition, the solution X will be of the same shape as the right hand side B, which can be a vector (single linear system) or a matrix (set of linear systems)
inverse(A)(f): computes the matrix inverse of A
sor(A,b,rel_tol,max_iter)(f): solves Ax=b via successive over-relaxation (an iterative method); terminates when relative error is less than rel_tol or when max_iter iterations have been performed
lssol(A,r)(f): solve least squares problem Ax \approx b
Used internally by other functions, but not the meat and potatoes of a linear algebra library.
All included in prim.f90 module currently.
is_square(A,n)(f): returns logical value on whether matrix A is square, if so sets n equal to dim(A)
is_triangular(A,type)(f): outputs a logical value, and assigns type to character corresponding to type of triangularity ("q"-not square, "n"-not triangular, "u"-upper triangular, "l"-lower triangular, "d"-diagonal)
get_triang(A,type)(f): returns lower (type="l") or upper (type="u") triangle (NOT including the diagonal) of matrix A as a matrix
Notes
I'm mainly using this library as a way to strengthen my programming in a standard scientific computing language and to better understand implementations of popular linear algebra methods.
This was originally supposed to be a C library, but implementing even vector vector addition of arbitrarily sized arrays was needlessly complicated (almost wholly because C does not allow the determination of array size from pointers to arrays). I think I have figured our the proper way to get around this now, so maybe I'll write a C version in the future.