From bdcf62a96350c9125bb2479d0f982f8fa1187fe6 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 3 Feb 2014 14:32:10 +0000 Subject: [PATCH] Introduce mmult and second factorization routine. git-svn-id: https://svn.dealii.org/trunk@32385 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 7 + .../deal.II/lac/full_matrix.templates.h | 3 +- .../include/deal.II/lac/lapack_full_matrix.h | 717 ++++++++---------- deal.II/source/lac/lapack_full_matrix.cc | 391 +++++++--- tests/lapack/full_matrix_08.cc | 68 ++ tests/lapack/full_matrix_08.output | 2 + tests/lapack/full_matrix_09.cc | 68 ++ tests/lapack/full_matrix_09.output | 2 + tests/lapack/full_matrix_10.cc | 68 ++ tests/lapack/full_matrix_10.output | 2 + tests/lapack/full_matrix_11.cc | 68 ++ tests/lapack/full_matrix_11.output | 2 + tests/lapack/full_matrix_12.cc | 71 ++ tests/lapack/full_matrix_12.output | 2 + 14 files changed, 979 insertions(+), 492 deletions(-) create mode 100644 tests/lapack/full_matrix_08.cc create mode 100644 tests/lapack/full_matrix_08.output create mode 100644 tests/lapack/full_matrix_09.cc create mode 100644 tests/lapack/full_matrix_09.output create mode 100644 tests/lapack/full_matrix_10.cc create mode 100644 tests/lapack/full_matrix_10.output create mode 100644 tests/lapack/full_matrix_11.cc create mode 100644 tests/lapack/full_matrix_11.output create mode 100644 tests/lapack/full_matrix_12.cc create mode 100644 tests/lapack/full_matrix_12.output diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index c09702f038..56b56bf4dc 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,13 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/02/03) +
  2. 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. diff --git a/deal.II/include/deal.II/lac/full_matrix.templates.h b/deal.II/include/deal.II/lac/full_matrix.templates.h index 44c3bb8933..2731918d25 100644 --- a/deal.II/include/deal.II/lac/full_matrix.templates.h +++ b/deal.II/include/deal.II/lac/full_matrix.templates.h @@ -575,8 +575,7 @@ void FullMatrix::mmult (FullMatrix &dst, 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); diff --git a/deal.II/include/deal.II/lac/lapack_full_matrix.h b/deal.II/include/deal.II/lac/lapack_full_matrix.h index b6000c48e5..92e0dd09e8 100644 --- a/deal.II/include/deal.II/lac/lapack_full_matrix.h +++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $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. // @@ -60,47 +60,31 @@ public: typedef types::global_dof_index size_type; /** - * Constructor. Initialize the - * matrix as a square matrix with - * dimension n. + * Constructor. Initialize the matrix as a square matrix with dimension + * n. * - * In order to avoid the implicit - * conversion of integers and - * other types to a matrix, this - * constructor is declared - * explicit. + * In order to avoid the implicit conversion of integers and other types to + * a matrix, this constructor is declared explicit. * - * 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 explicit, - * since that prevents the use of - * this class in containers, such - * as std::vector. 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 + * explicit, since that prevents the use of this class in + * containers, such as std::vector. The responsibility to check + * performance of programs must therefore remain with the user of this + * class. */ LAPACKFullMatrix (const LAPACKFullMatrix &); @@ -111,35 +95,26 @@ public: operator = (const LAPACKFullMatrix &); /** - * 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 LAPACKFullMatrix & operator = (const FullMatrix &); /** - * This operator assigns a scalar - * to a matrix. To avoid - * confusion with constructors, - * zero is the only value allowed - * for d + * This operator assigns a scalar to a matrix. To avoid confusion with + * constructors, zero is the only value allowed for d */ LAPACKFullMatrix & 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 void copy_from (const MATRIX &); @@ -147,24 +122,15 @@ public: /** * Fill rectangular block. * - * A rectangular block of the - * matrix src is copied into - * this. The upper left - * corner of the block being - * copied is - * (src_offset_i,src_offset_j). - * The upper left corner of the - * copied block is - * (dst_offset_i,dst_offset_j). - * The size of the rectangular - * block being copied is the - * maximum size possible, - * determined either by the size - * of this or src. - * - * The final two arguments allow - * to enter a multiple of the - * source or its transpose. + * A rectangular block of the matrix src is copied into + * this. The upper left corner of the block being copied is + * (src_offset_i,src_offset_j). The upper left corner of the + * copied block is (dst_offset_i,dst_offset_j). The size of the + * rectangular block being copied is the maximum size possible, determined + * either by the size of this or src. + * + * The final two arguments allow to enter a multiple of the source or its + * transpose. */ template void fill (const MATRIX &src, @@ -179,33 +145,19 @@ public: /** * Matrix-vector-multiplication. * - * Depending on previous - * transformations recorded in #state, the - * result of this function is one - * of - * - *
      - *
    • If #state is LAPACKSupport::matrix or - * LAPACKSupport::inverse_matrix, - * this will be a regular matrix - * vector product using LAPACK - * gemv(). - *
    • 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. - *
    - * - * The optional parameter - * adding determines, whether the - * result is stored in w or added - * to w. + * Depending on previous transformations recorded in #state, the result of + * this function is one of + *
      + *
    • If #state is LAPACKSupport::matrix or LAPACKSupport::inverse_matrix, + * this will be a regular matrix vector product using LAPACK gemv(). + *
    • 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. + *
    + * + * The optional parameter adding determines, whether the result is + * stored in w or added to w. * * if (adding) * w += A*v @@ -213,36 +165,28 @@ public: * if (!adding) * w = A*v * - * @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 VECTOR=Vector<number> + * @note This template only exists for compile-time compatibility with + * FullMatrix. Implementation is only available for + * VECTOR=Vector<number> */ template void vmult(VECTOR &dst, const VECTOR &src, const bool adding = false) const; + /** - * Adding Matrix-vector-multiplication. - * w += A*v + * Adding Matrix-vector-multiplication. w += A*v * - * See the documentation of - * vmult() for details on the - * implementation. + * See the documentation of vmult() for details on the implementation. */ template void vmult_add (VECTOR &w, const VECTOR &v) const; /** - * Transpose - * matrix-vector-multiplication. + * Transpose matrix-vector-multiplication. * - * The optional parameter - * adding determines, whether the - * result is stored in w or added - * to w. + * The optional parameter adding determines, whether the result is + * stored in w or added to w. * * if (adding) * w += AT*v @@ -250,22 +194,17 @@ public: * if (!adding) * w = AT*v * - * See the documentation of - * vmult() for details on the - * implementation. + * See the documentation of vmult() for details on the implementation. */ template void Tvmult (VECTOR &w, const VECTOR &v, const bool adding=false) const; /** - * Adding transpose - * matrix-vector-multiplication. - * w += AT*v + * Adding transpose matrix-vector-multiplication. w += + * AT*v * - * See the documentation of - * vmult() for details on the - * implementation. + * See the documentation of vmult() for details on the implementation. */ template void Tvmult_add (VECTOR &w, const VECTOR &v) const; @@ -280,103 +219,197 @@ public: const bool adding=false) const; void Tvmult_add (Vector &w, const Vector &v) const; + + + /** + * Matrix-matrix-multiplication. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) + * C += A*B + * + * if (!adding) + * C = A*B + * + * Assumes that A and B have compatible sizes and that + * C already has the right size. + * + * This function uses the BLAS function Xgemm. + */ + void mmult (LAPACKFullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Same as before, but stores the result in a FullMatrix, not in a + * LAPACKFullMatrix. + */ + void mmult (FullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of this. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) + * C += AT*B + * + * if (!adding) + * C = AT*B + * + * Assumes that A and B have compatible sizes and that + * C already has the right size. + * + * This function uses the BLAS function Xgemm. + */ + void Tmmult (LAPACKFullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Same as before, but stores the result in a FullMatrix, not in a + * LAPACKFullMatrix. + */ + void Tmmult (FullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of B. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) + * C += A*BT + * + * if (!adding) + * C = A*BT + * + * Assumes that A and B have compatible sizes and that + * C already has the right size. + * + * This function uses the BLAS function Xgemm. + */ + void mTmult (LAPACKFullMatrix &C, + const LAPACKFullMatrix &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 &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Matrix-matrix-multiplication using transpose of this and + * B. + * + * The optional parameter adding determines, whether the result is + * stored in C or added to C. + * + * if (adding) + * C += AT*BT + * + * if (!adding) + * C = AT*BT + * + * Assumes that A and B have compatible sizes and that + * C already has the right size. + * + * This function uses the BLAS function Xgemm. + */ + void TmTmult (LAPACKFullMatrix &C, + const LAPACKFullMatrix &B, + const bool adding=false) const; + + /** + * Same as before, but stores the result in a FullMatrix, not in a + * LAPACKFullMatrix. + */ + void TmTmult (FullMatrix &C, + const LAPACKFullMatrix &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 &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 &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, @@ -386,45 +419,25 @@ public: FullMatrix &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 &B, @@ -436,40 +449,20 @@ public: 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 &B, @@ -477,116 +470,70 @@ public: 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 - * threshold determines, - * when a singular value should - * be considered zero. It is the - * ratio of the smallest to the - * largest nonzero singular - * value - * smax. Thus, - * the inverses of all singular - * values less than - * smax/threshold - * 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 threshold determines, when a singular value should + * be considered zero. It is the ratio of the smallest to the largest + * nonzero singular value smax. Thus, the inverses of all + * singular values less than smax/threshold 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 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 precision - * denotes the number of trailing - * digits. - * - * @arg scientific is - * used to determine the number - * format, where - * scientific = - * false means fixed - * point notation. - * - * @arg width denotes - * the with of each column. A - * zero entry for width - * makes the function compute a - * width, but it may be changed - * to a positive value, if output - * is crude. - * - * @arg zero_string - * specifies a string printed for - * zero entries. - * - * @arg denominator - * Multiply the whole matrix by - * this common denominator to get - * nicer numbers. - * - * @arg threshold: 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 precision denotes the number of trailing digits. + * + * @arg scientific is used to determine the number format, where + * scientific = false means fixed point notation. + * + * @arg width denotes the with of each column. A zero entry for + * width makes the function compute a width, but it may be changed + * to a positive value, if output is crude. + * + * @arg zero_string specifies a string printed for zero entries. + * + * @arg denominator Multiply the whole matrix by this common + * denominator to get nicer numbers. + * + * @arg threshold: all entries with absolute value smaller than + * this are considered zero. */ void print_formatted (std::ostream &out, const unsigned int presicion=3, @@ -598,82 +545,64 @@ public: 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 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 ipiv; /** - * Workspace for calculating the - * inverse matrix from an LU - * factorization. + * Workspace for calculating the inverse matrix from an LU factorization. */ std::vector 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 wr; /** - * Imaginary parts of - * eigenvalues. Filled by - * compute_eigenvalues. + * Imaginary parts of eigenvalues. Filled by compute_eigenvalues. */ std::vector wi; /** - * Space where left eigenvectors - * can be stored. + * Space where left eigenvectors can be stored. */ std::vector vl; /** - * Space where right eigenvectors - * can be stored. + * Space where right eigenvectors can be stored. */ std::vector vr; /** - * The matrix U in the - * singular value decomposition + * The matrix U in the singular value decomposition * USVT. */ std_cxx1x::shared_ptr > svd_u; + /** - * The matrix - * VT in the - * singular value decomposition + * The matrix VT in the singular value decomposition * USVT. */ std_cxx1x::shared_ptr > svd_vt; diff --git a/deal.II/source/lac/lapack_full_matrix.cc b/deal.II/source/lac/lapack_full_matrix.cc index 215f9c07e1..c58afb3518 100644 --- a/deal.II/source/lac/lapack_full_matrix.cc +++ b/deal.II/source/lac/lapack_full_matrix.cc @@ -1,7 +1,7 @@ // --------------------------------------------------------------------- // $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. // @@ -37,6 +37,7 @@ LAPACKFullMatrix::LAPACKFullMatrix(const size_type n) {} + template LAPACKFullMatrix::LAPACKFullMatrix( const size_type m, @@ -47,6 +48,7 @@ LAPACKFullMatrix::LAPACKFullMatrix( {} + template LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) : @@ -55,6 +57,7 @@ LAPACKFullMatrix::LAPACKFullMatrix(const LAPACKFullMatrix &M) {} + template LAPACKFullMatrix & LAPACKFullMatrix::operator = (const LAPACKFullMatrix &M) @@ -65,12 +68,13 @@ LAPACKFullMatrix::operator = (const LAPACKFullMatrix &M) } + template template LAPACKFullMatrix & LAPACKFullMatrix::operator = (const FullMatrix &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; in_rows(); ++i) for (size_type j=0; jn_cols(); ++j) @@ -81,6 +85,7 @@ LAPACKFullMatrix::operator = (const FullMatrix &M) } + template LAPACKFullMatrix & LAPACKFullMatrix::operator = (const double d) @@ -95,6 +100,7 @@ LAPACKFullMatrix::operator = (const double d) } + template void LAPACKFullMatrix::vmult ( @@ -153,6 +159,7 @@ LAPACKFullMatrix::vmult ( } + template void LAPACKFullMatrix::Tvmult ( @@ -213,6 +220,211 @@ LAPACKFullMatrix::Tvmult ( } + +template +void +LAPACKFullMatrix::mmult(LAPACKFullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::mmult(FullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::Tmmult(LAPACKFullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::Tmmult(FullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::mTmult(LAPACKFullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::mTmult(FullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::TmTmult(LAPACKFullMatrix &C, + const LAPACKFullMatrix &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 +void +LAPACKFullMatrix::TmTmult(FullMatrix &C, + const LAPACKFullMatrix &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 void LAPACKFullMatrix::compute_lu_factorization() @@ -232,6 +444,7 @@ LAPACKFullMatrix::compute_lu_factorization() } + template void LAPACKFullMatrix::compute_svd() @@ -252,12 +465,10 @@ LAPACKFullMatrix::compute_svd() number *mvt = const_cast (&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, @@ -284,6 +495,7 @@ LAPACKFullMatrix::compute_svd() } + template void LAPACKFullMatrix::compute_inverse_svd(const double threshold) @@ -305,6 +517,7 @@ LAPACKFullMatrix::compute_inverse_svd(const double threshold) } + template void LAPACKFullMatrix::invert() @@ -337,6 +550,7 @@ LAPACKFullMatrix::invert() } + template void LAPACKFullMatrix::apply_lu_factorization(Vector &v, @@ -358,6 +572,29 @@ LAPACKFullMatrix::apply_lu_factorization(Vector &v, } + +template +void +LAPACKFullMatrix::apply_lu_factorization(LAPACKFullMatrix &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 void LAPACKFullMatrix::compute_eigenvalues( @@ -380,23 +617,20 @@ LAPACKFullMatrix::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); @@ -405,14 +639,10 @@ LAPACKFullMatrix::compute_eigenvalues( &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 @@ -425,9 +655,7 @@ LAPACKFullMatrix::compute_eigenvalues( &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? @@ -470,20 +698,17 @@ LAPACKFullMatrix::compute_eigenvalues_symmetric( // 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. @@ -498,14 +723,10 @@ LAPACKFullMatrix::compute_eigenvalues_symmetric( &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 @@ -522,9 +743,7 @@ LAPACKFullMatrix::compute_eigenvalues_symmetric( &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; @@ -583,23 +802,20 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( // 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); @@ -609,14 +825,10 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( 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 @@ -624,17 +836,14 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric( // resize workspace arrays work.resize(static_cast (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; @@ -686,23 +895,20 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( // 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); @@ -710,14 +916,10 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( 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 @@ -725,14 +927,11 @@ LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric ( // 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) diff --git a/tests/lapack/full_matrix_08.cc b/tests/lapack/full_matrix_08.cc new file mode 100644 index 0000000000..8d2a514269 --- /dev/null +++ b/tests/lapack/full_matrix_08.cc @@ -0,0 +1,68 @@ +// --------------------------------------------------------------------- +// $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 +#include +#include +#include + +#include +#include + + + +void test() +{ + const unsigned int m=2; + const unsigned int n=3; + const unsigned int k=4; + FullMatrix A(m, k), B(k, n), C(m, n), OC(m,n); + LAPACKFullMatrix AL(m, k), BL(k, n), CL(m, n); + for (unsigned int i=0; i +#include +#include +#include + +#include +#include + + + +void test() +{ + const unsigned int m=2; + const unsigned int n=3; + const unsigned int k=4; + FullMatrix A(k, m), B(k, n), C(m, n), OC(m,n); + LAPACKFullMatrix AL(k,m), BL(k, n), CL(m, n); + for (unsigned int i=0; i +#include +#include +#include + +#include +#include + + + +void test() +{ + const unsigned int m=2; + const unsigned int n=3; + const unsigned int k=4; + FullMatrix A(m, k), B(n, k), C(m, n), OC(m,n); + LAPACKFullMatrix AL(m, k), BL(n, k), CL(m, n); + for (unsigned int i=0; i +#include +#include +#include + +#include +#include + + + +void test() +{ + const unsigned int m=2; + const unsigned int n=3; + const unsigned int k=4; + FullMatrix A(k, m), B(n, k), C(m, n), OC(m,n); + LAPACKFullMatrix AL(k, m), BL(n, k), CL(m, n); + for (unsigned int i=0; i +#include +#include + +#include +#include + + + +void test() +{ + const unsigned int n=11; + LAPACKFullMatrix A(n,n); + for (unsigned int i=0; i rhs_orig(n, 3); + for (unsigned int i=0; i<3; ++i) + for (unsigned int j=0; j rhs(rhs_orig); + A.apply_lu_factorization(rhs, transpose); + for (unsigned int i=0; i<3; ++i) + { + Vector check(n); + for (unsigned int j=0; j