Matrix routines

A common matrix task is the solution of a set of simultaneous linear equations, since this can be conveniently expressed as

$$\m Ax=z$$

where $A$ is an $n$ by $n$ matrix and $x$ and $z$ are vectors of size $n$. CCATSL provides several routines for this: GaussElimCL solves the system above using Gaussian elimination with pivoting, while DecomposeCL and SolveCL respectively decompose a matrix into a more useful form and use the resulting decomposition to solve the equation. The DecomposeCL/SolveCL approach is preferable when you want to solve the equation above for several different $z$ since DecomposeCL does not involve $z$ and SolveCL is very quick.

Two special cases are where the matrix $A$ is tridiagonal

$$\m A= \left( \begin{array}{cccccc} b_1 & c_1 & 0 & 0 & \ldot... ... c_{n-1} \\ 0 & 0 & \ldots & 0 & a_n & b_n \end{array}\right)$$

for which the routine TridiagCL is better suited, or banded such as

$$\m A= \left( \begin{array}{cccccc} b_1 & c_1 & d_1 & 0 & \ld... ... c_{n-1} \\ 0 & 0 & \ldots & 0 & a_n & b_n \end{array}\right)$$

which are best handled with BandCL. The routine InvertCL inverts a matrix.

Another common problem is the determination of eigenvalues and eigenvectors. CCATSL has three routines for this: EigenMaxCL finds an eigenvalue of greatest modulus and an eigenvector with this eigenvalue, EigenValCL takes a parameter $\lambda$ and finds an eigenvalue $\theta$ minimising $\vert\theta-\lambda\vert,$ together with a corresponding eigenvector. For symmetric $A$, JacobiCL returns all the eigenvalues and eigenvectors, using the method of Jacobi.

If $A$ is an $m$ by $n$ matrix, a singular value decomposition (i.e., a decomposition of the form $A=USV^T$ where $U$ and $V$ are orthogonal and $S$ is diagonal) can be found using the routine SvdCL.

Gaussian Elimination with Pivoting, GaussElimCL

This routines use Gaussian elimination with pivoting to solve the system of simultaneous linear equations

$$\m Ax=z$$

where $A$ is an $n$ by $n$ matrix and $x$ and $z$ are vectors of size $n$.

	double GaussElimCL (	int n,
		int ncols,
		double *a,
		double *z,
		double *x);

	n		The order of the matrix $A.$
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		Array holding the matrix $A.$
	z		Array holding the vector $z.$
	x		Array to receive the solution vector $x.$

The return value is the determinant of the matrix $A$. A typical use of GaussElimCL might look like

/* /examples/chapter2/elim.c */ #include <catam.h> int MainCL(void) { double a[2][2]; double z[2]; double x[2]; double det; a[0][0]=1; a[0][1]=1; a[1][0]=0; a[1][1]=1; z[0]=2; z[1]=3; det=GaussElimCL(2,2,a,z,x); printf("determinant is %lf",det); return 0; }

A complete example of how to use GaussElimCL can be found in b09matrx.c.

DecomposeCL

DecomposeCL is designed to be used in combination with SolveCL to solve a set of simultaneous linear equations expressed in the form $Ax=z.$ Once the matrix $A$ has been decomposed with DecomposeCL, SolveCL can be called to solve for a given $z.$ Since SolveCL is much faster than GaussElimCL, if the same set of equations is to be solved for different vectors $z$, the DecomposeCL/SolveCL approach is to be preferred to GaussElimCL.

	double DecomposeCL (	int n,
		int ncols,
		double *a);

	n		The order of the matrix $A.$
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		Array holding the matrix $A.$ When DecomposeCL returns, $A$ will have been replaced by its decomposed form. If you want to retain the original matrix $A$, you must make a copy of it before calling DecomposeCL.

The return value is the determinant of the matrix $A$. See b09matrx.c for an example of how to use DecomposeCL.

ConditionCL

ConditionCL returns the last condition number of the last matrix decomposition (see DecomposeCL) or inversion (see InvertCL).

SolveCL

This routine complements DecomposeCL by taking a matrix in decomposed form and a vector $z$ and solving the system of linear equations $Ax=z.$ It is faster than Gaussian elimination (GaussElimCL) and is thus more efficient for solving for several different vectors $z.$

	void SolveCL (	int n,
		int ncols,
		double *a,
		double *z);

	n		The order of the matrix $A.$
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		The matrix $A.$
	z		Array to hold the right-hand side $z.$ When SolveCL returns, it will have been overwritten and will hold the solution vector $x.$

