/* Define to 1 if you have the `dgesvd_' function. */
#undef HAVE_DGESVD_
+/* Define to 1 if you have the `dgetrf_' function. */
+#undef HAVE_DGETRF_
+
+/* Define to 1 if you have the `dgetrs_' function. */
+#undef HAVE_DGETRS_
+
/* Define to 1 if you have the `dstev_' function. */
#undef HAVE_DSTEV_
/* Define to 1 if you have the `sgesvd_' function. */
#undef HAVE_SGESVD_
+/* Define to 1 if you have the `sgetrf_' function. */
+#undef HAVE_SGETRF_
+
+/* Define to 1 if you have the `sgetrs_' function. */
+#undef HAVE_SGETRS_
+
/* Define to 1 if you have the `sstev_' function. */
#undef HAVE_SSTEV_
+for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_
+do
+as_ac_var=`echo "ac_cv_func_$ac_func" | $as_tr_sh`
+echo "$as_me:$LINENO: checking for $ac_func" >&5
+echo $ECHO_N "checking for $ac_func... $ECHO_C" >&6
+if eval "test \"\${$as_ac_var+set}\" = set"; then
+ echo $ECHO_N "(cached) $ECHO_C" >&6
+else
+ cat >conftest.$ac_ext <<_ACEOF
+/* confdefs.h. */
+_ACEOF
+cat confdefs.h >>conftest.$ac_ext
+cat >>conftest.$ac_ext <<_ACEOF
+/* end confdefs.h. */
+/* Define $ac_func to an innocuous variant, in case <limits.h> declares $ac_func.
+ For example, HP-UX 11i <limits.h> declares gettimeofday. */
+#define $ac_func innocuous_$ac_func
+
+/* System header to define __stub macros and hopefully few prototypes,
+ which can conflict with char $ac_func (); below.
+ Prefer <limits.h> to <assert.h> if __STDC__ is defined, since
+ <limits.h> exists even on freestanding compilers. */
+
+#ifdef __STDC__
+# include <limits.h>
+#else
+# include <assert.h>
+#endif
+
+#undef $ac_func
+
+/* Override any gcc2 internal prototype to avoid an error. */
+#ifdef __cplusplus
+extern "C"
+{
+#endif
+/* We use char because int might match the return type of a gcc2
+ builtin and then its argument prototype would still apply. */
+char $ac_func ();
+/* The GNU C library defines this for functions which it implements
+ to always fail with ENOSYS. Some functions are actually named
+ something starting with __ and the normal name is an alias. */
+#if defined (__stub_$ac_func) || defined (__stub___$ac_func)
+choke me
+#else
+char (*f) () = $ac_func;
+#endif
+#ifdef __cplusplus
+}
+#endif
+
+int
+main ()
+{
+return f != $ac_func;
+ ;
+ return 0;
+}
+_ACEOF
+rm -f conftest.$ac_objext conftest$ac_exeext
+if { (eval echo "$as_me:$LINENO: \"$ac_link\"") >&5
+ (eval $ac_link) 2>conftest.er1
+ ac_status=$?
+ grep -v '^ *+' conftest.er1 >conftest.err
+ rm -f conftest.er1
+ cat conftest.err >&5
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); } &&
+ { ac_try='test -z "$ac_cxx_werror_flag"
+ || test ! -s conftest.err'
+ { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+ (eval $ac_try) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); }; } &&
+ { ac_try='test -s conftest$ac_exeext'
+ { (eval echo "$as_me:$LINENO: \"$ac_try\"") >&5
+ (eval $ac_try) 2>&5
+ ac_status=$?
+ echo "$as_me:$LINENO: \$? = $ac_status" >&5
+ (exit $ac_status); }; }; then
+ eval "$as_ac_var=yes"
+else
+ echo "$as_me: failed program was:" >&5
+sed 's/^/| /' conftest.$ac_ext >&5
+
+eval "$as_ac_var=no"
+fi
+rm -f conftest.err conftest.$ac_objext \
+ conftest$ac_exeext conftest.$ac_ext
+fi
+echo "$as_me:$LINENO: result: `eval echo '${'$as_ac_var'}'`" >&5
+echo "${ECHO_T}`eval echo '${'$as_ac_var'}'`" >&6
+if test `eval echo '${'$as_ac_var'}'` = yes; then
+ cat >>confdefs.h <<_ACEOF
+#define `echo "HAVE_$ac_func" | $as_tr_cpp` 1
+_ACEOF
+
+fi
+done
+
+
+
+
+
-for ac_func in dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_
+for ac_func in dgetrf_ sgetrf_ 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
argument is given, use -llapack. Default is not to
use -llapack.],
DEAL_II_WITH_LAPACK($withval))
-AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_ dstev_ sstev_])
+AC_CHECK_FUNCS([dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_ dgesvd_ sgesvd_])
+AC_CHECK_FUNCS([dgetrf_ sgetrf_ dgetrs_ sgetrs_ dstev_ sstev_])
dnl -------------------------------------------------------------
dnl set include paths of several libraries
const double* alpha, const double* A, const int* lda,
const double* x, const int* incx,
const double* b, double* y, const int* incy);
+
+// Compute LU factorization
+void dgetrf_ (const int* m, const int* n, double* A,
+ const int* lda, int* ipiv, 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,
+ double* b, const int* ldb, int* info);
+
// Compute eigenvalues and vectors
void dgeev_ (const char* jobvl, const char* jobvr,
const int* n, double* A, const int* lda,
// $Id$
// Version: $Name$
//
-// Copyright (C) 2004, 2005 by the deal.II authors
+// Copyright (C) 2004, 2005, 2006 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
*/
DeclException0 (ExcNotQuadratic);
+ /**
+ * The operation cannot be finished
+ * since the matrix is singular.
+ */
+ DeclException0 (ExcSingular);
+
/**
* Block indices of two block
* objects are different.
*/
- DeclException0(ExcDifferentBlockIndices);
+ DeclException0 (ExcDifferentBlockIndices);
/**
* An error of a PETSc function was
#include <base/config.h>
+#include <base/smartpointer.h>
#include <base/table.h>
#include <lac/lapack_support.h>
+#include <lac/vector_memory.h>
#include <vector>
#include <complex>
// forward declarations
template<typename number> class Vector;
+template<typename number> class BlockVector;
template<typename number> class FullMatrix;
-
-/*! @addtogroup Matrix1
- *@{
- */
+template<class VECTOR> class VectorMemory;
/**
- * A variant of FullMatrix using LAPACK functions whereever possible.
+ * A variant of FullMatrix using LAPACK functions whereever
+ * possible. In order to do this, the matrix is stored in transposed
+ * order. The element access functions hide this fact by reverting the
+ * transposition.
*
+ * @note In order to perform LAPACK functions, the class contains a lot of
+ * auxiliary data in the private section. The names of these data
+ * vectors are usually the names chosen for the arguments in the
+ * LAPACK documentation.
+ *
+ * @ingroup Matrix1
* @author Guido Kanschat, 2005
*/
template <typename number>
* maximum size possible,
* determined either by the size
* of <tt>this</tt> or <tt>src</tt>.
+ *
+ * The final two arguments allow
+ * to enter a multiple of the
+ * source or its transpose.
*/
template<class MATRIX>
void fill (const MATRIX &src,
const unsigned int dst_offset_i = 0,
const unsigned int dst_offset_j = 0,
const unsigned int src_offset_i = 0,
- const unsigned int src_offset_j = 0);
+ const unsigned int src_offset_j = 0,
+ const number factor = 1.,
+ const bool transpose = false);
/**
* Matrix-vector-multiplication.
void Tvmult_add (Vector<number> &w,
const Vector<number> &v) const;
+ /**
+ * Compute the LU factorization
+ * of the matrix using LAPACK
+ * function Xgetrf.
+ */
+ void compute_lu_factorization ();
+
+ /**
+ * Solve the linear system with
+ * right hand side given by
+ * applying forward/backward
+ * substitution to the previously
+ * computed LU
+ * factorization. Uses LAPACK
+ * function Xgetrs.
+ */
+ void apply_lu_factorization (Vector<number>& v,
+ const bool transposed) const;
+
/**
* Compute eigenvalues of the
* matrix. After this routine has
* some LAPACK functions.
*/
mutable std::vector<number> work;
+
+ /**
+ * The vector storing the
+ * permutations applied for
+ * pivoting in the
+ * LU-factorization.
+ */
+ std::vector<int> ipiv;
/**
* Real parts of
std::vector<number> vr;
};
-/**@}*/
+
+
+/**
+ * A preconditioner based on the LU-factorization of LAPACKFullMatrix.
+ *
+ * @ingroup Preconditioners
+ * @author Guido Kanschat, 2006
+ */
+template <typename number>
+class LUPrecondition
+ :
+ public Subscriptor
+{
+ public:
+ void initialize(const LAPACKFullMatrix<number>&);
+ void initialize(const LAPACKFullMatrix<number>&,
+ VectorMemory<Vector<number> >&);
+ void vmult(Vector<number>&, const Vector<number>&) const;
+ void Tvmult(Vector<number>&, const Vector<number>&) const;
+ void vmult(BlockVector<number>&,
+ const BlockVector<number>&) const;
+ void Tvmult(BlockVector<number>&,
+ const BlockVector<number>&) const;
+ private:
+// SmartPointer<const LAPACKFullMatrix<number> > matrix;
+// SmartPointer<VectorMemory<Vector<number> > mem;
+};
+
+
template <typename number>
template <class MATRIX>
const unsigned int dst_offset_i,
const unsigned int dst_offset_j,
const unsigned int src_offset_i,
- const unsigned int src_offset_j)
+ const unsigned int src_offset_j,
+ const number factor,
+ const bool transpose)
{
const typename MATRIX::const_iterator end = M.end();
for (typename MATRIX::const_iterator entry = M.begin(src_offset_i);
entry != end; ++entry)
{
- const unsigned int i = entry->row();
- const unsigned int j = entry->column();
-
+ const unsigned int i = transpose ? entry->column() : entry->row();
+ const unsigned int j = transpose ? entry->row() : entry->column();
+
const unsigned int dst_i=dst_offset_i+i-src_offset_i;
const unsigned int dst_j=dst_offset_j+j-src_offset_j;
if (dst_i<this->n_rows() && dst_j<this->n_cols())
- (*this)(dst_i, dst_j) = entry->value();
+ (*this)(dst_i, dst_j) = factor * entry->value();
}
state = LAPACKSupport::matrix;
const float* alpha, const float* A, const int* lda,
const float* x, const int* incx,
const float* b, float* y, const int* incy);
+// Compute LU factorization
+void dgetrf_ (const int* m, const int* n, double* A,
+ const int* lda, int* ipiv, int* info);
+void sgetrf_ (const int* m, const int* n, float* A,
+ const int* lda, int* ipiv, 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,
+ double* b, const int* ldb, int* info);
+void sgetrs_ (const char* trans, const int* n, const int* nrhs,
+ const float* A, const int* lda, const int* ipiv,
+ float* b, const int* ldb, int* info);
// Compute eigenvalues and vectors
void dgeev_ (const char* jobvl, const char* jobvr,
const int* n, double* A, const int* lda,
}
#else
inline void
-gemv (const char*, const int*, const int*, const float*, const float*, const int*, const float*, const int*, const float*, float*, const int*)
+gemv (const char*, const int*, const int*, const double*, const double*, const int*, const double*, const int*, const double*, double*, const int*)
{
LAPACKSupport::ExcMissing("sgemv");
}
#endif
+#ifdef HAVE_DGETRF_
+inline void
+getrf (const int* m, const int* n, double* A, const int* lda, int* ipiv, int* info)
+{
+ dgetrf_ (m,n,A,lda,ipiv,info);
+}
+#else
+inline void
+getrf (const int*, const int*, double*, const int*, int*, int*)
+{
+ LAPACKSupport::ExcMissing("dgetrf");
+}
+#endif
+
+
+#ifdef HAVE_SGETRF_
+inline void
+getrf (const int* m, const int* n, float* A, const int* lda, int* ipiv, int* info)
+{
+ sgetrf_ (m,n,A,lda,ipiv,info);
+}
+#else
+inline void
+getrf (const int*, const int*, double*, const int*, 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)
+{
+ dgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info);
+}
+#else
+inline void
+getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*)
+{
+ LAPACKSupport::ExcMissing("dgetrs");
+}
+#endif
+
+
+#ifdef HAVE_SGETRS_
+inline void
+getrs (const char* trans, const int* n, const int* nrhs, const float* A, const int* lda, const int* ipiv, float* b, const int* ldb, int* info)
+{
+ sgetrs_ (trans,n,nrhs,A,lda,ipiv,b,ldb,info);
+}
+#else
+inline void
+getrs (const char*, const int*, const int*, const double*, const int*, const int*, double*, const int*, int*)
+{
+ LAPACKSupport::ExcMissing("sgetrs");
+}
+#endif
+
+
#ifdef HAVE_DGEEV_
inline void
geev (const char* jobvl, const char* jobvr, const int* n, double* A, const int* lda, double* lambda_re, double* lambda_im, double* vl, const int* ldvl, double* vr, const int* ldva, double* work, const int* lwork, int* info)
}
#else
inline void
-geev (const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, float*, const int*, int*)
+geev (const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, double*, const int*, int*)
{
LAPACKSupport::ExcMissing("sgeev");
}
}
#else
inline void
-geevx (const char*, const char*, const char*, const char*, const int*, float*, const int*, float*, float*, float*, const int*, float*, const int*, int*, int*, float*, float*, float*, float*, float*, const int*, int*, int*)
+geevx (const char*, const char*, const char*, const char*, const int*, double*, const int*, double*, double*, double*, const int*, double*, const int*, int*, int*, double*, double*, double*, double*, double*, const int*, int*, int*)
{
LAPACKSupport::ExcMissing("sgeevx");
}
}
#else
inline void
-gesvd (int*, int*, const int*, const int*, float*, const int*, float*, float*, const int*, float*, const int*, float*, const int*, int*)
+gesvd (int*, int*, const int*, const int*, double*, const int*, double*, double*, const int*, double*, const int*, double*, const int*, int*)
{
LAPACKSupport::ExcMissing("sgesvd");
}
}
#else
inline void
-stev (const char*, const int*, float*, float*, float*, const int*, float*, int*)
+stev (const char*, const int*, double*, double*, double*, const int*, double*, int*)
{
LAPACKSupport::ExcMissing("sstev");
}
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const unsigned int n)
:
- TransposeTable<number> (n,n)
+ TransposeTable<number> (n,n),
+ state(matrix)
{}
const unsigned int m,
const unsigned int n)
:
- TransposeTable<number> (m,n)
+ TransposeTable<number> (m,n),
+ state(matrix)
{}
template <typename number>
LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
:
- TransposeTable<number> (M)
-{
- state = LAPACKSupport::matrix;
-}
+ TransposeTable<number> (M),
+ state(matrix)
+{}
template <typename number>
const Vector<number> &v,
const bool adding) const
{
- Assert (state == matrix, ExcInvalidState());
+ Assert (state == matrix, ExcState(state));
const int mm = this->n_rows();
const int nn = this->n_cols();
const Vector<number> &v,
const bool adding) const
{
- Assert (state == matrix, ExcInvalidState());
+ Assert (state == matrix, ExcState(state));
const int mm = this->n_rows();
const int nn = this->n_cols();
#ifdef HAVE_LIBLAPACK
+template <typename number>
+void
+LAPACKFullMatrix<number>::compute_lu_factorization()
+{
+ Assert(state == matrix, ExcState(state));
+ const int mm = this->n_rows();
+ const int nn = this->n_cols();
+ number* values = const_cast<number*> (this->data());
+ ipiv.resize(mm);
+ int info;
+ getrf(&mm, &nn, values, &mm, &ipiv[0], &info);
+
+ Assert(info >= 0, ExcInternalError());
+ Assert(info == 0, LACExceptions::ExcSingular());
+
+ state = lu;
+}
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>& v,
+ const bool transposed) const
+{
+ Assert(state == lu, ExcState(state));
+ Assert(this->n_rows() == this->n_cols(),
+ LACExceptions::ExcNotQuadratic());
+
+ const char* trans = transposed ? &T : &N;
+ const int nn = this->n_cols();
+ const number* values = this->data();
+ int info;
+
+ getrs(trans, &nn, &one, values, &nn, &ipiv[0],
+ v.begin(), &nn, &info);
+
+ Assert(info == 0, ExcInternalError());
+}
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_eigenvalues(
const bool right,
const bool left)
{
+ Assert(state == matrix, ExcState(state));
const int nn = this->n_cols();
wr.resize(nn);
wi.resize(nn);
#else
+template <typename number>
+LAPACKFullMatrix<number>::compute_lu_factorization()
+{
+Assert(false, ExcNeedsLAPACK());
+
+}
+
+
+template <typename number>
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number>&, bool)
+{
+ Assert(false, ExcNeedsLAPACK());
+}
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_eigenvalues(const bool /*right*/,
#endif
+
// template <typename number>
// LAPACKFullMatrix<number>::()
// {}