# LinearAlgebra¶

View LinearAlgebralib.chpl on GitHub

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)
```

```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 `LayoutCS`) 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;
```

```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 `LayoutCS` layout module.

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

```{
use LayoutCS;
const parentDom = {1..100, 1..100};
var csrDom: sparse subdomain(parentDom) dmapped CS();
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
}
```

### 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 LayoutCS;
var parentDom = {0..#3, 0..#3}; // Parent domain
var csrDom: sparse subdomain(parentDom) dmapped CS(); // 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

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

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

B = CSRMatrix(Bdom),
C = CSRMatrix(Cdom);

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

```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
```