From: Denis Davydov Date: Tue, 26 Sep 2017 13:40:07 +0000 (+0200) Subject: extend LAPACKFullMatrix to do Cholesky, inverse of SPD, norms and condition number... X-Git-Tag: v9.0.0-rc1~1015^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ad9ce0f491f32771c0173fd7ae6a8ecad7f577b2;p=dealii.git extend LAPACKFullMatrix to do Cholesky, inverse of SPD, norms and condition number estimate --- diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h index 34d9247560..593794810b 100644 --- a/include/deal.II/lac/lapack_full_matrix.h +++ b/include/deal.II/lac/lapack_full_matrix.h @@ -46,7 +46,7 @@ template class SparseMatrix; * usually the names chosen for the arguments in the LAPACK documentation. * * @ingroup Matrix1 - * @author Guido Kanschat, 2005 + * @author Guido Kanschat, 2005, Denis Davydov, 2017 */ template class LAPACKFullMatrix : public TransposeTable @@ -157,6 +157,11 @@ public: void reinit (const size_type rows, const size_type cols); + /** + * Assign @p property to this matrix. + */ + void set_property(const LAPACKSupport::Property property); + /** * Return the dimension of the codomain (or range) space. * @@ -404,6 +409,23 @@ public: */ void compute_lu_factorization (); + /** + * Compute the Cholesky factorization of the matrix using LAPACK function Xpotrf. + * + * @note The factorization is stored in the lower diagonal. + */ + void compute_cholesky_factorization (); + + /** + * Estimate reciprocal of the condition number (in 1-norm) of a SPD matrix + * using Cholesky factorization. This function can only be called if the matrix + * is already factorized. + * + * @param a_norm Is the 1-norm of the matrix before calling Cholesky + * factorization. + */ + number reciprocal_condition_number(const number a_norm) const; + /** * Compute the determinant of a matrix. As it requires the LU factorization of * the matrix, this function can only be called after @@ -412,8 +434,23 @@ public: number determinant () const; /** - * Invert the matrix by first computing an LU factorization with the LAPACK - * function Xgetrf and then building the actual inverse using Xgetri. + * Compute $L_1$ norm. + */ + number l1_norm() const; + + /** + * Compute $L_\infty$ norm. + */ + number linfty_norm() const; + + /** + * Compute Frobenius norm + */ + number frobenius_norm() const; + + /** + * Invert the matrix by first computing an LU/Cholesky factorization with the LAPACK + * function Xgetrf/Xpotrf and then building the actual inverse using Xgetri/Xpotri. */ void invert (); @@ -613,6 +650,10 @@ public: const double threshold = 0.) const; private: + /** + * Internal function to compute various norms. + */ + number norm(const char type) const; /** * Since LAPACK operations notoriously change the meaning of the matrix @@ -621,10 +662,10 @@ private: LAPACKSupport::State state; /** - * Additional properties of the matrix which may help to select more + * Additional property of the matrix which may help to select more * efficient LAPACK functions. */ - LAPACKSupport::Properties properties; + LAPACKSupport::Property property; /** * The working array used for some LAPACK functions. diff --git a/include/deal.II/lac/lapack_support.h b/include/deal.II/lac/lapack_support.h index 53bbb000a2..6b4032999a 100644 --- a/include/deal.II/lac/lapack_support.h +++ b/include/deal.II/lac/lapack_support.h @@ -85,7 +85,7 @@ namespace LAPACKSupport * A matrix can have certain features allowing for optimization, but hard to * test. These are listed here. */ - enum Properties + enum Property { /// No special properties general = 0, @@ -101,6 +101,28 @@ namespace LAPACKSupport hessenberg = 8 }; + /** + * Function printing the name of a Property. + */ + inline const char *property_name(const Property s) + { + switch (s) + { + case general: + return "general"; + case symmetric: + return "symmetric"; + case upper_triangular: + return "upper triangular"; + case lower_triangular: + return "lower triangular"; + case diagonal: + return "diagonal"; + case hessenberg: + return "Hessenberg"; + } + } + /** * Character constant. */ @@ -149,6 +171,15 @@ namespace LAPACKSupport << "The function cannot be called while the matrix is in state " << state_name(arg1)); + /** + * Exception thrown when a matrix does not have suitable properties for an + * operation. + */ + DeclException1(ExcProperty, Property, + << "The function cannot be called while the " + << property_name(arg1) + << " matrix."); + /** * This exception is thrown if a certain LAPACK function is not available * because no LAPACK installation was detected during configuration. diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h index 8f93986381..c110ff5a88 100644 --- a/include/deal.II/lac/lapack_templates.h +++ b/include/deal.II/lac/lapack_templates.h @@ -65,6 +65,33 @@ extern "C" int *ipiv, double *inv_work, const int *lwork, int *info); void sgetri_ (const int *n, float *A, const int *lda, int *ipiv, float *inv_work, const int *lwork, int *info); +// Compute Cholesky factorization of SPD + void dpotrf_ (const char *uplo, const int *n, double *A, + const int *lda, int *info); + void spotrf_ (const char *uplo, const int *n, float *A, + const int *lda, int *info); +// Estimate the reciprocal of the condition number in 1-norm from Cholesky + void dpocon_ (const char *uplo, const int *n, const double *A, + const int *lda, const double *anorm, double *rcond, + double *work, int *iwork, int *info); + void spocon_ (const char *uplo, const int *n, const float *A, + const int *lda, const float *anorm, float *rcond, + float *work, int *iwork, int *info); +// Computes the inverse from Cholesky + void dpotri_ (const char *uplo, const int *n, double *A, + const int *lda, int *info); + void spotri_ (const char *uplo, const int *n, float *A, + const int *lda, int *info); +// Norms + double dlansy_ (const char *norm, const char *uplo, const int *n, const double *A, + const int *lda, double *work); + float slansy_ (const char *norm, const char *uplo, const int *n, const float *A, + const int *lda, float *work); + double dlange_ (const char *norm, const int *m, const int *n, const double *A, + const int *lda, double *work); + float slange_ (const char *norm, const int *m, const int *n, const float *A, + const int *lda, float *work); + // Compute QR factorization (Householder) void dgeqrf_ (const int *m, const int *n, double *A, const int *lda, double *tau, double *work, @@ -365,6 +392,156 @@ gemm (const char *, const char *, const int *, const int *, const int *, const f #endif + +/// Template wrapper for potrf +template +inline void +potrf (const char *uplo, const int *n, number1 *A, const int *lda, int *info) +{ + Assert (false, ExcNotImplemented()); +} + +inline void +potrf (const char *uplo, const int *n, double *A, const int *lda, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + dpotrf_ (uplo,n,A,lda,info); +#else + Assert (false, LAPACKSupport::ExcMissing("dpotrf")); +#endif +} + +inline void +potrf (const char *uplo, const int *n, float *A, const int *lda, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + spotrf_ (uplo,n,A,lda,info); +#else + Assert (false, LAPACKSupport::ExcMissing("spotrf")); +#endif +} + + + +/// Template wrapper for pocon +template +inline void +pocon (const char *uplo, const int *n, const number1 *A, const int *lda, const number1 *anorm, number1 *rcond, number1 *work, int *iwork, int *info) +{ + Assert (false, ExcNotImplemented()); +} + +inline void +pocon (const char *uplo, const int *n, const double *A, const int *lda, const double *anorm, double *rcond, double *work, int *iwork, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + dpocon_ (uplo,n,A,lda,anorm,rcond,work,iwork,info); +#else + Assert (false, LAPACKSupport::ExcMissing("dpocon")); +#endif +} + +inline void +pocon (const char *uplo, const int *n, const float *A, const int *lda, const float *anorm, float *rcond, float *work, int *iwork, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + spocon_ (uplo,n,A,lda,anorm,rcond,work,iwork,info); +#else + Assert (false, LAPACKSupport::ExcMissing("dpocon")); +#endif +} + +/// Template wrapper for potri +template +inline void +potri(const char *uplo, const int *n, number1 *A, const int *lda, int *info) +{ + Assert (false, ExcNotImplemented()); +} + +inline void +potri(const char *uplo, const int *n, double *A, const int *lda, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + dpotri_(uplo,n,A,lda,info); +#else + Assert (false, LAPACKSupport::ExcMissing("dpotri")); +#endif +} + +inline void +potri(const char *uplo, const int *n, float *A, const int *lda, int *info) +{ +#ifdef DEAL_II_WITH_LAPACK + spotri_(uplo,n,A,lda,info); +#else + Assert (false, LAPACKSupport::ExcMissing("spotri")); +#endif +} + + +/// Template wrapper for lansy +template +inline +number lansy (const char *norm, const char *uplo, const int *n, const number *A, const int *lda, number *work) +{ + Assert (false, ExcNotImplemented()); +} + +inline +double lansy (const char *norm, const char *uplo, const int *n, const double *A, const int *lda, double *work) +{ +#ifdef DEAL_II_WITH_LAPACK + return dlansy_(norm,uplo,n,A,lda,work); +#else + Assert (false, LAPACKSupport::ExcMissing("lansy")); + return 0.; +#endif +} + +inline +float lansy (const char *norm, const char *uplo, const int *n, const float *A, const int *lda, float *work) +{ +#ifdef DEAL_II_WITH_LAPACK + return slansy_(norm,uplo,n,A,lda,work); +#else + Assert (false, LAPACKSupport::ExcMissing("lansy")); + return 0.; +#endif +} + + + +/// Template wrapper for lange +template +inline +number lange (const char *norm, const int *m, const int *n, const number *A, const int *lda, number *work) +{ + Assert (false, ExcNotImplemented()); +} + +inline +double lange (const char *norm, const int *m, const int *n, const double *A, const int *lda, double *work) +{ +#ifdef DEAL_II_WITH_LAPACK + return dlange_(norm,m,n,A,lda,work); +#else + Assert (false, LAPACKSupport::ExcMissing("lange")); + return 0.; +#endif +} + +inline +float lange (const char *norm, const int *m, const int *n, const float *A, const int *lda, float *work) +{ +#ifdef DEAL_II_WITH_LAPACK + return slange_(norm,m,n,A,lda,work); +#else + Assert (false, LAPACKSupport::ExcMissing("lange")); + return 0.; +#endif +} + /// Template wrapper for LAPACK functions dgetrf and sgetrf template inline void diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc index ea9a3b84b6..763daac138 100644 --- a/source/lac/lapack_full_matrix.cc +++ b/source/lac/lapack_full_matrix.cc @@ -33,7 +33,8 @@ template LAPACKFullMatrix::LAPACKFullMatrix (const size_type n) : TransposeTable (n,n), - state (matrix) + state (matrix), + property(general) {} @@ -42,7 +43,8 @@ LAPACKFullMatrix::LAPACKFullMatrix (const size_type m, const size_type n) : TransposeTable (m, n), - state (matrix) + state (matrix), + property(general) {} @@ -50,7 +52,8 @@ template LAPACKFullMatrix::LAPACKFullMatrix (const LAPACKFullMatrix &M) : TransposeTable (M), - state (matrix) + state (matrix), + property(general) {} @@ -59,7 +62,8 @@ LAPACKFullMatrix & LAPACKFullMatrix::operator = (const LAPACKFullMatrix &M) { TransposeTable::operator=(M); - state = LAPACKSupport::matrix; + state = M.state; + property = M.property; return *this; } @@ -95,6 +99,7 @@ LAPACKFullMatrix::operator = (const FullMatrix &M) (*this)(i,j) = M(i,j); state = LAPACKSupport::matrix; + property = LAPACKSupport::general; return *this; } @@ -111,6 +116,7 @@ LAPACKFullMatrix::operator = (const SparseMatrix &M) (*this)(i,j) = M.el(i,j); state = LAPACKSupport::matrix; + property = LAPACKSupport::general; return *this; } @@ -514,6 +520,127 @@ LAPACKFullMatrix::compute_lu_factorization() } + +template +void +LAPACKFullMatrix::set_property(const Property p) +{ + property = p; +} + + + +template +number LAPACKFullMatrix::l1_norm() const +{ + const char type('O'); + return norm(type); +} + + + +template +number LAPACKFullMatrix::linfty_norm() const +{ + const char type('I'); + return norm(type); +} + + + +template +number LAPACKFullMatrix::frobenius_norm() const +{ + const char type('F'); + return norm(type); +} + + + +template +number LAPACKFullMatrix::norm(const char type) const +{ + Assert (state == LAPACKSupport::matrix || + state == LAPACKSupport::inverse_matrix, + ExcMessage("norms can be called in matrix state only.")); + + const int N = this->n_cols(); + const int M = this->n_rows(); + const number *values = &this->values[0]; + std::vector w; + if (property == symmetric) + { + const int lda = std::max(1,N); + const int lwork = (type == 'I' || type == 'O') ? + std::max(1,N) : + 0; + w.resize(lwork); + return lansy (&type, &LAPACKSupport::L, &N, values, &lda, &w[0]); + } + else + { + const int lda = std::max(1,M); + const int lwork = (type == 'I') ? + std::max(1,M) : + 0; + w.resize(lwork); + return lange (&type, &M, &N, values, &lda, &w[0]); + } +} + + +template +void +LAPACKFullMatrix::compute_cholesky_factorization() +{ + Assert(state == matrix, ExcState(state)); + Assert(property == symmetric, ExcProperty(property)); + state = LAPACKSupport::unusable; + + const int mm = this->n_rows(); + const int nn = this->n_cols(); + Assert (mm == nn, ExcDimensionMismatch(mm,nn)); + + number *values = &this->values[0]; + int info = 0; + const int lda = std::max(1,nn); + potrf (&LAPACKSupport::L, &nn, values, &lda, &info); + + // info < 0 : the info-th argument had an illegal value + Assert(info >= 0, ExcInternalError()); + + state = cholesky; + AssertThrow(info == 0, LACExceptions::ExcSingular()); +} + + + +template +number +LAPACKFullMatrix::reciprocal_condition_number(const number a_norm) const +{ + Assert(state == cholesky, ExcState(state)); + number rcond = 0.; + + const int N = this->n_rows(); + const number *values = &this->values[0]; + int info = 0; + const int lda = std::max(1,N); + std::vector w(3*N); + std::vector iw(N); + + // use the same uplo as in Cholesky + pocon (&LAPACKSupport::L, &N, values, &lda, + &a_norm, &rcond, + &w[0], &iw[0], &info); + + Assert(info >= 0, ExcInternalError()); + + return rcond; +} + + + template void LAPACKFullMatrix::compute_svd() @@ -584,26 +711,36 @@ template void LAPACKFullMatrix::invert() { - Assert(state == matrix || state == lu, + Assert(state == matrix || state == lu || state == cholesky, ExcState(state)); const int mm = this->n_rows(); const int nn = this->n_cols(); Assert (nn == mm, ExcNotQuadratic()); number *values = const_cast (&this->values[0]); - ipiv.resize(mm); int info = 0; - if (state == matrix) + if (property != symmetric) { - getrf(&mm, &nn, values, &mm, &ipiv[0], &info); + if (state == matrix) + compute_lu_factorization(); - Assert(info >= 0, ExcInternalError()); - AssertThrow(info == 0, LACExceptions::ExcSingular()); + ipiv.resize(mm); + inv_work.resize (mm); + getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info); + } + else + { + if (state == matrix) + compute_cholesky_factorization(); + + const int lda = std::max(1,nn); + potri(&LAPACKSupport::L, &nn, values,&lda,&info); + // inverse is stored in lower diagonal, set the upper diagonal appropriately: + for (int i=0; i < nn; ++i) + for (int j=i+1; j < nn; ++j) + this->el(i,j) = this->el(j,i); } - - inv_work.resize (mm); - getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info); Assert(info >= 0, ExcInternalError()); AssertThrow(info == 0, LACExceptions::ExcSingular()); @@ -1016,7 +1153,8 @@ LAPACKFullMatrix::print_formatted ( Assert ((!this->empty()) || (this->n_cols()+this->n_rows()==0), ExcInternalError()); Assert (state == LAPACKSupport::matrix || - state == LAPACKSupport::inverse_matrix, + state == LAPACKSupport::inverse_matrix || + state == LAPACKSupport::cholesky, ExcState(state)); // set output format, but store old diff --git a/tests/lapack/create_matrix.h b/tests/lapack/create_matrix.h new file mode 100644 index 0000000000..2e147847f8 --- /dev/null +++ b/tests/lapack/create_matrix.h @@ -0,0 +1,52 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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. +// +// --------------------------------------------------------------------- + +// auxiliary function to create a SPD FullMatrix + +#include "../tests.h" + +template +void create_spd (FullMatrix &A) +{ + const unsigned int size = A.n(); + Assert (size == A.m(), ExcDimensionMismatch(size,A.m())); + + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j = i; j < size; ++j) + { + const double val = (double)rand()/RAND_MAX; + Assert (val >= 0. && val <= 1., + ExcInternalError()); + if (i==j) + // since A(i,j) < 1 and + // a symmetric diagonally dominant matrix is SPD + A(i,j) = val + size; + else + { + A(i,j) = val; + A(j,i) = val; + } + } +} + + + +template +void create_random (FullMatrix &A) +{ + for (unsigned int i = 0; i < A.m(); ++i) + for (unsigned int j = 0; j < A.n(); ++j) + A(i,j) = (double)rand()/RAND_MAX; +} diff --git a/tests/lapack/full_matrix_15.cc b/tests/lapack/full_matrix_15.cc new file mode 100644 index 0000000000..faec1b3193 --- /dev/null +++ b/tests/lapack/full_matrix_15.cc @@ -0,0 +1,78 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::compute_cholesky_factorization by comparing with FullMatrix + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size), C(size); + create_spd(F); + C.cholesky(F); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + M.set_property(LAPACKSupport::symmetric); + M.compute_cholesky_factorization(); + // factorization is stored in the lower diagonal part + for (unsigned int i = 0; i < size; ++i) + for (unsigned int j=i+1; j < size; ++j) + M(i,j) = 0.; + + FullMatrix diff(size); + diff = M; + diff.add(-1., C); + + const NumberType error = diff.frobenius_norm(); + deallog << size << " : " << diff.frobenius_norm() << std::endl; + if (false) + { + std::cout << "Lapack:" << std::endl; + M.print_formatted(std::cout); + std::cout << "FullMatrix:" << std::endl; + C.print_formatted(std::cout); + std::cout << std::flush; + AssertThrow(false, ExcInternalError()); + } +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{1,3,11,17,32,64,200,391}}; + for (const auto &s : sizes) + { + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_15.output b/tests/lapack/full_matrix_15.output new file mode 100644 index 0000000000..5e55595097 --- /dev/null +++ b/tests/lapack/full_matrix_15.output @@ -0,0 +1,9 @@ + +DEAL::1 : 0.00000 +DEAL::3 : 0.00000 +DEAL::11 : 0.00000 +DEAL::17 : 0.00000 +DEAL::32 : 0.00000 +DEAL::64 : 0.00000 +DEAL::200 : 0.00000 +DEAL::391 : 0.00000 diff --git a/tests/lapack/full_matrix_16.cc b/tests/lapack/full_matrix_16.cc new file mode 100644 index 0000000000..a8fcaaebcb --- /dev/null +++ b/tests/lapack/full_matrix_16.cc @@ -0,0 +1,61 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix:: norms + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size); + create_spd(F); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + M.set_property(LAPACKSupport::symmetric); + + deallog << "l1: " << (F.l1_norm() - M.l1_norm()) << std::endl + << "linfty: " << (F.linfty_norm() - M.linfty_norm()) << std::endl + << "frobenius: " << (F.frobenius_norm() - M.frobenius_norm()) << std::endl; +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{1,3,11,17,32,64,200,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_16.output b/tests/lapack/full_matrix_16.output new file mode 100644 index 0000000000..4fc3f3b1fe --- /dev/null +++ b/tests/lapack/full_matrix_16.output @@ -0,0 +1,33 @@ + +DEAL::size=1 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=3 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=11 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=17 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=32 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=64 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=200 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=391 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 diff --git a/tests/lapack/full_matrix_17.cc b/tests/lapack/full_matrix_17.cc new file mode 100644 index 0000000000..f1c9065dcc --- /dev/null +++ b/tests/lapack/full_matrix_17.cc @@ -0,0 +1,66 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::reciprocal_condition_number() by comparing estimated +// value to 1 / ( ||A||_1 * ||A^{-1}||_1 ) + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size), invF(size); + create_spd(F); + invF.invert(F); + const double l1 = F.l1_norm(); + const double inv_l1 = invF.l1_norm(); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + M.set_property(LAPACKSupport::symmetric); + const double la_l1 = M.l1_norm(); + M.compute_cholesky_factorization(); + const double rcond = M.reciprocal_condition_number(la_l1); + + deallog << 1./(l1*inv_l1) << " " << rcond << std::endl; +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{1,3,11,17,32,64,200,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_17.output b/tests/lapack/full_matrix_17.output new file mode 100644 index 0000000000..a0247b26be --- /dev/null +++ b/tests/lapack/full_matrix_17.output @@ -0,0 +1,17 @@ + +DEAL::size=1 +DEAL::1.00000 1.00000 +DEAL::size=3 +DEAL::0.482025 0.563918 +DEAL::size=11 +DEAL::0.425003 0.466425 +DEAL::size=17 +DEAL::0.426863 0.467023 +DEAL::size=32 +DEAL::0.441550 0.493985 +DEAL::size=64 +DEAL::0.452016 0.459357 +DEAL::size=200 +DEAL::0.467209 0.490028 +DEAL::size=391 +DEAL::0.468957 0.486322 diff --git a/tests/lapack/full_matrix_18.cc b/tests/lapack/full_matrix_18.cc new file mode 100644 index 0000000000..45a68bc9e9 --- /dev/null +++ b/tests/lapack/full_matrix_18.cc @@ -0,0 +1,76 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix::invert() wth SPD Cholesky by comparing to FullMatrix + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size), invF(size); + create_spd(F); + invF.invert(F); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + M.set_property(LAPACKSupport::symmetric); + M.compute_cholesky_factorization(); + M.invert(); + + FullMatrix diff(size); + diff = M; + diff.add(-1., invF); + const double error = diff.frobenius_norm(); + deallog << error << std::endl; + + if (error > 1e-10) + { + std::cout << "Lapack:" << std::endl; + M.print_formatted(std::cout); + std::cout << "FullMatrix:" << std::endl; + invF.print_formatted(std::cout); + std::cout << std::flush; + AssertThrow(false, ExcInternalError()); + } +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{1,3,11,17,32,64,200,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_18.output b/tests/lapack/full_matrix_18.output new file mode 100644 index 0000000000..7204301cdc --- /dev/null +++ b/tests/lapack/full_matrix_18.output @@ -0,0 +1,17 @@ + +DEAL::size=1 +DEAL::0.00000 +DEAL::size=3 +DEAL::0.00000 +DEAL::size=11 +DEAL::0.00000 +DEAL::size=17 +DEAL::0.00000 +DEAL::size=32 +DEAL::0.00000 +DEAL::size=64 +DEAL::0.00000 +DEAL::size=200 +DEAL::0.00000 +DEAL::size=391 +DEAL::0.00000 diff --git a/tests/lapack/full_matrix_19.cc b/tests/lapack/full_matrix_19.cc new file mode 100644 index 0000000000..c9df4fc780 --- /dev/null +++ b/tests/lapack/full_matrix_19.cc @@ -0,0 +1,60 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2005 - 2015 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. +// +// --------------------------------------------------------------------- + + +// test LAPACKFullMatrix:: norms for non-symmetric matrices + +#include "../tests.h" +#include "create_matrix.h" +#include +#include +#include + +#include + + +template +void +test(const unsigned int size) +{ + // Full matrix: + FullMatrix F(size); + create_random(F); + + // Lapack: + LAPACKFullMatrix M(size); + M = F; + + deallog << "l1: " << (F.l1_norm() - M.l1_norm()) << std::endl + << "linfty: " << (F.linfty_norm() - M.linfty_norm()) << std::endl + << "frobenius: " << (F.frobenius_norm() - M.frobenius_norm()) << std::endl; +} + + +int main() +{ + const std::string logname = "output"; + std::ofstream logfile(logname.c_str()); + logfile.precision(3); + deallog.attach(logfile); + + const std::vector sizes = {{1,3,11,17,32,64,200,391}}; + for (const auto &s : sizes) + { + deallog << "size=" << s << std::endl; + // test(s); + test(s); + } +} diff --git a/tests/lapack/full_matrix_19.output b/tests/lapack/full_matrix_19.output new file mode 100644 index 0000000000..4fc3f3b1fe --- /dev/null +++ b/tests/lapack/full_matrix_19.output @@ -0,0 +1,33 @@ + +DEAL::size=1 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=3 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=11 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=17 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=32 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=64 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=200 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000 +DEAL::size=391 +DEAL::l1: 0.00000 +DEAL::linfty: 0.00000 +DEAL::frobenius: 0.00000