Matrix Manipulation
*******************

There are a number of functions available for checking to see if the elements of a matrix meet some condition, and for rearranging the elements of a matrix.  For example, Octave can easily tell you if all the elements of a matrix are finite, or are less than some specified value.  Octave can also rotate the elements, extract the upper- or lower-triangular parts, or sort the columns of a matrix.

Finding Elements and Checking Conditions
========================================

The functions any and all are useful for determining whether any or all of the elements of a matrix satisfy some condition. The find function is also useful in determining which elements of a matrix meet a specified condition.

 -- TF =: any (X)
 -- TF =: any (X, DIM)
     For a vector argument, return true (logical 1) if any element of the vector is nonzero.

     For a matrix argument, return a row vector of logical ones and zeros with each element indicating whether any of the elements of the corresponding column of the matrix are nonzero.  For example:

          any (eye (2, 4))
           => [ 1, 1, 0, 0 ]

     If the optional argument DIM is supplied, work along dimension DIM. For example:

          any (eye (2, 4), 2)
           => [ 1; 1 ]

     all

 -- TF =: all (X)
 -- TF =: all (X, DIM)
     For a vector argument, return true (logical 1) if all elements of the vector are nonzero.

     For a matrix argument, return a row vector of logical ones and zeros with each element indicating whether all of the elements of the corresponding column of the matrix are nonzero.  For example:

          all ([2, 3; 1, 0])
              => [ 1, 0 ]

     If the optional argument DIM is supplied, work along dimension DIM. any

Since the comparison operators (see Comparison Ops) return matrices of ones and zeros, it is easy to test a matrix for many things, not just whether the elements are nonzero.  For example,

     all (all (rand (5) < 0.9))
          => 0

tests a random 5 by 5 matrix to see if all of its elements are less than 0.9.

Note that in conditional contexts (like the test clause of if and while statements) Octave treats the test as if you had typed all (all (condition)).

 -- Z =: xor (X, Y)
 -- Z =: xor (X1, X2, ...)
     Return the exclusive or of X and Y.

     For boolean expressions X and Y, xor (X, Y) is true if and only if one of X or Y is true.  Otherwise, if X and Y are both true or both false, xor returns false.

     The truth table for the xor operation is

       X   Y   Z
       -   -   -
       0   0   0
       1   0   1
       0   1   1
       1   1   0

     If more than two arguments are given the xor operation is applied cumulatively from left to right:

          (...((x1 XOR x2) XOR x3) XOR ...)

     and, or, not

 -- Y =: diff (X)
 -- Y =: diff (X, K)
 -- Y =: diff (X, K, DIM)
     If X is a vector of length n, diff (X) is the vector of first differences X(2) - X(1), ..., X(n) - X(n-1).

     If X is a matrix, diff (X) is the matrix of column differences along the first non-singleton dimension.

     The second argument is optional.  If supplied, diff (X, K), where K is a non-negative integer, returns the K-th differences.  It is possible that K is larger than the first non-singleton dimension of the matrix.  In this case, diff continues to take the differences along the next non-singleton dimension.

     The dimension along which to take the difference can be explicitly stated with the optional variable DIM.  In this case the K-th order differences are calculated along this dimension. In the case where K exceeds size (X, DIM) an empty matrix is returned. sort, merge

 -- TF =: isinf (X)
     Return a logical array which is true where the elements of X are infinite and false where they are not.

     For example:

          isinf ([13, Inf, NA, NaN])
                => [ 0, 1, 0, 0 ]

     isfinite, isnan, isna

 -- TF =: isnan (X)
     Return a logical array which is true where the elements of X are NaN values and false where they are not.

     NA values are also considered NaN values.  For example:

          isnan ([13, Inf, NA, NaN])
                => [ 0, 0, 1, 1 ]

     isna, isinf, isfinite

 -- TF =: isfinite (X)
     Return a logical array which is true where the elements of X are finite values and false where they are not.

     For example:

          isfinite ([13, Inf, NA, NaN])
               => [ 1, 0, 0, 0 ]

     isinf, isnan, isna

 -- [ERR, YI, ...] =: common_size (XI, ...)
     Determine if all input arguments are either scalar or of common size.

     If true, ERR is zero, and YI is a matrix of the common size with all entries equal to XI if this is a scalar or XI otherwise.  If the inputs cannot be brought to a common size, ERR is 1, and YI is XI.  For example:

          [err, a, b] = common_size ([1 2; 3 4], 5)
               => err = 0
               => a = [ 1, 2; 3, 4 ]
               => b = [ 5, 5; 5, 5 ]

     This is useful for implementing functions where arguments can either be scalars or of common size. size, size_equal, numel, ndims

 -- IDX =: find (X)
 -- IDX =: find (X, N)
 -- IDX =: find (X, N, DIRECTION)
 -- [i, j] =: find (...)
 -- [i, j, v] =: find (...)
     Return a vector of indices of nonzero elements of a matrix, as a row if X is a row vector or as a column otherwise.

     To obtain a single index for each matrix element, Octave pretends that the columns of a matrix form one long vector (like Fortran arrays are stored). For example:

          find (eye (2))
            => [ 1; 4 ]

     If two inputs are given, N indicates the maximum number of elements to find from the beginning of the matrix or vector.

     If three inputs are given, DIRECTION should be one of "first" or "last", requesting only the first or last N indices, respectively.  However, the indices are always returned in ascending order.

     If two outputs are requested, find returns the row and column indices of nonzero elements of a matrix.  For example:

          [i, j] = find (2 * eye (2))
              => i = [ 1; 2 ]
              => j = [ 1; 2 ]

     If three outputs are requested, find also returns a vector containing the nonzero values.  For example:

          [i, j, v] = find (3 * eye (2))
                 => i = [ 1; 2 ]
                 => j = [ 1; 2 ]
                 => v = [ 3; 3 ]

     If X is a multi-dimensional array of size m x n x p x ..., J contains the column locations as if X was flattened into a two-dimensional matrix of size m x (n + p + ...).

     Note that this function is particularly useful for sparse matrices, as it extracts the nonzero elements as vectors, which can then be used to create the original matrix.  For example:

          sz = size (a);
          [i, j, v] = find (a);
          b = sparse (i, j, v, sz(1), sz(2));

     nonzeros

 -- IDX =: lookup (TABLE, Y)
 -- IDX =: lookup (TABLE, Y, OPT)
     Lookup values in a sorted table.

     This function is usually used as a prelude to interpolation.

     If table is increasing, of length N and idx = lookup (table, y), then table(idx(i)) <= y(i) < table(idx(i+1)) for all y(i) within the table.  If y(i) < table(1) then idx(i) is 0.  If y(i) >= table(end) or isnan (y(i)) then idx(i) is N.

     If the table is decreasing, then the tests are reversed.  For non-strictly monotonic tables, empty intervals are always skipped.  The result is undefined if TABLE is not monotonic, or if TABLE contains a NaN.

     The complexity of the lookup is O(M*log(N)) where M is the size of Y. In the special case when Y is also sorted, the complexity is O(min (M*log(N), M+N)).

     TABLE and Y can also be cell arrays of strings (or Y can be a single string).  In this case, string lookup is performed using lexicographical comparison.

     If OPTS is specified, it must be a string with letters indicating additional options.

     m
          Match.  table(idx(i)) == y(i) if y(i) occurs in table; otherwise, idx(i) is zero.

     b
          Boolean.  idx(i) is a logical 1 or 0, indicating whether y(i) is contained in table or not.

     l
          Left.  For numeric lookups the leftmost subinterval shall be extended to minus infinity (i.e., all indices at least 1).

     r
          Right.  For numeric lookups the rightmost subinterval shall be extended to infinity (i.e., all indices at most N-1).

     Note: If TABLE is not sorted the results from lookup will be unpredictable.

