Example

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

The full example code is also included in the Download.

// 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);
    }
  }
  // fully threaded matrix finalization
  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;
  }

  // Read in problem size.
  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 
  // threading package, such as OpenMP. Please refer to the GaspiLS 
  // documentation for thread pinning options.
  GaspiLS::ThreadPool             pool(atoi(argv[4]));
  // 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;
}