forgeCol and forgeMat¶

forgeCol<T> and forgeMat<T> are FORGE's vector and matrix types. They look like Armadillo's Col and Mat, and they are what the solver and the operator inner loops use so that the same code can run on the CPU, through Apple Accelerate and Metal, or on CUDA.

Use float, double, forgeComplex<float> or forgeComplex<double> as the element type. forgeComplex<T> has the same memory layout as std::complex<T>.

Layout and sizing¶

Data is column-major: element (r, c) of a forgeMat is at index r + n_rows * c, the same as Armadillo, LAPACK and Fortran.

The forgeMat(nRows, nCols) constructor takes rows first, but forgeMat::set_size(nCols, nRows) takes columns first. This is the opposite of most matrix libraries, so check the order whenever you call set_size directly.

Converting to and from Armadillo¶

arma::Col<std::complex<float>> a = ...;
forgeCol<forgeComplex<float>> v(a);   // copies

arma::Col<std::complex<float>> b = v.getArma();   // copies

The const overload of getArma() returns an Armadillo view of the same memory instead of a copy. getArmaComplex() also returns a temporary view; don't keep it beyond the expression it's used in.

Operators that still take Armadillo arguments, such as Gnufft and Robject, are called from forgeCol code by converting at the call site:

Col<CxT1> result = (*G_obj) * weighted.getArma();
outData.set_col(ii, forgeCol<forgeComplex<T1>>(result));

Views¶

A view is a forgeCol that points at memory it doesn't own. forgeMat::col() and forgeCol::subvec() return views, so taking a column costs nothing and writing to it changes the matrix:

forgeCol<float> c = mat.col(2);         // view of column 2
c.is_view();                            // true
mat.col(1) %= scale;                    // scales column 1 of mat in place

forgeCol<float> s = vec.subvec(10, 19); // elements 10 to 19, inclusive
forgeCol<float> w = forgeCol<float>::view(ptr, length);  // wrap your own buffer

Things to keep in mind:

A view must not outlive the memory it points at.
Copy-constructing from a view gives you an owning deep copy.
Destroying a view doesn't free anything.
Calling set_size() on a view detaches it: it allocates new memory and stops being a view. The original memory is left alone.

When you need a column that lives longer than the matrix, use mat.col_copy(ii). To write a vector into a column, use mat.set_col(ii, src).

Functions¶

Element-wise arithmetic works as in Armadillo: +, -, % (element-wise multiply), /, the compound forms +=, -=, %=, /=, and scalar versions. A real vector times a complex vector (W % X) gives a complex result, and s / A divides a scalar by each element.

Free functions on forgeCol:

Function	Result
`sum(A)`, `accu(A)`	Sum of elements
`cdot(A, B)`	Complex dot product
`norm(A)`	L2 norm
`accu_schur(A, B)`	Sum of element-wise products, without conjugation
`norm_grad(g, yi, W)`	Scale-invariant gradient norm used by the PCG solver
`conj`, `abs`, `real`, `imag`	Element-wise
`sqrt`, `square`, `log1p`	Element-wise
`to_complex(A)`	Real vector to complex, with zero imaginary part

Free functions on forgeMat:

Function	Result
`sum(M, dim)`	Sum down columns (`dim = 0`) or across rows (`dim = 1`)
`vectorise(M)`	Flatten to a `forgeCol`
`conj(M)`	Element-wise conjugate
`matmul(A, B)`	Matrix-matrix or matrix-vector product (via Armadillo)
`trans(A)`, `strans(A)`	Conjugate transpose and plain transpose
`reshape(A, rows, cols)`	Reshape, keeping column-major order
`cols(A, first, last)`	Deep copy of columns `first` to `last`, inclusive

Where the work runs¶

On macOS, float and forgeComplex<float> arithmetic goes to Apple Accelerate (vDSP) for vectors of 32 elements or more. In a Metal build, complex multiply, divide, scalar multiply, and real-times-complex multiply try the Metal GPU first for vectors of at least 4096 elements. You can change that threshold with the PG_METAL_MIN_SIZE environment variable (minimum 256). double types always use the CPU loops.

In a CUDA build, a forgeCol can be moved to the GPU with putOnGPU() and back with getFromGPU(). Operations on vectors that are already on the GPU run there.

In a Metal build, forgeCol and forgeMat allocate with std::aligned_alloc(16384, ...), matching the Apple Silicon page size. Page-aligned buffers let Metal use the memory directly instead of copying it.

Checking for NaN¶

has_nan() inspects the bit pattern, so it still works when the code is compiled with -ffast-math:

if (v.has_nan()) {
    // handle the error
}