1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
|
.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)
|
