Fundamental Operations

Operations in LAPACK are expressed to the greatest extent possible in terms of BLAS operations. Expressing the linear algebra operations in terms of core operations reduces errors and allows optimizations implemented in the core operations to be propogated to the higher-level operations. The fundamental linear algebra operations summarized in the following table are based on a review of the BLAS Technical Forum documentation. These operations provide the fundamental building blocks necessary for higher level linear algebra operations.

Operation	Description
sumsq	Return the scaling parameter and sum of the squares.
sump	Return the scaling parameter and sum of the P powers.
norm	Vector and matrix norms.
transpose	Transpose the vector or matrix.
ntranspose	Destructive transpose of the vector or matrix.
scale	Scale each element of the vector or matrix.
nscale	Destructively scale each element of the vector or matrix.
add	Vector or matrix binary addition.
nadd	Destructive vector or matrix binary addition.
substract	Vector or matrix binary subtraction.
nsubstract	Destructive vector or matrix binary subtraction.
product	Vector-vector, matrix-vector, or matrix-matrix product.

Scaled Vector and Matrix Sums

(defgeneric sumsq (vector-or-matrix &key scale sumsq)
  (:documentation
   "Return the scaling parameter and the sum of the squares."))

Scaled vector and matrix sums are primarily required to support vector and matrix norms. sumsq returns 2 values, the scaling parameter and the sum of the squares. The function iterates through the elements of the vector or matrix updating scale and sumsq so that they satisfy the following equations.

[ (scale, sumsq) \leftarrow \sum_i x_i^2 ]

[ (scale, sumsq) \leftarrow \sum_i \sum_j A_{ij}^2 ]

such that

[ scale^2 \times sumsq = \sum_i x_i^2 ]

[ scale^2 \times sumsq = \sum_i \sum_j A_{ij}^2 ]

where

[ scale \leftarrow \max ( scale, |x|_i ) ]

for vectors and

[ scale \leftarrow \max ( scale, |A|_{ij} ) ]

for matrices.

(defgeneric sump (vector-or-matrix p &key scale sump)
  (:documentation
   "Return the scaling parameter and the sum of the P powers."))

Similarly, sump returns 2 values, the scaling parameter and the sum of the P powers of the absolute values of the vector or matrix. The function iterates through the vector or matrix updating the scaling parameter and the sum so that they satisfy the following equation.

[ (scale, sump) \leftarrow \sum_i x_i^p ]

[ (scale, sump) \leftarrow \sum_i \sum_j A_{ij}^p ]

such that

[ scale^p \times sump = \sum_i x_i^p ]

[ scale^p \times sump = \sum_i \sum_j A_{ij}^p ]

where

[ scale \leftarrow \max ( scale, |x|_i ) ]

for vectors and

[ scale \leftarrow \max ( scale, |A|_{ij} ) ] for matrices.

Note: The keyword arguments for SUMSQ and SUMP are present to mirror the argument list of LASSQ and friends subroutines in LAPACK. It appears that these are present to facilitate using LASSQ with arrays. If this can be verified, these keywords should be removed from the lambda list.

Vector and Matrix Norms

(defgeneric norm (vector-or-matrix &key measure)
  (:documentation
   "Return the norm according to measure."))

The norm function returns the norm of the vector or matrix according to the measure. There is one required argument, either a vector or a matrix, and one keyword argument, the measure of the norm.

[ \mbox{result} \leftarrow || x ||_p ]

[ \mbox{result} \leftarrow || A ||_p ]

Three specific vector norms are implemented, listed in the following table. The default norm is the taxicab, or one, norm. It is the sum of the absolute values of the vector elements. The next norm is the Euclidean, or two, norm. It is the square root of the sum of the squares of the vector elements. The last specific norm is the infinity norm. It is simply the maximum of the absolute values of the vector elements. All other measures are handled by the general P-norm.

Vector Norm
Taxicab	\( ||x||_1 = \sum_i |x|_i\)
Euclidean	\(||x||_2 = \Bigl( \sum_i |x_i|^2 \Bigr)^{1/2}\)
P-norm	\(||x||_p = \Bigl( \sum_i |x_i|^p \Bigr)^{1/p}\)
Infinity	\(||x||_\infty = \max_i |x_i|\)

