*/
void invert ();
+ /**
+ * Solve the linear system with right hand side. The matrix should
+ * be either triangular or LU factorization should be previously computed.
+ *
+ * The flag transposed indicates whether the solution of the transposed
+ * system is to be performed.
+ */
+ void solve(Vector<number> &v,
+ const bool transposed) const;
+
+ /**
+ * Same as above but for multiple right hand sides (as many as there
+ * are columns in the matrix @p B).
+ */
+ void solve(LAPACKFullMatrix<number> &B,
+ const bool transposed) const;
+
/**
* Solve the linear system with right hand side given by applying
* forward/backward substitution to the previously computed LU
*
* The flag transposed indicates whether the solution of the transposed
* system is to be performed.
+ *
+ * @deprecated use solve() instead.
*/
+ DEAL_II_DEPRECATED
void apply_lu_factorization (Vector<number> &v,
const bool transposed) const;
*
* The flag transposed indicates whether the solution of the transposed
* system is to be performed.
+ *
+ * @deprecated use solve() instead.
*/
+ DEAL_II_DEPRECATED
void apply_lu_factorization (LAPACKFullMatrix<number> &B,
const bool transposed) const;
const int *lda, int *info);
void spotrf_ (const char *uplo, const int *n, float *A,
const int *lda, int *info);
+// Apply forward/backward substitution to Cholesky factorization
+ void dpotrs_ (const char *uplo, const int *n, const int *nrhs,
+ const double *A, const int *lda,
+ double *B, const int *ldb,
+ int *info);
+ void spotrs_ (const char *uplo, const int *n, const int *nrhs,
+ const float *A, const int *lda,
+ float *B, const int *ldb,
+ 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,
#endif
+
+/// Template wrapper for LAPACK functions dpotrs and spotrs
+template <typename number>
+inline void
+potrs (const char *, const int *, const int *,
+ const number *, const int *,
+ number *, const int *,
+ int *)
+{
+ Assert (false, ExcNotImplemented());
+}
+
+#ifdef DEAL_II_WITH_LAPACK
+inline void
+potrs (const char *uplo, const int *n, const int *nrhs,
+ const double *A, const int *lda,
+ double *B, const int *ldb,
+ int *info)
+{
+ dpotrs_(uplo,n,nrhs,A,lda,B,ldb,info);
+}
+inline void
+potrs (const char *uplo, const int *n, const int *nrhs,
+ const float *A, const int *lda,
+ float *B, const int *ldb,
+ int *info)
+{
+ spotrs_(uplo,n,nrhs,A,lda,B,ldb,info);
+}
+#else
+inline void
+potrs (const char *, const int *, const int *,
+ const double *, const int *,
+ double *, const int *,
+ int *)
+{
+ Assert (false, LAPACKSupport::ExcMissing("dpotrs"));
+}
+inline void
+potrs (const char *, const int *, const int *,
+ const float *, const int *,
+ flaot *, const int *,
+ int *)
+{
+ Assert (false, LAPACKSupport::ExcMissing("spotrs"));
+}
+#endif
+
+
+
+
/// Template wrapper for LAPACK functions dgetri and sgetri
template <typename number1, typename number2>
inline void
}
+
template <typename number>
void
-LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number> &v,
- const bool transposed) const
+LAPACKFullMatrix<number>::solve(Vector<number> &v,
+ const bool transposed) const
{
- Assert(state == lu, ExcState(state));
Assert(this->n_rows() == this->n_cols(),
LACExceptions::ExcNotQuadratic());
AssertDimension(this->n_rows(), v.size());
-
const char *trans = transposed ? &T : &N;
const int nn = this->n_cols();
const number *values = &this->values[0];
+ const int n_rhs = 1;
int info = 0;
- getrs(trans, &nn, &one, values, &nn, ipiv.data(),
- v.begin(), &nn, &info);
+ if (state == lu)
+ {
+ getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(),
+ v.begin(), &nn, &info);
+ }
+ else if (state == cholesky)
+ {
+ potrs(&LAPACKSupport::L, &nn, &n_rhs,
+ values, &nn, v.begin(), &nn, &info);
+ }
+ else if (property == upper_triangular ||
+ property == lower_triangular)
+ {
+ const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
+
+ const int lda = nn;
+ const int ldb = nn;
+ trtrs (&uplo, trans, "N",
+ &nn, &n_rhs,
+ values, &lda,
+ v.begin(), &ldb, &info);
+ }
+ else
+ {
+ Assert(false,
+ ExcMessage("The matrix has to be either factorized or triangular."));
+ }
Assert(info == 0, ExcInternalError());
}
+
template <typename number>
void
-LAPACKFullMatrix<number>::apply_lu_factorization(LAPACKFullMatrix<number> &B,
- const bool transposed) const
+LAPACKFullMatrix<number>::solve(LAPACKFullMatrix<number> &B,
+ const bool transposed) const
{
- Assert(state == lu, ExcState(state));
- Assert(B.state == matrix, ExcState(state));
- Assert(this->n_rows() == this->n_cols(), LACExceptions::ExcNotQuadratic());
- AssertDimension(this->n_rows(), B.n_rows());
+ Assert(B.state == matrix, ExcState(B.state));
+ Assert(this->n_rows() == this->n_cols(),
+ LACExceptions::ExcNotQuadratic());
+ AssertDimension(this->n_rows(), B.n_rows());
const char *trans = transposed ? &T : &N;
const int nn = this->n_cols();
- const int kk = B.n_cols();
const number *values = &this->values[0];
+ const int n_rhs = B.n_cols();
int info = 0;
- getrs(trans, &nn, &kk, values, &nn, ipiv.data(), &B.values[0], &nn, &info);
+ if (state == lu)
+ {
+ getrs(trans, &nn, &n_rhs, values, &nn, ipiv.data(),
+ &B.values[0], &nn, &info);
+ }
+ else if (state == cholesky)
+ {
+ potrs(&LAPACKSupport::L, &nn, &n_rhs,
+ values, &nn, &B.values[0], &nn, &info);
+ }
+ else if (property == upper_triangular ||
+ property == lower_triangular)
+ {
+ const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
+
+ const int lda = nn;
+ const int ldb = nn;
+ trtrs (&uplo, trans, "N",
+ &nn, &n_rhs,
+ values, &lda,
+ &B.values[0], &ldb, &info);
+ }
+ else
+ {
+ Assert(false,
+ ExcMessage("The matrix has to be either factorized or triangular."));
+ }
Assert(info == 0, ExcInternalError());
}
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(Vector<number> &v,
+ const bool transposed) const
+{
+ solve(v,transposed);
+}
+
+
+
+template <typename number>
+void
+LAPACKFullMatrix<number>::apply_lu_factorization(LAPACKFullMatrix<number> &B,
+ const bool transposed) const
+{
+ solve(B,transposed);
+}
+
+
+
template <typename number>
number
LAPACKFullMatrix<number>::determinant() const
--- /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::solve() for triangular matrices
+
+/* MWE for size=3 in Octave:
+R = [1,2,3; 0, 5, 6; 0, 0, 9]
+x = [2; -7; 1]
+
+> R\x
+ans =
+
+ 4.73333
+ -1.53333
+ 0.11111
+
+> R'\x
+ans =
+
+ 2.00000
+ -2.20000
+ 0.91111
+*/
+
+#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 = 3;
+ LAPACKFullMatrix<NumberType> M(size);
+ M.set_property(LAPACKSupport::upper_triangular);
+
+ M = 0.;
+ unsigned int counter = 1;
+ for (unsigned int i = 0; i < size; ++i)
+ for (unsigned int j = 0; j < size; ++j)
+ {
+ if (j >= i)
+ M(i,j) = counter;
+
+ counter++;
+ }
+
+ Vector<NumberType> x(size), y(size);
+ x[0] = 2;
+ x[1] = -7;
+ x[2] = 1;
+
+ y = x;
+ M.solve(y, false);
+ y.print(deallog.get_file_stream(), 6, false);
+
+ y = x;
+ M.solve(y, true);
+ y.print(deallog.get_file_stream(), 6, false);
+}
+
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+
+ test<double>();
+
+}
--- /dev/null
+
+4.733333 -1.533333 0.111111
+2.000000 -2.200000 0.911111