// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* well-behaved for positive
* definite matrices, but be
* aware of round-off errors in
- * the indefinite case.
+ * the indefinite case.
+ *
+ * In case deal.II was configured with
+ * LAPACK, the functions Xgetrf and
+ * Xgetri build an LU factorization and
+ * invert the matrix upon that
+ * factorization, providing best
+ * performance up to matrices with a
+ * few hundreds rows and columns.
*
* The numerical effort to invert
* an <tt>n x n</tt> matrix is of the
//@}
friend class Accessor;
+
+ private:
+
+#ifdef HAVE_LIBLAPACK
+ /**
+ * The vector storing the
+ * permutations applied for
+ * pivoting in the
+ * LU-factorization.
+ */
+ std::vector<int> ipiv;
+
+ /**
+ * Workspace for calculating the
+ * inverse matrix from an LU
+ * factorization.
+ */
+ std::vector<number> inv_work;
+#endif
+
};
/**@}*/
// $Id$
// Version: $Name$
//
-// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 by the deal.II authors
+// Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
* inverse matrix from an LU
* factorization.
*/
- std::vector<number> iwork;
+ std::vector<number> inv_work;
/**
* Real parts of
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);
+ int* ipiv, double* inv_work, const int* lwork, int* info);
void sgetri_ (const int* n, float* A, const int* lda,
- int* ipiv, float* iwork, const int* lwork, int* info);
+ int* ipiv, float* inv_work, 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_
+#ifdef HAVE_DGETRF_
inline void
-getri (const int* n, double* A, const int* lda, int* ipiv, double* iwork, const int* lwork, int* info)
+getri (const int* n, double* A, const int* lda, int* ipiv, double* inv_work, const int* lwork, int* info)
{
- dgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+ dgetri_ (n,A,lda,ipiv,inv_work,lwork,info);
}
#else
inline void
getri (const int*, double*, const int*, int*, double*, const int*, int*)
{
- LAPACKSupport::ExcMissing("dgetrf");
+ LAPACKSupport::ExcMissing("dgetri");
}
#endif
-#ifdef HAVE_SGETRI_
+#ifdef HAVE_SGETRF_
inline void
-getri (const int* n, float* A, const int* lda, int* ipiv, float* iwork, const int* lwork, int* info)
+getri (const int* n, float* A, const int* lda, int* ipiv, float* inv_work, const int* lwork, int* info)
{
- sgetri_ (n,A,lda,ipiv,iwork,lwork,info);
+ sgetri_ (n,A,lda,ipiv,inv_work,lwork,info);
}
#else
inline void
getri (const int*, float*, const int*, int*, float*, const int*, int*)
{
- LAPACKSupport::ExcMissing("sgetrf");
+ LAPACKSupport::ExcMissing("sgetri");
}
#endif
:
BlockVectorBase<Vector > ()
{
- reinit(v);
+ this->components.resize (v.n_blocks());
+ this->block_indices = v.block_indices;
+
+ for (unsigned int i=0; i<this->n_blocks(); ++i)
+ this->components[i] = v.components[i];
}
#include <lac/full_matrix.templates.h>
+#include <lac/lapack_templates.h>
#include <base/logstream.h>
DEAL_II_NAMESPACE_OPEN
+
+ // Need to explicitly state the Lapack
+ // inversion since it only works with
+ // floats and doubles in case LAPACK was
+ // detected by configure.
+#ifdef HAVE_LIBLAPACK
+
+template <>
+void
+FullMatrix<float>::gauss_jordan ()
+{
+ Assert (!this->empty(), ExcEmptyMatrix());
+ Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
+
+ // In case we have the LAPACK functions
+ // getrf and getri detected at configure,
+ // we use these algorithms for inversion
+ // since they provide better performance
+ // than the deal.II native functions.
+ //
+ // Note that BLAS/LAPACK stores matrix
+ // elements column-wise (i.e., all values in
+ // one column, then all in the next, etc.),
+ // whereas the FullMatrix stores them
+ // row-wise.
+ // We ignore that difference, and give our
+ // row-wise data to LAPACK,
+ // let LAPACK build the inverse of the
+ // transpose matrix, and read the result as
+ // if it were row-wise again. In other words,
+ // we just got ((A^T)^{-1})^T, which is
+ // A^{-1}.
+
+ const int nn = this->n();
+ float* values = const_cast<float*> (this->data());
+ ipiv.resize(nn);
+ int info;
+
+ // Use the LAPACK function getrf for
+ // calculating the LU factorization.
+ getrf(&nn, &nn, values, &nn, &ipiv[0], &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+
+ inv_work.resize (nn);
+ // Use the LAPACK function getri for
+ // calculating the actual inverse using
+ // the LU factorization.
+ getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+}
+
+template <>
+void
+FullMatrix<double>::gauss_jordan ()
+{
+ Assert (!this->empty(), ExcEmptyMatrix());
+ Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
+
+ // In case we have the LAPACK functions
+ // getrf and getri detected at configure,
+ // we use these algorithms for inversion
+ // since they provide better performance
+ // than the deal.II native functions.
+ //
+ // Note that BLAS/LAPACK stores matrix
+ // elements column-wise (i.e., all values in
+ // one column, then all in the next, etc.),
+ // whereas the FullMatrix stores them
+ // row-wise.
+ // We ignore that difference, and give our
+ // row-wise data to LAPACK,
+ // let LAPACK build the inverse of the
+ // transpose matrix, and read the result as
+ // if it were row-wise again. In other words,
+ // we just got ((A^T)^{-1})^T, which is
+ // A^{-1}.
+
+ const int nn = this->n();
+ double* values = const_cast<double*> (this->data());
+ ipiv.resize(nn);
+ int info;
+
+ // Use the LAPACK function getrf for
+ // calculating the LU factorization.
+ getrf(&nn, &nn, values, &nn, &ipiv[0], &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+
+ inv_work.resize (nn);
+ // Use the LAPACK function getri for
+ // calculating the actual inverse using
+ // the LU factorization.
+ std::cout << " blas";
+ getri(&nn, values, &nn, &ipiv[0], &inv_work[0], &nn, &info);
+ std::cout << " 2 mal" << std::endl;
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+}
+
+ // ... and now the usual instantiations
+ // of gauss_jordan() and all the rest.
+template void FullMatrix<long double>::gauss_jordan ();
+template void FullMatrix<std::complex<float> >::gauss_jordan ();
+template void FullMatrix<std::complex<double> >::gauss_jordan ();
+template void FullMatrix<std::complex<long double> >::gauss_jordan ();
+
+#else
+
+template void FullMatrix<float>::gauss_jordan ();
+template void FullMatrix<double>::gauss_jordan ();
+template void FullMatrix<long double>::gauss_jordan ();
+template void FullMatrix<std::complex<float> >::gauss_jordan ();
+template void FullMatrix<std::complex<double> >::gauss_jordan ();
+template void FullMatrix<std::complex<long double> >::gauss_jordan ();
+
+#endif
+
+
#include "full_matrix.inst"
#undef TEMPL_OP_EQ
+
DEAL_II_NAMESPACE_CLOSE
Assert(info == 0, LACExceptions::ExcSingular());
}
- iwork.resize (mm);
- getri(&mm, values, &mm, &ipiv[0], &iwork[0], &mm, &info);
+ inv_work.resize (mm);
+ getri(&mm, values, &mm, &ipiv[0], &inv_work[0], &mm, &info);
Assert(info >= 0, ExcInternalError());
Assert(info == 0, LACExceptions::ExcSingular());
components.resize(n_blocks());
for (unsigned int i=0;i<n_blocks();++i)
- components[i].reinit(v.block(i), true, false);
-
- collect_sizes();
+ components[i] = v.block(i);
}