Four specific matrix norms are implemented, listed in the following table. The default norm is the 1 norm. It is the maximum of the sums of the matrix columns. The second norm is the max norm that is the maximum of the absolute values of the elements. Next is the Frobenius norm that is the matrix equivalent of the vector Euclidean norm. Finally, the infinity norm of a matrix is the maximum of the sums of the matrix rows.

Matrix Norm
1-norm	\(||A||_1 = \max_j \sum_i |a_{ij}|\)
Max	\(||A||_{max} = \max\Bigl\{|a_{ij}|\Bigr\}\)
Frobenius	\(||A||_F = \Bigl( \sum_{ij} |a_{ij}|^2 \Bigr)^{1/2}\)
Infinity	\(||A||_\infty = \max_i \sum_j |a_{ij}|\)

Vector and Matrix Transpose

(defgeneric transpose (vector-or-matrix &key conjugate)
  (:documentation
   "Transpose the vector or matrix."))

The transpose function returns the transpose or conjugate transpose of a vector or matrix. It has one required argument that is either the vector or matrix. There is a conjugate keyword that is only applicable to complex valued vectors or matrices. Setting conjugate to true returns the conjugate transpose. The conjugate keyword has no effect on real valued vectors or matrices.A new vector or matrix is returned from transpose.

[ \mbox{result} \leftarrow x^T ]

[ \mbox{result} \leftarrow x^* ]

[ \mbox{result} \leftarrow A^T ]

[ \mbox{result} \leftarrow A^* ]

(defgeneric ntranspose (vector-or-matrix &key conjugate)
  (:documentation
   "Destructively transpose the vector or matrix."))

The destructive transpose, ntranspose, updates the vector or matrix.

[ x \leftarrow x^T ]

[ x \leftarrow x^* ]

[ A \leftarrow A^T ]

[ A \leftarrow A^* ]

Vector and Matrix Permutation

(defgeneric permute (vector-or-matrix vector-or-matrix)
  (:documentation
   "Permute the vector or matrix."))

The permute function applies a permutation matrix to a vector or matrix and accepts 2 arguments. The permutation can be applied to either the rows or the columns depending on the order of the arguments. If the first argument is the permutation matrix and the second argument the vector or matrix, the rows are permuted. If the first argument is the vector or matrix and the second argument is the permutation matrix, the columns are permuted. A new vector or matrix is returned from this function.

[ \mbox{result} \leftarrow Px ]

[ \mbox{result} \leftarrow xP ]

[ \mbox{result} \leftarrow PA ]

[ \mbox{result} \leftarrow AP ]

(defgeneric npermute (vector-or-matrix vector-or-matrix)
  (:documentation
   "Destructively permute the vector or matrix."))

npermute destructively permutes the vector or matrix.

[ x \leftarrow Px ]

[ x \leftarrow xP ]

[ A \leftarrow PA ]

[ A \leftarrow AP ]

Vector and Matrix Scaling

(defgeneric scale (scalar vector-or-matrix)
  (:documentation
   "Scale each element by the scalar."))

The scale function scales each element of the vector or matrix by a scalar. Two arguments are required, the scalar multiple and the vector or matrix to be scaled. A new vector or matrix is returned from this function.

[ \mbox{result} \leftarrow scalar \times x ]

[ \mbox{result} \leftarrow scalar \times A ]

(defgeneric nscale (scalar vector-or-matrix)
  (:documentation
   "Destructively scale each element by the scalar."))

The destructive version, nscale, updates the vector or matrix.

[ x \leftarrow scalar \times x ]

[ A \leftarrow scalar \times A ]

Vector and Matrix Addition

(defgeneric add (vector-or-matrix-1 vector-or-matrix-2
                 &key scalar1 scalar2)
  (:documentation
   "Vector or matrix binary addition."))

The add function returns the addition of 2 vectors or matrices, optionally scaled. Two arguments are required, the 2 vectors or matrices to perform addition on. The scalar values are specified using keywords. A new vector or matrix is returned from add.

