# Example

In the code example on the right, you can see how easy and straightforward the use of GaspiLS can be.

```// Solve Poisson's equation in 3D with finite differences.

#include <GASPI.h>
#include <GASPI_Ext.h>
#include <GaspiLS/GaspiLS.hpp>
#include <GaspiLS/Solvers.hpp>
#include <cmath>
#include <cstdlib>
#include <cstdio>
#include <ctime>

// The problem settings
int    n1d;
int    n;
double h;
double h2inv;

// Matrix and vector classes used in this example
typedef GaspiLS::Vector<double, int, int>    vec;
typedef GaspiLS::CSRMatrix<double, int, int> mat;

// The right-hand side of Poisson's equation.
double func_f(double x, double y, double z)
{
// Approximate a Dirac Delta.
x -= 0.5;
y -= 0.5;
z -= 0.5;
if(sqrt(x*x + y*y + z*z) < h) return 1.0/(h*h*h);
return 0.0;
}

// Set up the right-hand side of the linear system.
void setRHS(vec &f)
{
double x, y, z;
double *fEntries = f.localEntryPointer();
int    k         = f.distribution().minGlobalIndex();

for(int i = 0; i < f.localSize(); i++, k++) {
x = (k%n1d + 1)*h;
y = (k%(n1d*n1d)/n1d + 1)*h;
z = (k/(n1d*n1d) + 1)*h;
fEntries[i] = -func_f(x, y, z);
}
}

// Discretize the Laplace operator using finite differences.
void setOperator(mat &A)
{
for(int k  = A.rowDistribution().minGlobalIndex();
k <= A.rowDistribution().maxGlobalIndex();
k++) {

A.insertEntry(k, k, 6.0*h2inv);

if(k%n1d > 0) {
A.insertEntry(k, k - 1, -1.0*h2inv);
}
if(k%n1d < n1d - 1) {
A.insertEntry(k, k + 1, -1.0*h2inv);
}
if(k%(n1d*n1d) >= n1d) {
A.insertEntry(k, k - n1d, -1.0*h2inv);
}
if(k%(n1d*n1d) < n1d*n1d - n1d) {
A.insertEntry(k, k + n1d, -1.0*h2inv);
}
if(k >= n1d*n1d) {
A.insertEntry(k, k - n1d*n1d, -1.0*h2inv);
}
if(k < n1d*n1d*n1d - n1d*n1d) {
A.insertEntry(k, k + n1d*n1d, -1.0*h2inv);
}
}
A.finalize();
}

int main(int argc, char *argv[])
{
if(argc < 5 || argc > 6) {
printf("Usage: poisson3d N_1D rtol maxit N_THREADS [OUTPUT_FILE]\n");
return -1;
}

n1d   = atoi(argv[1]);
n     = n1d*n1d*n1d;
h     = 1.0/(n1d + 1);
h2inv = 1.0/(h*h);

const size_t segmentSize = static_cast<size_t>(1) << 30;
// Create a GASPI interface for the group GASPI_GROUP_ALL
// with the segment 0 and queue 0.
gaspi_proc_init(GASPI_BLOCK);
gaspi_segment_create(0,
segmentSize,
GASPI_GROUP_ALL,
GASPI_BLOCK,
GASPI_MEM_INITIALIZED);

// Next, create a thread pool with the number of threads given by
// argv[4]. Thread pools are used by many GaspiLS objects to per-
// form parallel computations. They are independent from any user-
// generated threads, and the user might even use a different
// documentation for thread pinning options.
// Create a GASPI interface for the group GASPI_GROUP_ALL
// using segment 0 and queue 0.
GaspiLS::GASPIInterface         interface(GASPI_GROUP_ALL, 0, 0);
// Now create a uniform distribution of n elements on the GASPI
// interface. Distributions have two template parameters which
// specify the local and global index type respectively. Using int
// here will improve performance in later computations, but re-
// stricts the distributon to 2^31 - 1 elements.
// If the distribution has more elements, long must be used for
// the global index type, while the local index type should still
// be int as long as possible.
GaspiLS::Distribution<int, int> distribution(interface, n);
// Vectors also have three template parameters. The first para-
// meter specifies the scalar type, such as double, or
// std::complex<float>. The second and third parameter again
// specify the local and global index type. Each vector requires
// a distribution to determine which vector entries are stored
// on which rank, as well as a thread pool for parallel vector
// computations.
vec                             u(distribution, pool);
vec                             f(distribution, pool);
// The parameters for creating a matrix are very similar to those
// of the vector class. The main difference is that matrices have
// two distributions. The first one is called the row distribution
// and it specifies which matrix rows are owned by which rank.
// The second one is called the domain distribution and it
// corresponds to the distribution of the vector x during the
// matrix-vector product y = Ax. It also determines the number
// of columns for the matrix A.
mat                             A(distribution, distribution, pool);
GaspiLS::CGSolver<mat, vec>     cg;

// Discretize the PDE.
setOperator(A);
setRHS(f);

// Set up linear solver.
cg.setProblem(A, u, f);
cg.absoluteTolerance(1.0e-20);
cg.relativeTolerance(atof(argv[2]));
cg.divergenceTolerance(1.0e9);
cg.maxIterations(atoi(argv[3]));
cg.zeroInitialGuess(true);

// Solve linear system.
cg.solve();

// Write solution to disk.
if(argc == 4) {
u.exportToFile(argv[3]);
}

if(interface.groupRank() == 0) {
printf("FINISHED (%d iterations, rel. residual: %e, "
"abs. residual: %e\n",
cg.iterations(),
cg.relativeResidual(),
cg.absoluteResidual());
}

gaspi_segment_delete(0);
gaspi_proc_term(GASPI_BLOCK);
return 0;
}```