# LAPACK¶

Example usage of the LAPACK module in Chapel.

This particular file demonstrates how to use the gesv procedure with Chapel’s
arrays. The gesv procedure solves `A*X = B`

for `X`

given both `A`

and
`B`

, where `A`

is a square matrix.

Here, we create `A`

and `X`

, then compute `B`

with a matrix
multiplication, and show that the result of the gesv procedure is the same as
our synthesized `X`

array (within a small margin of error for floating
point).

To compile a program with LAPACK, you may need to add some additional flags depending on how LAPACK was installed on your system:

```
chpl -I$PATH_TO_LAPACK_INCLUDE_DIR \
-L$PATH_TO_LAPACK_BINARIES \
-lgfortran
```

See the `LAPACK`

module documentation for more information on compiling.

Start by using the `LAPACK`

module to gain access to the
`gesv`

function. The `Random`

module will be used to fill
arrays with random values.

```
use LAPACK;
use Random;
```

- Here we set up several config consts that represent, in order:
`N`

, the dimension for the square array`A`

`K`

, the second dimension for arrays`X`

and`B`

`epsilon`

, the margin of error for success`seed`

, a seed for random number generation

```
config const N = 2;
config const K = 1;
config const epsilon = 1e-13;
config const seed = 41;
```

Create the arrays `A`

, `X`

, and `B`

. Fill `A`

and `X`

with random
values.

```
var A : [1..N, 1..N] real;
fillRandom(A, seed);
var X : [1..N, 1..K] real;
fillRandom(X, seed);
var B : [1..N, 1..K] real;
```

Matrix multiply `A*X`

, store result in `B`

```
for i in 1..N do
for j in 1..K do
for k in 1..N do
B[i,j] += A[i,k] * X[k,j];
writeln("Matrix A:\n", A, "\n");
writeln("Matrix X:\n", X, "\n");
writeln("Matrix B:\n", B, "\n");
```

Copy original arrays into temporary arrays.

Input and work array. Becomes garbage.

```
var WorkA = A;
```

Input and output. Becomes result of solution (`X`

)

```
var WorkBX = B;
```

Output array. Stores pivot indices

```
var ipiv : [1..N] c_int;
```

Call the `gesv`

function to solve for `X`

. Note that Chapel arrays are
row-major order by default.

```
var info = gesv(lapack_memory_order.row_major, WorkA, ipiv, WorkBX);
```

LAPACK returns an error code to indicate a failure.

```
if info != 0 {
writeln("There was an error!");
if info < 0 {
writeln("Argument ", -info, " was incorrect.");
} else {
writeln("The matrix A is a singular matrix. U", (info,info), " is zero");
}
}
writeln("gesv result for X:\n", WorkBX, "\n");
```

The arrays may not be identical due to floating point errors. Use a small value as a margin of error to measure success.

```
const closeEnough = && reduce [d in (WorkBX - X)] abs(d) < epsilon;
if closeEnough then
writeln("SUCCESS");
else
writeln("FAILURE");
```