ccmath_grass_wrapper.c File Reference

#include <grass/ccmath_grass.h>

Include dependency graph for ccmath_grass_wrapper.c:

Go to the source code of this file.

Functions
int	G_math_solv (double **a, double *b, int n)
	Solve a general linear system A*x = b. More...
int	G_math_solvps (double **a, double *b, int n)
	Solve a symmetric positive definite linear system S*x = b. More...
void	G_math_solvtd (double *a, double *b, double *c, double *x, int m)
	Solve a tridiagonal linear system M*x = y. More...
int	G_math_solvru (double **a, double *b, int n)
int	G_math_minv (double **a, int n)
	Invert (in place) a general real matrix A -> Inv(A). More...
int	G_math_psinv (double **a, int n)
	Invert (in place) a symmetric real matrix, V -> Inv(V). More...
int	G_math_ruinv (double **a, int n)
	Invert an upper right triangular matrix T -> Inv(T). More...
void	G_math_eigval (double **a, double *ev, int n)
	Compute the eigenvalues of a real symmetric matrix A. More...
void	G_math_eigen (double **a, double *ev, int n)
	Compute the eigenvalues and eigenvectors of a real symmetric matrix A. More...
double	G_math_evmax (double **a, double *u, int n)
int	G_math_svdval (double *d, double **a, int m, int n)
	Compute the singular values of a real m by n matrix A. More...
int	G_math_sv2val (double *d, double **a, int m, int n)
	Compute singular values when m >> n. More...
int	G_math_svduv (double *d, double **a, double **u, int m, double **v, int n)
int	G_math_sv2uv (double *d, double **a, double **u, int m, double **v, int n)
	Compute the singular value transformation when m >> n. More...
int	G_math_svdu1v (double *d, double **a, int m, double **v, int n)
	Compute the singular value transformation with A overloaded by the partial U-matrix. More...

Function Documentation

G_math_eigen()

void G_math_eigen	(	double **	a,
		double *	ev,
		int	n
	)

Compute the eigenvalues and eigenvectors of a real symmetric matrix A.

The input and output matrices are related by

A = E*D*E~ where D is the diagonal matrix of eigenvalues
D[i,j] = ev[i] if i=j and 0 otherwise.

The columns of E are the eigenvectors.

Parameters

a	= pointer to store for symmetric n by n input matrix A. The computation overloads this with an orthogonal matrix of eigenvectors E.
ev	= pointer to the array of the output eigenvalues
n	= dimension parameter (dim(a)= n*n, dim(ev)= n)

Definition at line 310 of file ccmath_grass_wrapper.c.

G_math_eigval()

void G_math_eigval	(	double **	a,
		double *	ev,
		int	n
	)

Compute the eigenvalues of a real symmetric matrix A.

Parameters

a	= pointer to array of symmetric n by n input matrix A. The computation alters these values.
ev	= pointer to array of the output eigenvalues
n	= dimension parameter (dim(a)= n*n, dim(ev)= n)

Definition at line 289 of file ccmath_grass_wrapper.c.

G_math_evmax()

double G_math_evmax	(	double **	a,
		double *	u,
		int	n
	)

Definition at line 326 of file ccmath_grass_wrapper.c.

G_math_minv()

int G_math_minv	(	double **	a,
		int	n
	)

Invert (in place) a general real matrix A -> Inv(A).

Parameters

a	= array containing the input matrix A. This is converted to the inverse matrix.
n	= dimension of the system (i.e. A is n x n )

Returns

: 0 -> normal exit, 1 -> singular input matrix

Definition at line 242 of file ccmath_grass_wrapper.c.

G_math_psinv()

int G_math_psinv	(	double **	a,
		int	n
	)

Invert (in place) a symmetric real matrix, V -> Inv(V).

The input matrix V is symmetric (V[i,j] = V[j,i]).

Parameters

a	= array containing a symmetric input matrix. This is converted to the inverse matrix.
n	= dimension of the system (dim(v)=n*n)

Returns

: 0 -> normal exit 1 -> input matrix not positive definite

Definition at line 256 of file ccmath_grass_wrapper.c.

G_math_ruinv()

int G_math_ruinv	(	double **	a,
		int	n
	)

Invert an upper right triangular matrix T -> Inv(T).

Parameters

a	= pointer to array of upper right triangular matrix, This is replaced by the inverse matrix.
n	= dimension (dim(a)=n*n)

Returns

value: status flag, with 0 -> matrix inverted -1 -> matrix singular

Definition at line 269 of file ccmath_grass_wrapper.c.

G_math_solv()

int G_math_solv	(	double **	a,
		double *	b,
		int	n
	)

