From 9d27faee5b1e6f392b15b952cce6af5bd9920460 Mon Sep 17 00:00:00 2001
From: "Roy H. Stogner" <roystgnr@ices.utexas.edu>
Date: Wed, 27 Mar 2013 14:00:07 -0500
Subject: [PATCH] Patch based on Manav Bhatia's type-mixing methods

---
 include/numerics/dense_matrix.h | 63 +++++++++++++++++++++++++++++----
 include/numerics/dense_vector.h | 12 +++----
 2 files changed, 63 insertions(+), 12 deletions(-)

diff --git a/include/numerics/dense_matrix.h b/include/numerics/dense_matrix.h
index 49b44923db6..5399ef1e156 100644
--- a/include/numerics/dense_matrix.h
+++ b/include/numerics/dense_matrix.h
@@ -68,7 +68,7 @@ class DenseMatrix : public DenseMatrixBase<T>
    */
   virtual ~DenseMatrix() {}
 
-
+   
   /**
    * Set every element in the matrix to 0.
    */
@@ -143,6 +143,16 @@ class DenseMatrix : public DenseMatrixBase<T>
    */
   DenseMatrix<T>& operator = (const DenseMatrix<T>& other_matrix);
 
+  /**
+   * Assignment-from-other-matrix-type operator
+   *
+   * Copies the dense matrix of type T2 into the present matrix.  This
+   * is useful for copying real matrices into complex ones for further
+   * operations
+   */
+  template <typename T2>
+  DenseMatrix<T>& operator = (const DenseMatrix<T2>& other_matrix);
+ 
   /**
    * STL-like swap method
    */
@@ -160,6 +170,12 @@ class DenseMatrix : public DenseMatrixBase<T>
    */
   void scale (const T factor);
 
+    
+    /**
+     * Multiplies every element in the column \p col matrix by \p factor.
+     */
+    void scale_column (const unsigned int col, const T factor);
+
   /**
    * Multiplies every element in the matrix by \p factor.
    */
@@ -168,8 +184,10 @@ class DenseMatrix : public DenseMatrixBase<T>
   /**
    * Adds \p factor times \p mat to this matrix.
    */
-  void add (const T factor,
-            const DenseMatrix<T>& mat);
+    template<typename T2, typename T3>
+    typename boostcopy::enable_if_c<
+    ScalarTraits<T2>::value, void >::type add (const T2 factor,
+                                               const DenseMatrix<T3>& mat);
 
   /**
    * Tests if \p mat is exactly equal to this matrix.
@@ -584,6 +602,22 @@ void DenseMatrix<T>::swap(DenseMatrix<T>& other_matrix)
 }
 
 
+template <typename T>
+template <typename T2>
+inline
+DenseMatrix<T>&
+DenseMatrix<T>::operator=(const DenseMatrix<T2>& mat)
+{
+  unsigned int mat_m = mat.m(), mat_n = mat.n();
+  this->resize(mat_m, mat_n);
+  for (unsigned int i=0; i<mat_m; i++)
+    for (unsigned int j=0; j<mat_n; j++)
+      (*this)(i,j) = mat(i,j);
+
+  return *this;
+}
+
+    
 
 template<typename T>
 inline
@@ -674,6 +708,15 @@ void DenseMatrix<T>::scale (const T factor)
     _val[i] *= factor;
 }
 
+    
+template<typename T>
+inline
+void DenseMatrix<T>::scale_column (const unsigned int col, const T factor)
+{
+    for (unsigned int i=0; i<this->m(); i++)
+        (*this)(i, col) *= factor;
+}
+
 
 
 template<typename T>
@@ -687,11 +730,19 @@ DenseMatrix<T>& DenseMatrix<T>::operator *= (const T factor)
 
 
 template<typename T>
+template<typename T2, typename T3>
 inline
-void DenseMatrix<T>::add (const T factor, const DenseMatrix<T>& mat)
+typename boostcopy::enable_if_c<
+ScalarTraits<T2>::value, void >::type
+DenseMatrix<T>::add (const T2 factor,
+                     const DenseMatrix<T3>& mat)
 {
-  for (unsigned int i=0; i<_val.size(); i++)
-    _val[i] += factor * mat._val[i];
+    libmesh_assert_equal_to (this->m(), mat.m());
+    libmesh_assert_equal_to (this->n(), mat.n());
+    
+    for (unsigned int i=0; i<this->m(); i++)
+        for (unsigned int j=0; j<this->n(); j++)
+            (*this)(i,j) += factor * mat(i,j);
 }
 
 
diff --git a/include/numerics/dense_vector.h b/include/numerics/dense_vector.h
index 715eff17d27..76fb309244c 100644
--- a/include/numerics/dense_vector.h
+++ b/include/numerics/dense_vector.h
@@ -149,14 +149,14 @@ class DenseVector : public DenseVectorBase<T>
    * complex conjugate of vec.
    */
   template <typename T2>
-  Number dot (const DenseVector<T2> &vec) const;
+  typename CompareTypes<T, T2>::supertype dot (const DenseVector<T2> &vec) const;
 
   /**
    * Evaluate dot product with \p vec. In the complex-valued case, do not
    * use the complex conjugate of vec.
    */
   template <typename T2>
-  Number indefinite_dot (const DenseVector<T2> &vec) const;
+    typename CompareTypes<T, T2>::supertype indefinite_dot (const DenseVector<T2> &vec) const;
 
   /**
    * Tests if \p vec is exactly equal to this vector.
@@ -397,11 +397,11 @@ DenseVector<T>::add (const T2 factor,
 template<typename T>
 template<typename T2>
 inline
-Number DenseVector<T>::dot (const DenseVector<T2>& vec) const
+  typename CompareTypes<T, T2>::supertype DenseVector<T>::dot (const DenseVector<T2>& vec) const
 {
   libmesh_assert_equal_to (this->size(), vec.size());
 
-  Number val = 0.;
+  typename CompareTypes<T, T2>::supertype val = 0.;
 
   for (unsigned int i=0; i<this->size(); i++)
     val += (*this)(i)*libmesh_conj(vec(i));
@@ -412,11 +412,11 @@ Number DenseVector<T>::dot (const DenseVector<T2>& vec) const
 template<typename T>
 template<typename T2>
 inline
-Number DenseVector<T>::indefinite_dot (const DenseVector<T2>& vec) const
+  typename CompareTypes<T, T2>::supertype DenseVector<T>::indefinite_dot (const DenseVector<T2>& vec) const
 {
   libmesh_assert_equal_to (this->size(), vec.size());
 
-  Number val = 0.;
+  typename CompareTypes<T, T2>::supertype val = 0.;
 
   for (unsigned int i=0; i<this->size(); i++)
     val += (*this)(i)*(vec(i));