## Functions in library SURVOMAT.LIB

``````int mat_add(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
``````
computes `T`=`X`+`Y`, where `X` and `Y` are `m`*`n` matrices.
##### mat_sub

``````int mat_sub(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
``````
computes `T`=`X`-`Y`, where `X` and `Y` are `m`*`n` matrices.
mat_sub always returns 1.
##### mat_mlt

``````int mat_mlt(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
``````
computes `T`=`X`*`Y`, where `X` is an `m`*`n` and `Y` is an `n`*`r` matrix.
mat_mlt always returns 1.
##### mat_inv

``````int mat_inv(T,X,n,pdet)
double *T,*X;
int n;
double *pdet;
``````
computes matrix `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.
##### mat_transp

``````int mat_transp(T,X,m,n)
double *T,*X;
int m,n;
``````
transposes an `m`*`n` matrix `X` to an `n`*`m` matrix `T`.
mat_transp always returns 1.
##### mat_mtm

``````int mat_mtm(T,X,m,n)
double *T,*X;
int m,n;
``````
computes `T`=`X`'`X`, where `X` is an `m`*`n` matrix.
mat_mtm always returns 1.
##### mat_mmt

``````int mat_mmt(T,X,m,n)
double *T,*X;
int m,n;
``````
computes `T`=`XX`', where `X` is an `m`*`n` matrix.
mat_mmt always returns 1.
##### mat_2mtm

``````int mat_2mtm(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
``````
computes `T`=`X`'`Y`, where `X` is an `m`*`n` matrix and `Y` is an `m`*`r` matrix.
mat_2mtm always returns 1.
##### mat_2mmt

``````int mat_2mmt(T,X,Y,m,n,r)
double *T,*X,*Y;
int m,n,r;
``````
computes `T`=`XY`' where `X` is an `m`*`n` matrix and `Y` is an `r`*`n` matrix.
mat_2mmt always returns 1.
##### mat_dmlt

``````int mat_dmlt(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
``````
computes `T`=`X`*`Y`, where `X` is an `m`*`m` diagonal matrix and `Y` is an `m`*`n` matrix.
mat_dmlt always returns 1.
##### mat_mltd

``````int mat_mltd(T,X,Y,m,n)
double *T,*X,*Y;
int m,n;
``````
computes `T`=`X`*`Y`, where `X` is an `m`*`n` matrix and `Y` is an `n`*`n` diagonal matrix.
mat_mltd always returns 1.
##### mat_center

``````int mat_center(T,X,m,n)
double *T,*X;
int m,n;
``````
centers an `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.
mat_center always returns 1.
##### mat_nrm

``````int mat_nrm(T,X,m,n)
double *T,*X;
int m,n;
``````
normalizes the columns of an `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.
mat_nrm always returns 1.
##### mat_sum

``````int mat_sum(T,X,m,n)
double *T,*X;
int m,n;
``````
computes the column sums of an `m`*`n` matrix `X` as an `n` element vector `T`.
mat_sum always returns 1.
##### mat_p

``````int mat_p(X,n,k)
double *X;
int n,k;
``````
transforms the `n`*`n` matrix `X` by the pivotal operation by using the diagonal element `k` `(k=0,1,...,n-1)` as the pivot.
##### mat_chol

``````int mat_chol(T,X,n)
double *T,*X;
int n;
``````
performs the Cholesky decomposition of an `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.
##### mat_cholinv

``````int mat_cholinv(A,n)
double *A;
int n;
``````
inverts an `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   *
``````
After mat_cholinv:
``````     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.

##### mat_cholmove

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
``````
##### mat_gram_schmidt

``````int mat_gram_schmidt(S,U,X,m,n,tol)
double *S,*U,*X;
int m,n;
double tol;
``````
computes the Gram-Schmidt decomposition `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.

##### mat_svd

``````int mat_svd(X,D,V,m,n,eps,tol)
double *X,*D,*V;
int m,n;
double eps,tol;
``````
makes the singular value decomposition `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).

##### mat_tred2

``````int mat_tred2(d,e,A,n,tol)
double *d,*e,*A;
int n;
double tol;
``````
reduces an `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).

##### mat_tql2

``````int mat_tql2(d,e,A,n,eps,maxiter)
double *d,*e,*A;
int n;
double eps;
int maxiter;
``````
finds the eigenvalues and vectors of the `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).

##### solve_upper, solve_lower, solve_diag

``````int solve_upper(X,A,B,m,k,eps)
double *X,*A,*B;
int m,k;
double eps;
``````
solves the system of linear equations `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`.

##### solve_symm

``````int solve_symm(X,A,B,m,k,eps)
double *X,*A,*B;
int m,k;
double eps;
``````
solves the system of linear equations `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).

##### ortholin1

``````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 */
``````
gives least squares solutions for `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).

##### sis_tulo

``````double sis_tulo(a,b,sa,sb,n)
double *a,*b;
int sa,sb,n;
``````
is an assembler routine (written by T.Patovaara) for computation of the inner product
``````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.

Front page of Survo C libraries