The example program b09matrx.c demonstrates the use of SolveCL.

Tridiagonal Linear System Solver, TridiagCL

This routine solves the linear system $Ax=z$ for the special case where $A$ is tridiagonal

$$\m A=\left( \begin{array}{cccccc} b_1 & c_1 & 0 & 0 & \ldots... ...c_{n-1} \\ 0 & 0 & \ldots & 0 & a_n & b_n \end{array}\right).$$

	void TridiagCL (	int n,
		double *a,
		double *b,
		double *c,
		double *z,
		double *x);

	n		The order of the matrix $A.$
	a		Array holding the values of  $a_1,\ldots{,}\,a_n.$ (The value of $a_1$ is irrelevant but it must be present.)
	b		Array holding  $b_1,\ldots{,}\, b_n.$
	c		Array holding  $c_1,\ldots{,}\, c_n.$ (The value of $c_n$ is irrelevant.)
	z		Array holding the vector $z.$
	x		Array to hold the solution vector $x.$

A typical use of TridiagCL is the following:

/* /example/chapter2/tridiag.c */ #include <catam.h> int MainCL(void) { double a[5],b[5],c[5],z[5],x[5]; c[0]= 1; c[1]= 1; c[2]= 1; c[3]=1; b[0]=-2; b[1]=-2; b[2]=-2; b[3]=-2; b[4]=-2; a[1]= 1; a[2]= 1; a[3]= 1; a[4]= 1; z[0]=0.5; z[1]=0.5; z[2]=0.5; z[3]=0.5; z[4]=0.5; TridiagCL(5,a,b,c,z,x); }

Banded Linear System Solver, BandCL

This routine solves the linear system $Ax=z$ when $A$ is banded: $A_{ij}=0$ if $i-j>l$ or $i-j<-u$ where $0\le l, u\le n.$ The numbers $l$ and $u$ are called the lower bandwidth and the upper bandwidth respectively. For example, the matrix

$$\m A= \left( \begin{array}{ccccc} b_1 & c_1 & d_1 & 0 & 0 \\... ...& a_4 & b_4 & c_4 \\ 0 & 0 & 0 & a_5 & b_5 \end{array}\right)$$

has a lower bandwidth of $1$ and an upper bandwidth of 2. The advantage of BandCL over a standard Gaussian elimination routine such as GaussElimCL is that it exploits the sparse structure to speed up the computation, and that it can handle much larger matrices, since the full matrix $A$ is never actually stored. Instead a compact form of $A$ must be setup and passed to BandCL. The compact form of the matrix above is

$$\m \left( \begin{array}{ccccc} \cdot & b_1 & c_1 & d_1 \\ a... ...& c_4 & \cdot \\ a_5 & b_5 & \cdot & \cdot \end{array}\right)$$

where dotted entries are irrelevant and ignored by BandCL. Note that if $A$ is $n$ by $n$ then the compact form of $A$ is $n$ by $(1+l+u)$ where $l$ and $u$ are the lower and upper bandwidths.

	double BandCL (	int n,
		int ncols,
		int lbw,
		int ubw,
		double *cpact,
		double *z,
		double *x);

	n		The order of the matrix $A.$
	ncols		Normally ncols=1+lbw+ubw. Specifically, ncols is the size of the second dimension in the declaration of the array cpact.
	lbw		The lower bandwidth of the matrix.
	ubw		The upper bandwidth.
	cpact		Array holding the matrix $A$ stored in compact form. (See above for details of the compact form.)
	z		Array holding the vector $z.$
	x		Array to hold the solution vector $x.$

The return value is the determinant of the matrix. b09matrx.c contains an example of how to use BandCL.

Matrix Inversion, InvertCL

InvertCL inverts a square matrix and returns its determinant. The inverse overwrites the original matrix; thus if you want to retain the original matrix you must make a copy of it prior to calling InvertCL.

	double InvertCL (	int n,
		int ncols,
		void *a);

	n		The order of the matrix.
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		Array holding the matrix $A.$ When InvertCL returns, it will have been overwritten with the inverse matrix.

The return value is the determinant of the matrix $A$. b09matrx.c has an example of how to use InvertCL.

Singular Value Decomposition, SvdCL

