Added: Option to return the sum of a specified matrix column

Author: kmokstad

Diff for: src/ASM/LR/ASMu2D.C

@@ -503,12 +503,12 @@ bool ASMu2D::evaluateBasis (int iel, FiniteElement& fe, int derivs) const
   dNdu.fillColumn(2,C*Bv);
 
 #ifdef SP_DEBUG
-  if (fabs(dNdu.getColumn(1).sum()) > 1.0e-10) {
+  if (fabs(dNdu.sum(-1)) > 1.0e-10) {
     std::cerr <<"dNdu do not sum to zero at integration point #"
               << fe.iGP << std::endl;
     return false;
   }
-  else if (fabs(dNdu.getColumn(2).sum()) > 1.0e-10) {
+  else if (fabs(dNdu.sum(-2)) > 1.0e-10) {
     std::cerr <<"dNdv do not sums to zero at integration point #"
               << fe.iGP << std::endl;
     return false;

Diff for: src/ASM/LR/ASMu2Dnurbs.C

@@ -57,7 +57,7 @@ bool ASMu2D::evaluateBasisNurbs (int iel, FiniteElement& fe,
     if (fabs(fe.N.sum()-1.0) > 1.0e-10)
       std::cerr <<"fe.N do not sum to one at integration point #"
                 << fe.iGP << std::endl;
-    else if (fabs(static_cast<const Vector&>(B).sum()-1.0) > 1.0e-10)
+    else if (fabs(B.sum()-1.0) > 1.0e-10)
       std::cerr <<"Bezier basis do not sum to one at integration point #"
                 << fe.iGP << std::endl;
     else
@@ -107,10 +107,10 @@ bool ASMu2D::evaluateBasisNurbs (int iel, FiniteElement& fe,
     else if (fabs(B.sum()-1.0) > 1.0e-10)
       std::cerr <<"Bezier basis do not sum to one at integration point #"
                 << fe.iGP << std::endl;
-    else if (fabs(dNdu.getColumn(1).sum()) > 1.0e-10)
+    else if (fabs(dNdu.sum(-1)) > 1.0e-10)
       std::cerr <<"dNdu not sums to zero at integration point #"
                 << fe.iGP << std::endl;
-    else if (fabs(dNdu.getColumn(2).sum()) > 1.0e-10)
+    else if (fabs(dNdu.sum(-2)) > 1.0e-10)
       std::cerr <<"dNdv not sums to zero at integration point #"
                 << fe.iGP << std::endl;
     else if (fabs(Bu.sum()) > 1.0e-10 || fabs(Bv.sum()) > 1.0e-10)

Diff for: src/LinAlg/Test/TestMatrix.C

@@ -122,6 +122,17 @@ TEST(TestMatrix, AugmentCols)
 }
 
 
+TEST(TestMatrix, SumCols)
+{
+  utl::matrix<int> a(5,3);
+  std::iota(a.begin(),a.end(),1);
+  std::cout <<"A:"<< a;
+  EXPECT_EQ(a.sum(-1), 15);
+  EXPECT_EQ(a.sum(-2), 40);
+  EXPECT_EQ(a.sum(-3), 65);
+}
+
+
 TEST(TestMatrix, Multiply)
 {
   utl::vector<double> u(14), v(9), x, y;

Diff for: src/LinAlg/matrix.h

@@ -414,8 +414,21 @@ namespace utl //! General utility classes and functions.
     //! \param[in] inc Increment in the matrix element vector indices
     T asum(int inc = 1) const { return elem.asum(0,inc); }
     //! \brief Return the sum of the matrix elements.
-    //! \param[in] inc Increment in the matrix element vector indices
-    T sum(int inc = 1) const { return elem.sum(0,inc); }
+    //! \param[in] inc Increment in the matrix element vector indices.
+    //! If negative, the sum of matrix column \a -inc will be returned instead.
+    T sum(int inc = 1) const
+    {
+      if (inc > 0)
+        return elem.sum(0,inc);
+      else if (inc == 0 || -inc > static_cast<int>(n[1]))
+        return T(0);
+
+      T colsum = T(0);
+      size_t ofs = n[0]*(-inc-1);
+      for (size_t i = 0; i < n[0]; i++, ofs++)
+        colsum += elem[ofs];
+      return colsum;
+    }
 
   protected:
     size_t     n[4]; //!< Dimension of the matrix