Solve a general linear system A*x = b.

Documentation and ccmath library version 2.2.1 by Daniel A. Atkinson

                         Chapter 1

                       LINEAR ALGEBRA

                          Summary

        The matrix algebra library contains functions that
        perform the standard computations of linear algebra.
        General areas covered are:

                  o Solution of Linear Systems
                  o Matrix Inversion
                  o Eigensystem Analysis
                  o Matrix Utility Operations
                  o Singular Value Decomposition

        The operations covered here are fundamental to many
        areas of mathematics and statistics. Thus, functions
        in this library segment are called by other library
        functions. Both real and complex valued matrices
        are covered by functions in the first four of these
        categories.

Notes on Contents

Functions in this library segment provide the basic operations of

numerical linear algebra and some useful utility functions for operations on vectors and matrices. The following list describes the functions available for operations with real-valued matrices.

o Solving and Inverting Linear Systems:

solv ---—— solve a general system of real linear equations. solvps -—— solve a real symmetric linear system. solvru -—— solve a real right upper triangular linear system. solvtd -—— solve a tridiagonal real linear system.

minv ---—— invert a general real square matrix. psinv --—— invert a real symmetric matrix. ruinv --—— invert a right upper triangular matrix.

 The solution of a general linear system and efficient algorithms for

solving special systems with symmetric and tridiagonal matrices are provided by these functions. The general solution function employs a LU factorization with partial pivoting and it is very robust. It will work efficiently on any problem that is not ill-conditioned. The symmetric matrix solution is based on a modified Cholesky factorization. It is best used on positive definite matrices that do not require pivoting for numeric stability. Tridiagonal solvers require order-N operations (N = dimension). Thus, they are highly recommended for this important class of sparse systems. Two matrix inversion routines are provided. The general inversion function is again LU based. It is suitable for use on any stable (ie. well-conditioned) problem. The Cholesky based symmetric matrix inversion is efficient and safe for use on matrices known to be positive definite, such as the variance matrices encountered in statistical computations. Both the solver and the inverse functions are designed to enhance data locality. They are very effective on modern microprocessors.

o Eigensystem Analysis:

eigen —— extract all eigen values and vectors of a real symmetric matrix. eigval –— extract the eigen values of a real symmetric matrix. evmax —— compute the eigen value of maximum absolute magnitude and its corresponding vector for a symmetric matrix.

 Eigensystem functions operate on real symmetric matrices. Two forms of

the general eigen routine are provided because the computation of eigen values only is much faster when vectors are not required. The basic algorithms use a Householder reduction to tridiagonal form followed by QR iterations with shifts to enhance convergence. This has become the accepted standard for symmetric eigensystem computation. The evmax function uses an efficient iterative power method algorithm to extract the eigen value of maximum absolute size and the corresponding eigenvector.

o Singular Value Decomposition:

svdval –— compute the singular values of a m by n real matrix. sv2val –— compute the singular values of a real matrix efficiently for m >> n. svduv —— compute the singular values and the transformation matrices u and v for a real m by n matrix. sv2uv —— compute the singular values and transformation matrices efficiently for m >> n. svdu1v –— compute the singular values and transformation matrices u1 and v, where u1 overloads the input with the first n column vectors of u. sv2u1v –— compute the singular values and the transformation matrices u1 and v efficiently for m >> n.

 Singular value decomposition is extremely useful when dealing with linear

systems that may be singular. Singular values with values near zero are flags of a potential rank deficiency in the system matrix. They can be used to identify the presence of an ill-conditioned problem and, in some cases, to deal with the potential instability. They are applied to the linear least squares problem in this library. Singular values also define some important matrix norm parameters such as the 2-norm and the condition value. A complete decomposition provides both singular values and an orthogonal decomposition of vector spaces related to the matrix identifying the range and null-space. Fortunately, a highly stable algorithm based on Householder reduction to bidiagonal form and QR rotations can be used to implement the decomposition. The library provides two forms with one more efficient when the dimensions satisfy m > (3/2)n.

General Technical Comments

Efficient computation with matrices on modern processors must be

adapted to the storage scheme employed for matrix elements. The functions of this library segment do not employ the multidimensional array intrinsic of the C language. Access to elements employs the simple row-major scheme described here.

Matrices are modeled by the library functions as arrays with elements

stored in row order. Thus, the element in the jth row and kth column of the n by n matrix M, stored in the array mat[], is addressed by

      M[j,k] = mat[n*j+k]  , with   0 =< j,k <= n-1 .

(Remember that C employs zero as the starting index.) The storage order has important implications for data locality.

The algorithms employed here all have excellent numerical stability, and

