```
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.

```
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.

```
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.

```
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
`X`

will
be `*pdet`

. If any of the pivot elements are less than 1e-15, `X`

is
considered singular and no `T`

is computed; `X`

to be inverted
is not preserved during the

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

.

```
int mat_mtm(T,X,m,n)
double *T,*X;
int m,n;
```

computes `T`

=`X`

'`X`

, where `X`

is an `m`

*`n`

matrix.

```
int mat_mmt(T,X,m,n)
double *T,*X;
int m,n;
```

computes `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;
```

computes `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;
```

computes `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;
```

computes `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;
```

computes `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;
```

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.

```
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.

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

.

```
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.

```
int mat_chol(T,X,n)
double *T,*X;
int n;
```

performs the `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,

```
int mat_cholinv(A,n)
double *A;
int n;
```

inverts an `n`

*`n`

positive definite matrix `A`

by the `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 ```
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;
```

computes the `X`

=`S`

*`U`

for an `m`

*`n`

matrix `X`

(with `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;
```

makes the singular value decomposition `X`

=`U`

*`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;
```

reduces an `n`

*`n`

symmetric matrix `A`

to tridiagonal form using
`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
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;
```

finds the eigenvalues and vectors of the `n`

*`n`

tridiagonal matrix `A`

obtained by `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;
```

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 `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;
```

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 `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 */
```

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).

```
double sis_tulo(a,b,sa,sb,n)
double *a,*b;
int sa,sb,n;
```

is an assembler routine ```
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