[ \mbox{result} \leftarrow scalar1 \times x_1 + scalar2 \times x_2 ]

[ \mbox{result} \leftarrow scalar1 \times A_1 + scalar2 \times A_2 ]

(defgeneric nadd (vector-or-matrix-1 vector-or-matrix-2
                  &key scalar1 scalar2)
  (:documentation
   "Destructive vector or matrix addition."))

The destructive version, nadd, updates the first vector or matrix.

[ x_1 \leftarrow scalar1 \times x_1 + scalar2 \times x_2 ]

[ A_1 \leftarrow scalar1 \times A_1 + scalar2 \times A_2 ]

Vector and Matrix Subtraction

(defgeneric subtract (vector-or-matrix-1 vector-or-matrix-2
                      &key scalar1 scalar2)
  (:documentation
   "Vector or matrix binary subtraction."))

The subtract function returns the subtraction of a vector or matrix from a vector or matrix, respectively. Two arguments are required, the 2 vectors or matrices to operate on. The vectors or matrices are optionally scaled by scalars specified with keywords. The subtraction is performed by subtracting the elements of the second argument from the elements of the first argument. A new vector or matrix is returned from subtract.

[ \mbox{result} \leftarrow scalar1 \times x_1 - scalar2 \times x_2 ]

[ \mbox{result} \leftarrow scalar1 \times A_1 - scalar2 \times A_2 ]

(defgeneric nsubtract (vector-or-matrix-1 vector-or-matrix2
                       &key scalar1 scalar2)
  (:documentation
   "Destructive vector or matrix subtraction."))

The destructive version, nsubtract, updates the first vector or matrix.

[ x_1 \leftarrow scalar1 \times x_1 - scalar2 \times x_2 ]

[ A_1 \leftarrow scalar1 \times A_1 - scalar2 \times A_2 ]

Note: The subtraction function is in a sense redundant to the addition function. This could be resolved in a couple ways. The first way would be to only provide the addition function and perform subtraction by passing a negative value for the second scalar. The second way would be to abstract binary combination out to a general function that is wrapped by addition and subtraction.

Vector and Matrix Product

(defgeneric product (vector-or-matrix-1 vector-or-matrix-2 &key scalar)
  (:documentation
   "Return the vector-vector, matrix-vector or matrix-matrix product."))

product is the interface to the vector dot product, matrix-vector product and matrix-matrix product operations. There are 2 required arguments and 2 keyword arguments, scalar and conjugate. Matrix-vector product implementations require a matrix and a vector argument. If the matrix is pre-multiplied by the vector, the first argument must be a row vector and the second argument a matrix. Otherwise, the first argument is a matrix and the second argument is a column vector. A new vector is returned from matrix-vector products. A new matrix is returned from matrix-matrix products.

[ \mbox{result} \leftarrow scalar \times x^T y = scalar \times \sum_i x_i y_i ]

[ \mbox{result} \leftarrow scalar \times x^* y = scalar \times \sum_i \bar{a}_i y_i ]

[ \mbox{result} \leftarrow scalar \times x^T A = scalar \times \sum_i x_i A_{ij} ]

[ \mbox{result} \leftarrow scalar \times x^* A = scalar \times \sum_i \bar{x}i A{ij} ]

[ \mbox{result} \leftarrow scalar \times A x = scalar \times \sum_j A_{ij} x_j ]

[ \mbox{result} \leftarrow scalar \times A B = scalar \times \sum_k A_{ik} B_{kj} ]

Note: Need to add the conjugate keyword to the lambda list.

Note: The BLAS xxMV routines and GEMM routine have optional transpose arguments. In the case of matrix-vector products, it is presumed that the transpose argument indicates whether it is a vector-matrix or matrix-vector product. In this library, this is handled using generic dispatch. For example, it is valid to only pre-multiply a matrix by a row vector. In the case of the GEMM routine, the transpose arguments account for matrix-matrix product with non-square matrices. This is automatically handled in this library. Consequently, there is no need for a transpose keyword that is analogous to BLAS optional argument. The presumption for this will be verified once the entire library is implemented. If this presumption is wrong, a transpose keyword will be added.