If you wish to check if a variable exists at all, instead of properties its elements may have, consult Status of Variables.

Rearranging Matrices
====================

 -- B =: fliplr (A)
     Flip array left to right.

     Return a copy of A with the order of the columns reversed.  In other words, A is flipped left-to-right about a vertical axis.  For example:

          fliplr ([1, 2; 3, 4])
               =>  2  1
                   4  3

     flipud, flip, rot90, rotdim

 -- B =: flipud (A)
     Flip array upside down.

     Return a copy of A with the order of the rows reversed.  In other words, A is flipped upside-down about a horizontal axis.  For example:

          flipud ([1, 2; 3, 4])
               =>  3  4
                   1  2

     fliplr, flip, rot90, rotdim

 -- B =: flip (A)
 -- B =: flip (A, DIM)
     Return a copy of array A flipped across dimension DIM.

     If DIM is unspecified it defaults to the first non-singleton dimension.

     Examples:

          ## row vector
          flip ([1  2  3  4])
                =>  4  3  2  1

          ## column vector
          flip ([1; 2; 3; 4])
                =>  4
                    3
                    2
                    1

          ## 2-D matrix along dimension 1
          flip ([1 2; 3 4])
                =>  3  4
                    1  2

          ## 2-D matrix along dimension 2
          flip ([1 2; 3 4], 2)
                =>  2  1
                    4  3

     fliplr, flipud, rot90, rotdim, permute, transpose

 -- B =: rot90 (A)
 -- B =: rot90 (A, K)
     Rotate array by 90 degree increments.

     Return a copy of A with the elements rotated counterclockwise in 90-degree increments.

     The second argument is optional, and specifies how many 90-degree rotations are to be applied (the default value is 1).  Negative values of K rotate the matrix in a clockwise direction. For example,

          rot90 ([1, 2; 3, 4], -1)
              =>  3  1
                  4  2

     rotates the given matrix clockwise by 90 degrees.  The following are all equivalent statements:

          rot90 ([1, 2; 3, 4], -1)
          rot90 ([1, 2; 3, 4], 3)
          rot90 ([1, 2; 3, 4], 7)

     The rotation is always performed on the plane of the first two dimensions, i.e., rows and columns.  To perform a rotation on any other plane, use rotdim.

     rotdim, fliplr, flipud, flip

 -- B =: rotdim (A)
 -- B =: rotdim (A, N)
 -- B =: rotdim (A, N, PLANE)
     Return a copy of A with the elements rotated counterclockwise in 90-degree increments.

     The second argument N is optional, and specifies how many 90-degree rotations are to be applied (the default value is 1).  Negative values of N rotate the matrix in a clockwise direction.

     The third argument is also optional and defines the plane of the rotation. If present, PLANE is a two element vector containing two different valid dimensions of the matrix.  When PLANE is not given the first two non-singleton dimensions are used.

     For example,

          rotdim ([1, 2; 3, 4], -1, [1, 2])
               =>  3  1
                   4  2

     rotates the given matrix clockwise by 90 degrees.  The following are all equivalent statements:

          rotdim ([1, 2; 3, 4], -1, [1, 2])
          rotdim ([1, 2; 3, 4], 3, [1, 2])
          rotdim ([1, 2; 3, 4], 7, [1, 2])

     rot90, fliplr, flipud, flip

 -- A =: cat (DIM, ARRAY1, ARRAY2, ..., ARRAYN)
     Return the concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension DIM.

          A = ones (2, 2);
          B = zeros (2, 2);
          cat (2, A, B)
            => 1 1 0 0
               1 1 0 0

     Alternatively, we can concatenate A and B along the second dimension in the following way:

          [A, B]

     DIM can be larger than the dimensions of the N-D array objects and the result will thus have DIM dimensions as the following example shows:

          cat (4, ones (2, 2), zeros (2, 2))
            => ans(:,:,1,1) =

                 1 1
                 1 1

               ans(:,:,1,2) =

                 0 0
                 0 0

     horzcat, vertcat

 -- A =: horzcat (ARRAY1, ARRAY2, ..., ARRAYN)
     Return the horizontal concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension 2.

     Arrays may also be concatenated horizontally using the syntax for creating new matrices.  For example:

          A = [ ARRAY1, ARRAY2, ... ]

     This syntax is slightly more efficient because the Octave parser can concatenate the arrays without the overhead of a function call. cat, vertcat

 -- A =: vertcat (ARRAY1, ARRAY2, ..., ARRAYN)
     Return the vertical concatenation of N-D array objects, ARRAY1, ARRAY2, ..., ARRAYN along dimension 1.

     Arrays may also be concatenated vertically using the syntax for creating new matrices.  For example:

          A = [ ARRAY1; ARRAY2; ... ]

     This syntax is slightly more efficient because the Octave parser can concatenate the arrays without the overhead of a function call. cat, horzcat

 -- B =: permute (A, PERM)
     Return the generalized transpose for an N-D array object A.

     The permutation vector PERM must contain the elements 1:ndims (A) (in any order, but each element must appear only once).  The Nth dimension of A gets remapped to dimension PERM(N).  For example:

          X = zeros ([2, 3, 5, 7]);
          size (X)
             =>  2   3   5   7

          size (permute (X, [2, 1, 3, 4]))
             =>  3   2   5   7

          size (permute (X, [1, 3, 4, 2]))
             =>  2   5   7   3

          ## The identity permutation
          size (permute (X, [1, 2, 3, 4]))
             =>  2   3   5   7

     ipermute

 -- A =: ipermute (B, IPERM)
     The inverse of the permute function.

     The expression

          ipermute (permute (A, perm), perm)

     returns the original array A. permute

 -- B =: reshape (A, M, N, ...)
 -- B =: reshape (A, [M N ...])
 -- B =: reshape (A, ..., [], ...)
 -- B =: reshape (A, SIZE)
     Return a matrix with the specified dimensions (M, N, ...) whose elements are taken from the matrix A.

     The elements of the matrix are accessed in column-major order (like Fortran arrays are stored).

     The following code demonstrates reshaping a 1x4 row vector into a 2x2 square matrix.

          reshape ([1, 2, 3, 4], 2, 2)
                =>  1  3
                    2  4

     Note that the total number of elements in the original matrix (prod (size (A))) must match the total number of elements in the new matrix (prod ([M N ...])).

     A single dimension of the return matrix may be left unspecified and Octave will determine its size automatically.  An empty matrix ([]) is used to flag the unspecified dimension. resize, vec, postpad, cat, squeeze

 -- B =: resize (A, M)
 -- B =: resize (A, M, N, ...)
 -- B =: resize (A, [M N ...])
     Resize A cutting off elements as necessary.

     In the result, element with certain indices is equal to the corresponding element of A if the indices are within the bounds of A; otherwise, the element is set to zero.

     In other words, the statement

          B = resize (A, dv)

     is equivalent to the following code:

          B = zeros (dv, class (A));
          sz = min (dv, size (A));
          for i = 1:length (sz)
            idx{i} = 1:sz(i);
          endfor
          B(idx{:}) = A(idx{:});

     but is performed more efficiently.

     If only M is supplied, and it is a scalar, the dimension of the result is M-by-M.  If M, N, ... are all scalars, then the dimensions of the result are M-by-N-by-....  If given a vector as input, then the dimensions of the result are given by the elements of that vector.

     An object can be resized to more dimensions than it has; in such case the missing dimensions are assumed to be 1.  Resizing an object to fewer dimensions is not possible. reshape, postpad, prepad, cat

 -- Y =: circshift (X, N)
 -- Y =: circshift (X, N, DIM)
     Circularly shift the values of the array X.

     N must be a vector of integers no longer than the number of dimensions in X.  The values of N can be either positive or negative, which determines the direction in which the values of X are shifted.  If an element of N is zero, then the corresponding dimension of X will not be shifted.  If N is a scalar and no DIM is specified then the shift is applied to the first non-singular dimension.

     If a scalar DIM is given then operate along the specified dimension. In this case N must be a scalar as well.

     Examples:

          x = [1, 2, 3;
               4, 5, 6;
               7, 8, 9];
          ## positive shift on rows (1st non-singular dim)
          circshift (x, 1)
            =>
                 7   8   9
                 1   2   3
                 4   5   6
          ## negative shift on rows (1st non-singular dim)
          circshift (x, -2)
            =>
                 7   8   9
                 1   2   3
                 4   5   6
          ## no shift of rows, shift columns by 1 (2nd dimension)
          circshift (x, [0,1])
            =>
                 3   1   2
                 6   4   5
                 9   7   8
          ## shift columns (2nd dimension)
          circshift (x, 1, 2)
            =>
                 3   1   2
                 6   4   5
                 9   7   8

     permute, ipermute, shiftdim

 -- Y =: shiftdim (X, N)
 -- [Y, NS] =: shiftdim (X)
     Shift the dimensions of X by N, where N must be an integer scalar.

     When N is positive, the dimensions of X are shifted to the left, with the leading dimensions circulated to the end.  If N is negative, then the dimensions of X are shifted to the right, with N leading singleton dimensions added.

     Called with a single argument, shiftdim, removes the leading singleton dimensions, returning the number of dimensions removed in the second output argument NS.

     For example:

          x = ones (1, 2, 3);
          size (shiftdim (x, -1))
            =>   1   1   2   3
          size (shiftdim (x, 1))
            =>   2   3
          [b, ns] = shiftdim (x)
            => b =
                  1   1   1
                  1   1   1
            => ns = 1

     reshape, permute, ipermute, circshift, squeeze

 -- [S, I] =: sort (X)
 -- [S, I] =: sort (X, DIM)
 -- [S, I] =: sort (X, MODE)
 -- [S, I] =: sort (X, DIM, MODE)
     Return a copy of X with the elements arranged in increasing order.

     For matrices, sort orders the elements within columns

     For example:

          sort ([1, 2; 2, 3; 3, 1])
             =>  1  1
                 2  2
                 3  3

     If the optional argument DIM is given, then the matrix is sorted along the dimension defined by DIM.  The optional argument MODE defines the order in which the values will be sorted.  Valid values of MODE are "ascend" or "descend".

     The sort function may also be used to produce a matrix containing the original row indices of the elements in the sorted matrix.  For example:

          [s, i] = sort ([1, 2; 2, 3; 3, 1])
            => s = 1  1
                   2  2
                   3  3
            => i = 1  3
                   2  1
                   3  2

     For equal elements, the indices are such that equal elements are listed in the order in which they appeared in the original list.

     Sorting of complex entries is done first by magnitude (abs (Z)) and for any ties by phase angle (angle (z)).  For example:

          sort ([1+i; 1; 1-i])
              => 1 + 0i
                 1 - 1i
                 1 + 1i

     NaN values are treated as being greater than any other value and are sorted to the end of the list.

     The sort function may also be used to sort strings and cell arrays of strings, in which case ASCII dictionary order (uppercase 'A' precedes lowercase 'a') of the strings is used.

     The algorithm used in sort is optimized for the sorting of partially ordered lists. sortrows, issorted

 -- [S, I] =: sortrows (A)
 -- [S, I] =: sortrows (A, C)
     Sort the rows of the matrix A according to the order of the columns specified in C.

     By default (C omitted, or a particular column unspecified in C) an ascending sort order is used.  However, if elements of C are negative then the corresponding column is sorted in descending order.  If the elements of A are strings then a lexicographical sort is used.

     Example: sort by column 2 in descending order, then 3 in ascending order

          x = [ 7, 1, 4;
                8, 3, 5;
                9, 3, 6 ];
          sortrows (x, [-2, 3])
             => 8  3  5
                9  3  6
                7  1  4

     sort

 -- TF =: issorted (A)
 -- TF =: issorted (A, MODE)
 -- TF =: issorted (A, "rows", MODE)
     Return true if the vector A is sorted according to MODE, which may be either "ascend", "descend", or "either".

     By default, MODE is "ascend".  NaNs are treated in the same manner as sort.

     If the optional argument "rows" is supplied, check whether the matrix is sorted by rows as output by the function sortrows (with no options).

     This function does not support sparse matrices. sort, sortrows

 -- NEL =: nth_element (X, N)
 -- NEL =: nth_element (X, N, DIM)
     Select the n-th smallest element of a vector, using the ordering defined by sort.

     The result is equivalent to sort(X)(N).

     N can also be a contiguous range, either ascending l:u or descending u:-1:l, in which case a range of elements is returned.

     If X is an array, nth_element operates along the dimension defined by DIM, or the first non-singleton dimension if DIM is not given.

     Programming Note: nth_element encapsulates the C++ standard library algorithms nth_element and partial_sort.  On average, the complexity of the operation is O(M*log(K)), where M = size (X, DIM) and K = length (N).  This function is intended for cases where the ratio K/M is small; otherwise, it may be better to use sort. sort, min, max

 -- A_LO =: tril (A)
 -- A_LO =: tril (A, K)
 -- A_LO =: tril (A, K, PACK)
     Return a new matrix formed by extracting the lower triangular part of the matrix A, and setting all other elements to zero.

     The optional second argument specifies how many diagonals above or below the main diagonal should also be set to zero.  The default value of K is zero which includes the main diagonal as part of the result.  If the value of K is a nonzero integer then the selection of elements starts at an offset of K diagonals above the main diagonal for positive K or below the main diagonal for negative K.  The absolute value of K may not be greater than the number of subdiagonals or superdiagonals.

     Example 1 : exclude main diagonal

          tril (ones (3), -1)
               =>  0  0  0
                   1  0  0
                   1  1  0

     Example 2 : include first superdiagonal

          tril (ones (3), 1)
               =>  1  1  0
                   1  1  1
                   1  1  1

     If the optional third argument "pack" is given then the extracted elements are not inserted into a matrix, but instead stacked column-wise one above another, and returned as a column vector. triu, istril, diag

 -- A_UP =: triu (A)
 -- A_UP =: triu (A, K)
 -- A_UP =: triu (A, K, PACK)
     Return a new matrix formed by extracting the upper triangular part of the matrix A, and setting all other elements to zero.

     The optional second argument specifies how many diagonals above or below the main diagonal should also be set to zero.  The default value of K is zero which includes the main diagonal as part of the result.  If the value of K is a nonzero integer then the selection of elements starts at an offset of K diagonals above the main diagonal for positive K or below the main diagonal for negative K.  The absolute value of K may not be greater than the number of subdiagonals or superdiagonals.

     Example 1 : exclude main diagonal

          triu (ones (3), 1)
               =>  0  1  1
                   0  0  1
                   0  0  0

     Example 2 : include first subdiagonal

          triu (ones (3), -1)
               =>  1  1  1
                   1  1  1
                   0  1  1

     If the optional third argument "pack" is given then the extracted elements are not inserted into a matrix, but instead stacked column-wise one above another, and returned as a column vector. tril, istriu, diag

 -- V =: vec (X)
 -- V =: vec (X, DIM)
     Return the vector obtained by stacking the columns of the matrix X one above the other.

     Without DIM this is equivalent to X(:).

     If DIM is supplied, the dimensions of V are set to DIM with all elements along the last dimension.  This is equivalent to shiftdim (X(:), 1-DIM). vech, resize, cat

 -- V =: vech (X)
     Return the vector obtained by eliminating all superdiagonal elements of the square matrix X and stacking the result one column above the other.

     This has uses in matrix calculus where the underlying matrix is symmetric and it would be pointless to keep values above the main diagonal. vec

 -- B =: prepad (A, L)
 -- B =: prepad (A, L, C)
 -- B =: prepad (A, L, C, DIM)
     Prepend the scalar value C to the vector A until it is of length L.  If C is not given, a value of 0 is used.

     If length (A) > L, elements from the beginning of A are removed until a vector of length L is obtained.

     If A is a matrix, elements are prepended or removed from each row.

     If the optional argument DIM is given, operate along this dimension.

     If DIM is larger than the dimensions of A, the result will have DIM dimensions. postpad, cat, resize

 -- B =: postpad (A, L)
 -- B =: postpad (A, L, C)
 -- B =: postpad (A, L, C, DIM)
     Append the scalar value C to the vector A until it is of length L.  If C is not given, a value of 0 is used.

     If length (A) > L, elements from the end of A are removed until a vector of length L is obtained.

     If A is a matrix, elements are appended or removed from each row.

     If the optional argument DIM is given, operate along this dimension.

     If DIM is larger than the dimensions of A, the result will have DIM dimensions. prepad, cat, resize

 -- M =: diag (V)
 -- M =: diag (V, K)
 -- M =: diag (V, M, N)
 -- V =: diag (M)
 -- V =: diag (M, K)
     Return a diagonal matrix with vector V on diagonal K.

     The second argument is optional.  If it is positive, the vector is placed on the K-th superdiagonal.  If it is negative, it is placed on the -K-th subdiagonal.  The default value of K is 0, and the vector is placed on the main diagonal.  For example:

          diag ([1, 2, 3], 1)
             =>  0  1  0  0
                 0  0  2  0
                 0  0  0  3
                 0  0  0  0

     The 3-input form returns a diagonal matrix with vector V on the main diagonal and the resulting matrix being of size M rows x N columns.

     Given a matrix argument, instead of a vector, diag extracts the K-th diagonal of the matrix.

 -- M =: blkdiag (A, B, C, ...)
     Build a block diagonal matrix from A, B, C, ...

     All arguments must be numeric and either two-dimensional matrices or scalars.  If any argument is of type sparse, the output will also be sparse. diag, horzcat, vertcat, sparse

