Matrices and Linear Algebra

Matrices arise in many, many, many different contexts. Most generally a matrix is simply a rectangular array of entities also called the components of the matrix. Depending on the context in which the matrix comes into existence, the entities themselves may be elements of number field, such as the field of real numbers, complex numbers, or the finite fields or a ring, such as the integers, ring of functions, or polynomials over a field. When specifying a particular component of a matrix A, one usually specifies the row first then the column, so a_ij is the component of row i column j, the matrix A then is sometimes denoted (a_ij).
The C source and NASM source for various matrix functions can be accessed by following the links below:

Various Collections of Matrix Routines

General Purpose Matrix Routines
The generic nature of the routines in this collection of matrix routines consists of routines which are applicable to matrices regardless of whether or not the components of the matrix are members of any algebraic object. For instance, it is possible to set the matrix A equal to the matrix B, A = B, as long as A and B have the same dimensions, without regard for the nature of the components, i.e. this operation depends only on the nature of something being a matrix and the nature of the operation =. Although the generic nature of the routines are applicable to matrices in general, all the implementations available at this website are for real matrices whose components are of type double.
Arithmetic Routines
The routines in this collection of matrix routines consists of those routines for which the components are themselves elements of a ring or in some cases a ring with identity. The algebraic operations of the ring define algebraic operations on matrices, e.g. matrix addition may be defined component-wise using the addition operation of the ring or matrix multiplication may be defined using the addition operation and multiplication operation defined on the ring. In particular, matrix inversion is not included in this collection, but rather in the link below. Again, although the generic nature of the routines are applicable to matrices over a ring (with identity), all the implementations available at this website are for real matrices whose components are of type double.
Systems of Linear Equations
For matrices defined over a field it is possible to find a solution x to the matrix equation Ax = B for a fixed n×n matrix A and a fixed n×1 matrix B providing that the matrix A is non-singular, x = A^-1B. In general, it is better to solve the equation Ax = B directly for x than to find the inverse A^-1. Again, although the generic nature of the routines are applicable to matrices over a field, all the implementations available at this website are for real matrices whose components are of type double.
Normed and Inner Product Spaces
The routines in this collection supplement those of the above collections and are applicable to matrices which arise from linear transformations defined on an inner product space or a normed linear space. Again, although the generic nature of the routines are applicable to matrices over a field, all the implementations available at this website are for real matrices whose components are of type double.
Eigenvalues and Eigenvectors
The routines in this collection of matrix routines consists of routines which calculate the eigenvalues and eigenvectors of a square matrix. Generally in order to calculate the eigenvalues for square matrix the ground field needs to be closed i.e. the characteristic polynomial splits into linear factors with constants belonging to the field. This means that arbitrary matrices should be defined over the field of complex numbers. Note that real symmetric matrices, however, have real eigenvalues. Again, although the generic nature of the routines are applicable to matrices over a closed field, all the implementations available at this website are for real matrices whose components are of type double. For this reason some of the routines may fail if it is not known apriori that the eigenvalues are real.

Various Examples of the Birth of Matrices

Linear Transformations:
Given an n-dimensional real vector space V an m-dimensional real vector space W and a linear transformation L:V→W with domain V and range W, then given basis vectors {v₀, . . ., v_n-1} of V and basis vectors {w₀, . . ., w_m-1} of W, the linear transformation L can be represented as a matrix of real numbers where l_ij where L(v_j) = Σ_il_ijw_i. Moveover if W is an inner product space and {w₀, . . ., w_m-1} is an orthonormal basis then l_ij = <w_i, L(v_j)>, where < > denotes the inner product. Further, the addition of linear transformations corresponds to addition of the corresponding matrices and composition of linear transformations corresponds to multiplication of the corresponding matrices. If V and W are complex vector spaces, then L can be represented as matrix of complex numbers. In general if V and W are vector spaces over a field F, then the linear transformation L can be represented as a matrix of elements of the field F.

Bilinear Functions:
Given an n-dimensional real vector space V, an m-dimensional real vector space W and a bilinear function B:V×W→R, then given basis vectors {v₀, . . , v_n-1} of V and basis vectors {w₀, . . , w_m-1} of W, the bilinear function B can be represented as a matrix of real numbers where b_ij = B(v_i,w_j). If v =Σ_ia_iv_i and w =Σ_jc_jw_j, then B(v,w) = Σ_iΣ_ja_ib_ijc_j.

Probability Transition Matrices:
Given a stationary first-order discrete Markov process on states X_i, the probability transition matrix P is a matrix whose components are non-negative real numbers such that the sum of the components along a row equals 1. The component p_ij is the probability of transition from state X_i to X_j during one unit of time. For non-stationary first-order discrete Markov processes, the probability transition matrices depend on the time parameter, i.e. at time n, P⁽ⁿ⁾ is a probability transition matrix whose i,j-th component p⁽ⁿ⁾_ij is the probability of transition from state X_i to X_j during the time interval n to n+1.

Systems of Linear Equations:
Given a system of linear equations
a_i,0 x₀ + . . . + a_i,n-1 x_n-1 = b_i, for i = 0, . . ., n-1. The coefficients can be represented as the n×n matrix A = (a_ij), the right-hand side can be written as the n × 1 matrix B = (b_i) and the unknown quantities x_i as the n × 1 matrix x = (x_i). The system of linear equations in matrix form is then Ax = B.
Graphs and Networks:
In graph theory and network theory many different matrices are defined the most common being the incidence matrix, the circuit matrix, the path matrix, the adjacency matrix, and in switching theory the switching matrix, connection matrix, the transmission matrix etc. Some of these matrices are boolean matrices and some are defined over the ring of boolean polynomials. Currently, at this website, the matrix routines for which the C source code is provided deals only with real-valued matrices in which the components are of type double.

Machine Storage of Matrices

In C there are two ways to declare a matrix. The first is double A[M][N] which declares A to be a matrix with M rows and N columns and the second is double **A which declares A to be a pointer to an array of pointers. For the second method the user must call malloc() or an equivalent function to dynamically allocate the memory for the arrays. Given a statement involving A[i][j], the way the compiler accesses a memory location depends on the way A was declared. If A was declared as double A[M][N]; then the memory location at &A[0][0] + i*N + j is used, (C stores a matrix in consecutive memory locations sequentially by the row of the matrix.) And if A was declared as double **A; then the memory location at *(A+i)+j is used.

All matrices for which the C source code and ASM source code is provided at this website assume that the matrices are stored sequentially as if having been declared in the calling routine as double A[M][N].

For statically defined matrices declared as double A[M][N] the C source code at this website would be called using either the cast (double*) A or as &A[0][0]. Dynamic memory allocation of matrix then can be achieved as follows:

#include <stdlib> // required for malloc()
int i,m,n;
double *A;
double **pA;
   .
   .
   .
   A = (double*) malloc(sizeof(double) * m * n);
   pA = (double**) malloc(sizeof(double) * m);
   for(i = 1, *pA = A; i < m; i++) *(pA+i) = *(pA + i - 1) + n;

Note that if *(pA+i) is changed for some i, then the matrix A and the matrix pA will no longer agree. The C source code at this website then would be called using A as the matrix argument.