number reciprocal_condition_number(const number l1_norm) const;
/**
- * Same as above but for triangular matrices. The matrix has to have
- * the LAPACKSupport::Property set, see set_property().
+ * Estimate the reciprocal of the condition number $1/k(A)$ in $L_1$ norm
+ * for triangular matrices. The matrix has to have
+ * the LAPACKSupport::Property set to either LAPACKSupport::Property::upper_triangular
+ * or LAPACKSupport::Property::lower_triangular, see set_property().
*/
number reciprocal_condition_number() const;
void invert ();
/**
- * Solve the linear system with right hand side. The matrix should
- * be either triangular or LU factorization should be previously computed.
+ * Solve the linear system with right hand side @p v and put the solution
+ * back to @p v. The matrix should be either triangular or LU/Cholesky
+ * 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;
+ const bool transposed = false) 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;
+ const bool transposed = false) const;
/**
* Solve the linear system with right hand side given by applying
--- /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 Cholesky factorization
+
+/* MWE for size=3 in Octave:
+R = [10,2,3; 2, 20, 6; 3, 6, 90]
+x = [2; -7; 1]
+
+> R\x
+ans =
+
+ 0.268693
+ -0.385220
+ 0.027836
+*/
+
+#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::symmetric);
+
+ M(0,0) = 10;
+ M(0,1) = 2;
+ M(1,0) = 2;
+ M(0,2) = 3;
+ M(2,0) = 3;
+ M(1,1) = 20;
+ M(1,2) = 6;
+ M(2,1) = 6;
+ M(2,2) = 90;
+
+ Vector<NumberType> x(size), y(size);
+ x[0] = 2;
+ x[1] = -7;
+ x[2] = 1;
+
+ y = x;
+ M.compute_cholesky_factorization();
+ M.solve(y);
+ 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>();
+
+}