First stab at a class structure.

Author: memmett

README.md

@@ -1,4 +1,4 @@
 PFASST
 ======
 
-C++ PFASST
+Proposed class structure for a C++ PFASST implementation.

src/pfasst-encapsulated.hpp

@@ -0,0 +1,117 @@
+/*
+ * Host based encapsulated base sweeper.
+ */
+
+#ifndef _PFASST_ENCAPSULATED_HPP_
+#define _PFASST_ENCAPSULATED_HPP_
+
+#include <vector>
+
+#include "pfasst.hpp"
+
+using namespace std;
+
+namespace pfasst {
+
+  template<typename T>
+  vector<T> compute_nodes(unsigned int nnodes, string qtype) {
+    vector<T> nodes(nnodes);
+
+    // ...
+
+    return nodes;
+  }
+
+  typedef enum encaptype { solution, function } encaptype;
+
+  template<typename T>
+  class matrix : public vector<T> {
+
+  public:
+    unsigned int n, m;
+    matrix() { }
+    matrix(unsigned int n, unsigned int m) {
+      zeros(n, m);
+    }
+    void zeros(unsigned int n, unsigned int m) {
+      this->n = n; this->m = m;
+      this->resize(n*m);
+      // ...
+    }
+    T& operator()(unsigned int i, unsigned int j) {
+      return (*this)[i*m+j];
+    }
+  };
+
+  //
+  // encapsulation
+  //
+
+  struct encapsulation {
+    virtual ~encapsulation() { }
+
+    // required for time-parallel communications
+    virtual unsigned int nbytes() { }
+    virtual void pack(char *buf) { }
+    virtual void unpack(char *buf) { }
+
+    // required for interp/restrict helpers
+    virtual void interpolate(const encapsulation *) { }
+    virtual void restrict(const encapsulation *) { }
+
+    // required for host based encap helpers
+    virtual void setval(double) { }
+    virtual void copy(const encapsulation *) { }
+    virtual void mat_apply(encapsulation dst[], double a, matrix m, const encapsulation src[]) { }
+  };
+
+  struct encapsulation_factory {
+    virtual encapsulation* create(const encaptype) = 0;
+  };
+
+
+  template<typename T>
+  class encapsulated_sweeper_mixin : public isweeper {
+    shared_ptr<vector<T>>             nodes;
+    shared_ptr<encapsulation_factory> encap;
+
+  public:
+    vector<encapsulation*> q;
+    vector<T>* get_nodes() { return nodes.get(); }
+
+    virtual void set_q0(const encapsulation* q0) { }
+    virtual encapsulation* get_qend() { }
+  };
+
+  template<class T>
+  class poly_interp_mixin : public T {
+    virtual void interpolate(const isweeper*) { }
+    virtual void restrict(const isweeper*) { }
+  };
+
+  template<typename T>
+  struct vector_encapsulation : public vector<T>, public encapsulation {
+    vector_encapsulation(int size) : vector<T>(size) { }
+    virtual unsigned int nbytes() const {
+      return sizeof(T) * this->size();
+    }
+    void setval(double v) {
+      for (int i=0; i<this->size(); i++)
+	this->data()[i] = v;
+    }
+    // ...
+  };
+
+  template<typename T>
+  class vector_factory : public pfasst::encapsulation_factory {
+    int size;
+  public:
+    vector_factory(const int size) : size(size) { }
+    encapsulation* create(const pfasst::encap_type) {
+      return new vector_encapsulation<T>(size);
+    }
+  };
+
+}
+
+#endif

src/pfasst.hpp

@@ -0,0 +1,111 @@
+/*
+ * Interfaces for SDC/MLSDC/PFASST algorithms.
+ */
+
+#ifndef _PFASST_HPP_
+#define _PFASST_HPP_
+
+#include <deque>
+#include <exception>
+#include <iostream>
+#include <iterator>
+#include <memory>
+#include <string>
+
+using namespace std;
+
+namespace pfasst {
+
+  //
+  // sdc sweeper interface
+  //
+
+  struct isweeper {
+    virtual void setup() { }
+    virtual ~isweeper() { }
+
+    // required for all sdc schemes
+    virtual void sweep(double t, double dt) = 0;
+    virtual void predict(double t, double dt) = 0;
+
+    // required for multi-level sdc and pfasst schemes
+    virtual void interpolate(const isweeper*) { }
+    virtual void restrict(const isweeper*) { }
+
+    // required for pfasst schemes
+    virtual void post() { }
+    virtual void send() { }
+    virtual void recv() { }
+  };
+
+  //
+  // pfasst controller
+  //
+
+  class pfasst {
+    deque<shared_ptr<isweeper>> levels;
+
+  public:
+    int nstep, niter;
+    double dt, t;
+
+    void add_level(isweeper *sweeper, bool coarse) {
+      if (coarse)
+	levels.push_front(shared_ptr<isweeper>(sweeper));
+      else
+	levels.push_back(shared_ptr<isweeper>(sweeper));
+    }
+
+    template<typename R=isweeper> R* get_level(int level) {
+      return dynamic_cast<R*>(levels[level].get());
+    }
+
+    int nlevels() { return levels.size(); }
+    void setup();
+    void run();
+
+    void predictor();
+    void iteration();
+
+    struct leveliter {
+      int level;
+      pfasst *pf;
+
+      leveliter(int level, pfasst *pf) : level(level), pf(pf) {}
+
+      template<typename R=isweeper> R* current() {
+	return pf->get_level<R>(level);
+      }
+      template<typename R=isweeper> R* fine() {
+	return pf->get_level<R>(level+1);
+      }
+      template<typename R=isweeper> R* coarse() {
+	return pf->get_level<R>(level-1);
+      }
+
+      isweeper *operator*() { return current(); }
+      bool operator==(leveliter i) { return level == i.level; }
+      bool operator!=(leveliter i) { return level != i.level; }
+      bool operator<=(leveliter i) { return level <= i.level; }
+      bool operator>=(leveliter i) { return level >= i.level; }
+      bool operator< (leveliter i) { return level <  i.level; }
+      bool operator> (leveliter i) { return level >  i.level; }
+      leveliter operator- (int i) { return leveliter(level-1, pf); }
+      leveliter operator+ (int i) { return leveliter(level+1, pf); }
+      void operator++() { level++; }
+      void operator--() { level--; }
+    };
+
+    leveliter finest()   { return leveliter(nlevels()-1, this); }
+    leveliter coarsest() { return leveliter(0, this); }
+
+    leveliter cycle_down(leveliter levels, double t, double dt);
+    leveliter cycle_up(leveliter levels, double t, double dt);
+    leveliter cycle_bottom(leveliter levels, double t, double dt);
+    leveliter cycle_top(leveliter levels, double t, double dt);
+    leveliter cycle_v(leveliter levels, double t, double dt);
+  };
+
+}
+
+#endif