/**
- * Rectangular/quadratic full matrix.
+ * @brief Rectangular/quadratic full matrix.
*
* Implementation of a classical rectangular scheme of numbers. The
* data type of the entries is provided in the template argument
*/
Accessor accessor;
};
-
+/**
+ * @name Constructors Constructors and initalization.
+ * See also the base class Table.
+ */
+//@{
+
/**
* Constructor. Initialize the
* matrix as a square matrix with
template <class MATRIX>
void copy_from (const MATRIX&);
- /**
- * Comparison operator. Be
- * careful with this thing, it
- * may eat up huge amounts of
- * computing time! It is most
- * commonly used for internal
- * consistency checks of
- * programs.
- */
- bool operator == (const FullMatrix<number> &) const;
-
/**
* Fill rectangular block.
*
void fill_permutation (const FullMatrix<number2> &src,
const std::vector<unsigned int> &p_rows,
const std::vector<unsigned int> &p_cols);
+
+//@}
+///@name Non-modifying operators
+//@{
+ /**
+ * Comparison operator. Be
+ * careful with this thing, it
+ * may eat up huge amounts of
+ * computing time! It is most
+ * commonly used for internal
+ * consistency checks of
+ * programs.
+ */
+ bool operator == (const FullMatrix<number> &) const;
+
/**
* Number of rows of this matrix.
* To remember: this matrix is an
*/
bool all_zero () const;
+ /**
+ * Return the square of the norm
+ * of the vector <tt>v</tt> with
+ * respect to the norm induced by
+ * this matrix,
+ * i.e. <i>(v,Mv)</i>. This is
+ * useful, e.g. in the finite
+ * element context, where the
+ * <i>L<sup>2</sup></i> 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.
+ *
+ * Obviously, the matrix needs to
+ * be quadratic for this operation.
+ */
+ template<typename number2>
+ number2 matrix_norm_square (const Vector<number2> &v) const;
+
+ /**
+ * Build the matrix scalar product
+ * <tt>u^T M v</tt>. This function is mostly
+ * useful when building the cellwise
+ * scalar product of two functions in
+ * the finite element context.
+ */
+ template<typename number2>
+ number2 matrix_scalar_product (const Vector<number2> &u,
+ const Vector<number2> &v) const;
+
+ /**
+ * Return the $l_1$-norm of the matrix, i.e.
+ * $|M|_1=max_{all columns j}\sum_{all
+ * rows i} |M_ij|$,
+ * (max. sum of columns). This is the
+ * natural matrix norm that is compatible
+ * to the $l_1$-norm for vectors, i.e.
+ * $|Mv|_1\leq |M|_1 |v|_1$.
+ * (cf. Rannacher Numerik0)
+ */
+ number l1_norm () const;
+
+ /**
+ * Return the $l_\infty$-norm of the
+ * matrix, i.e.
+ * $|M|_\infty=\max_{all rows i}\sum_{all
+ * columns j} |M_{ij}|$,
+ * (max. sum of rows).
+ * This is the
+ * natural matrix norm that is compatible
+ * to the $l_\infty$-norm of vectors, i.e.
+ * $|Mv|_\infty \leq |M|_\infty |v|_\infty$.
+ * (cf. Rannacher Numerik0)
+ */
+ number linfty_norm () const;
+
+ /**
+ * Compute the quadratic matrix norm.
+ * Return value is the root of the square
+ * sum of all matrix entries. Also called
+ * Frobenius norm.
+ *
+ * This norm is compatible with the $l_2$
+ * vector norm. But it is not a natural
+ * matrix norm (cf Rannacher Numeric0),
+ * therefore it is not called $l_2$-norm.
+ */
+ number norm2 () const;
+
+ /**
+ * Compute the relative norm of
+ * the skew-symmetric part. The
+ * return value is the Frobenius
+ * norm of the skew-symmetric
+ * part of the matrix divided by
+ * that of the matrix.
+ *
+ * Main purpose of this function
+ * is to check, if a matrix is
+ * symmetric within a certain
+ * accuracy, or not.
+ */
+ number relative_symmetry_norm2 () const;
+
+ /**
+ * Computes the determinant of a
+ * matrix. This is only
+ * implemented for one, two, and
+ * three dimensions, since for
+ * higher dimensions the
+ * numerical work explodes.
+ * Obviously, the matrix needs to
+ * be quadratic for this function.
+ */
+ double determinant () const;
+
+ /**
+ * Output of the matrix in
+ * user-defined format.
+ */
+ void print (std::ostream &s,
+ const unsigned int width=5,
+ const unsigned int precision=2) const;
+
+ /**
+ * Print the matrix in the usual
+ * format, i.e. as a matrix and
+ * not as a list of nonzero
+ * elements. For better
+ * readability, elements not in
+ * the matrix are displayed as
+ * empty space, while matrix
+ * elements which are explicitly
+ * set to zero are displayed as
+ * such.
+ *
+ * The parameters allow for a
+ * flexible setting of the output
+ * format: <tt>precision</tt> and
+ * <tt>scientific</tt> are used to
+ * determine the number format,
+ * where <tt>scientific</tt> = <tt>false</tt>
+ * means fixed point notation. A
+ * zero entry for <tt>width</tt> makes
+ * the function compute a width,
+ * but it may be changed to a
+ * positive value, if output is
+ * crude.
+ *
+ * Additionally, a character for
+ * an empty value may be
+ * specified.
+ *
+ * Finally, the whole matrix can
+ * be multiplied with a common
+ * denominator to produce more
+ * readable output, even
+ * integers.
+ *
+ * @attention This function
+ * may produce <b>large</b> amounts of
+ * output if applied to a large matrix!
+ */
+ void print_formatted (std::ostream &out,
+ const unsigned int presicion=3,
+ const bool scientific = true,
+ const unsigned int width = 0,
+ const char *zero_string = " ",
+ const double denominator = 1.) const;
+
+ /**
+ * Determine an estimate for the
+ * memory consumption (in bytes)
+ * of this object.
+ */
+ unsigned int memory_consumption () const;
+
+//@}
+///@name Iterator functions
+//@{
+
/**
* STL-like iterator with the
* first entry.
* Final iterator of row <tt>r</tt>.
*/
const_iterator end (const unsigned int r) const;
+
+//@}
+///@name Modifying operators
+//@{
/**
* Scale the entire matrix by a
const unsigned int dst_offset_j = 0,
const unsigned int src_offset_i = 0,
const unsigned int src_offset_j = 0);
-
+
/**
- * Matrix-matrix-multiplication.
- *
- * The optional parameter
- * <tt>adding</tt> determines, whether the
- * result is stored in <tt>C</tt> or added
- * to <tt>C</tt>.
- *
- * if (adding)
- * $C += A*B$
- *
- * if (!adding)
- * $C = A*B$
- *
- * Assumes that <tt>A</tt> and <tt>B</tt> have
- * compatible sizes and that <tt>C</tt>
- * already has the right size.
+ * $A(i,1...n)+=s*A(j,1...n)$.
+ * Simple addition of rows of this
*/
- template<typename number2>
- void mmult (FullMatrix<number2> &C,
- const FullMatrix<number2> &B,
- const bool adding=false) const;
-
+ void add_row (const unsigned int i,
+ const number s,
+ const unsigned int j);
+
/**
- * Matrix-matrix-multiplication using
- * transpose of <tt>this</tt>.
- *
- * The optional parameter
- * <tt>adding</tt> determines, whether the
- * result is stored in <tt>C</tt> or added
- * to <tt>C</tt>.
- *
- * if (adding)
- * $C += A^T*B$
- *
- * if (!adding)
- * $C = A^T*B$
- *
- * Assumes that <tt>A</tt> and <tt>B</tt> have
- * compatible sizes and that <tt>C</tt>
- * already has the right size.
+ * $A(i,1...n)+=s*A(j,1...n)+t*A(k,1...n)$.
+ * Multiple addition of rows of this.
*/
- template<typename number2>
- void Tmmult (FullMatrix<number2> &C,
- const FullMatrix<number2> &B,
- const bool adding=false) const;
-
+ void add_row (const unsigned int i,
+ const number s, const unsigned int j,
+ const number t, const unsigned int k);
+
/**
- * Matrix-vector-multiplication.
- *
- * The optional parameter
- * <tt>adding</tt> determines, whether the
- * result is stored in <tt>w</tt> or added
- * to <tt>w</tt>.
- *
- * if (adding)
- * $w += A*v$
- *
- * if (!adding)
- * $w = A*v$
- *
- * Source and destination must
- * not be the same vector.
+ * $A(1...n,i)+=s*A(1...n,j)$.
+ * Simple addition of columns of this.
*/
- template<typename number2>
- void vmult (Vector<number2> &w,
- const Vector<number2> &v,
- const bool adding=false) const;
-
+ void add_col (const unsigned int i,
+ const number s,
+ const unsigned int j);
+
/**
- * Transpose
- * matrix-vector-multiplication.
- * See vmult() above.
- *
- * Source and destination must
- * not be the same vector.
+ * $A(1...n,i)+=s*A(1...n,j)+t*A(1...n,k)$.
+ * Multiple addition of columns of this.
*/
- template<typename number2>
- void Tvmult (Vector<number2> &w,
- const Vector<number2> &v,
- const bool adding=false) const;
+ void add_col (const unsigned int i,
+ const number s, const unsigned int j,
+ const number t, const unsigned int k);
/**
- * Return the square of the norm
- * of the vector <tt>v</tt> with
- * respect to the norm induced by
- * this matrix,
- * i.e. <i>(v,Mv)</i>. This is
- * useful, e.g. in the finite
- * element context, where the
- * <i>L<sup>2</sup></i> 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.
- *
- * Obviously, the matrix needs to
- * be quadratic for this operation.
+ * Swap A(i,1...n) <-> A(j,1...n).
+ * Swap rows i and j of this
*/
- template<typename number2>
- number2 matrix_norm_square (const Vector<number2> &v) const;
+ void swap_row (const unsigned int i,
+ const unsigned int j);
/**
- * Build the matrix scalar product
- * <tt>u^T M v</tt>. This function is mostly
- * useful when building the cellwise
- * scalar product of two functions in
- * the finite element context.
+ * Swap A(1...n,i) <-> A(1...n,j).
+ * Swap columns i and j of this
+ */
+ void swap_col (const unsigned int i,
+ const unsigned int j);
+
+ /**
+ * A(i,i)+=B(i,1...n). Addition of complete
+ * rows of B to diagonal-elements of this ; <p>
+ * ( i = 1 ... m )
*/
template<typename number2>
- number2 matrix_scalar_product (const Vector<number2> &u,
- const Vector<number2> &v) const;
+ void add_diag (const number s,
+ const FullMatrix<number2> &B);
+
+ /**
+ * Add constant to diagonal
+ * elements of this, i.e. add a
+ * multiple of the identity
+ * matrix.
+ */
+ void diagadd (const number s);
/**
* Symmetrize the matrix by
* quadratic for this operation.
*/
void symmetrize ();
-
- /**
- * Return the $l_1$-norm of the matrix, i.e.
- * $|M|_1=max_{all columns j}\sum_{all
- * rows i} |M_ij|$,
- * (max. sum of columns). This is the
- * natural matrix norm that is compatible
- * to the $l_1$-norm for vectors, i.e.
- * $|Mv|_1\leq |M|_1 |v|_1$.
- * (cf. Rannacher Numerik0)
- */
- number l1_norm () const;
/**
- * Return the $l_\infty$-norm of the
- * matrix, i.e.
- * $|M|_\infty=\max_{all rows i}\sum_{all
- * columns j} |M_{ij}|$,
- * (max. sum of rows).
- * This is the
- * natural matrix norm that is compatible
- * to the $l_\infty$-norm of vectors, i.e.
- * $|Mv|_\infty \leq |M|_\infty |v|_\infty$.
- * (cf. Rannacher Numerik0)
- */
- number linfty_norm () const;
-
- /**
- * Compute the quadratic matrix norm.
- * Return value is the root of the square
- * sum of all matrix entries. Also called
- * Frobenius norm.
- *
- * This norm is compatible with the $l_2$
- * vector norm. But it is not a natural
- * matrix norm (cf Rannacher Numeric0),
- * therefore it is not called $l_2$-norm.
- */
- number norm2 () const;
-
- /**
- * Compute the relative norm of
- * the skew-symmetric part. The
- * return value is the Frobenius
- * norm of the skew-symmetric
- * part of the matrix divided by
- * that of the matrix.
- *
- * Main purpose of this function
- * is to check, if a matrix is
- * symmetric within a certain
- * accuracy, or not.
- */
- number relative_symmetry_norm2 () const;
-
- /**
* A=Inverse(A). Inversion of
* this matrix by Gauss-Jordan
* algorithm with partial
*/
void gauss_jordan ();
- /**
- * Computes the determinant of a
- * matrix. This is only
- * implemented for one, two, and
- * three dimensions, since for
- * higher dimensions the
- * numerical work explodes.
- * Obviously, the matrix needs to
- * be quadratic for this function.
- */
- double determinant () const;
-
/**
* Assign the inverse of the
* given matrix to
template <typename number2>
void invert (const FullMatrix<number2> &M);
- /**
- * Apply the Jacobi
- * preconditioner, which
- * multiplies every element of
- * the <tt>src</tt> vector by the
- * inverse of the respective
- * diagonal element and
- * multiplies the result with the
- * damping factor <tt>omega</tt>.
- */
- template <typename somenumber>
- void precondition_Jacobi (Vector<somenumber> &dst,
- const Vector<somenumber> &src,
- const number omega = 1.) const;
-
- /**
- * $A(i,1...n)+=s*A(j,1...n)$.
- * Simple addition of rows of this
- */
- void add_row (const unsigned int i,
- const number s,
- const unsigned int j);
-
- /**
- * $A(i,1...n)+=s*A(j,1...n)+t*A(k,1...n)$.
- * Multiple addition of rows of this.
- */
- void add_row (const unsigned int i,
- const number s, const unsigned int j,
- const number t, const unsigned int k);
/**
- * $A(1...n,i)+=s*A(1...n,j)$.
- * Simple addition of columns of this.
+ * QR-factorization of a matrix.
+ * The orthogonal transformation
+ * Q is applied to the vector y
+ * and this matrix.
+ *
+ * After execution of
+ * householder, the upper
+ * triangle contains the
+ * resulting matrix R, the lower
+ * the incomplete factorization
+ * matrices.
*/
- void add_col (const unsigned int i,
- const number s,
- const unsigned int j);
+ template<typename number2>
+ void householder (Vector<number2> &y);
+//@}
+///@name Multiplications
+//@{
+
/**
- * $A(1...n,i)+=s*A(1...n,j)+t*A(1...n,k)$.
- * Multiple addition of columns of this.
+ * Matrix-matrix-multiplication.
+ *
+ * The optional parameter
+ * <tt>adding</tt> determines, whether the
+ * result is stored in <tt>C</tt> or added
+ * to <tt>C</tt>.
+ *
+ * if (adding)
+ * $C += A*B$
+ *
+ * if (!adding)
+ * $C = A*B$
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have
+ * compatible sizes and that <tt>C</tt>
+ * already has the right size.
*/
- void add_col (const unsigned int i,
- const number s, const unsigned int j,
- const number t, const unsigned int k);
-
+ template<typename number2>
+ void mmult (FullMatrix<number2> &C,
+ const FullMatrix<number2> &B,
+ const bool adding=false) const;
+
/**
- * Swap A(i,1...n) <-> A(j,1...n).
- * Swap rows i and j of this
+ * Matrix-matrix-multiplication using
+ * transpose of <tt>this</tt>.
+ *
+ * The optional parameter
+ * <tt>adding</tt> determines, whether the
+ * result is stored in <tt>C</tt> or added
+ * to <tt>C</tt>.
+ *
+ * if (adding)
+ * $C += A^T*B$
+ *
+ * if (!adding)
+ * $C = A^T*B$
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have
+ * compatible sizes and that <tt>C</tt>
+ * already has the right size.
*/
- void swap_row (const unsigned int i,
- const unsigned int j);
-
+ template<typename number2>
+ void Tmmult (FullMatrix<number2> &C,
+ const FullMatrix<number2> &B,
+ const bool adding=false) const;
+
/**
- * Swap A(1...n,i) <-> A(1...n,j).
- * Swap columns i and j of this
+ * Matrix-vector-multiplication.
+ *
+ * The optional parameter
+ * <tt>adding</tt> determines, whether the
+ * result is stored in <tt>w</tt> or added
+ * to <tt>w</tt>.
+ *
+ * if (adding)
+ * $w += A*v$
+ *
+ * if (!adding)
+ * $w = A*v$
+ *
+ * Source and destination must
+ * not be the same vector.
*/
- void swap_col (const unsigned int i,
- const unsigned int j);
-
+ template<typename number2>
+ void vmult (Vector<number2> &w,
+ const Vector<number2> &v,
+ const bool adding=false) const;
+
/**
- * A(i,i)+=B(i,1...n). Addition of complete
- * rows of B to diagonal-elements of this ; <p>
- * ( i = 1 ... m )
+ * Transpose
+ * matrix-vector-multiplication.
+ * See vmult() above.
+ *
+ * Source and destination must
+ * not be the same vector.
*/
template<typename number2>
- void add_diag (const number s,
- const FullMatrix<number2> &B);
+ void Tvmult (Vector<number2> &w,
+ const Vector<number2> &v,
+ const bool adding=false) const;
/**
- * Add constant to diagonal
- * elements of this, i.e. add a
- * multiple of the identity
- * matrix.
+ * Apply the Jacobi
+ * preconditioner, which
+ * multiplies every element of
+ * the <tt>src</tt> vector by the
+ * inverse of the respective
+ * diagonal element and
+ * multiplies the result with the
+ * damping factor <tt>omega</tt>.
*/
- void diagadd (const number s);
-
+ template <typename somenumber>
+ void precondition_Jacobi (Vector<somenumber> &dst,
+ const Vector<somenumber> &src,
+ const number omega = 1.) const;
+
/**
* <i>dst=b-A*x</i>. Residual calculation,
* returns the <i>l<sub>2</sub></i>-norm
void backward (Vector<number2> &dst,
const Vector<number2> &src) const;
- /**
- * QR-factorization of a matrix.
- * The orthogonal transformation
- * Q is applied to the vector y
- * and this matrix.
- *
- * After execution of
- * householder, the upper
- * triangle contains the
- * resulting matrix R, the lower
- * the incomplete factorization
- * matrices.
- */
- template<typename number2>
- void householder (Vector<number2> &y);
-
/**
* Least-Squares-Approximation by
* QR-factorization. The return
template<typename number2>
double least_squares (Vector<number2> &dst,
Vector<number2> &src);
+//@}
- /**
- * Output of the matrix in
- * user-defined format.
- */
- void print (std::ostream &s,
- const unsigned int width=5,
- const unsigned int precision=2) const;
-
- /**
- * Print the matrix in the usual
- * format, i.e. as a matrix and
- * not as a list of nonzero
- * elements. For better
- * readability, elements not in
- * the matrix are displayed as
- * empty space, while matrix
- * elements which are explicitly
- * set to zero are displayed as
- * such.
- *
- * The parameters allow for a
- * flexible setting of the output
- * format: <tt>precision</tt> and
- * <tt>scientific</tt> are used to
- * determine the number format,
- * where <tt>scientific</tt> = <tt>false</tt>
- * means fixed point notation. A
- * zero entry for <tt>width</tt> makes
- * the function compute a width,
- * but it may be changed to a
- * positive value, if output is
- * crude.
- *
- * Additionally, a character for
- * an empty value may be
- * specified.
- *
- * Finally, the whole matrix can
- * be multiplied with a common
- * denominator to produce more
- * readable output, even
- * integers.
- *
- * @attention This function
- * may produce <b>large</b> amounts of
- * output if applied to a large matrix!
- */
- void print_formatted (std::ostream &out,
- const unsigned int presicion=3,
- const bool scientific = true,
- const unsigned int width = 0,
- const char *zero_string = " ",
- const double denominator = 1.) const;
-
- /**
- * Determine an estimate for the
- * memory consumption (in bytes)
- * of this object.
- */
- unsigned int memory_consumption () const;
-
/**
* Exception
*/
/**
* Exception
*/
- DeclException2 (ExcDimensionMismatch,
- int, int,
- << "The two dimensions " << arg1 << " and " << arg2
- << " do not match here.");
- /**
- * Exception
- */
DeclException0 (ExcNotQuadratic);
/**
* Exception
projects / dealii-svn.git / commitdiff