*/
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
}
}*/
+
+
+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";