LinearAlgebra ¶

Example usage of the LinearAlgebra module in Chapel.

use LinearAlgebra;

Compiling with Linear Algebra ¶

Some of the linear algebra module procedures rely on the BLAS and LAPACK modules. If using routines that rely on these modules, be sure to have a BLAS and LAPACK implementation available on your system. See the BLAS and LAPACK documentation for further details.

To explicitly opt out of using the BLAS and LAPACK procedures, compile your Chapel program with the flags:

chpl --set blasImpl=off --set lapackImpl=off myProgram.chpl

This will result in a cleaner compiler error when using a procedure that is only available with BLAS or LAPACK.

Basics ¶

Linear Algebra Types ¶

In Chapel, there are no Vector or Matrix types. Instead, vectors and matrices are represented as 1D and 2D native Chapel arrays, respectively. As a result, the terms vector and 1D array are used interchangeably, as well as matrix and 2D array.

The following example demonstrates this:

var vector1: [0..#5] int;             // Using Chapel array syntax
var vector2 = Vector(5, eltType=int); // Using LinearAlgebra Vector function
assert(vector1.type == vector2.type); // ... produces same result

var matrix1: [0..#3, 0..#5] int;         // Using Chapel array syntax
var matrix2 = Matrix(3, 5, eltType=int); // Using LinearAlgebra Matrix function
assert(matrix1.type == matrix2.type);    // ... produces same result

Array operations ¶

This section demonstrates some of the functionality built into Chapel arrays that can be useful in a linear algebra context. None of the operations in this section require usage of LinearAlgebra. Learn more about Chapel arrays in the Arrays chapter of the language specification and the Array Primer.

Create a 3 x 5 matrix of reals, default-initialized to values of 0.0

var a: [0..#3, 0..#5] real; // Equivalent to: var a = Matrix(3, 5);

Set all elements to 1.0

a = 1.0;

Print matrix

writeln(a);
  // Prints:
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0

Number of dimensions

a.rank; // 2

Size

a.size; // 15

Shape

a.shape; // (3, 5)

Element type

writeln(a.eltType: string); // real(64)

Array type in format of: [domain(rank, index-type, strides)] element-type

writeln(a.type: string); // [domain(2,int(64),one)] real(64)

Element-wise addition (and subtraction)

a = a + 1; // or, a += 1
writeln(a);
  // Prints:
  //   2.0 2.0 2.0 2.0 2.0
  //   2.0 2.0 2.0 2.0 2.0
  //   2.0 2.0 2.0 2.0 2.0

Element-wise multiplication

a = 2 * a;
writeln(a);
  // Prints:
  //   4.0 4.0 4.0 4.0 4.0
  //   4.0 4.0 4.0 4.0 4.0
  //   4.0 4.0 4.0 4.0 4.0

Element-wise division

a = a / 4.0;
writeln(a);
  // Prints:
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0

Array slicing

a[0, ..] = 0.0; // Sets first row to 0.0
writeln(a);
  // Prints:
  //   0.0 0.0 0.0 0.0 0.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0

a[.., 1] = 3.0; // Sets second column to 3.0
writeln(a);
  // Prints:
  //   0.0 3.0 0.0 0.0 0.0
  //   1.0 3.0 1.0 1.0 1.0
  //   1.0 3.0 1.0 1.0 1.0

Factory Functions ¶

The Linear Algebra module contains several factory functions for creating 1D (vector) and 2D (matrix) arrays. These functions are purely a convenience that provides a less verbose syntax for array declarations as well as an interface similar to other common linear algebra libraries. Below are some examples. See the documentation starting at LinearAlgebra.Vector for a comprehensive list of the available factory functions.

Vectors ¶

Shown below are a variety of way to create the same 5-element array of real type.

Note

The LinearAlgebra module represents row-vectors the same way as column-vectors. They are both 1D arrays.

Creating 1D arrays (vectors):

var v0: [0..#5] real;     // Chapel array syntax
writeln(v0);
  // Prints:
  // 0.0 0.0 0.0 0.0 0.0
var v1 = Vector(5);       // from a length
var v2 = Vector(0..#5);   // from a range
var v3 = Vector({0..#5}); // from a domain
var v4 = Vector([0, 0, 0, 0, 0], eltType=real); // from a 1D array
var v5 = Vector(0, 0, 0, 0, 0, eltType=real);   // from a list of values as arguments

If the eltType is omitted for the argument case, it will be inferred from the first argument. For example, this creates a vector of reals:

var v6 = Vector(0.0, 0, 0, 0, 0);
assert(v6.eltType == real);

Matrices ¶

Shown below are a variety of way to create the same 3 x 3 array of real type.

Note

The base LinearAlgebra module strictly supports DefaultRectangular arrays for matrix operations. Other layouts (such as CompressedSparseLayout) and distributions (such as BlockDist) will be supported through LinearAlgebra submodules.

Creating 2D arrays (matrices):

var M0: [0..#3, 0..#3] real; // Chapel array syntax
writeln(M0);
  // Prints:
  //   0.0 0.0 0.0
  //   0.0 0.0 0.0
  //   0.0 0.0 0.0
var M1 = Matrix(3);                 // from a length
var M2 = Matrix(3, 3);              // from a length and width
var M3 = Matrix(0..#3);             // from 1 range
var M4 = Matrix(0..#3, 0..#3);      // from 2 ranges
var M5 = Matrix({0..#3, 0..#3});    // from a 2D domain
var M6 = Matrix(M5);                // from a 2D array
var M7 = Matrix([0,0,0],[0,0,0],[0,0,0], eltType=real); // from a list of 1D arrays as arguments

Creating identity matrices:

var I1 = eye(3); // From a length
writeln(I1);
  // Prints:
  //   1.0 0.0 0.0
  //   0.0 1.0 0.0
  //   0.0 0.0 1.0
var I2 = eye(3,3); // From a length and width
var I3 = eye({0..#3, 0..#3}); // from a 2D domain

Creating diagonal matrix from a vector:

var vec = Vector(1, 2, 3, eltType=real);
var diagMatrix = diag(vec);
writeln(diagMatrix);
  // Prints:
  //   1.0 0.0 0.0
  //   0.0 2.0 0.0
  //   0.0 0.0 3.0

Operations ¶

Setup some matrices

{ // Operations scope

var A = Matrix(3,5),
    B = Matrix(3,5);

A = 1.0;
B = 2.0;

Element-wise addition

var ApB = A + B;

Element-wise subtraction

var AmB = A - B;

Element-wise multiplication

var AtB = A * B;

Element-wise division

var AdB = A / B;

Taking the transpose of a matrix:

var M0T = transpose(M0);
M0T = M0.T; // short-hand syntax

Taking the transpose of a vector simply returns the vector because row and column vectors are equivalent:

var t1 = v0.T;
assert(t1.type == v0.type);
assert(t1 == v0);

Creating a matrix and some vectors to use:

var X = Matrix(3, 5),
    y = Vector(3),
    z = Vector(5);

X = 1;
X -= eye(3,5);
writeln(X);
  // Prints:
  //   0.0 1.0 1.0 1.0 1.0
  //   1.0 0.0 1.0 1.0 1.0
  //   1.0 1.0 0.0 1.0 1.0
y = 2;
z = 1;

Matrix-matrix multiplication is computed by dot():

var MM = dot(X, X.T);

… as well as matrix-vector multiplication:

var Mv = dot(X, z);

… or vector-matrix multiplication:

var vM = dot(y, X);

… and vector-vector multiplication:

var vv = dot(y, y);

Because row and column vectors are indistinguishable, an dot() always computes the inner product of 2 vectors.

To compute an outer product, one must explicitly call outer():

var yz = outer(y, z);
writeln('outer:', yz);

Matrix powers can be computed on square matrices with matPow(). Only integer powers are currently supported:

var MM4 = matPow(MM, 4);

Cross products of 3-element vectors can be computed with cross():

var crossProduct = cross(y, y);

} // Operations scope

Properties ¶

There are several functions for computing matrix properties.

Trace of a matrix

var tr = trace(diagMatrix); // 6.0

Norm of a vector

var N1 = norm(vec); // 3.74166

Norm of a matrix

var N2 = norm(diagMatrix); // 3.74166

Structure ¶

There are several functions for inspecting and extracting matrix structure.

Below are a few examples. For the comprehensive list of available structural functions, see the documentation starting at LinearAlgebra.diag.

Extract diagonal of matrix

var diagVec = diag(diagMatrix); // 1.0 2.0 3.0

Confirm that a matrix is diagonal

writeln(isDiag(diagMatrix)); // true

Extract upper triangle of matrix

var onesMatrix = Matrix(5,5);
onesMatrix = 1.0;
writeln(onesMatrix);
  // Prints:
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
  //   1.0 1.0 1.0 1.0 1.0
var upper = triu(onesMatrix, k=0); // k=0 includes diagonal
writeln(upper);
  // Prints:
  //   1.0 1.0 1.0 1.0 1.0
  //   0.0 1.0 1.0 1.0 1.0
  //   0.0 0.0 1.0 1.0 1.0
  //   0.0 0.0 0.0 1.0 1.0
  //   0.0 0.0 0.0 0.0 1.0

Confirm that a matrix is upper triangular

writeln(isTriu(upper));         // true
writeln(isTriu(upper, k=1));    // false (k=1 does not include diagonal)

LinearAlgebra.Sparse ¶

The Sparse submodule supports operations on sparse matrices.

When using the LinearAlgebra module, this submodule can be accessed like so:

{
  var A = Sparse.CSRDomain(100, 100);
}

Alternatively, one can expose all the submodule symbols with an additional use statement:

use LinearAlgebra.Sparse;

Now we can do:

{
  var A = CSRDomain(100, 100);
}

Supported Sparse Layouts ¶

Currently only CSR matrices are supported. A CSR matrix is composed of a 2D array that is domain-mapped to a CS() layout, from the CompressedSparseLayout layout module.

For example, the following CSR matrix (CS(compressRows=true)) is supported by this submodule:

{
  use CompressedSparseLayout;
  const parentDom = {1..100, 1..100};
  var csrDom: sparse subdomain(parentDom) dmapped new csrLayout();
  var csrMatrix: [csrDom] real; // Supported by LinearAlgebra.Sparse
}

And as a counter-example, the following COO matrix (DefaultSparse) is not supported by this module:

{
  const parentDom = {1..100, 1..100};
  var cooDom: sparse subdomain(parentDom);
  var cooMatrix: [cooDom] real; // NOT supported by LinearAlgebra.Sparse
}

For more information about working with sparse arrays in Chapel, see the Sparse Primer.

Factory Functions ¶

The Sparse submodule provides factory functions for both sparse domains and sparse matrices.

In Chapel, sparse domains can only be modified directly, so it is important to maintain access to both the sparse array and sparse domain.

Consider the following example:

{
  // Bad example of LinearAlgebra.Sparse usage:
  var A = CSRMatrix({1..100, 1..100}); // Create a CSR matrix from dense domain
  writeln(A.domain); // empty sparse domain, with no way to add elements..
}

The above sparse matrix is initialized with an empty domain, and the user has no way of accessing that domain. Therefore, it is not possible to add indices to the sparse domain.

Instead, users should interact with both the sparse domain and sparse array:

{
  // Good example of LinearAlgebra.Sparse usage:
  var D = CSRDomain(100, 100); // empty sparse domain
  var A = CSRMatrix(D); // Create a CSR matrix from a sparse domain
  D += [(10,20), (30, 24), (64, 42)]; // Add some indices to the sparse domain
}

Below are a list of the available factory functions:

{ // Factory functions scope

Creating CSR domains & arrays in Chapel without LinearAlgebra.Sparse:

use CompressedSparseLayout;
var parentDom = {0..#3, 0..#3}; // Parent domain
var csrDom: sparse subdomain(parentDom) dmapped new csrLayout(); // CSR domain
var csrArr: [csrDom] real; // CSR array

Adding indices to a sparse domain:

csrDom += (0,0);
csrDom += [(1,1), (2,2)];

Creating CSR domains with LinearAlgebra.Sparse:

var D1 = CSRDomain(3);                 // from a length
var D2 = CSRDomain(3, 3);              // from a length and width
var D3 = CSRDomain(0..#3);             // from 1 range
var D4 = CSRDomain(0..#3, 0..#3);      // from 2 ranges
var D5 = CSRDomain({0..#3, 0..#3});    // from a 2D domain

Creating CSR matrices:

var M1 = CSRMatrix(D1);             // from a CSR domain
var M2 = CSRMatrix({0..#3, 0..#3}); // from a dense domain
var M3 = CSRMatrix(M1);             // From a CSR matrix
const I = eye(3,3);
var M4 = CSRMatrix(I);              // From a dense matrix
                                    // Indices holding non-zero elements are
                                    // added to M4's sparse domain


} // Factory functions scope

Operations ¶

Below is a list of the currently supported operations.

Note

For sparse matrices, the element-wise operations such as +, -, *, and / are promoted, and will result in an error if there is a sparse domain mismatch. To avoid this, it is preferred to use the element-wise methods: .plus(), .minus(), .times(), and .elementDiv().

Setup some sparse domains and arrays

{ // Operations scope

var Adom = CSRDomain(1..100, 1..100),
    Bdom = CSRDomain(1..100, 1..100),
    Cdom = CSRDomain(1..100, 1..100);

Adom += (1,1);
Bdom += (2,2);
Cdom += [(1,1), (2,2)];

var A = CSRMatrix(Adom),
    B = CSRMatrix(Bdom),
    C = CSRMatrix(Cdom);

A = 1.0;
B = 2.0;
C = 4.0;

Element-wise addition.

A.plus(B); // 1.0, 2.0
A.plus(C); // 5.0, 4.0

Element-wise subtraction

A.minus(B); // 1.0, -2.0

Element-wise multiplication

A.times(B); // 0.0, 0.0
A.times(C); // 4.0, 4.0

Element-wise division

A.elementDiv(C); // 0.25, 0.0
B.elementDiv(C); // 0.0, 0.5

Matrix-matrix multiplication

C.dot(C); // 16.0, 16.0

Matrix-vector multiplication

var a = Vector(100);
a = 2.0;
C.dot(a);

Matrix-scalar multiplication

C.dot(2.0); // equivalent to C.times(2.0);

Matrix transpose

transpose(C);
C.T; // short-hand syntax

} // Operations scope