the default double precision arithmetic of C enhances this. Thus, any problems encountered in using the matrix algebra functions will almost certainly be due to an ill-conditioned matrix. (The Hilbert matrices,

            H[i,j] = 1/(1+i+j)  for i,j < n

form a good example of such ill-conditioned systems.) We remind the reader that the appropriate response to such ill-conditioning is to seek an alternative approach to the problem. The option of increasing precision has already been exploited. Modification of the linear algebra algorithm code is not normally effective in an ill-conditioned problem.

---

FUNCTION SYNOPSES

Linear System Solutions:

Parameters

a	= array containing system matrix A in row order (altered to L-U factored form by computation)
b	= array containing system vector b at entry and solution vector x at exit
n	= dimension of system

Returns

0 -> normal exit; -1 -> singular input

Definition at line 184 of file ccmath_grass_wrapper.c.

G_math_solvps()

int G_math_solvps	(	double **	a,
		double *	b,
		int	n
	)

Solve a symmetric positive definite linear system S*x = b.

Parameters

a	= array containing system matrix S (altered to Cholesky upper right factor by computation)
b	= array containing system vector b as input and solution vector x as output
n	= dimension of system

Returns

: 0 -> normal exit; -1 -> input matrix not positive definite

Definition at line 198 of file ccmath_grass_wrapper.c.

G_math_solvru()

int G_math_solvru	(	double **	a,
		double *	b,
		int	n
	)

Definition at line 229 of file ccmath_grass_wrapper.c.

G_math_solvtd()

void G_math_solvtd	(	double *	a,
		double *	b,
		double *	c,
		double *	x,
		int	m
	)

Solve a tridiagonal linear system M*x = y.

Parameters

a	= array containing m+1 diagonal elements of M
b	= array of m elements below the main diagonal of M
c	= array of m elements above the main diagonal
x	= array containing the system vector y initially, and the solution vector at exit (m+1 elements)
m	= dimension parameter ( M is (m+1)x(m+1) )

Definition at line 214 of file ccmath_grass_wrapper.c.

G_math_sv2uv()

int G_math_sv2uv	(	double *	d,
		double **	a,
		double **	u,
		int	m,
		double **	v,
		int	n
	)

Compute the singular value transformation when m >> n.

Parameters

d	= pointer to double array of dimension n (output = singular values of A)
a	= pointer to store of the m by n input matrix A (A is altered by the computation)
u	= pointer to store for m by m orthogonal matrix U
v	= pointer to store for n by n orthogonal matrix V
m	= number of rows in A
n	= number of columns in A (m>=n required)

Returns

value: status flag with: 0 -> success -1 -> input error m < n

Definition at line 457 of file ccmath_grass_wrapper.c.

G_math_sv2val()

int G_math_sv2val	(	double *	d,
		double **	a,
		int	m,
		int	n
	)

Compute singular values when m >> n.

Parameters

d	= pointer to double array of dimension n (output = singular values of A)
a	= pointer to store of the m by n input matrix A (A is altered by the computation)
m	= number of rows in A
n	= number of columns in A (m>=n required)

Returns

value: status flag with: 0 -> success -1 -> input error m < n

Definition at line 423 of file ccmath_grass_wrapper.c.

G_math_svdu1v()

int G_math_svdu1v	(	double *	d,
		double **	a,
		int	m,
		double **	v,
		int	n
	)

Compute the singular value transformation with A overloaded by the partial U-matrix.

Parameters

d	= pointer to double array of dimension n (output = singular values of A)
a	= pointer to store of the m by n input matrix A (At output a is overloaded by the matrix U1 whose n columns are orthogonal vectors equal to the first n columns of U.)
v	= pointer to store for n by n orthogonal matrix V
m	= number of rows in A
n	= number of columns in A (m>=n required)

Returns

value: status flag with: 0 -> success -1 -> input error m < n

Definition at line 476 of file ccmath_grass_wrapper.c.

G_math_svduv()

int G_math_svduv	(	double *	d,
		double **	a,
		double **	u,
		int	m,
		double **	v,
		int	n
	)

Definition at line 440 of file ccmath_grass_wrapper.c.

G_math_svdval()

int G_math_svdval	(	double *	d,
		double **	a,
		int	m,
		int	n
	)

Compute the singular values of a real m by n matrix A.

Parameters

d	= pointer to double array of dimension n (output = singular values of A)
a	= pointer to store of the m by n input matrix A (A is altered by the computation)
m	= number of rows in A
n	= number of columns in A (m>=n required)

Returns

value: status flag with: 0 -> success -1 -> input error m < n

Definition at line 407 of file ccmath_grass_wrapper.c.