Clean of dgshell

Integration on bulk and virtual interface terms are performed by a
matrix vector product. Add 2 functions in fullMatrix.h and cpp to
realize this operation

Numeric/fullMatrix.cpp

@@ -156,6 +156,27 @@ void fullMatrix<std::complex<double> >::multAddy(const fullVector<std::complex<d
                  &beta, y._data, &INCY);
 }
 
+
+template<>
+void fullMatrix<double>::multOnBlock(const fullMatrix<double> &b, const int ncol, const int fcol, const int alpha_, const int beta_, fullVector<double> &c,const int row) const
+{
+  int M = 1, N = ncol, K = b.size1() ;
+  int LDA = _r, LDB = b.size1(), LDC = 1;
+  double alpha = alpha_, beta = beta_;
+  F77NAME(dgemm)("N", "N", &M, &N, &K, &alpha, _data, &LDA, &(b._data[fcol*K]), &LDB,
+                 &beta, &(c._data[fcol]), &LDC);
+}
+
+template<>
+void fullMatrix<double>::multWithATranspose(const fullVector<double> &x, const int alpha_, const int beta_,fullVector<double> &y) const
+{
+  int M = _r, N = _c, LDA = _r, INCX = 1, INCY = 1;
+  double alpha = alpha_, beta = beta_;
+  F77NAME(dgemv)("T", &M, &N, &alpha, _data, &LDA, x._data, &INCX,
+                 &beta, y._data, &INCY);
+
+}
+
 #endif
 
 

Numeric/fullMatrix.h

@@ -517,6 +517,35 @@ class fullMatrix
     }
   }
 
+// specific functions for dgshell
+  void mult_naiveBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol, const int alpha, const int beta, fullVector<scalar> &c, const int row=0) const
+  {
+    if(beta != 1)
+      c.scale(beta);
+    for(int j = fcol; j < fcol+ncol; j++)
+      for(int k = 0; k < _c ; k++)
+          c._data[j] += alpha*(*this)(row, k) * b(k, j);
+  }
+  void multOnBlock(const fullMatrix<scalar> &b, const int ncol, const int fcol, const int alpha, const int beta, fullVector<scalar> &c,const int row=0) const
+#if !defined(HAVE_BLAS)
+  {
+    mult_naiveBlock(b,ncol,fcol,alpha,beta,c);
+  }
+#endif
+  ;
+
+  void multWithATranspose(const fullVector<scalar> &x, const int alpha_, const int beta_, fullVector<scalar> &y) const
+#if !defined(HAVE_BLAS)
+  {
+    y.scale(beta_);
+    for(int j = 0; j < _c; j++)
+      for(int i = 0; i < _r; i++)
+        y._data[j] += (*this)(i, j) * x(i);
+  }
+#endif
+  ;
+
+
   static void registerBindings(binding *b);
 };