Special Utility Matrices
========================

 -- I =: eye (N)
 -- I =: eye (M, N)
 -- I =: eye ([M N])
 -- I =: eye (..., CLASS)
     Return an identity matrix.

     If invoked with a single scalar argument N, return a square NxN identity matrix.

     If supplied two scalar arguments (M, N), eye takes them to be the number of rows and columns.  If given a vector with two elements, eye uses the values of the elements as the number of rows and columns, respectively.  For example:

          eye (3)
           =>  1  0  0
               0  1  0
               0  0  1

     The following expressions all produce the same result:

          eye (2)
          ==
          eye (2, 2)
          ==
          eye (size ([1, 2; 3, 4]))

     The optional argument CLASS, allows eye to return an array of the specified type, like

          val = zeros (n,m, "uint8")

     Calling eye with no arguments is equivalent to calling it with an argument of 1.  Any negative dimensions are treated as zero.  These odd definitions are for compatibility with MATLAB. speye, ones, zeros

 -- VAL =: ones (N)
 -- VAL =: ones (M, N)
 -- VAL =: ones (M, N, K, ...)
 -- VAL =: ones ([M N ...])
 -- VAL =: ones (..., "like", VAR)
 -- VAL =: ones (..., CLASS)
     Return a matrix or N-dimensional array whose elements are all 1.

     If invoked with a single scalar integer argument N, return a square NxN matrix.

     If invoked with two or more scalar integer arguments, or a vector of integer values, return an array with the given dimensions.

     To create a constant matrix whose values are all the same use an expression such as

          val_matrix = val * ones (m, n)

     If a variable VAR is specified after "like", the output VAL will have the same data type, complexity, and sparsity as VAR.

     The optional argument CLASS specifies the class of the return array and defaults to double.  For example:

          val = ones (m,n, "uint8")

     zeros

 -- VAL =: zeros (N)
 -- VAL =: zeros (M, N)
 -- VAL =: zeros (M, N, K, ...)
 -- VAL =: zeros ([M N ...])
 -- VAL =: zeros (..., "like", VAR)
 -- VAL =: zeros (..., CLASS)
     Return a matrix or N-dimensional array whose elements are all 0.

     If invoked with a single scalar integer argument, return a square NxN matrix.

     If invoked with two or more scalar integer arguments, or a vector of integer values, return an array with the given dimensions.

     If a variable VAR is specified after "like", the output VAL will have the same data type, complexity, and sparsity as VAR.

     The optional argument CLASS specifies the class of the return array and defaults to double.  For example:

          val = zeros (m,n, "uint8")

     ones

 -- B =: repmat (A, M)
 -- B =: repmat (A, M, N)
 -- B =: repmat (A, M, N, P ...)
 -- B =: repmat (A, [M N])
 -- B =: repmat (A, [M N P ...])
     Repeat matrix or N-D array.

     Form a block matrix of size M by N, with a copy of matrix A as each element.

     If N is not specified, form an M by M block matrix.  For copying along more than two dimensions, specify the number of times to copy across each dimension M, N, P, ..., in a vector in the second argument.

     bsxfun, kron, repelems

 -- Y =: repelems (X, R)
     Construct a vector of repeated elements from X.

     R is a 2xN integer matrix specifying which elements to repeat and how often to repeat each element.  Entries in the first row, R(1,j), select an element to repeat.  The corresponding entry in the second row, R(2,j), specifies the repeat count.  If X is a matrix then the columns of X are imagined to be stacked on top of each other for purposes of the selection index.  A row vector is always returned.

     Conceptually the result is calculated as follows:

          y = [];
          for i = 1:columns (R)
            y = [y, X(R(1,i)*ones(1, R(2,i)))];
          endfor

     repmat, cat

 -- XXX =: repelem (X, R)
 -- XXX =: repelem (X, R_1, ..., R_N)
     Construct an array of repeated elements from X and repeat instructions R_1, ....

     X must be a scalar, vector, or N-dimensional array.

     A repeat instruction R_J must either be a scalar or a vector.  If the instruction is a scalar then each component of X in dimension J is repeated R_J times.  If the instruction is a vector then it must have the same number of elements as the corresponding dimension J of X.  In this case, the Kth component of dimension J is repeated R_J(K) times.

     If X is a scalar or vector then repelem may be called with just a single repeat instruction R and repelem will return a vector with the same orientation as the input.

     If X is a matrix then at least two R_Js must be specified.

     Note: Using repelem with a vector X and a vector for R_J is equivalent to Run Length Decoding.

     Examples:

          A = [1 2 3 4 5];
          B = [2 1 0 1 2];
          repelem (A, B)
            =>   1   1   2   4   5   5

          A = magic (3)
            => A =
                 8   1   6
                 3   5   7
                 4   9   2
          B1 = [1 2 3];
          B2 = 2;
          repelem (A, B1, B2)
            =>     8   8   1   1   6   6
                   3   3   5   5   7   7
                   3   3   5   5   7   7
                   4   4   9   9   2   2
                   4   4   9   9   2   2
                   4   4   9   9   2   2

     More R_J may be specified than the number of dimensions of X. Any excess R_J must be scalars (because X's size in those dimensions is only 1), and X will be replicated in those dimensions accordingly.

          A = [1 2 3 4 5];
          B1 = 2;
          B2 = [2 1 3 0 2];
          B3 = 3;
          repelem (A, B1, B2, B3)
            =>    ans(:,:,1) =
                     1   1   2   3   3   3   5   5
                     1   1   2   3   3   3   5   5

                  ans(:,:,2) =

                     1   1   2   3   3   3   5   5
                     1   1   2   3   3   3   5   5

                  ans(:,:,3) =
                     1   1   2   3   3   3   5   5
                     1   1   2   3   3   3   5   5

     R_J must be specified in order.  A placeholder of 1 may be used for dimensions which do not need replication.

          repelem ([-1, 0; 0, 1], 1, 2, 1, 2)
            =>  ans(:,:,1,1) =
                  -1  -1   0   0
                   0   0   1   1

                ans(:,:,1,2) =
                  -1  -1   0   0
                   0   0   1   1

     If fewer R_J are given than the number of dimensions in X, repelem will assume R_J is 1 for those dimensions.

          A = cat (3, [-1 0; 0 1], [-1 0; 0 1])
            =>  ans(:,:,1) =
                  -1   0
                   0   1

                ans(:,:,2) =
                  -1   0
                   0   1

          repelem (A,2,3)
            =>  ans(:,:,1) =
                  -1  -1  -1   0   0   0
                  -1  -1  -1   0   0   0
                   0   0   0   1   1   1
                   0   0   0   1   1   1

                ans(:,:,2) =
                  -1  -1  -1   0   0   0
                  -1  -1  -1   0   0   0
                   0   0   0   1   1   1
                   0   0   0   1   1   1

     repelem preserves the class of X, and works with strings, cell arrays, NA, and NAN inputs.  If any R_J is 0 the output will be an empty array.

          repelem ("Octave", 2, 3)
            =>    OOOccctttaaavvveee
                  OOOccctttaaavvveee

          repelem ([1 2 3; 1 2 3], 2, 0)
            =>    [](4x0)

     cat, kron, repmat

The functions linspace and logspace make it very easy to create vectors with evenly or logarithmically spaced elements. See Ranges.

 -- Y =: linspace (START, END)
 -- Y =: linspace (START, END, N)
     Return a row vector with N linearly spaced elements between START and END.

     If the number of elements N is greater than one, then the endpoints START and END are always included in the range.  If START is greater than END, the elements are stored in decreasing order.  If the number of points N is not specified, a value of 100 is used.

     The linspace function returns a row vector when both START and END are scalars.  If one, or both, inputs are vectors, then linspace transforms them to column vectors and returns a matrix where each row is an independent sequence between START(ROW_N), END(ROW_N).

     Programming Notes: For compatibility with MATLAB, return the second argument (END) when a single value (N = 1) is requested.  If N is not an integer then floor (N) is used to round the number of elements.  If N is zero or negative then an empty 1x0 matrix is returned. colon, logspace

 -- Y =: logspace (A, B)
 -- Y =: logspace (A, B, N)
 -- Y =: logspace (A, pi)
 -- Y =: logspace (A, pi, N)
     Return a row vector with N elements logarithmically spaced from 10^A to 10^B.

     If the number of elements N is unspecified it defaults to 50.

     If B is equal to pi, the points are between 10^A and pi, not 10^A and 10^pi, which is useful in digital signal processing.

     Programming Notes: For compatibility with MATLAB, return the right-hand side of the range (10^B) when a single value (N = 1) is requested. If N is not an integer then floor (N) is used to round the number of elements.  If N is zero or negative then an empty 1x0 matrix is returned. linspace

 -- X =: rand (N)
 -- X =: rand (M, N, ...)
 -- X =: rand ([M N ...])
 -- X =: rand (..., "single")
 -- X =: rand (..., "double")
 -- V =: rand ("state")
 -- rand: ("state", V)
 -- rand: ("state", "reset")
 -- V =: rand ("seed")
 -- rand: ("seed", V)
 -- rand: ("seed", "reset")
     Return a matrix with random elements uniformly distributed on the interval (0, 1).

     The arguments are handled the same as the arguments for eye.

     You can query the state of the random number generator using the form

          v = rand ("state")

     This returns a column vector V of length 625.  Later, you can restore the random number generator to the state V using the form

          rand ("state", v)

     You may also initialize the state vector from an arbitrary vector of length ≤ 625 for V.  This new state will be a hash based on the value of V, not V itself.

     By default, the generator is initialized by contributing entropy from the wall clock time, the CPU time, the current fraction of a second, the process ID and---if available---up to 1024 bits from the C++ random numbers source random_device, which might be non-deterministic (implementation specific).  Note that this differs from MATLAB, which always initializes the state to the same state at startup.  To obtain behavior comparable to MATLAB, initialize with a deterministic state vector in Octave's startup files (see Startup Files).

     To compute the pseudo-random sequence, rand uses the Mersenne Twister with a period of 2^19937-1 (See M. Matsumoto and T. Nishimura, Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator, ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, pp. 3--30, January 1998, http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html). Do not use for cryptography without securely hashing several returned values together, otherwise the generator state can be learned after reading 624 consecutive values.

     Older versions of Octave used a different random number generator. The new generator is used by default as it is significantly faster than the old generator, and produces random numbers with a significantly longer cycle time.  However, in some circumstances it might be desirable to obtain the same random sequences as produced by the old generators.  To do this the keyword "seed" is used to specify that the old generators should be used, as in

          rand ("seed", val)

     which sets the seed of the generator to VAL.  The seed of the generator can be queried with

          s = rand ("seed")

     However, it should be noted that querying the seed will not cause rand to use the old generators, only setting the seed will.  To cause rand to once again use the new generators, the keyword "state" should be used to reset the state of the rand.

     The state or seed of the generator can be reset to a new random value using the "reset" keyword.

     The class of the value returned can be controlled by a trailing "double" or "single" argument.  These are the only valid classes. randn, rande, randg, randp

 -- R =: randi (IMAX)
 -- R =: randi (IMAX, N)
 -- R =: randi (IMAX, M, N, ...)
 -- R =: randi ([IMIN IMAX], ...)
 -- R =: randi (..., "CLASS")
     Return random integers in the range 1:IMAX.

     Additional arguments determine the shape of the return matrix.  When no arguments are specified a single random integer is returned.  If one argument N is specified then a square matrix (N x N) is returned.  Two or more arguments will return a multi-dimensional matrix (M x N x ...).

     The integer range may optionally be described by a two-element matrix with a lower and upper bound in which case the returned integers will be on the interval [IMIN, IMAX].

     The optional argument CLASS will return a matrix of the requested type.  The default is "double".

     The following example returns 150 integers in the range 1--10.

          ri = randi (10, 150, 1)

     Implementation Note: randi relies internally on rand which uses class "double" to represent numbers.  This limits the maximum integer (IMAX) and range (IMAX - IMIN) to the value returned by the flintmax function.  For IEEE floating point numbers this value is 2^53 - 1.

     rand, randn

 -- X =: randn (N)
 -- X =: randn (M, N, ...)
 -- X =: randn ([M N ...])
 -- X =: randn (..., "single")
 -- X =: randn (..., "double")
 -- V =: randn ("state")
 -- randn: ("state", V)
 -- randn: ("state", "reset")
 -- V =: randn ("seed")
 -- randn: ("seed", V)
 -- randn: ("seed", "reset")
     Return a matrix with normally distributed random elements having zero mean and variance one.

     The arguments are handled the same as the arguments for rand.

     By default, randn uses the Marsaglia and Tsang ``Ziggurat technique'' to transform from a uniform to a normal distribution.

     The class of the value returned can be controlled by a trailing "double" or "single" argument.  These are the only valid classes.

     Reference: G. Marsaglia and W.W. Tsang, Ziggurat Method for Generating Random Variables, J. Statistical Software, vol 5, 2000, https://www.jstatsoft.org/v05/i08/

     rand, rande, randg, randp

 -- X =: rande (N)
 -- X =: rande (M, N, ...)
 -- X =: rande ([M N ...])
 -- X =: rande (..., "single")
 -- X =: rande (..., "double")
 -- V =: rande ("state")
 -- rande: ("state", V)
 -- rande: ("state", "reset")
 -- V =: rande ("seed")
 -- rande: ("seed", V)
 -- rande: ("seed", "reset")
     Return a matrix with exponentially distributed random elements.

     The arguments are handled the same as the arguments for rand.

     By default, rande uses the Marsaglia and Tsang ``Ziggurat technique'' to transform from a uniform to an exponential distribution.

     The class of the value returned can be controlled by a trailing "double" or "single" argument.  These are the only valid classes.

     Reference: G. Marsaglia and W.W. Tsang, Ziggurat Method for Generating Random Variables, J. Statistical Software, vol 5, 2000, https://www.jstatsoft.org/v05/i08/

     rand, randn, randg, randp

 -- X =: randp (L, N)
 -- X =: randp (L, M, N, ...)
 -- X =: randp (L, [M N ...])
 -- X =: randp (..., "single")
 -- X =: randp (..., "double")
 -- V =: randp ("state")
 -- randp: ("state", V)
 -- randp: ("state", "reset")
 -- V =: randp ("seed")
 -- randp: ("seed", V)
 -- randp: ("seed", "reset")
     Return a matrix with Poisson distributed random elements with mean value parameter given by the first argument, L.

     The arguments are handled the same as the arguments for rand, except for the argument L.

     Five different algorithms are used depending on the range of L and whether or not L is a scalar or a matrix.

     For scalar L ≤ 12, use direct method.
          W.H. Press, et al., Numerical Recipes in C, Cambridge University Press, 1992.

     For scalar L > 12, use rejection method.[1]
          W.H. Press, et al., Numerical Recipes in C, Cambridge University Press, 1992.

     For matrix L ≤ 10, use inversion method.[2]
          E. Stadlober, et al., WinRand source code, available via FTP.

     For matrix L > 10, use patchwork rejection method.
          E. Stadlober, et al., WinRand source code, available via FTP, or H. Zechner, Efficient sampling from continuous and discrete unimodal distributions, Doctoral Dissertation, 156pp., Technical University Graz, Austria, 1994.

     For L > 1e8, use normal approximation.
          L. Montanet, et al., Review of Particle Properties, Physical Review D 50 p1284, 1994.

     The class of the value returned can be controlled by a trailing "double" or "single" argument.  These are the only valid classes. rand, randn, rande, randg

 -- X =: randg (A, N)
 -- X =: randg (A, M, N, ...)
 -- X =: randg (A, [M N ...])
 -- X =: randg (..., "single")
 -- X =: randg (..., "double")
 -- V =: randg ("state")
 -- randg: ("state", V)
 -- randg: ("state", "reset")
 -- V =: randg ("seed")
 -- randg: ("seed", V)
 -- randg: ("seed", "reset")

     Return a matrix with gamma (A,1) distributed random elements.

     The arguments are handled the same as the arguments for rand, except for the argument A.

     This can be used to generate many distributions:

     gamma (a, b) for a > -1, b > 0

               r = b * randg (a)

     beta (a, b) for a > -1, b > -1

               r1 = randg (a, 1)
               r = r1 / (r1 + randg (b, 1))

     Erlang (a, n)

               r = a * randg (n)

     chisq (df) for df > 0

               r = 2 * randg (df / 2)

     t (df) for 0 < df < inf (use randn if df is infinite)

               r = randn () / sqrt (2 * randg (df / 2) / df)

     F (n1, n2) for 0 < n1, 0 < n2

               ## r1 equals 1 if n1 is infinite
               r1 = 2 * randg (n1 / 2) / n1
               ## r2 equals 1 if n2 is infinite
               r2 = 2 * randg (n2 / 2) / n2
               r = r1 / r2

     negative binomial (n, p) for n > 0, 0 < p <= 1

               r = randp ((1 - p) / p * randg (n))

     non-central chisq (df, L), for df >= 0 and L > 0
          (use chisq if L = 0)

               r = randp (L / 2)
               r(r > 0) = 2 * randg (r(r > 0))
               r(df > 0) += 2 * randg (df(df > 0)/2)

     Dirichlet (a1, ... ak)

               r = (randg (a1), ..., randg (ak))
               r = r / sum (r)

     The class of the value returned can be controlled by a trailing "double" or "single" argument.  These are the only valid classes. rand, randn, rande, randp

 -- rng: (SEED)
 -- rng: (SEED, "GENERATOR")
 -- rng: ("shuffle")
 -- rng: ("shuffle", "GENERATOR")
 -- rng: ("default")
 -- S =: rng ()
 -- rng: (S)
 -- S =: rng (...)
     Set or query the seed of the random number generator used by rand and randn.

     The input SEED is a scalar numeric value used to initialize the state vector of the random number generator.

     The optional string GENERATOR specifies the type of random number generator to be used.  Its value can be "twister", "v5uniform", or "v5normal".  The "twister" keyword is described below.  "v5uniform" and "v5normal" refer to older versions of Octave that used to use a different random number generator.

     The state or seed of the random number generator can be reset to a new random value using the "shuffle" keyword.

     The random number generator can be reset to default values using the "default" keyword.  The default values are to use the Mersenne Twister generator with a seed of 0.

     The optional return value S contains the state of the random number generator at the time the function is called (i.e., before it might be modified according to the input arguments).  It is encoded as a structure variable with three fields: "Type", "Seed", and "State".  The random number generator can be restored to the state S using rng (S).  This is useful when the identical sequence of pseudo-random numbers is required for an algorithm.

     By default, and with the "twister" option, pseudo-random sequences are computed using the Mersenne Twister with a period of 2^19937-1 (See M. Matsumoto and T. Nishimura, Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator, ACM Trans. on Modeling and Computer Simulation Vol. 8, No. 1, pp. 3--30, January 1998, http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/emt.html). Do not use for cryptography without securely hashing several returned values together, otherwise the generator state can be learned after reading 624 consecutive values.

     rand, randn

The generators operate in the new or old style together, it is not possible to mix the two.  Initializing any generator with "state" or "seed" causes the others to switch to the same style for future calls.

The state of each generator is independent and calls to different generators can be interleaved without affecting the final result.  For example,

     rand ("state", [11, 22, 33]);
     randn ("state", [44, 55, 66]);
     u = rand (100, 1);
     n = randn (100, 1);

and

     rand ("state", [11, 22, 33]);
     randn ("state", [44, 55, 66]);
     u = zeros (100, 1);
     n = zeros (100, 1);
     for i = 1:100
       u(i) = rand ();
       n(i) = randn ();
     end

produce equivalent results.  When the generators are initialized in the old style with "seed" only rand and randn are independent, because the old rande, randg and randp generators make calls to rand and randn.

The generators are initialized with random states at start-up, so that the sequences of random numbers are not the same each time you run Octave. (The old versions of rand and randn obtain their initial seeds from the system clock.) If you really do need to reproduce a sequence of numbers exactly, you can set the state or seed to a specific value.

If invoked without arguments, rand and randn return a single element of a random sequence.

The original rand and randn functions use Fortran code from RANLIB, a library of Fortran routines for random number generation, compiled by Barry W. Brown and James Lovato of the Department of Biomathematics at The University of Texas, M.D. Anderson Cancer Center, Houston, TX 77030.

 -- V =: randperm (N)
 -- V =: randperm (N, M)
     Return a row vector containing a random permutation of 1:N.

     If M is supplied, return M unique entries, sampled without replacement from 1:N.

     The complexity is O(N) in memory and O(M) in time, unless M < N/5, in which case O(M) memory is used as well.  The randomization is performed using rand().  All permutations are equally likely. perms

Famous Matrices
===============

The following functions return famous matrix forms.

 -- gallery: (NAME)
 -- gallery: (NAME, ARGS)
     Create interesting matrices for testing.

 -- C =: gallery ("cauchy", X)
 -- C =: gallery ("cauchy", X, Y)
     Create a Cauchy matrix.

 -- C =: gallery ("chebspec", N)
 -- C =: gallery ("chebspec", N, K)
     Create a Chebyshev spectral differentiation matrix.

 -- C =: gallery ("chebvand", P)
 -- C =: gallery ("chebvand", M, P)
     Create a Vandermonde-like matrix for the Chebyshev polynomials.

 -- A =: gallery ("chow", N)
 -- A =: gallery ("chow", N, ALPHA)
 -- A =: gallery ("chow", N, ALPHA, DELTA)
     Create a Chow matrix -- a singular Toeplitz lower Hessenberg matrix.

 -- C =: gallery ("circul", V)
     Create a circulant matrix.

 -- A =: gallery ("clement", N)
 -- A =: gallery ("clement", N, K)
     Create a tridiagonal matrix with zero diagonal entries.

 -- C =: gallery ("compar", A)
 -- C =: gallery ("compar", A, K)
     Create a comparison matrix.

 -- A =: gallery ("condex", N)
 -- A =: gallery ("condex", N, K)
 -- A =: gallery ("condex", N, K, THETA)
     Create a "counterexample" matrix to a condition estimator.

 -- A =: gallery ("cycol", [M N])
 -- A =: gallery ("cycol", N)
 -- A =: gallery (..., K)
     Create a matrix whose columns repeat cyclically.

 -- [C, D, E] =: gallery ("dorr", N)
 -- [C, D, E] =: gallery ("dorr", N, THETA)
 -- A =: gallery ("dorr", ...)
     Create a diagonally dominant, ill-conditioned, tridiagonal matrix.

 -- A =: gallery ("dramadah", N)
 -- A =: gallery ("dramadah", N, K)
     Create a (0, 1) matrix whose inverse has large integer entries.

 -- A =: gallery ("fiedler", C)
     Create a symmetric Fiedler matrix.

 -- A =: gallery ("forsythe", N)
 -- A =: gallery ("forsythe", N, ALPHA)
 -- A =: gallery ("forsythe", N, ALPHA, LAMBDA)
     Create a Forsythe matrix (a perturbed Jordan block).

 -- F =: gallery ("frank", N)
 -- F =: gallery ("frank", N, K)
     Create a Frank matrix (ill-conditioned eigenvalues).

 -- C =: gallery ("gcdmat", N)
     Create a greatest common divisor matrix.

     C is an N-by-N matrix whose values correspond to the greatest common divisor of its coordinate values, i.e., C(i,j) correspond gcd (i, j).

 -- A =: gallery ("gearmat", N)
 -- A =: gallery ("gearmat", N, I)
 -- A =: gallery ("gearmat", N, I, J)
     Create a Gear matrix.

 -- G =: gallery ("grcar", N)
 -- G =: gallery ("grcar", N, K)
     Create a Toeplitz matrix with sensitive eigenvalues.

 -- A =: gallery ("hanowa", N)
 -- A =: gallery ("hanowa", N, D)
     Create a matrix whose eigenvalues lie on a vertical line in the complex plane.

 -- V =: gallery ("house", X)
 -- [V, BETA] =: gallery ("house", X)
     Create a householder matrix.

 -- A =: gallery ("integerdata", IMAX, [M N ...], J)
 -- A =: gallery ("integerdata", IMAX, M, N, ..., J)
 -- A =: gallery ("integerdata", [IMIN, IMAX], [M N ...], J)
 -- A =: gallery ("integerdata", [IMIN, IMAX], M, N, ..., J)
 -- A =: gallery ("integerdata", ..., "CLASS")
     Create a matrix with random integers in the range [1, IMAX]. If IMIN is given then the integers are in the range [IMIN, IMAX].

     The second input is a matrix of dimensions describing the size of the output.  The dimensions can also be input as comma-separated arguments.

     The input J is an integer index in the range [0, 2^32-1].  The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index.

     The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "uint8", "uint16", "uint32", "int8", "int16", int32", "single", "double".  The default is "double".

 -- A =: gallery ("invhess", X)
 -- A =: gallery ("invhess", X, Y)
     Create the inverse of an upper Hessenberg matrix.

 -- A =: gallery ("invol", N)
     Create an involutory matrix.

 -- A =: gallery ("ipjfact", N)
 -- A =: gallery ("ipjfact", N, K)
     Create a Hankel matrix with factorial elements.

 -- A =: gallery ("jordbloc", N)
 -- A =: gallery ("jordbloc", N, LAMBDA)
     Create a Jordan block.

 -- U =: gallery ("kahan", N)
 -- U =: gallery ("kahan", N, THETA)
 -- U =: gallery ("kahan", N, THETA, PERT)
     Create a Kahan matrix (upper trapezoidal).

 -- A =: gallery ("kms", N)
 -- A =: gallery ("kms", N, RHO)
     Create a Kac-Murdock-Szego Toeplitz matrix.

 -- B =: gallery ("krylov", A)
 -- B =: gallery ("krylov", A, X)
 -- B =: gallery ("krylov", A, X, J)
     Create a Krylov matrix.

 -- A =: gallery ("lauchli", N)
 -- A =: gallery ("lauchli", N, MU)
     Create a Lauchli matrix (rectangular).

 -- A =: gallery ("lehmer", N)
     Create a Lehmer matrix (symmetric positive definite).

 -- T =: gallery ("lesp", N)
     Create a tridiagonal matrix with real, sensitive eigenvalues.

 -- A =: gallery ("lotkin", N)
     Create a Lotkin matrix.

 -- A =: gallery ("minij", N)
     Create a symmetric positive definite matrix MIN(i,j).

 -- A =: gallery ("moler", N)
 -- A =: gallery ("moler", N, ALPHA)
     Create a Moler matrix (symmetric positive definite).

 -- [A, T] =: gallery ("neumann", N)
     Create a singular matrix from the discrete Neumann problem (sparse).

 -- A =: gallery ("normaldata", [M N ...], J)
 -- A =: gallery ("normaldata", M, N, ..., J)
 -- A =: gallery ("normaldata", ..., "CLASS")
     Create a matrix with random samples from the standard normal distribution (mean = 0, std = 1).

     The first input is a matrix of dimensions describing the size of the output. The dimensions can also be input as comma-separated arguments.

     The input J is an integer index in the range [0, 2^32-1].  The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index.

     The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "single", "double". The default is "double".

 -- Q =: gallery ("orthog", N)
 -- Q =: gallery ("orthog", N, K)
     Create orthogonal and nearly orthogonal matrices.

 -- A =: gallery ("parter", N)
     Create a Parter matrix (a Toeplitz matrix with singular values near pi).

 -- P =: gallery ("pei", N)
 -- P =: gallery ("pei", N, ALPHA)
     Create a Pei matrix.

 -- A =: gallery ("poisson", N)
     Create a block tridiagonal matrix from Poisson's equation (sparse).

 -- A =: gallery ("prolate", N)
 -- A =: gallery ("prolate", N, W)
     Create a prolate matrix (symmetric, ill-conditioned Toeplitz matrix).

 -- H =: gallery ("randhess", X)
     Create a random, orthogonal upper Hessenberg matrix.

 -- A =: gallery ("rando", N)
 -- A =: gallery ("rando", N, K)
     Create a random matrix with elements -1, 0 or 1.

 -- A =: gallery ("randsvd", N)
 -- A =: gallery ("randsvd", N, KAPPA)
 -- A =: gallery ("randsvd", N, KAPPA, MODE)
 -- A =: gallery ("randsvd", N, KAPPA, MODE, KL)
 -- A =: gallery ("randsvd", N, KAPPA, MODE, KL, KU)
     Create a random matrix with pre-assigned singular values.

 -- A =: gallery ("redheff", N)
     Create a zero and ones matrix of Redheffer associated with the Riemann hypothesis.

 -- A =: gallery ("riemann", N)
     Create a matrix associated with the Riemann hypothesis.

 -- A =: gallery ("ris", N)
     Create a symmetric Hankel matrix.

 -- A =: gallery ("smoke", N)
 -- A =: gallery ("smoke", N, K)
     Create a complex matrix, with a "smoke ring" pseudospectrum.

 -- T =: gallery ("toeppd", N)
 -- T =: gallery ("toeppd", N, M)
 -- T =: gallery ("toeppd", N, M, W)
 -- T =: gallery ("toeppd", N, M, W, THETA)
     Create a symmetric positive definite Toeplitz matrix.

 -- P =: gallery ("toeppen", N)
 -- P =: gallery ("toeppen", N, A)
 -- P =: gallery ("toeppen", N, A, B)
 -- P =: gallery ("toeppen", N, A, B, C)
 -- P =: gallery ("toeppen", N, A, B, C, D)
 -- P =: gallery ("toeppen", N, A, B, C, D, E)
     Create a pentadiagonal Toeplitz matrix (sparse).

 -- A =: gallery ("tridiag", X, Y, Z)
 -- A =: gallery ("tridiag", N)
 -- A =: gallery ("tridiag", N, C, D, E)
     Create a tridiagonal matrix (sparse).

 -- T =: gallery ("triw", N)
 -- T =: gallery ("triw", N, ALPHA)
 -- T =: gallery ("triw", N, ALPHA, K)
     Create an upper triangular matrix discussed by Kahan, Golub, and Wilkinson.

 -- A =: gallery ("uniformdata", [M N ...], J)
 -- A =: gallery ("uniformdata", M, N, ..., J)
 -- A =: gallery ("uniformdata", ..., "CLASS")
     Create a matrix with random samples from the standard uniform distribution (range [0,1]).

     The first input is a matrix of dimensions describing the size of the output. The dimensions can also be input as comma-separated arguments.

     The input J is an integer index in the range [0, 2^32-1].  The values of the output matrix are always exactly the same (reproducibility) for a given size input and J index.

     The final optional argument determines the class of the resulting matrix. Possible values for CLASS: "single", "double". The default is "double".

 -- A =: gallery ("wathen", NX, NY)
 -- A =: gallery ("wathen", NX, NY, K)
     Create the Wathen matrix.

 -- [A, B] =: gallery ("wilk", N)
     Create various specific matrices devised/discussed by Wilkinson.

 -- H =: hadamard (N)
     Construct a Hadamard matrix (Hn) of size N-by-N.

     The size N must be of the form 2^k * p in which p is one of 1, 12, 20 or 28.  The returned matrix is normalized, meaning Hn(:,1) == 1 and Hn(1,:) == 1.

     Some of the properties of Hadamard matrices are:

       * kron (Hm, Hn) is a Hadamard matrix of size M-by-N.

       * Hn * Hn' = N * eye (N).

       * The rows of Hn are orthogonal.

       * det (A) <= abs (det (Hn)) for all A with abs (A(i, j)) <= 1.

       * Multiplying any row or column by -1 and the matrix will remain a Hadamard matrix.

     compan, hankel, toeplitz

 -- H =: hankel (C)
 -- H =: hankel (C, R)
     Return the Hankel matrix constructed from the first column C, and (optionally) the last row R.

     If the last element of C is not the same as the first element of R, the last element of C is used.  If the second argument is omitted, it is assumed to be a vector of zeros with the same size as C.

     A Hankel matrix formed from an m-vector C, and an n-vector R, has the elements

          H(i,j) = c(i+j-1),  i+j-1 <= m;
          H(i,j) = r(i+j-m),  otherwise

     hadamard, toeplitz

 -- H =: hilb (N)
     Return the Hilbert matrix of order N.

     The i,j element of a Hilbert matrix is defined as

          H(i, j) = 1 / (i + j - 1)

     Hilbert matrices are close to being singular which make them difficult to invert with numerical routines.  Comparing the condition number of a random matrix 5x5 matrix with that of a Hilbert matrix of order 5 reveals just how difficult the problem is.

          cond (rand (5))
             => 14.392
          cond (hilb (5))
             => 4.7661e+05

     invhilb

 -- HINV =: invhilb (N)
     Return the inverse of the Hilbert matrix of order N.

     This can be computed exactly using

                     (i+j)         /n+i-1\  /n+j-1\   /i+j-2\ 2
          A(i,j) = -1      (i+j-1)(       )(       ) (       )
                                   \ n-j /  \ n-i /   \ i-2 /

                 = p(i) p(j) / (i+j-1)

     where

                   k  /k+n-1\   /n\
          p(k) = -1  (       ) (   )
                      \ k-1 /   \k/

     The validity of this formula can easily be checked by expanding the binomial coefficients in both formulas as factorials.  It can be derived more directly via the theory of Cauchy matrices.  See J. W. Demmel, Applied Numerical Linear Algebra, p. 92.

     Compare this with the numerical calculation of inv (hilb (n)), which suffers from the ill-conditioning of the Hilbert matrix, and the finite precision of your computer's floating point arithmetic. hilb

 -- M =: magic (N)

     Create an N-by-N magic square.

     A magic square is an arrangement of the integers 1:n^2 such that the row sums, column sums, and diagonal sums are all equal to the same value.

     Note: N must be a scalar greater than or equal to 3.  If you supply N less than 3, magic returns either a nonmagic square, or else the degenerate magic squares 1 and [].

 -- P =: pascal (N)
 -- P =: pascal (N, T)
     Return the Pascal matrix of order N if T = 0.

     The default value of T is 0.

     When T = 1, return the pseudo-lower triangular Cholesky factor of the Pascal matrix (The sign of some columns may be negative).  This matrix is its own inverse, that is pascal (N, 1) ^ 2 == eye (N).

     If T = -1, return the true Cholesky factor with strictly positive values on the diagonal.

     If T = 2, return a transposed and permuted version of pascal (N, 1), which is the cube root of the identity matrix. That is, pascal (N, 2) ^ 3 == eye (N).

     chol

 -- R =: rosser ()
     Return the Rosser matrix.

     This is a difficult test case used to evaluate eigenvalue algorithms. wilkinson, eig

 -- T =: toeplitz (C)
 -- T =: toeplitz (C, R)
     Return the Toeplitz matrix constructed from the first column C, and optionally the first row R.

     If the second argument is omitted, the first row is taken to be the same as the first column.  If the first element of R is not the same as the first element of C, the first element of C is used.

     A Toeplitz, or diagonal-constant, matrix has the same value along each diagonal.  Although it need not be square, it often is.  An MxN Toeplitz matrix has the form:

          c(1)  r(2)   r(3)  ...  r(n)
          c(2)  c(1)   r(2)  ... r(n-1)
          c(3)  c(2)   c(1)  ... r(n-2)
           .     .      .   .      .
           .     .      .     .    .
           .     .      .       .  .
          c(m) c(m-1) c(m-2) ... c(m-n+1)

     hankel

 -- V =: vander (C)
 -- V =: vander (C, N)
     Return the Vandermonde matrix whose next to last column is C.

     If N is specified, it determines the number of columns; otherwise, N is taken to be equal to the length of C.

     A Vandermonde matrix has the form:

          c(1)^(n-1) ... c(1)^2  c(1)  1
          c(2)^(n-1) ... c(2)^2  c(2)  1
              .     .      .      .    .
              .       .    .      .    .
              .         .  .      .    .
          c(n)^(n-1) ... c(n)^2  c(n)  1

     polyfit

 -- W =: wilkinson (N)
     Return the Wilkinson matrix of order N.

     Wilkinson matrices are symmetric and tridiagonal with pairs of nearly, but not exactly, equal eigenvalues.  They are useful in testing the behavior and performance of eigenvalue solvers.

     rosser, eig