int mat_add(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
T=X+Y, where X and Y are m*n matrices.
int mat_sub(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
T=X-Y, where X and Y are m*n matrices.
int mat_mlt(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
T=X*Y, where X is an m*n and Y is an n*r matrix.
int mat_inv(T,X,n,pdet)
double *T,*X;
int n;
double *pdet;
T as the inverse matrix of an n*n matrix X by the
Gauss-Jordan elimination method. As a by-product, determinant of X will
be *pdet. If any of the pivot elements are less than 1e-15, X is
considered singular and no T is computed; -i will then be returned where
i is the current row index of the pivot element (i=0,1,...,n-1). In
non-singular cases, 1 is returned. Warning: The matrix X to be inverted
is not preserved during the mat_inv call.
int mat_transp(T,X,m,n)
double *T,*X;
int m,n;
m*n matrix X to an n*m matrix T.
int mat_mtm(T,X,m,n)
double *T,*X;
int m,n;
T=X'X, where X is an m*n matrix.
int mat_mmt(T,X,m,n)
double *T,*X;
int m,n;
T=XX', where X is an m*n matrix.
int mat_2mtm(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
T=X'Y, where X is an m*n matrix and Y is an m*r matrix.
int mat_2mmt(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
T=XY' where X is an m*n matrix and Y is an r*n matrix.
int mat_dmlt(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
T=X*Y, where X is an m*m diagonal matrix and Y is an m*n matrix.
int mat_mltd(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
T=X*Y, where X is an m*n matrix and Y is an n*n diagonal matrix.
int mat_center(T,X,m,n)
double *T,*X;
int m,n;
m*n matrix X by computing the means of the X columns as an
n element T vector and subtracting them from the corresponding columns.
int mat_nrm(T,X,m,n)
double *T,*X;
int m,n;
m*n matrix X to length 1. The original
column lengths (square root of sum of squares) will be stored as an
n element vector T. Columns of length=0 are not changed.
int mat_sum(T,X,m,n)
double *T,*X;
int m,n;
m*n matrix X as an n element vector T.
int mat_p(X,n,k)
double *X;
int n,k;
n*n matrix X by the pivotal operation by using
the diagonal element k (k=0,1,...,n-1) as the pivot.
int mat_chol(T,X,n)
double *T,*X;
int n;
n*n positive definite matrix
X. Hence an n*n lower triangular matrix T satisfying X=TT' will be
computed. If X is not positive definite, mat_chol returns -i, where i
(i=0,1,...,n-1) represents the column index where this assumption fails.
If decomposition is successful, 1 is returned.
int mat_cholinv(A,n)
double *A;
int n;
n*n positive definite matrix A by the Cholesky method and
writes the inverted matrix B partially on A according to the following
scheme:
Before mat_cholinv: (Here n=5 assumed)
0 1 2 n-1 n
0 a00 a01 a02 a03 a04 *
1 a10 a11 a12 a13 a14 *
2 a20 a21 a22 a23 a24 *
a30 a31 a32 a33 a34 *
n-1 a40 a41 a42 a43 a44 *
0 1 2 n-1 n
0 a00 b00 b01 b02 b03 b04
1 a10 a11 b11 b12 b13 b14
2 a20 a21 a22 b22 b23 b24
a30 a31 a32 a33 b33 b34
n-1 a40 a41 a42 a43 a44 b44
Please note that the elements are assumed to be saved columnwise.
To have enough space for B, at least n*(n+1) elements (of type double)
must have been allocated for A before the mat_cholinv call.
If A is not positive definite, -i (where i is the first column dependent
on previous ones) is returned. In a successful case 1 is returned.
To overwrite A by its inverse completely, use mat_cholmove(A,n) after
mat_cholinv(A,n) to obtain
0 1 2 n-1 n
0 b00 b01 b02 b03 b04 b04
1 b10 b11 b12 b13 b14 b14
2 b20 b21 b22 b23 b24 b24
b30 b31 b32 b33 b34 b34
n-1 b40 b41 b42 b43 b44 b44
int mat_gram_schmidt(S,U,X,m,n,tol)
double *S,*U,*X;
int m,n;
double tol;
X=S*U for an m*n matrix X
(with rank(X)=n<=m), where S is an m*n matrix with orthonormal
columns and U is n*n upper triangular.
The accuracy in checking the linear independency of columns is given by
tol. The value tol=1e-15 is recommended.
Return value -i indicates that column i (i=0,1,...,n-1) is linearly dependent on previous ones. After a successful decomposition, 1 is returned.
int mat_svd(X,D,V,m,n,eps,tol)
double *X,*D,*V;
int m,n;
double eps,tol;
X=U*diag(D)*V' for an m*n
matrix X with m>=n. After the mat_svd call X will be overwritten by an
m*n matrix U which is columnwise orthogonal. D will be an n element
vector consisting of singular values and V an n*n orthogonal matrix.
eps and tol are tolerance constants (see the source cited below).
Suitable values are eps=1e-16 and tol=(1e-300)/eps.
mat_svd has been written using the ALGOL procedure by G.H.Golub and C.Reinsch as the basis. See Handbook for Automatic Computation, Volume II, edited by J.H.Wilkinson and C.Reinsch, pp. 134-151 (Springer 1971).
int mat_tred2(d,e,A,n,tol)
double *d,*e,*A;
int n;
double tol;
n*n symmetric matrix A to tridiagonal form using
Householder's reduction. The diagonal of the result is stored as an n
element vector d and the sub-diagonal as the last n-1 elements of an n
element vector e. A will be overwritten by the transformation matrices.
tol is an accuracy constant (see mat_svd).
Space for d and e (n elements each of type double) must be allocated
before mat_tred2 is called.
To get the eigenvalues and vectors after mat_tred2(d,e,A,n,tol),
function mat_tql2 has to be called.
mat_tred2 and mat_tql2 have been written using the ALGOL procedures tred2 and tql2 as the basis. See Handbook for Automatic Computation, Volume II, edited by J.H.Wilkinson and C.Reinsch, (Springer 1971).
int mat_tql2(d,e,A,n,eps,maxiter)
double *d,*e,*A;
int n;
double eps;
int maxiter;
n*n tridiagonal matrix A
obtained by mat_tred2. Matrix A will be overwritten by the eigenvectors
and the eigenvalues will be saved in descending order as an n element
vector d.
eps in an accuracy constant (see mat_svd). Maximum number of iterations
for one eigenvalue is maxiter. maxiter=30 is recommended. In case of no
convergence within maxiter iterations, -1 is returned. If the
eigenvalues and vectors are obtained, 1 is the return value.
mat_tred2 and mat_tql2 have been written using the ALGOL procedures tred2 and tql2 as the basis. See Handbook for Automatic Computation, Volume II, edited by J.H.Wilkinson and C.Reinsch, (Springer 1971).
int solve_upper(X,A,B,m,k,eps)
double *X,*A,*B;
int m,k;
double eps;
AX=B where A is an m*m upper
triangular matrix and B is an m*k matrix. Before calling solve_upper,
space must also be allocated to the m*k solution matrix X.
If any of the pivot elements is smaller than eps, solve_upper returns -i
where i=0,1,...,m-1 is the current column. After a successful solution,
1 is returned.
solve_lower works as solve_upper but with an m*m lower triangular
matrix A.
solve_diag works as solve_upper but with an m*m diagonal matrix A.
int solve_symm(X,A,B,m,k,eps)
double *X,*A,*B;
int m,k;
double eps;
AX=B where A is an m*m positive
definite matrix and B is an m*k matrix. Before calling solve_symm,
space must also be allocated to the m*k solution matrix X.
If any of the pivot elements is smaller than eps, solve_symm returns -i,
where i=0,1,...,m-1 is the current column. After a successful solution,
1 is returned. If A is not positive definite, solve_symm calls
ortholin1.
solve_symm is based on the ALGOL procedures choldet1 and cholsol1 in Handbook for Automatic Computation, Volume II, edited by J.H.Wilkinson and C.Reinsch, (Springer 1971).
int ortholin1(A,n,m,B,k,eps,X,improvement)
double *A;
int n,m;
double *B;
int k;
double eps;
double *X;
int improvement;
/* iterative improvement 1=yes 0=no */
AX=B, where A is an n*m matrix, B an
n*k matrix and n>=m.
eps is the maximal relative rounding error (typically eps=1e-15).
ortholin1 is based on the ALGOL procedure ortholin1 in
Handbook for Automatic Computation, Volume II, edited by
J.H.Wilkinson and C.Reinsch, (Springer 1971).
double sis_tulo(a,b,sa,sb,n)
double *a,*b;
int sa,sb,n;
a[0]*b[0]+a[sa]*b[sb]+a[2*sa]*b[2*sb]+ ... +a[(n-1)*sa]*b[(n-1)*sb]
To speed up computations, many of the SURVOMAT.LIB functions use sis_tulo for scalar products.