# linalg

linalg - linear algebra primitives

`include "math.m";
math:= load Math Math->PATH;
dot : fn(x, y: array of real) : real;
norm1, norm2: fn(x: array of real) : real;
iamax : fn(x: array of real) : int;
gemm : fn(transa, transb: int, # upper case N or T
m, n, k: int, alpha : real,
a: array of real, lda : int,
b: array of real, ldb : int, beta: real,
c: array of real, ldc : int);
sort : fn(x: array of real, p: array of int);
`

## Description

These routines implement the basic functions of linear algebra. The standard vector inner product and norms are defined:

` dot(`*x*,* y*) = sum(*x* [*i*]** y* [*i*])
norm1(*x*) = sum(fabs(*x* [*i*]))
norm2(*x*) = sqrt(dot(*x*, *x*))

For the infinity norm, the function iamax(*x*) computes an integer *i* such that `fabs(x[i])` is maximal. These are all standard BLAS (basic linear algebra subroutines) except that in Inferno they apply only to contiguous (unit stride) vectors.
We assume the convention that matrices are represented as singly-subscripted arrays with Fortran storage order. So for an *m* by *n* matrix *A*, we use loops with 0<=* i* < *m* and 0<= *j* <* n* and access *A [i + m *j]*.

Rather than provide the huge number of entry points in BLAS2 and BLAS3, Inferno provides the matrix multiply primitive, gemm, defined by

* ***A = ***** A'*B' + *****C**

where the apostrophes indicate an optional transposition. As shown by the work of Kagstrom, Ling, and Van Loan, the other BLAS functionality can be built efficiently on top of gemm.

Matrix *a* ' is *m* by *k*, matrix *b* ' is *k* by *n*, and matrix *c* is *m* by *n*. Here *a* ' means *a* *(*a ') if *transa* is the constant 'N' ('T'), and similarly for *b* '. The sizes *m*, *n*, and *k* denote the 'active' part of the matrix. The parameters *lda*, *ldb*, and *ldc* denote the 'leading dimension' or size for purposes of indexing. So to loop over *c* use loops with 0<= *i* < *m* and 0<= *j* < *j* and access *a* [*i* +*ldc* **j*].

The sorting function `sort` (*x, p*) updates a 0-origin permutation *p* so that *x* [*p*[*i*]] <= x [p[i+1]], leaving *x* unchanged.

*Limbo Math Modules*
Bo Kagstrom, Per Ling, and Charles Van Loan. Portable high performance GEMM-based level 3 BLAS. In R.F. Sincovec et al., editor, parallel Processing for Scientific Computing, pages 339-346. SIAM Publications, 1993.

infernosupport@lucent.com
*Copyright © 1996,Lucent Technologies, Inc. All rights
reserved.*