The CCATSL routine SvdCL computes a decomposition of the form $A=USV^T$ where $U$ and $V$ are orthogonal and $S$ is diagonal. If $A$ is $m$ by $n$, then $U$ is $m$ by $m$, $S$ is $m$ by $n$ and $V$ is $n$ by $n$, thus if $m>n$, the last $m-n$ rows of $S$ will be zero and the last $m-n$ columns of $U$ will not have any special significance. Similarly if $n>m$, the last $n-m$ columns of $S$ will be zero and the last $n-m$ columns of $V$ are not significant. Since $S$ is diagonal, SvdCL just returns the diagonal elements (the singular values) of $S$ as a one-dimensional vector.

	int SvdCL (	int m,
		int n,
		int ncols,
		double *a,
		double *u,
		double *v,
		double *w);

	m		The number of rows in the matrix $A$.
	n		The number of columns in the matrix $A$.
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		Array holding the matrix $A$.
	u		Array to hold the matrix $U$.
	v		Array to hold the matrix $V$.
	w		Array to hold the singular values (a one-dimensional vector).

The return value is zero if the decomposition was successful. A non-zero value, $k$ say, indicates that the iteration for the $k$th singular value failed. A typical use of SvdCL might look like

/* /examples/chapter2/svd.c */ #include <catam.h> int MainCL(void) { double a[2][3]; double u[2][2]; double v[3][3]; double w[3]; int ok; a[0][0]=1; a[0][1]=1; a[0][2]=2; a[1][0]=0; a[1][1]=2; a[1][2]=2; ok=SvdCL(2,3,3,a,u,v,w); return 0; }

Eigenvalues and Eigenvectors: the Power method, EigenMaxCL

EigenMaxCL finds an eigenvalue of largest modulus of a matrix $A$, together with a corresponding eigenvector.

	double EigenMaxCL (	int n,
		int ncols,
		double acc,
		double *a,
		double *evec);

	n		The order of the matrix.
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	acc		The absolute accuracy required.
	Ap		Array holding the matrix $A$.
	evec		Array holding an initial `guess' $v$ at the desired eigenvector. This guess does not have to be good, but if it is way out (and you are also very unlucky) you might have $A^n v=0$ for some $n$, at which point EigenMaxCL will fall over (by setting ErrorFlagCD=true and terminating). Assuming this does not happen, when EigenMaxCL returns, it will hold a genuine eigenvector.

The return value is an eigenvalue of largest modulus A typical use of EigenMaxCL might look like

/* /examples/chapter2/power.c */ #include <catam.h> int MainCL(void) { double a[2][2]; double evect[2]; double eval; a[0][0]=1; a[0][1]=1; a[1][0]=1; a[1][1]=0; evect[0]=1; evect[1]=1; eval=EigenMaxCL(2,2,0.00001,a,evect); return 0; }

Eigenvalues and Eigenvectors: the Modified Power method, EigenValCL

The drawback of EigenMaxCL is that it can only give an eigenvector corresponding to an eigenvalue of largest modulus. EigenValCL takes a parameter $\theta$ and returns an eigenvalue $\lambda$ and a corresponding eigenvector minimising $\vert\lambda-\theta\vert$ for a given matrix $A$.

	double EigenValCL (	int n,
		int ncols,
		double acc,
		double *a,
		double *evec,
		double theta);

	n		The order of the matrix $A$.
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	acc		The absolute accuracy required.
	a		Array holding the matrix $A$.
	evec		Array holding an initial `guess' $v$ at the desired eigenvector. (This guess does not have to be very good.)
	theta		The parameter $\theta$. If EigenValCL encounters a problem (as it will, for example, if $\theta$ actually is an eigenvalue) it will set ErrorFlagCD=true and terminate.

See EigenMaxCL for an illustration of how to use a very similar function. The return value is the desired eigenvalue, $\lambda$.

Eigenvalues and Eigenvectors: the Jacobi method for real symmetric matrices, JacobiCL

The Jacobi method is a powerful technique for finding the eigenvalues and eigenvectors of a real symmetric matrix $A$.

	void JacobiCL (	int n,
		int ncols,
		double *a,
		double *evals);

	n		The order of the matrix $A$.
	ncols		Normally ncols=n. Specifically, ncols is the size of the second dimension in the declaration of the array a.
	a		Array holding the matrix $A$. This array will be corrupted by the call to JacobiCL. The eigenvectors will be returned as the rows of the array a.
	evals		Array to hold the eigenvalues.

An example of how to use JacobiCL is shown below:

/* /examples/chapter2/jacobi.c */ #include <catam.h> int MainCL(void) { double a[2][2]; double evals[2]; a[0][0]=1; a[0][1]=-1; a[1][0]=1; a[1][1]= 1; JacobiCL(2,2,false,a,evals); return 0; }