<h3>Specific improvements</h3>
<ol>
+ <li>Added: The class LAPACKFullMatrix now implements interfaces to
+ matrix-matrix multiplication. Also, LAPACKFullMatrix::apply_lu_factorization
+ now also operates on multiple right hand sides in form of another
+ LAPACKFullMatrix.
+ <br>
+ (Martin Kronbichler, 2014/02/03)
+
<li>Added: A sanity check for the full link interface at configure time.
Hopefully this prevents some people from compiling the whole library just
to hit a link error.
const number alpha = 1.;
const number beta = (adding == true) ? 1. : 0.;
- // Use the BLAS function gemm for
- // calculating the matrix-matrix
+ // Use the BLAS function gemm for calculating the matrix-matrix
// product.
gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
&this->values[0], &k, &beta, &dst(0,0), &m);
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2005 - 2013 by the deal.II authors
+// Copyright (C) 2005 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
typedef types::global_dof_index size_type;
/**
- * Constructor. Initialize the
- * matrix as a square matrix with
- * dimension <tt>n</tt>.
+ * Constructor. Initialize the matrix as a square matrix with dimension
+ * <tt>n</tt>.
*
- * In order to avoid the implicit
- * conversion of integers and
- * other types to a matrix, this
- * constructor is declared
- * <tt>explicit</tt>.
+ * In order to avoid the implicit conversion of integers and other types to
+ * a matrix, this constructor is declared <tt>explicit</tt>.
*
- * By default, no memory is
- * allocated.
+ * By default, no memory is allocated.
*/
explicit LAPACKFullMatrix (const size_type n = 0);
/**
- * Constructor. Initialize the
- * matrix as a rectangular
- * matrix.
+ * Constructor. Initialize the matrix as a rectangular matrix.
*/
LAPACKFullMatrix (const size_type rows,
const size_type cols);
/**
- * Copy constructor. This
- * constructor does a deep copy
- * of the matrix. Therefore, it
- * poses a possible efficiency
- * problem, if for example,
- * function arguments are passed
- * by value rather than by
- * reference. Unfortunately, we
- * can't mark this copy
- * constructor <tt>explicit</tt>,
- * since that prevents the use of
- * this class in containers, such
- * as <tt>std::vector</tt>. The
- * responsibility to check
- * performance of programs must
- * therefore remain with the
- * user of this class.
+ * Copy constructor. This constructor does a deep copy of the
+ * matrix. Therefore, it poses a possible efficiency problem, if for
+ * example, function arguments are passed by value rather than by
+ * reference. Unfortunately, we can't mark this copy constructor
+ * <tt>explicit</tt>, since that prevents the use of this class in
+ * containers, such as <tt>std::vector</tt>. The responsibility to check
+ * performance of programs must therefore remain with the user of this
+ * class.
*/
LAPACKFullMatrix (const LAPACKFullMatrix &);
operator = (const LAPACKFullMatrix<number> &);
/**
- * Assignment operator for a
- * regular FullMatrix. Note that
- * since LAPACK expects matrices
- * in transposed order, this
- * transposition is included here.
+ * Assignment operator for a regular FullMatrix. Note that since LAPACK
+ * expects matrices in transposed order, this transposition is included
+ * here.
*/
template <typename number2>
LAPACKFullMatrix<number> &
operator = (const FullMatrix<number2> &);
/**
- * This operator assigns a scalar
- * to a matrix. To avoid
- * confusion with constructors,
- * zero is the only value allowed
- * for <tt>d</tt>
+ * This operator assigns a scalar to a matrix. To avoid confusion with
+ * constructors, zero is the only value allowed for <tt>d</tt>
*/
LAPACKFullMatrix<number> &
operator = (const double d);
/**
- * Assignment from different
- * matrix classes, performing the
- * usual conversion to the
- * transposed format expected by LAPACK. This
- * assignment operator uses
- * iterators of the class
- * MATRIX. Therefore, sparse
- * matrices are possible sources.
+ * Assignment from different matrix classes, performing the usual conversion
+ * to the transposed format expected by LAPACK. This assignment operator
+ * uses iterators of the class MATRIX. Therefore, sparse matrices are
+ * possible sources.
*/
template <class MATRIX>
void copy_from (const MATRIX &);
/**
* Fill rectangular block.
*
- * A rectangular block of the
- * matrix <tt>src</tt> is copied into
- * <tt>this</tt>. The upper left
- * corner of the block being
- * copied is
- * <tt>(src_offset_i,src_offset_j)</tt>.
- * The upper left corner of the
- * copied block is
- * <tt>(dst_offset_i,dst_offset_j)</tt>.
- * The size of the rectangular
- * block being copied is the
- * maximum size possible,
- * determined either by the size
- * of <tt>this</tt> or <tt>src</tt>.
- *
- * The final two arguments allow
- * to enter a multiple of the
- * source or its transpose.
+ * A rectangular block of the matrix <tt>src</tt> is copied into
+ * <tt>this</tt>. The upper left corner of the block being copied is
+ * <tt>(src_offset_i,src_offset_j)</tt>. The upper left corner of the
+ * copied block is <tt>(dst_offset_i,dst_offset_j)</tt>. The size of the
+ * rectangular block being copied is the maximum size possible, determined
+ * either by the size of <tt>this</tt> or <tt>src</tt>.
+ *
+ * The final two arguments allow to enter a multiple of the source or its
+ * transpose.
*/
template<class MATRIX>
void fill (const MATRIX &src,
/**
* Matrix-vector-multiplication.
*
- * Depending on previous
- * transformations recorded in #state, the
- * result of this function is one
- * of
- *
- * <ul>
- * <li> If #state is LAPACKSupport::matrix or
- * LAPACKSupport::inverse_matrix,
- * this will be a regular matrix
- * vector product using LAPACK
- * gemv().
- * <li> If #state is
- * LAPACKSupport::svd or
- * LAPACKSupport::inverse_svd,
- * this function first multiplies
- * with the right transformation
- * matrix, then with the diagonal
- * matrix of singular values or
- * their reciprocal values, and
- * finally with the left
- * trandformation matrix.
- * </ul>
- *
- * The optional parameter
- * <tt>adding</tt> determines, whether the
- * result is stored in <tt>w</tt> or added
- * to <tt>w</tt>.
+ * Depending on previous transformations recorded in #state, the result of
+ * this function is one of
+ * <ul>
+ * <li> If #state is LAPACKSupport::matrix or LAPACKSupport::inverse_matrix,
+ * this will be a regular matrix vector product using LAPACK gemv().
+ * <li> If #state is LAPACKSupport::svd or LAPACKSupport::inverse_svd, this
+ * function first multiplies with the right transformation matrix, then with
+ * the diagonal matrix of singular values or their reciprocal values, and
+ * finally with the left trandformation matrix.
+ * </ul>
+ *
+ * The optional parameter <tt>adding</tt> determines, whether the result is
+ * stored in <tt>w</tt> or added to <tt>w</tt>.
*
* if (adding)
* <i>w += A*v</i>
* if (!adding)
* <i>w = A*v</i>
*
- * @note Source and destination must
- * not be the same vector.
+ * @note Source and destination must not be the same vector.
*
- * @note This template only
- * exists for compile-time
- * compatibility with
- * FullMatrix. Implementation is
- * only available for <tt>VECTOR=Vector<number></tt>
+ * @note This template only exists for compile-time compatibility with
+ * FullMatrix. Implementation is only available for
+ * <tt>VECTOR=Vector<number></tt>
*/
template <class VECTOR>
void vmult(VECTOR &dst, const VECTOR &src, const bool adding = false) const;
+
/**
- * Adding Matrix-vector-multiplication.
- * <i>w += A*v</i>
+ * Adding Matrix-vector-multiplication. <i>w += A*v</i>
*
- * See the documentation of
- * vmult() for details on the
- * implementation.
+ * See the documentation of vmult() for details on the implementation.
*/
template <class VECTOR>
void vmult_add (VECTOR &w, const VECTOR &v) const;
/**
- * Transpose
- * matrix-vector-multiplication.
+ * Transpose 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>.
+ * The optional parameter <tt>adding</tt> determines, whether the result is
+ * stored in <tt>w</tt> or added to <tt>w</tt>.
*
* if (adding)
* <i>w += A<sup>T</sup>*v</i>
* if (!adding)
* <i>w = A<sup>T</sup>*v</i>
*
- * See the documentation of
- * vmult() for details on the
- * implementation.
+ * See the documentation of vmult() for details on the implementation.
*/
template <class VECTOR>
void Tvmult (VECTOR &w, const VECTOR &v,
const bool adding=false) const;
/**
- * Adding transpose
- * matrix-vector-multiplication.
- * <i>w += A<sup>T</sup>*v</i>
+ * Adding transpose matrix-vector-multiplication. <i>w +=
+ * A<sup>T</sup>*v</i>
*
- * See the documentation of
- * vmult() for details on the
- * implementation.
+ * See the documentation of vmult() for details on the implementation.
*/
template <class VECTOR>
void Tvmult_add (VECTOR &w, const VECTOR &v) const;
const bool adding=false) const;
void Tvmult_add (Vector<number> &w,
const Vector<number> &v) const;
+
+
+ /**
+ * 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)
+ * <i>C += A*B</i>
+ *
+ * if (!adding)
+ * <i>C = A*B</i>
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+ * <tt>C</tt> already has the right size.
+ *
+ * This function uses the BLAS function Xgemm.
+ */
+ void mmult (LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Same as before, but stores the result in a FullMatrix, not in a
+ * LAPACKFullMatrix.
+ */
+ void mmult (FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * 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)
+ * <i>C += A<sup>T</sup>*B</i>
+ *
+ * if (!adding)
+ * <i>C = A<sup>T</sup>*B</i>
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+ * <tt>C</tt> already has the right size.
+ *
+ * This function uses the BLAS function Xgemm.
+ */
+ void Tmmult (LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Same as before, but stores the result in a FullMatrix, not in a
+ * LAPACKFullMatrix.
+ */
+ void Tmmult (FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Matrix-matrix-multiplication using transpose of <tt>B</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)
+ * <i>C += A*B<sup>T</sup></i>
+ *
+ * if (!adding)
+ * <i>C = A*B<sup>T</sup></i>
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+ * <tt>C</tt> already has the right size.
+ *
+ * This function uses the BLAS function Xgemm.
+ */
+ void mTmult (LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
/**
- * Compute the LU factorization
- * of the matrix using LAPACK
- * function Xgetrf.
+ * Same as before, but stores the result in a FullMatrix, not in a
+ * LAPACKFullMatrix.
+ */
+ void mTmult (FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Matrix-matrix-multiplication using transpose of <tt>this</tt> and
+ * <tt>B</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)
+ * <i>C += A<sup>T</sup>*B<sup>T</sup></i>
+ *
+ * if (!adding)
+ * <i>C = A<sup>T</sup>*B<sup>T</sup></i>
+ *
+ * Assumes that <tt>A</tt> and <tt>B</tt> have compatible sizes and that
+ * <tt>C</tt> already has the right size.
+ *
+ * This function uses the BLAS function Xgemm.
+ */
+ void TmTmult (LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Same as before, but stores the result in a FullMatrix, not in a
+ * LAPACKFullMatrix.
+ */
+ void TmTmult (FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding=false) const;
+
+ /**
+ * Compute the LU factorization of the matrix using LAPACK function Xgetrf.
*/
void compute_lu_factorization ();
/**
- * Invert the matrix by first computing
- * an LU factorization with the LAPACK
- * function Xgetrf and then building
- * the actual inverse using Xgetri.
+ * Invert the matrix by first computing an LU factorization with the LAPACK
+ * function Xgetrf and then building the actual inverse using Xgetri.
*/
void invert ();
/**
- * Solve the linear system with
- * right hand side given by
- * applying forward/backward
- * substitution to the previously
- * computed LU
- * factorization. Uses LAPACK
- * function Xgetrs.
+ * Solve the linear system with right hand side given by applying
+ * forward/backward substitution to the previously computed LU
+ * factorization. Uses LAPACK function Xgetrs.
+ *
+ * The flag transposed indicates whether the solution of the transposed
+ * system is to be performed.
*/
void apply_lu_factorization (Vector<number> &v,
const bool transposed) const;
/**
- * Compute eigenvalues of the
- * matrix. After this routine has
- * been called, eigenvalues can
- * be retrieved using the
- * eigenvalue() function. The
- * matrix itself will be
- * LAPACKSupport::unusable after
- * this operation.
- *
- * The optional arguments allow
- * to compute left and right
- * eigenvectors as well.
- *
- * Note that the function does
- * not return the computed
- * eigenvalues right away since
- * that involves copying data
- * around between the output
- * arrays of the LAPACK functions
- * and any return array. This is
- * often unnecessary since one
- * may not be interested in all
- * eigenvalues at once, but for
- * example only the extreme
- * ones. In that case, it is
- * cheaper to just have this
- * function compute the
- * eigenvalues and have a
- * separate function that returns
- * whatever eigenvalue is
- * requested.
- *
- * @note Calls the LAPACK
- * function Xgeev.
+ * Solve the linear system with multiple right hand sides (as many as there
+ * are columns in the matrix b) given by applying forward/backward
+ * substitution to the previously computed LU factorization. Uses LAPACK
+ * function Xgetrs.
+ *
+ * The flag transposed indicates whether the solution of the transposed
+ * system is to be performed.
+ */
+ void apply_lu_factorization (LAPACKFullMatrix<number> &B,
+ const bool transposed) const;
+
+ /**
+ * Compute eigenvalues of the matrix. After this routine has been called,
+ * eigenvalues can be retrieved using the eigenvalue() function. The matrix
+ * itself will be LAPACKSupport::unusable after this operation.
+ *
+ * The optional arguments allow to compute left and right eigenvectors as
+ * well.
+ *
+ * Note that the function does not return the computed eigenvalues right
+ * away since that involves copying data around between the output arrays of
+ * the LAPACK functions and any return array. This is often unnecessary
+ * since one may not be interested in all eigenvalues at once, but for
+ * example only the extreme ones. In that case, it is cheaper to just have
+ * this function compute the eigenvalues and have a separate function that
+ * returns whatever eigenvalue is requested.
+ *
+ * @note Calls the LAPACK function Xgeev.
*/
void compute_eigenvalues (const bool right_eigenvectors = false,
const bool left_eigenvectors = false);
/**
- * Compute eigenvalues and
- * eigenvectors of a real symmetric
- * matrix. Only eigenvalues in the
- * interval (lower_bound, upper_bound]
- * are computed with the absolute
- * tolerance abs_accuracy. An approximate
- * eigenvalue is accepted as converged
- * when it is determined to lie in an
- * interval [a,b] of width less than or
- * equal to abs_accuracy + eps * max( |a|,|b| ),
- * where eps is the machine precision.
- * If abs_accuracy is less than
- * or equal to zero, then eps*|t| will
- * be used in its place, where |t| is the
- * 1-norm of the tridiagonal matrix obtained
- * by reducing A to tridiagonal form.
- * Eigenvalues will be computed most accurately
- * when abs_accuracy is set to twice the
- * underflow threshold, not zero.
- * After this routine has
- * been called, all eigenvalues in
- * (lower_bound, upper_bound] will be
- * stored in eigenvalues and the
- * corresponding eigenvectors will be stored
- * in the columns of eigenvectors, whose
- * dimension is set accordingly.
- *
- * @note Calls the LAPACK function
- * Xsyevx. For this to work, ./configure
- * has to be told to use LAPACK.
+ * Compute eigenvalues and eigenvectors of a real symmetric matrix. Only
+ * eigenvalues in the interval (lower_bound, upper_bound] are computed with
+ * the absolute tolerance abs_accuracy. An approximate eigenvalue is
+ * accepted as converged when it is determined to lie in an interval [a,b]
+ * of width less than or equal to abs_accuracy + eps * max( |a|,|b| ), where
+ * eps is the machine precision. If abs_accuracy is less than or equal to
+ * zero, then eps*|t| will be used in its place, where |t| is the 1-norm of
+ * the tridiagonal matrix obtained by reducing A to tridiagonal form.
+ * Eigenvalues will be computed most accurately when abs_accuracy is set to
+ * twice the underflow threshold, not zero. After this routine has been
+ * called, all eigenvalues in (lower_bound, upper_bound] will be stored in
+ * eigenvalues and the corresponding eigenvectors will be stored in the
+ * columns of eigenvectors, whose dimension is set accordingly.
+ *
+ * @note Calls the LAPACK function Xsyevx. For this to work, ./configure has
+ * to be told to use LAPACK.
*/
void compute_eigenvalues_symmetric(
const number lower_bound,
FullMatrix<number> &eigenvectors);
/**
- * Compute generalized eigenvalues
- * and eigenvectors of
- * a real generalized symmetric
- * eigenproblem of the form
- * itype = 1: $Ax=\lambda B x$
- * itype = 2: $ABx=\lambda x$
- * itype = 3: $BAx=\lambda x$,
- * where A is this matrix.
- * A and B are assumed to be symmetric,
- * and B has to be positive definite.
- * Only eigenvalues in the interval
- * (lower_bound, upper_bound] are
- * computed with the absolute tolerance
- * abs_accuracy.
- * An approximate eigenvalue is accepted
- * as converged when it is determined to
- * lie in an interval [a,b] of width less
- * than or equal to abs_accuracy + eps * max( |a|,|b| ),
- * where eps is the machine precision.
- * If abs_accuracy is less than
- * or equal to zero, then eps*|t| will
- * be used in its place, where |t| is the
- * 1-norm of the tridiagonal matrix obtained
- * by reducing A to tridiagonal form.
- * Eigenvalues will be computed most accurately
- * when abs_accuracy is set to twice the
- * underflow threshold, not zero.
- * After this routine has
- * been called, all eigenvalues in
- * (lower_bound, upper_bound] will be
- * stored in eigenvalues and the
- * corresponding eigenvectors will be stored
- * in eigenvectors, whose dimension is set
- * accordingly.
- *
- * @note Calls the LAPACK
- * function Xsygvx. For this to
- * work, ./configure has to
- * be told to use LAPACK.
+ * Compute generalized eigenvalues and eigenvectors of a real generalized
+ * symmetric eigenproblem of the form itype = 1: $Ax=\lambda B x$ itype = 2:
+ * $ABx=\lambda x$ itype = 3: $BAx=\lambda x$, where A is this matrix. A
+ * and B are assumed to be symmetric, and B has to be positive definite.
+ * Only eigenvalues in the interval (lower_bound, upper_bound] are computed
+ * with the absolute tolerance abs_accuracy. An approximate eigenvalue is
+ * accepted as converged when it is determined to lie in an interval [a,b]
+ * of width less than or equal to abs_accuracy + eps * max( |a|,|b| ), where
+ * eps is the machine precision. If abs_accuracy is less than or equal to
+ * zero, then eps*|t| will be used in its place, where |t| is the 1-norm of
+ * the tridiagonal matrix obtained by reducing A to tridiagonal form.
+ * Eigenvalues will be computed most accurately when abs_accuracy is set to
+ * twice the underflow threshold, not zero. After this routine has been
+ * called, all eigenvalues in (lower_bound, upper_bound] will be stored in
+ * eigenvalues and the corresponding eigenvectors will be stored in
+ * eigenvectors, whose dimension is set accordingly.
+ *
+ * @note Calls the LAPACK function Xsygvx. For this to work, ./configure has
+ * to be told to use LAPACK.
*/
void compute_generalized_eigenvalues_symmetric(
LAPACKFullMatrix<number> &B,
const int itype = 1);
/**
- * Same as the other
- * compute_generalized_eigenvalues_symmetric
- * function except that all
- * eigenvalues are computed
- * and the tolerance is set
- * automatically.
- * Note that this function does
- * not return the computed
- * eigenvalues right away since
- * that involves copying data
- * around between the output
- * arrays of the LAPACK functions
- * and any return array. This is
- * often unnecessary since one
- * may not be interested in all
- * eigenvalues at once, but for
- * example only the extreme
- * ones. In that case, it is
- * cheaper to just have this
- * function compute the
- * eigenvalues and have a
- * separate function that returns
- * whatever eigenvalue is
- * requested. Eigenvalues can
- * be retrieved using the
- * eigenvalue() function.
- * The number of computed
- * eigenvectors is equal
- * to eigenvectors.size()
- *
- * @note Calls the LAPACK
- * function Xsygv. For this to
- * work, ./configure has to
- * be told to use LAPACK.
+ * Same as the other compute_generalized_eigenvalues_symmetric function
+ * except that all eigenvalues are computed and the tolerance is set
+ * automatically. Note that this function does not return the computed
+ * eigenvalues right away since that involves copying data around between
+ * the output arrays of the LAPACK functions and any return array. This is
+ * often unnecessary since one may not be interested in all eigenvalues at
+ * once, but for example only the extreme ones. In that case, it is cheaper
+ * to just have this function compute the eigenvalues and have a separate
+ * function that returns whatever eigenvalue is requested. Eigenvalues can
+ * be retrieved using the eigenvalue() function. The number of computed
+ * eigenvectors is equal to eigenvectors.size()
+ *
+ * @note Calls the LAPACK function Xsygv. For this to work, ./configure has
+ * to be told to use LAPACK.
*/
void compute_generalized_eigenvalues_symmetric (
LAPACKFullMatrix<number> &B,
const int itype = 1);
/**
- * Compute the singular value
- * decomposition of the
- * matrix using LAPACK function
- * Xgesdd.
+ * Compute the singular value decomposition of the matrix using LAPACK
+ * function Xgesdd.
*
- * Requires that the #state is
- * LAPACKSupport::matrix, fills
- * the data members #wr, #svd_u,
- * and #svd_vt, and leaves the
- * object in the #state
+ * Requires that the #state is LAPACKSupport::matrix, fills the data members
+ * #wr, #svd_u, and #svd_vt, and leaves the object in the #state
* LAPACKSupport::svd.
*/
void compute_svd();
+
/**
- * Compute the inverse of the
- * matrix by singular value
- * decomposition.
- *
- * Requires that #state is either
- * LAPACKSupport::matrix or
- * LAPACKSupport::svd. In the
- * first case, this function
- * calls compute_svd(). After
- * this function, the object will
- * have the #state
+ * Compute the inverse of the matrix by singular value decomposition.
+ *
+ * Requires that #state is either LAPACKSupport::matrix or
+ * LAPACKSupport::svd. In the first case, this function calls
+ * compute_svd(). After this function, the object will have the #state
* LAPACKSupport::inverse_svd.
*
- * For a singular value
- * decomposition, the inverse is
- * simply computed by replacing
- * all singular values by their
- * reciprocal values. If the
- * matrix does not have maximal
- * rank, singular values 0 are
- * not touched, thus computing
- * the minimal norm right inverse
- * of the matrix.
- *
- * The parameter
- * <tt>threshold</tt> determines,
- * when a singular value should
- * be considered zero. It is the
- * ratio of the smallest to the
- * largest nonzero singular
- * value
- * <i>s</i><sub>max</sub>. Thus,
- * the inverses of all singular
- * values less than
- * <i>s</i><sub>max</sub>/<tt>threshold</tt>
- * will be set to zero.
+ * For a singular value decomposition, the inverse is simply computed by
+ * replacing all singular values by their reciprocal values. If the matrix
+ * does not have maximal rank, singular values 0 are not touched, thus
+ * computing the minimal norm right inverse of the matrix.
+ *
+ * The parameter <tt>threshold</tt> determines, when a singular value should
+ * be considered zero. It is the ratio of the smallest to the largest
+ * nonzero singular value <i>s</i><sub>max</sub>. Thus, the inverses of all
+ * singular values less than <i>s</i><sub>max</sub>/<tt>threshold</tt> will
+ * be set to zero.
*/
void compute_inverse_svd (const double threshold = 0.);
/**
- * Retrieve eigenvalue after
- * compute_eigenvalues() was
- * called.
+ * Retrieve eigenvalue after compute_eigenvalues() was called.
*/
std::complex<number>
eigenvalue (const size_type i) const;
/**
- * Retrieve singular values after
- * compute_svd() or
- * compute_inverse_svd() was
+ * Retrieve singular values after compute_svd() or compute_inverse_svd() was
* called.
*/
number
singular_value (const size_type i) const;
/**
- * Print the matrix and allow
- * formatting of entries.
- *
- * The parameters allow for a
- * flexible setting of the output
- * format:
- *
- * @arg <tt>precision</tt>
- * denotes the number of trailing
- * digits.
- *
- * @arg <tt>scientific</tt> is
- * used to determine the number
- * format, where
- * <tt>scientific</tt> =
- * <tt>false</tt> means fixed
- * point notation.
- *
- * @arg <tt>width</tt> denotes
- * the with of each column. 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.
- *
- * @arg <tt>zero_string</tt>
- * specifies a string printed for
- * zero entries.
- *
- * @arg <tt>denominator</tt>
- * Multiply the whole matrix by
- * this common denominator to get
- * nicer numbers.
- *
- * @arg <tt>threshold</tt>: all
- * entries with absolute value
- * smaller than this are
- * considered zero.
+ * Print the matrix and allow formatting of entries.
+ *
+ * The parameters allow for a flexible setting of the output format:
+ *
+ * @arg <tt>precision</tt> denotes the number of trailing digits.
+ *
+ * @arg <tt>scientific</tt> is used to determine the number format, where
+ * <tt>scientific</tt> = <tt>false</tt> means fixed point notation.
+ *
+ * @arg <tt>width</tt> denotes the with of each column. 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.
+ *
+ * @arg <tt>zero_string</tt> specifies a string printed for zero entries.
+ *
+ * @arg <tt>denominator</tt> Multiply the whole matrix by this common
+ * denominator to get nicer numbers.
+ *
+ * @arg <tt>threshold</tt>: all entries with absolute value smaller than
+ * this are considered zero.
*/
void print_formatted (std::ostream &out,
const unsigned int presicion=3,
private:
/**
- * Since LAPACK operations
- * notoriously change the meaning
- * of the matrix entries, we
- * record the current state after
- * the last operation here.
+ * Since LAPACK operations notoriously change the meaning of the matrix
+ * entries, we record the current state after the last operation here.
*/
LAPACKSupport::State state;
+
/**
- * Additional properties of the
- * matrix which may help to
- * select more efficient LAPACK
- * functions.
+ * Additional properties of the matrix which may help to select more
+ * efficient LAPACK functions.
*/
LAPACKSupport::Properties properties;
/**
- * The working array used for
- * some LAPACK functions.
+ * The working array used for some LAPACK functions.
*/
mutable std::vector<number> work;
/**
- * The vector storing the
- * permutations applied for
- * pivoting in the
+ * The vector storing the permutations applied for pivoting in the
* LU-factorization.
*
- * Also used as the scratch array
- * IWORK for LAPACK functions
- * needing it.
+ * Also used as the scratch array IWORK for LAPACK functions needing it.
*/
std::vector<int> ipiv;
/**
- * Workspace for calculating the
- * inverse matrix from an LU
- * factorization.
+ * Workspace for calculating the inverse matrix from an LU factorization.
*/
std::vector<number> inv_work;
/**
- * Real parts of eigenvalues or
- * the singular values. Filled by
+ * Real parts of eigenvalues or the singular values. Filled by
* compute_eigenvalues() or compute_svd().
*/
std::vector<number> wr;
/**
- * Imaginary parts of
- * eigenvalues. Filled by
- * compute_eigenvalues.
+ * Imaginary parts of eigenvalues. Filled by compute_eigenvalues.
*/
std::vector<number> wi;
/**
- * Space where left eigenvectors
- * can be stored.
+ * Space where left eigenvectors can be stored.
*/
std::vector<number> vl;
/**
- * Space where right eigenvectors
- * can be stored.
+ * Space where right eigenvectors can be stored.
*/
std::vector<number> vr;
/**
- * The matrix <i>U</i> in the
- * singular value decomposition
+ * The matrix <i>U</i> in the singular value decomposition
* <i>USV<sup>T</sup></i>.
*/
std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_u;
+
/**
- * The matrix
- * <i>V<sup>T</sup></i> in the
- * singular value decomposition
+ * The matrix <i>V<sup>T</sup></i> in the singular value decomposition
* <i>USV<sup>T</sup></i>.
*/
std_cxx1x::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
// ---------------------------------------------------------------------
// $Id$
//
-// Copyright (C) 2005 - 2013 by the deal.II authors
+// Copyright (C) 2005 - 2014 by the deal.II authors
//
// This file is part of the deal.II library.
//
{}
+
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(
const size_type m,
{}
+
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
:
{}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number> &M)
}
+
template <typename number>
template <typename number2>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator = (const FullMatrix<number2> &M)
{
- Assert (this->n_rows() == M.m(), ExcDimensionMismatch(this->n_rows(), M.m()));
+ Assert (this->n_rows() == M.n_rows(), ExcDimensionMismatch(this->n_rows(), M.n_rows()));
Assert (this->n_cols() == M.n(), ExcDimensionMismatch(this->n_cols(), M.n()));
for (size_type i=0; i<this->n_rows(); ++i)
for (size_type j=0; j<this->n_cols(); ++j)
}
+
template <typename number>
LAPACKFullMatrix<number> &
LAPACKFullMatrix<number>::operator = (const double d)
}
+
template <typename number>
void
LAPACKFullMatrix<number>::vmult (
}
+
template <typename number>
void
LAPACKFullMatrix<number>::Tvmult (
}
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert(C.state == matrix, ExcState(state));
+ Assert (this->n_cols() == B.n_rows(), ExcDimensionMismatch(this->n_cols(), B.n_rows()));
+ Assert (C.n_cols() == B.n_cols(), ExcDimensionMismatch(C.n_cols(), B.n_cols()));
+ Assert (C.n_rows() == this->n_rows(), ExcDimensionMismatch(this->n_rows(), C.n_rows()));
+ const int mm = this->n_rows();
+ const int nn = B.n_cols();
+ const int kk = this->n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ gemm("N", "N", &mm, &nn, &kk, &alpha, &this->values[0], &mm, &B.values[0],
+ &kk, &beta, &C.values[0], &mm);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::mmult(FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert (this->n_cols() == B.n_rows(), ExcDimensionMismatch(this->n_cols(), B.n_rows()));
+ Assert (C.n_cols() == B.n_cols(), ExcDimensionMismatch(C.n_cols(), B.n_cols()));
+ Assert (C.n_rows() == this->n_rows(), ExcDimensionMismatch(this->n_rows(), C.n_rows()));
+ const int mm = this->n_rows();
+ const int nn = B.n_cols();
+ const int kk = this->n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ // since FullMatrix stores the matrix in transposed order compared to this
+ // matrix, compute B^T * A^T = (A * B)^T
+ gemm("T", "T", &nn, &mm, &kk, &alpha, &B.values[0], &kk, &this->values[0],
+ &mm, &beta, &C(0,0), &nn);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert(C.state == matrix, ExcState(state));
+ Assert (this->n_rows() == B.n_rows(), ExcDimensionMismatch(this->n_rows(), B.n_rows()));
+ Assert (C.n_cols() == B.n_cols(), ExcDimensionMismatch(C.n_cols(), B.n_cols()));
+ Assert (C.n_rows() == this->n_cols(), ExcDimensionMismatch(this->n_cols(), C.n_rows()));
+ const int mm = this->n_cols();
+ const int nn = B.n_cols();
+ const int kk = B.n_rows();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ gemm("T", "N", &mm, &nn, &kk, &alpha, &this->values[0], &kk, &B.values[0],
+ &kk, &beta, &C.values[0], &mm);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::Tmmult(FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert (this->n_rows() == B.n_rows(), ExcDimensionMismatch(this->n_rows(), B.n_rows()));
+ Assert (C.n_cols() == B.n_cols(), ExcDimensionMismatch(C.n_cols(), B.n_cols()));
+ Assert (C.n_rows() == this->n_cols(), ExcDimensionMismatch(this->n_cols(), C.n_rows()));
+ const int mm = this->n_cols();
+ const int nn = B.n_cols();
+ const int kk = B.n_rows();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ // since FullMatrix stores the matrix in transposed order compared to this
+ // matrix, compute B^T * A = (A^T * B)^T
+ gemm("T", "N", &nn, &mm, &kk, &alpha, &B.values[0], &kk, &this->values[0],
+ &kk, &beta, &C(0,0), &nn);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::mTmult(LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert(C.state == matrix, ExcState(state));
+ Assert (this->n_cols() == B.n_cols(), ExcDimensionMismatch(this->n_cols(), B.n_cols()));
+ Assert (C.n_cols() == B.n_rows(), ExcDimensionMismatch(C.n_cols(), B.n_rows()));
+ Assert (C.n_rows() == this->n_rows(), ExcDimensionMismatch(this->n_rows(), C.n_rows()));
+ const int mm = this->n_rows();
+ const int nn = B.n_rows();
+ const int kk = B.n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ gemm("N", "T", &mm, &nn, &kk, &alpha, &this->values[0], &mm, &B.values[0],
+ &nn, &beta, &C.values[0], &mm);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::mTmult(FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert (this->n_cols() == B.n_cols(), ExcDimensionMismatch(this->n_cols(), B.n_cols()));
+ Assert (C.n_cols() == B.n_rows(), ExcDimensionMismatch(C.n_cols(), B.n_rows()));
+ Assert (C.n_rows() == this->n_rows(), ExcDimensionMismatch(this->n_rows(), C.n_rows()));
+ const int mm = this->n_rows();
+ const int nn = B.n_rows();
+ const int kk = B.n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ // since FullMatrix stores the matrix in transposed order compared to this
+ // matrix, compute B * A^T = (A * B^T)^T
+ gemm("N", "T", &nn, &mm, &kk, &alpha, &B.values[0], &nn, &this->values[0],
+ &mm, &beta, &C(0,0), &nn);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::TmTmult(LAPACKFullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert(C.state == matrix, ExcState(state));
+ Assert (this->n_rows() == B.n_cols(), ExcDimensionMismatch(this->n_rows(), B.n_cols()));
+ Assert (C.n_cols() == B.n_rows(), ExcDimensionMismatch(C.n_cols(), B.n_rows()));
+ Assert (C.n_rows() == this->n_cols(), ExcDimensionMismatch(this->n_cols(), C.n_rows()));
+ const int mm = this->n_cols();
+ const int nn = B.n_rows();
+ const int kk = B.n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ gemm("T", "T", &mm, &nn, &kk, &alpha, &this->values[0], &kk, &B.values[0],
+ &nn, &beta, &C.values[0], &mm);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::TmTmult(FullMatrix<number> &C,
+ const LAPACKFullMatrix<number> &B,
+ const bool adding) const
+{
+ Assert(state == matrix, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert (this->n_rows() == B.n_cols(), ExcDimensionMismatch(this->n_rows(), B.n_cols()));
+ Assert (C.n_cols() == B.n_rows(), ExcDimensionMismatch(C.n_cols(), B.n_rows()));
+ Assert (C.n_rows() == this->n_cols(), ExcDimensionMismatch(this->n_cols(), C.n_rows()));
+ const int mm = this->n_cols();
+ const int nn = B.n_rows();
+ const int kk = B.n_cols();
+ const number alpha = 1.;
+ const number beta = (adding ? 1. : 0.);
+ const number null = 0.;
+
+ // since FullMatrix stores the matrix in transposed order compared to this
+ // matrix, compute B * A = (A^T * B^T)^T
+ gemm("N", "N", &nn, &mm, &kk, &alpha, &B.values[0], &nn, &this->values[0],
+ &kk, &beta, &C(0,0), &nn);
+}
+
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_lu_factorization()
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_svd()
number *mvt = const_cast<number *> (&svd_vt->values[0]);
int info = 0;
- // see comment on this #if
- // below. Another reason to love Petsc
+ // see comment on this #if below. Another reason to love Petsc
#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
- // First determine optimal
- // workspace size
+ // First determine optimal workspace size
work.resize(1);
int lwork = -1;
gesdd(&LAPACKSupport::A, &mm, &nn, values, &mm,
}
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
}
+
template <typename number>
void
LAPACKFullMatrix<number>::invert()
}
+
template <typename number>
void
LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number> &v,
}
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(LAPACKFullMatrix<number> &B,
+ const bool transposed) const
+{
+ Assert(state == lu, ExcState(state));
+ Assert(B.state == matrix, ExcState(state));
+ Assert(this->n_rows() == this->n_cols(), LACExceptions::ExcNotQuadratic());
+
+ const char *trans = transposed ? &T : &N;
+ const int nn = this->n_cols();
+ const int kk = B.n_cols();
+ const number *values = &this->values[0];
+ int info = 0;
+
+ getrs(trans, &nn, &kk, values, &nn, &ipiv[0], &B.values[0], &nn, &info);
+
+ AssertThrow(info == 0, ExcInternalError());
+}
+
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_eigenvalues(
// Optimal workspace query:
- // The LAPACK routine DGEEV requires
- // a sufficient large workspace variable,
+ // The LAPACK routine DGEEV requires a sufficient large workspace variable,
// minimum requirement is
// work.size>=4*nn.
- // However, to improve performance, a
- // somewhat larger workspace may be needed.
+ // However, to improve performance, a somewhat larger workspace may be
+ // needed.
- // SOME implementations of the LAPACK routine
- // provide a workspace query call,
+ // SOME implementations of the LAPACK routine provide a workspace query
+ // call,
// info:=0, lwork:=-1
- // which returns an optimal value for the
- // size of the workspace array
- // (the PETSc 2.3.0 implementation does NOT
- // provide this functionality!).
+ // which returns an optimal value for the size of the workspace array (the
+ // PETSc 2.3.0 implementation does NOT provide this functionality!).
- // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
- // disable the workspace query.
+ // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to disable the workspace
+ // query.
#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
lwork = -1;
work.resize(1);
&wr[0], &wi[0],
&vl[0], &nn, &vr[0], &nn,
&work[0], &lwork, &info);
- // geev returns info=0 on
- // success. Since we only queried
- // the optimal size for work,
- // everything else would not be
- // acceptable.
+ // geev returns info=0 on success. Since we only queried the optimal size
+ // for work, everything else would not be acceptable.
Assert (info == 0, ExcInternalError());
- // Allocate working array according
- // to suggestion.
+ // Allocate working array according to suggestion.
lwork = (int) (work[0]+.1);
#else
lwork = 4*nn; // no query mode
&wr[0], &wi[0],
&vl[0], &nn, &vr[0], &nn,
&work[0], &lwork, &info);
- // Negative return value implies a
- // wrong argument. This should be
- // internal.
+ // Negative return value implies a wrong argument. This should be internal.
Assert (info >=0, ExcInternalError());
//TODO:[GK] What if the QR method fails?
// Optimal workspace query:
- // The LAPACK routine ?SYEVX requires
- // a sufficient large workspace variable,
+ // The LAPACK routine ?SYEVX requires a sufficient large workspace variable,
// minimum requirement is
// work.size>=3*nn-1.
- // However, to improve performance, a
- // somewhat larger workspace may be needed.
+ // However, to improve performance, a somewhat larger workspace may be
+ // needed.
- // SOME implementations of the LAPACK routine
- // provide a workspace query call,
+ // SOME implementations of the LAPACK routine provide a workspace query
+ // call,
// info:=0, lwork:=-1
- // which returns an optimal value for the
- // size of the workspace array
- // (the PETSc 2.3.0 implementation does NOT
- // provide this functionality!).
+ // which returns an optimal value for the size of the workspace array (the
+ // PETSc 2.3.0 implementation does NOT provide this functionality!).
// define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
// disable the workspace query.
&n_eigenpairs, &wr[0], values_eigenvectors,
&nn, &work[0], &lwork, &iwork[0],
&ifail[0], &info);
- // syevx returns info=0 on
- // success. Since we only queried
- // the optimal size for work,
- // everything else would not be
- // acceptable.
+ // syevx returns info=0 on success. Since we only queried the optimal size
+ // for work, everything else would not be acceptable.
Assert (info == 0, ExcInternalError());
- // Allocate working array according
- // to suggestion.
+ // Allocate working array according to suggestion.
lwork = (int) (work[0]+.1);
#else
lwork = 8*nn > 1 ? 8*nn : 1; // no query mode
&nn, &work[0], &lwork, &iwork[0],
&ifail[0], &info);
- // Negative return value implies a
- // wrong argument. This should be
- // internal.
+ // Negative return value implies a wrong argument. This should be internal.
Assert (info >=0, ExcInternalError());
if (info != 0)
std::cerr << "LAPACK error in syevx" << std::endl;
// Optimal workspace query:
- // The LAPACK routine ?SYGVX requires
- // a sufficient large workspace variable,
+ // The LAPACK routine ?SYGVX requires a sufficient large workspace variable,
// minimum requirement is
// work.size>=3*nn-1.
- // However, to improve performance, a
- // somewhat larger workspace may be needed.
+ // However, to improve performance, a somewhat larger workspace may be
+ // needed.
- // SOME implementations of the LAPACK routine
- // provide a workspace query call,
+ // SOME implementations of the LAPACK routine provide a workspace query
+ // call,
// info:=0, lwork:=-1
- // which returns an optimal value for the
- // size of the workspace array
- // (the PETSc 2.3.0 implementation does NOT
- // provide this functionality!).
+ // which returns an optimal value for the size of the workspace array (the
+ // PETSc 2.3.0 implementation does NOT provide this functionality!).
- // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
- // disable the workspace query.
+ // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to disable the workspace
+ // query.
#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
lwork = -1;
work.resize(1);
dummy, dummy, &abs_accuracy, &n_eigenpairs,
&wr[0], values_eigenvectors, &nn, &work[0],
&lwork, &iwork[0], &ifail[0], &info);
- // sygvx returns info=0 on
- // success. Since we only queried
- // the optimal size for work,
- // everything else would not be
- // acceptable.
+ // sygvx returns info=0 on success. Since we only queried the optimal size
+ // for work, everything else would not be acceptable.
Assert (info == 0, ExcInternalError());
- // Allocate working array according
- // to suggestion.
+ // Allocate working array according to suggestion.
lwork = (int) (work[0]+.1);
#else
lwork = 8*nn > 1 ? 8*nn : 1; // no query mode
// resize workspace arrays
work.resize(static_cast<size_type> (lwork));
- // Finally compute the generalized
- // eigenvalues.
+ // Finally compute the generalized eigenvalues.
sygvx (&itype, jobz, range, uplo, &nn, values_A, &nn,
values_B, &nn, &lower_bound, &upper_bound,
dummy, dummy, &abs_accuracy, &n_eigenpairs,
&wr[0], values_eigenvectors, &nn, &work[0],
&lwork, &iwork[0], &ifail[0], &info);
- // Negative return value implies a
- // wrong argument. This should be
- // internal.
+ // Negative return value implies a wrong argument. This should be internal.
Assert (info >=0, ExcInternalError());
if (info != 0)
std::cerr << "LAPACK error in sygvx" << std::endl;
// Optimal workspace query:
- // The LAPACK routine DSYGV requires
- // a sufficient large workspace variable,
+ // The LAPACK routine DSYGV requires a sufficient large workspace variable,
// minimum requirement is
// work.size>=3*nn-1.
- // However, to improve performance, a
- // somewhat larger workspace may be needed.
+ // However, to improve performance, a somewhat larger workspace may be
+ // needed.
- // SOME implementations of the LAPACK routine
- // provide a workspace query call,
+ // SOME implementations of the LAPACK routine provide a workspace query
+ // call,
// info:=0, lwork:=-1
- // which returns an optimal value for the
- // size of the workspace array
- // (the PETSc 2.3.0 implementation does NOT
- // provide this functionality!).
+ // which returns an optimal value for the size of the workspace array (the
+ // PETSc 2.3.0 implementation does NOT provide this functionality!).
- // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to
- // disable the workspace query.
+ // define the DEAL_II_LIBLAPACK_NOQUERYMODE flag to disable the workspace
+ // query.
#ifndef DEAL_II_LIBLAPACK_NOQUERYMODE
lwork = -1;
work.resize(1);
sygv (&itype, jobz, uplo, &nn, values_A, &nn,
values_B, &nn,
&wr[0], &work[0], &lwork, &info);
- // sygv returns info=0 on
- // success. Since we only queried
- // the optimal size for work,
- // everything else would not be
- // acceptable.
+ // sygv returns info=0 on success. Since we only queried the optimal size
+ // for work, everything else would not be acceptable.
Assert (info == 0, ExcInternalError());
- // Allocate working array according
- // to suggestion.
+ // Allocate working array according to suggestion.
lwork = (int) (work[0]+.1);
#else
lwork = 3*nn-1 > 1 ? 3*nn-1 : 1; // no query mode
// resize workspace array
work.resize((size_type) lwork);
- // Finally compute the generalized
- // eigenvalues.
+ // Finally compute the generalized eigenvalues.
sygv (&itype, jobz, uplo, &nn, values_A, &nn,
values_B, &nn,
&wr[0], &work[0], &lwork, &info);
- // Negative return value implies a
- // wrong argument. This should be
- // internal.
+ // Negative return value implies a wrong argument. This should be internal.
Assert (info >=0, ExcInternalError());
if (info != 0)
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::mmult
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+ const unsigned int m=2;
+ const unsigned int n=3;
+ const unsigned int k=4;
+ FullMatrix<double> A(m, k), B(k, n), C(m, n), OC(m,n);
+ LAPACKFullMatrix<double> AL(m, k), BL(k, n), CL(m, n);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<k; ++j)
+ A(i,j) = AL(i,j) = (double)rand()/RAND_MAX;
+ for (unsigned int i=0; i<k; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ B(i,j) = BL(i,j) = (double)rand()/RAND_MAX;
+
+ A.mmult(C, B);
+ AL.mmult(CL, BL);
+ AL.mmult(OC, BL);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ {
+ Assert(std::abs(C(i,j)-CL(i,j)) < 1e-13, ExcInternalError());
+ Assert(std::abs(C(i,j)-OC(i,j)) < 1e-13, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::Tmmult
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+ const unsigned int m=2;
+ const unsigned int n=3;
+ const unsigned int k=4;
+ FullMatrix<double> A(k, m), B(k, n), C(m, n), OC(m,n);
+ LAPACKFullMatrix<double> AL(k,m), BL(k, n), CL(m, n);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<k; ++j)
+ A(j,i) = AL(j,i) = (double)rand()/RAND_MAX;
+ for (unsigned int i=0; i<k; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ B(i,j) = BL(i,j) = (double)rand()/RAND_MAX;
+
+ A.Tmmult(C, B);
+ AL.Tmmult(CL, BL);
+ AL.Tmmult(OC, BL);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ {
+ Assert(std::abs(C(i,j)-CL(i,j)) < 1e-13, ExcInternalError());
+ Assert(std::abs(C(i,j)-OC(i,j)) < 1e-13, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::mTmult
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+ const unsigned int m=2;
+ const unsigned int n=3;
+ const unsigned int k=4;
+ FullMatrix<double> A(m, k), B(n, k), C(m, n), OC(m,n);
+ LAPACKFullMatrix<double> AL(m, k), BL(n, k), CL(m, n);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<k; ++j)
+ A(i,j) = AL(i,j) = (double)rand()/RAND_MAX;
+ for (unsigned int i=0; i<n; ++i)
+ for (unsigned int j=0; j<k; ++j)
+ B(i,j) = BL(i,j) = (double)rand()/RAND_MAX;
+
+ A.mTmult(C, B);
+ AL.mTmult(CL, BL);
+ AL.mTmult(OC, BL);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ {
+ Assert(std::abs(C(i,j)-CL(i,j)) < 1e-13, ExcInternalError());
+ Assert(std::abs(C(i,j)-OC(i,j)) < 1e-13, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::TmTmult
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+ const unsigned int m=2;
+ const unsigned int n=3;
+ const unsigned int k=4;
+ FullMatrix<double> A(k, m), B(n, k), C(m, n), OC(m,n);
+ LAPACKFullMatrix<double> AL(k, m), BL(n, k), CL(m, n);
+ for (unsigned int i=0; i<k; ++i)
+ for (unsigned int j=0; j<m; ++j)
+ A(i,j) = AL(i,j) = (double)rand()/RAND_MAX;
+ for (unsigned int i=0; i<n; ++i)
+ for (unsigned int j=0; j<k; ++j)
+ B(i,j) = BL(i,j) = (double)rand()/RAND_MAX;
+
+ A.TmTmult(C, B);
+ AL.TmTmult(CL, BL);
+ AL.TmTmult(OC, BL);
+ for (unsigned int i=0; i<m; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ {
+ Assert(std::abs(C(i,j)-CL(i,j)) < 1e-13, ExcInternalError());
+ Assert(std::abs(C(i,j)-OC(i,j)) < 1e-13, ExcInternalError());
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+}
--- /dev/null
+
+DEAL::OK
--- /dev/null
+// ---------------------------------------------------------------------
+// $Id$
+//
+// Copyright (C) 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// Tests LAPACKFullMatrix::apply_lu_factorization in two different variants
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+void test()
+{
+ const unsigned int n=11;
+ LAPACKFullMatrix<double> A(n,n);
+ for (unsigned int i=0; i<n; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ A(i,j) = (double)rand()/RAND_MAX;
+ A.compute_lu_factorization();
+ LAPACKFullMatrix<double> rhs_orig(n, 3);
+ for (unsigned int i=0; i<3; ++i)
+ for (unsigned int j=0; j<n; ++j)
+ rhs_orig(j,i) = 2.32 * rand() / RAND_MAX - 0.9923;
+
+ for (unsigned int transpose=0; transpose < 2; ++transpose)
+ {
+ LAPACKFullMatrix<double> rhs(rhs_orig);
+ A.apply_lu_factorization(rhs, transpose);
+ for (unsigned int i=0; i<3; ++i)
+ {
+ Vector<double> check(n);
+ for (unsigned int j=0; j<n; ++j)
+ check(j) = rhs_orig(j, i);
+ A.apply_lu_factorization(check, transpose);
+ for (unsigned int j=0; j<n; ++j)
+ Assert(std::abs(check(j) - rhs(j,i)) < 1e-13, ExcInternalError());
+ }
+ }
+
+ deallog << "OK" << std::endl;
+}
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ test();
+}
--- /dev/null
+
+DEAL::OK
projects / dealii-svn.git / commitdiff