From d23ba50998419f852e7d81c4a4a685eab53eeb9f Mon Sep 17 00:00:00 2001
From: Wolfgang Bangerth <bangerth@math.tamu.edu>
Date: Tue, 21 Jul 1998 08:53:37 +0000
Subject: [PATCH] Implement a matrix scalar product for the dFMatrix class.

git-svn-id: https://svn.dealii.org/trunk@450 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/lac/include/lac/dfmatrix.h | 32 ++++++++++++++++++
 deal.II/lac/source/dfmatrix.cc     | 53 ++++++++++++++++++++++++++++++
 2 files changed, 85 insertions(+)

diff --git a/deal.II/lac/include/lac/dfmatrix.h b/deal.II/lac/include/lac/dfmatrix.h
index c0254a2ed6..104513ee17 100644
--- a/deal.II/lac/include/lac/dfmatrix.h
+++ b/deal.II/lac/include/lac/dfmatrix.h
@@ -237,6 +237,38 @@ class dFMatrix
 				      */
     void Tvmult (dVector& w, const dVector& v, const bool adding=false) const;
 
+				     /**
+				      * Return the norm of the vector #v# with
+				      * respect to the norm induced by this
+				      * matrix, i.e. $\left<v,Mv\right>$. This
+				      * is useful, e.g. in the finite element
+				      * context, where the $L_2$ norm of a
+				      * function equals the matrix norm with
+				      * respect to the mass matrix of the vector
+				      * representing the nodal values of the
+				      * finite element function.
+				      *
+				      * Note the order in which the matrix
+				      * appears. For non-symmetric matrices
+				      * there is a difference whether the
+				      * matrix operates on the first
+				      * or on the second operand of the
+				      * scalar product.
+				      *
+				      * Obviously, the matrix needs to be square
+				      * for this operation.
+				      */
+    double matrix_norm (const dVector &v) const;
+
+				     /**
+				      * Build the matrix scalar product
+				      * #u^T M v#. This function is mostly
+				      * useful when building the cellwise
+				      * scalar product of two functions in
+				      * the finite element context.
+				      */
+    double matrix_scalar_product (const dVector &u, const dVector &v) const;
+    
 				     /**
 				      * A=Inverse(A). Inversion of this by
 				      * Gauss-Jordan-algorithm
diff --git a/deal.II/lac/source/dfmatrix.cc b/deal.II/lac/source/dfmatrix.cc
index 2f50134288..2b95838794 100644
--- a/deal.II/lac/source/dfmatrix.cc
+++ b/deal.II/lac/source/dfmatrix.cc
@@ -571,6 +571,59 @@ void dFMatrix::Tmmult (dFMatrix& dst, const dFMatrix& src) const
       }
 }*/
 
+
+
+double dFMatrix::matrix_norm (const dVector &v) const {
+  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
+  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
+
+  double sum = 0.;
+  const unsigned int n_rows = m();
+  const double *val_ptr = &val[0];
+  const double *v_ptr;
+  
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+      double s = 0.;
+      const double * const val_end_of_row = val_ptr+n_rows;
+      v_ptr = v.begin();
+      while (val_ptr != val_end_of_row)
+	s += *val_ptr++ * *v_ptr++;
+
+      sum += s* v(row);
+    };
+
+  return sum;
+};
+
+
+
+double dFMatrix::matrix_scalar_product (const dVector &u, const dVector &v) const {
+  Assert(m() == u.size(), ExcDimensionMismatch(m(),v.size()));
+  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
+
+  double sum = 0.;
+  const unsigned int n_rows = m();
+  const unsigned int n_cols = n();
+  const double *val_ptr = &val[0];
+  const double *v_ptr;
+  
+  for (unsigned int row=0; row<n_rows; ++row)
+    {
+      double s = 0.;
+      const double * const val_end_of_row = val_ptr+n_cols;
+      v_ptr = v.begin();
+      while (val_ptr != val_end_of_row)
+	s += *val_ptr++ * *v_ptr++;
+
+      sum += s* u(row);
+    };
+
+  return sum;
+};
+
+
+
 void dFMatrix::print (FILE* f, const char* format) const
 {
   if (!format) format = " %5.2f";
-- 
2.39.5