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

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