-for ac_func in dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_
+for ac_func in dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_
do
as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
{ echo "$as_me:$LINENO: checking for $ac_func" >&5
fi
AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_])
-AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_])
+AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetri_ sgetri_ dgetrs_ sgetrs_ dstev_ sstev_])
dnl -------------------------------------------------------------
dnl set include paths of several libraries
<h3>lac</h3>
<ol>
+ <li>
+ <p>
+ New: The class LAPACKFullMatrix can now invert full matrices using
+ the (optimized) LAPACK functions getrf and getri. The speedup over
+ the FullMatrix::gauss_jordan() function is a factor of two for matrices
+ with 100 rows and columns, and grows with matrix size.
+ <br>
+ (Martin Kronbichler 2008/11/11)
+ </p>
+
<li>
<p>
Fixed: The BlockMatrixBase::clear() function that is used by all other
#include <base/table.h>
#include <lac/exceptions.h>
#include <lac/identity_matrix.h>
+#include <lac/lapack_full_matrix.h>
#include <vector>
#include <iomanip>
FullMatrix<number> &
operator = (const IdentityMatrix &id);
+ /**
+ * Assignment operator for a
+ * LapackFullMatrix. The calling matrix
+ * must be of the same size as the
+ * LAPACK matrix.
+ */
+ template <typename number2>
+ FullMatrix<number> &
+ operator = (const LAPACKFullMatrix<number2>&);
+
/**
* Assignment from different
* matrix classes. This
+template <typename number>
+template <typename number2>
+FullMatrix<number> &
+FullMatrix<number>::operator = (const LAPACKFullMatrix<number2>& M)
+{
+ Assert (this->m() == M.n_rows(), ExcDimensionMismatch(this->m(), M.n_rows()));
+ Assert (this->n() == M.n_cols(), ExcDimensionMismatch(this->n(), M.n_rows()));
+ for (unsigned int i=0;i<this->m();++i)
+ for (unsigned int j=0;j<this->n();++j)
+ (*this)(i,j) = M(i,j);
+
+ return *this;
+}
+
+
+
template <typename number>
bool
FullMatrix<number>::all_zero () const
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005, 2006 by the deal.II authors
+// Copyright (C) 2005, 2006, 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
*/
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.
+ */
+ void invert ();
+
/**
* Solve the linear system with
* right hand side given by
*/
std::vector<int> ipiv;
+ /**
+ * Workspace for calculating the
+ * inverse matrix from an LU
+ * factorization.
+ */
+ std::vector<number> iwork;
+
/**
* Real parts of
* eigenvalues. Filled by
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005, 2006 by the deal.II authors
+// Copyright (C) 2005, 2006, 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
{
/// Contents is actually a matrix.
matrix,
+ /// Contents is the inverse of a matrix.
+ inverse_matrix,
/// Contents is an LU decomposition.
lu,
/// Eigenvalue vector is filled
{
case matrix:
return "matrix";
+ case inverse_matrix:
+ return "inverse matrix";
case lu:
return "lu decomposition";
case eigenvalues:
// This file was automatically generated from blas.h.in
// See blastemplates in the deal.II contrib directory
//
-// Copyright (C) 2005, 2006, 2007 by the deal authors
+// Copyright (C) 2005, 2006, 2007, 2008 by the deal authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
const int* lda, int* ipiv, int* info);
void sgetrf_ (const int* m, const int* n, float* A,
const int* lda, int* ipiv, int* info);
+// Invert matrix from LU factorization
+void dgetri_ (const int* n, double* A, const int* lda,
+ int* ipiv, double* iwork, const int* lwork, int* info);
+void sgetri_ (const int* n, float* A, const int* lda,
+ int* ipiv, float* iwork, const int* lwork, int* info);
// Apply forward/backward substitution to LU factorization
void dgetrs_ (const char* trans, const int* n, const int* nrhs,
const double* A, const int* lda, const int* ipiv,
#endif
+#ifdef HAVE_DGETRI_
+inline void
+getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info)
+{
+ dgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+}
+#else
+inline void
+getri (const int*, double*, const int*, int*, double*, const int*, int*)
+{
+ LAPACKSupport::ExcMissing("dgetrf");
+}
+#endif
+
+
+#ifdef HAVE_SGETRI_
+inline void
+getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info)
+{
+ sgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+}
+#else
+inline void
+getri (const int*, float*, const int*, int*, float*, const int*, int*)
+{
+ LAPACKSupport::ExcMissing("sgetrf");
+}
+#endif
+
+
#ifdef HAVE_DGETRS_
inline void
getrs (const char* trans, const int* n, const int* nrhs, const double* A, const int* lda, const int* ipiv, double* b, const int* ldb, int* info)
for (S1, S2 : REAL_SCALARS)
{
+ template
+ FullMatrix<S1>& FullMatrix<S1>::operator = (const LAPACKFullMatrix<S2>&);
+
template
void FullMatrix<S1>::fill<S2> (
const FullMatrix<S2>&, unsigned, unsigned, unsigned, unsigned);
// $Id$
// Version: $Name$
//
-// Copyright (C) 2005, 2006 by the deal.II authors
+// Copyright (C) 2005, 2006, 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
}
+template <typename number>
+void
+LAPACKFullMatrix<number>::invert()
+{
+ Assert(state == matrix || state == lu,
+ ExcState(state));
+ const int mm = this->n_rows();
+ const int nn = this->n_cols();
+ Assert (nn == mm, ExcNotQuadratic());
+
+ number* values = const_cast<number*> (this->data());
+ ipiv.resize(mm);
+ int info;
+
+ if (state == matrix)
+ {
+ getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+ }
+
+ iwork.resize (mm);
+ getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+
+ state = inverse_matrix;
+}
+
+
template <typename number>
void
LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>& v,