*/
number frobenius_norm() const;
+ /**
+ * Compute trace of the matrix, i.e. the sum of the diagonal values.
+ * Obviously, the matrix needs to be quadratic for this function.
+ */
+ number trace() const;
+
/**
* Invert the matrix by first computing an LU/Cholesky factorization with the LAPACK
* function Xgetrf/Xpotrf and then building the actual inverse using Xgetri/Xpotri.
}
+
+template <typename number>
+number LAPACKFullMatrix<number>::trace() const
+{
+ Assert (state == LAPACKSupport::matrix ||
+ state == LAPACKSupport::inverse_matrix,
+ ExcMessage("Trace can be called in matrix state only."));
+ Assert (this->n_cols() == this->n_rows(),
+ ExcDimensionMismatch(this->n_cols(), this->n_rows()));
+
+ number tr = 0;
+ for (size_type i=0; i<this->n_rows(); ++i)
+ tr += (*this)(i,i);
+
+ return tr;
+}
+
+
+
template <typename number>
void
LAPACKFullMatrix<number>::compute_cholesky_factorization()
--- /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::trace() by comparing to FullMatrix
+
+#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)
+{
+ // Full matrix:
+ FullMatrix<NumberType> F(size);
+ create_random(F);
+
+ // Lapack:
+ LAPACKFullMatrix<NumberType> M(size);
+ M = F;
+
+ deallog << "trace difference: " << (F.trace() - M.trace())<< std::endl;
+}
+
+
+int main()
+{
+ const std::string logname = "output";
+ std::ofstream logfile(logname.c_str());
+ logfile.precision(3);
+ deallog.attach(logfile);
+
+ const std::vector<unsigned int> sizes = {{17,391}};
+ for (const auto &s : sizes)
+ {
+ deallog << "size=" << s << std::endl;
+ test<double>(s);
+ }
+}