//////////
void compress ();
//////////
- int operator() (const unsigned int i, const unsigned int j);
+ int operator() (const unsigned int i, const unsigned int j) const;
//////////
void add (const unsigned int i, const unsigned int j);
//////////
* access to the rowstart array, but
* readonly.
*
+ * Though the return value is declared
+ * #const#, you shoudl be aware that it
+ * may change if you call any nonconstant
+ * function of objects which operate on
+ * it.
+ *
* You should use this interface very
* carefully and only if you are absolutely
* sure to know what you do. You should
* access to the colnums array, but
* readonly.
*
+ * Though the return value is declared
+ * #const#, you shoudl be aware that it
+ * may change if you call any nonconstant
+ * function of objects which operate on
+ * it.
+ *
* You should use this interface very
* carefully and only if you are absolutely
* sure to know what you do. You should
*/
class dSMatrix
{
- dSMatrixStruct * cols;
+ const dSMatrixStruct * cols;
double* val;
unsigned int max_len;
public:
//
void Tvmult (dVector& dst,const dVector& src) 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;
+
//
double residual (dVector& dst, const dVector& x, const dVector& b);
//
* Return a (constant) reference to the
* underlying sparsity pattern of this
* matrix.
+ *
+ * Though the return value is declared
+ * #const#, you shoudl be aware that it
+ * may change if you call any nonconstant
+ * function of objects which operate on
+ * it.
*/
const dSMatrixStruct & get_sparsity_pattern () const;
int
-dSMatrixStruct::operator () (const unsigned int i, const unsigned int j)
+dSMatrixStruct::operator () (const unsigned int i, const unsigned int j) const
{
+ Assert (i<rows, ExcInvalidIndex(i,rows));
+ Assert (j<cols, ExcInvalidIndex(j,cols));
+
for (unsigned int k=rowstart[i] ; k<rowstart[i+1] ; k++)
if (colnums[k] == (signed int)j) return k;
return -1;
}
}
+
+
+double
+dSMatrix::matrix_norm (const dVector& v) const
+{
+ Assert (cols != 0, ExcMatrixNotInitialized());
+ Assert(m() == v.n(), ExcDimensionsDontMatch(m(),v.n()));
+ Assert(n() == v.n(), ExcDimensionsDontMatch(n(),v.n()));
+
+ double sum = 0.;
+ for (unsigned int i=0;i<m();i++)
+ {
+ double s = 0.;
+ for (unsigned int j=cols->rowstart[i]; j<cols->rowstart[i+1]; ++j)
+ s += val[j] * v(cols->colnums[j]);
+ sum += s*v(i);
+ };
+
+ return sum;
+}
+
+
+
double
dSMatrix::residual (dVector& dst, const dVector& u, const dVector& b)
{