* 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 <typename number>
class LAPACKFullMatrix : public TransposeTable<number>
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.
*
*/
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
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 ();
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
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.
* 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,
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.
*/
<< "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.
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,
#endif
+
+/// Template wrapper for potrf
+template <typename number1>
+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 <typename number1>
+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 <typename number1>
+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 <typename number>
+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 <typename number>
+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 <typename number1>
inline void
LAPACKFullMatrix<number>::LAPACKFullMatrix (const size_type n)
:
TransposeTable<number> (n,n),
- state (matrix)
+ state (matrix),
+ property(general)
{}
const size_type n)
:
TransposeTable<number> (m, n),
- state (matrix)
+ state (matrix),
+ property(general)
{}
LAPACKFullMatrix<number>::LAPACKFullMatrix (const LAPACKFullMatrix &M)
:
TransposeTable<number> (M),
- state (matrix)
+ state (matrix),
+ property(general)
{}
LAPACKFullMatrix<number>::operator = (const LAPACKFullMatrix<number> &M)
{
TransposeTable<number>::operator=(M);
- state = LAPACKSupport::matrix;
+ state = M.state;
+ property = M.property;
return *this;
}
(*this)(i,j) = M(i,j);
state = LAPACKSupport::matrix;
+ property = LAPACKSupport::general;
return *this;
}
(*this)(i,j) = M.el(i,j);
state = LAPACKSupport::matrix;
+ property = LAPACKSupport::general;
return *this;
}
}
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::set_property(const Property p)
+{
+ property = p;
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::l1_norm() const
+{
+ const char type('O');
+ return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::linfty_norm() const
+{
+ const char type('I');
+ return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::frobenius_norm() const
+{
+ const char type('F');
+ return norm(type);
+}
+
+
+
+template <typename number>
+number LAPACKFullMatrix<number>::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<number> 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 <typename number>
+void
+LAPACKFullMatrix<number>::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 <typename number>
+number
+LAPACKFullMatrix<number>::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<number> w(3*N);
+ std::vector<int> 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 <typename number>
void
LAPACKFullMatrix<number>::compute_svd()
void
LAPACKFullMatrix<number>::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<number *> (&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());
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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <typename FullMatrix>
+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 <typename FullMatrix>
+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;
+}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+ // Full matrix:
+ FullMatrix<NumberType> F(size), C(size);
+ create_spd(F);
+ C.cholesky(F);
+
+ // Lapack:
+ LAPACKFullMatrix<NumberType> 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<NumberType> 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<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+ for (const auto &s : sizes)
+ {
+ // test<float>(s);
+ test<double>(s);
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+ // Full matrix:
+ FullMatrix<NumberType> F(size);
+ create_spd(F);
+
+ // Lapack:
+ LAPACKFullMatrix<NumberType> 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<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+ for (const auto &s : sizes)
+ {
+ deallog << "size=" << s << std::endl;
+ // test<float>(s);
+ test<double>(s);
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+ // Full matrix:
+ FullMatrix<NumberType> 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<NumberType> 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<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+ for (const auto &s : sizes)
+ {
+ deallog << "size=" << s << std::endl;
+ // test<float>(s);
+ test<double>(s);
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+ // Full matrix:
+ FullMatrix<NumberType> F(size), invF(size);
+ create_spd(F);
+ invF.invert(F);
+
+ // Lapack:
+ LAPACKFullMatrix<NumberType> M(size);
+ M = F;
+ M.set_property(LAPACKSupport::symmetric);
+ M.compute_cholesky_factorization();
+ M.invert();
+
+ FullMatrix<NumberType> 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<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+ for (const auto &s : sizes)
+ {
+ deallog << "size=" << s << std::endl;
+ // test<float>(s);
+ test<double>(s);
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/lac/lapack_full_matrix.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <iostream>
+
+
+template <typename NumberType>
+void
+test(const unsigned int size)
+{
+ // Full matrix:
+ FullMatrix<NumberType> F(size);
+ create_random(F);
+
+ // Lapack:
+ LAPACKFullMatrix<NumberType> 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<unsigned int> sizes = {{1,3,11,17,32,64,200,391}};
+ for (const auto &s : sizes)
+ {
+ deallog << "size=" << s << std::endl;
+ // test<float>(s);
+ test<double>(s);
+ }
+}
--- /dev/null
+
+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