.TH MATH-LINALG 2.SH NAMEMath: dot, norm1, norm2, iamax, gemm, sort \- linear algebra primitives.SH SYNOPSIS.EXinclude "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);.EE.SH DESCRIPTIONThese routines implement the basic functions of linear algebra.The standard vector inner product and norms are defined as follows:.IP.BI dot( "x , y" )=.BI sum( x [ i ]* y\fP[\fPi\fP]).IP.BI norm1( x )=.BI sum(fabs( x [ i\ \f5])).IP.BI norm2( x )=.BI sqrt(dot( "x , x" )).PPFor the infinity norm, the functioniamax(x)computes an integer.I isuch that.BI 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..PPWe assume the convention that matrices are represented assingly-subscripted arrays with Fortran storage order.So for an.I mby.I nmatrix.I Awe use loopswith.BI 0\(<= i < mand.BI 0\(<= j < nand access\f2A\f5[\f2i\f5+\f2m\f5*\f2j\f5]\f1..PPRather than provide the huge number of entry points in BLAS2 and BLAS3,Inferno provides the matrix multiply primitive,.BR gemm ,defined by.PP.EX \f2A\fP = \f1\(*a\fP*\f2A\fP\f1'\fP*\f2B\fP\f1'\fP + \f1\(*b\fP*\f2C\fP.EE.PPwhere the apostrophes indicate an optional transposition.As shown by thework of Kagstrom, Ling, and Van Loan, the other BLAS functionality canbe built efficiently on top of.BR gemm ..PPMatrix.IR a 'is.I mby.IR k ,matrix.IR b 'is.I kby.IR n ,and matrix.I cis.I mby.IR n .Here.IR a 'means.I a.RI ( a ')if.I transais the constant.B 'N'.RB ( 'T' ),and similarly for.IR b '.The sizes.IR m ,.IR n ,and.I kdenote the `active' part ofthe matrix.The parameters.IR lda ,.IR ldb ,and.I ldcdenote the `leading dimension'or size for purposes of indexing.So to loop over.I cuse loopswith.BI 0\(<= i < mand.BI 0\(<= j < nand access\f2a\f5[\f2i\f5+\f2ldc\f5*\f2j\f5]\f1..PPThe sorting function.BI sort( x , p )updates a 0-origin permutation.I pso that.IB x [ p [ i ]]\(<=.IB x [ p [ i +1]]\f1,leaving.I xunchanged..SH SOURCE.B /libinterp/math.c.SH SEE ALSO.IR math-intro (2)
.TH MATH-LINALG 2
.SH NAME
Math: dot, norm1, norm2, iamax, gemm, sort \- linear algebra primitives
.SH SYNOPSIS
.EX
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);
.EE
.SH DESCRIPTION
These routines implement the basic functions of linear algebra.
The standard vector inner product and norms are defined as follows:
.IP
.BI dot( "x , y" )
=
.BI sum( x [ i ]* y\fP[\fPi\fP])
.IP
.BI norm1( x )
=
.BI sum(fabs( x [ i\ \f5]))
.IP
.BI norm2( x )
=
.BI sqrt(dot( "x , x" ))
.PP
For the infinity norm, the function
iamax(x)
computes an integer
.I i
such that
.BI 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.
.PP
We assume the convention that matrices are represented as
singly-subscripted arrays with Fortran storage order.
So for an
.I m
by
.I n
matrix
.I A
we use loops
with
.BI 0\(<= i < m
and
.BI 0\(<= j < n
and access
\f2A\f5[\f2i\f5+\f2m\f5*\f2j\f5]\f1.
.PP
Rather than provide the huge number of entry points in BLAS2 and BLAS3,
Inferno provides the matrix multiply primitive,
.BR gemm ,
defined by
.PP
.EX
\f2A\fP = \f1\(*a\fP*\f2A\fP\f1'\fP*\f2B\fP\f1'\fP + \f1\(*b\fP*\f2C\fP
.EE
.PP
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
.BR gemm .
.PP
Matrix
.IR a '
is
.I m
by
.IR k ,
matrix
.IR b '
is
.I k
by
.IR n ,
and matrix
.I c
is
.I m
by
.IR n .
Here
.IR a '
means
.I a
.RI ( a ')
if
.I transa
is the constant
.B 'N'
.RB ( 'T' ),
and similarly for
.IR b '.
The sizes
.IR m ,
.IR n ,
and
.I k
denote the `active' part of
the matrix.
The parameters
.IR lda ,
.IR ldb ,
and
.I ldc
denote the `leading dimension'
or size for purposes of indexing.
So to loop over
.I c
use loops
with
.BI 0\(<= i < m
and
.BI 0\(<= j < n
and access
\f2a\f5[\f2i\f5+\f2ldc\f5*\f2j\f5]\f1.
.PP
The sorting function
.BI sort( x , p )
updates a 0-origin permutation
.I p
so that
.IB x [ p [ i ]]
\(<=
.IB x [ p [ i +1]]\f1,
leaving
.I x
unchanged.
.SH SOURCE
.B /libinterp/math.c
.SH SEE ALSO
.IR math-intro (2)