```
MATH-LINALG(2)                                     MATH-LINALG(2)

NAME
Math: dot, norm1, norm2, iamax, gemm, sort - linear algebra
primitives

SYNOPSIS
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 as follows:

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

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*A'*B' + 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

Page 1                       Plan 9             (printed 10/2/22)

MATH-LINALG(2)                                     MATH-LINALG(2)

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 permuta-
tion p so that x[p[i]] _< x[p[i+1]], leaving x unchanged.

SOURCE
/libinterp/math.c