quad: Use Boost matrix.

Signed-off-by: Matthew Emmett <[email protected]>

Author: memmett

src/pfasst-encapsulated.hpp

@@ -114,7 +114,7 @@ namespace pfasst {
 	auto* dst = dynamic_cast<EncapsulatedSweeperMixin<scalar>*>(DST);
 	auto* src = dynamic_cast<const EncapsulatedSweeperMixin<scalar>*>(SRC);
 
-	if (tmat.size() == 0)
+	if (tmat.size1() == 0)
 	  tmat = pfasst::compute_interp<scalar>(dst->get_nodes(), src->get_nodes());
 
 	int ndst = dst->get_nodes().size();

src/pfasst-quadrature.hpp

@@ -8,30 +8,15 @@
 #include <limits>
 #include <vector>
 
-using namespace std;
+#include <boost/numeric/ublas/matrix.hpp>
+
+using boost::numeric::ublas::matrix;
+using std::vector;
 
 namespace pfasst {
 
   typedef unsigned int uint;
 
-  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);
-      fill(this->begin(), this->end(), 0.0);
-    }
-    T& operator()(unsigned int i, unsigned int j) {
-      return (*this)[i*m+j];
-    }
-  };
-
   template<typename coeffT>
   class polynomial {
     vector<coeffT> c;
@@ -156,23 +141,23 @@ namespace pfasst {
 
   //#define pi 3.1415926535897932384626433832795028841971693993751
 
-  template<typename nodeT>
-  vector<nodeT> compute_nodes(int nnodes, string qtype)
+  template<typename time>
+  vector<time> compute_nodes(int nnodes, string qtype)
   {
-    vector<nodeT> nodes(nnodes);
+    vector<time> nodes(nnodes);
 
     if (qtype == "gauss-legendre") {
-      auto roots = polynomial<nodeT>::legendre(nnodes).roots();
+      auto roots = polynomial<time>::legendre(nnodes).roots();
       for (int j=0; j<nnodes; j++)
       	nodes[j] = 0.5 * (1.0 + roots[j]);
     } else if (qtype == "gauss-lobatto") {
-      auto roots = polynomial<nodeT>::legendre(nnodes-1).differentiate().roots();
+      auto roots = polynomial<time>::legendre(nnodes-1).differentiate().roots();
       for (int j=0; j<nnodes-2; j++)
 	nodes[j+1] = 0.5 * (1.0 + roots[j]);
       nodes[0] = 0.0; nodes[nnodes-1] = 1.0;
     } else if (qtype == "gauss-radau") {
-      auto l   = polynomial<nodeT>::legendre(nnodes);
-      auto lm1 = polynomial<nodeT>::legendre(nnodes-1);
+      auto l   = polynomial<time>::legendre(nnodes);
+      auto lm1 = polynomial<time>::legendre(nnodes-1);
       for (int i=0; i<nnodes; i++)
 	l[i] += lm1[i];
       auto roots = l.roots();
@@ -184,18 +169,18 @@ namespace pfasst {
     return nodes;
   }
 
-  template<typename nodeT>
-  matrix<nodeT> compute_quadrature(vector<nodeT> dst, vector<nodeT> src, char type)
+  template<typename time>
+  matrix<time> compute_quadrature(vector<time> dst, vector<time> src, char type)
   {
     const int ndst = dst.size();
     const int nsrc = src.size();
 
-    matrix<nodeT> mat(ndst-1, nsrc);
+    matrix<time> mat(ndst-1, nsrc);
 
     //   /* for (int n=0; n<(ndst-1)*nsrc; n++) */
     //   /*   smat[n] = 0.0; */
 
-    polynomial<nodeT> p(nsrc+1), p1(nsrc+1);
+    polynomial<time> p(nsrc+1), p1(nsrc+1);
 
     for (int i=0; i<nsrc; i++) {
       //      if ((flags[i] & SDC_NODE_PROPER) == 0) continue;
@@ -216,7 +201,7 @@ namespace pfasst {
       auto den = p.evaluate(src[i]);
       auto P = p.integrate();
       for (int j=1; j<ndst; j++) {
-	nodeT q = 0.0;
+	time q = 0.0;
 	if (type == 's')
 	  q = P.evaluate(dst[j]) - P.evaluate(dst[j-1]);
 	else

src/pfasst-vector.hpp

@@ -50,19 +50,23 @@ namespace pfasst {
       void mat_apply(vector<Encapsulation<scalar>*> DST, time a, matrix<time> mat,
 		     vector<Encapsulation<scalar>*> SRC, bool zero=true) {
 
-	vector<VectorEncapsulation<scalar,time>*> dst(mat.n), src(mat.m);
-	for (int n=0; n<mat.n; n++)
+	int ndst = DST.size();
+	int nsrc = SRC.size();
+
+	vector<VectorEncapsulation<scalar,time>*> dst(ndst), src(nsrc);
+	for (int n=0; n<ndst; n++)
 	  dst[n] = dynamic_cast<VectorEncapsulation*>(DST[n]);
-	for (int m=0; m<mat.m; m++)
+	for (int m=0; m<nsrc; m++)
 	  src[m] = dynamic_cast<VectorEncapsulation*>(SRC[m]);
 
 	if (zero)
-	  for (int n=0; n<mat.n; n++)
+	  for (int n=0; n<ndst; n++)
 	    dst[n]->setval(0.0);
 
-	for (int i=0; i<(*dst[0]).size(); i++)
-	  for (int n=0; n<mat.n; n++)
-	    for (int m=0; m<mat.m; m++)
+	int ndofs = (*dst[0]).size();
+	for (int i=0; i<ndofs; i++)
+	  for (int n=0; n<ndst; n++)
+	    for (int m=0; m<nsrc; m++)
 	      (*dst[n])[i] += a * mat(n,m) * (*src[m])[i];
       }