lapack: add reciprocal_condition_number() for triangular matrices

diff --git a/include/deal.II/lac/lapack_full_matrix.h b/include/deal.II/lac/lapack_full_matrix.h
index 7d75ba796b3f6d4353019f2e2ad2e1a10cb8ac55..83660f69b990108de365684ef7fc9020099cd07b 100644 (file)
--- a/include/deal.II/lac/lapack_full_matrix.h
+++ b/include/deal.II/lac/lapack_full_matrix.h
@@ -444,6 +444,12 @@ public:
    */
   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().
+   */
+  number reciprocal_condition_number() const;
+
   /**
    * Compute the determinant of a matrix. As it requires the LU factorization of
    * the matrix, this function can only be called after

diff --git a/include/deal.II/lac/lapack_templates.h b/include/deal.II/lac/lapack_templates.h
index 956e232e82a0780c401b98a19afff69100781607..a48ee3739af10429e6507047e463ee8f23b1d3ce 100644 (file)
--- a/include/deal.II/lac/lapack_templates.h
+++ b/include/deal.II/lac/lapack_templates.h
@@ -83,6 +83,16 @@ extern "C"
   void spocon_ (const char *uplo, const int *n, const float *A,
                 const int *lda, const float *anorm, float *rcond,
                 float *work, int *iwork, int *info);
+// Estimate the reciprocal of the condition number of triangular matrices
+// http://www.netlib.org/lapack/explore-html/da/dba/group__double_o_t_h_e_rcomputational_gaff914510b1673e90752c095f5b9dcedf.html#gaff914510b1673e90752c095f5b9dcedf
+  void dtrcon_ (const char *norm, const char *uplo, const char *diag,
+                const int *n, const double *A, const int *lda,
+                double *rcond,
+                double *work, int *iwork, int *info);
+  void strcon_ (const char *norm, const char *uplo, const char *diag,
+                const int *n, const float *A, const int *lda,
+                float *rcond,
+                float *work, int *iwork, int *info);
 // Computes the inverse from Cholesky
   void dpotri_ (const char *uplo, const int *n, double *A,
                 const int *lda, int *info);
@@ -491,6 +501,67 @@ potrf (const char *uplo, const int *n, float *A, const int *lda, int *info)
 
 
 
+/// Tempalte wrapper for trcon
+template <typename number>
+inline void
+trcon(const char *norm, const char *uplo, const char *diag,
+      const int *n, const number *A, const int *lda,
+      number *rcond,
+      number *work, int *iwork, int *info)
+{
+  Assert (false, ExcNotImplemented());
+}
+
+inline void
+trcon  (const char *norm, const char *uplo, const char *diag,
+        const int *n, const double *A, const int *lda,
+        double *rcond,
+        double *work, int *iwork, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  dtrcon_ (norm, uplo, diag, n, A, lda, rcond, work, iwork, info);
+#else
+  (void) norm;
+  (void) uplo;
+  (void) diag;
+  (void) n;
+  (void) A;
+  (void) lda;
+  (void) rcond;
+  (void) work;
+  (void) iwork;
+  (void) info;
+
+  Assert (false, LAPACKSupport::ExcMissing("dtrcon"));
+#endif
+}
+
+inline void
+trcon  (const char *norm, const char *uplo, const char *diag,
+        const int *n, const float *A, const int *lda,
+        float *rcond,
+        float *work, int *iwork, int *info)
+{
+#ifdef DEAL_II_WITH_LAPACK
+  strcon_ (norm, uplo, diag, n, A, lda, rcond, work, iwork, info);
+#else
+  (void) norm;
+  (void) uplo;
+  (void) diag;
+  (void) n;
+  (void) A;
+  (void) lda;
+  (void) rcond;
+  (void) work;
+  (void) iwork;
+  (void) info;
+
+  Assert (false, LAPACKSupport::ExcMissing("dtrcon"));
+#endif
+}
+
+
+
 /// Template wrapper for pocon
 template <typename number1>
 inline void

diff --git a/source/lac/lapack_full_matrix.cc b/source/lac/lapack_full_matrix.cc
index f53a0bb94dcd7c0decce91c714eeddd86753aead..344f09ec5711ffbc5bda79d70d6f967b89572497 100644 (file)
--- a/source/lac/lapack_full_matrix.cc
+++ b/source/lac/lapack_full_matrix.cc
@@ -689,6 +689,38 @@ LAPACKFullMatrix<number>::reciprocal_condition_number(const number a_norm) const
 
 
 
+template <typename number>
+number
+LAPACKFullMatrix<number>::reciprocal_condition_number() const
+{
+  Threads::Mutex::ScopedLock lock (mutex);
+  Assert(property == upper_triangular ||
+         property == lower_triangular,
+         ExcProperty(property));
+  number rcond = 0.;
+
+  const int N = this->n_rows();
+  const number *values = &this->values[0];
+  int info = 0;
+  const int lda = std::max(1,N);
+  work.resize(3*N);
+  iwork.resize(N);
+
+  const char norm = '1';
+  const char diag = 'N';
+  const char uplo = (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
+  trcon(&norm, &uplo, &diag,
+        &N, values, &lda,
+        &rcond,
+        work.data(), iwork.data(), &info);
+
+  Assert(info >= 0, ExcInternalError());
+
+  return rcond;
+}
+
+
+
 template <typename number>
 void
 LAPACKFullMatrix<number>::compute_svd()

diff --git a/tests/lapack/full_matrix_22.cc b/tests/lapack/full_matrix_22.cc
new file mode 100644 (file)
index 0000000..8f94ff1
--- /dev/null
+++ b/tests/lapack/full_matrix_22.cc
@@ -0,0 +1,73 @@
+// ---------------------------------------------------------------------
+//
+// 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::reciprocal_condition_number() for triangular matrices
+
+/* MWE for size=3 in Octave:
+R = [1,2,3; 0, 5, 6; 0, 0, 9]
+rcond(R)
+ans =  0.055556
+*/
+
+#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)
+{
+  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++;
+      }
+
+  // M.print_formatted(deallog.get_file_stream(), 3, false, 8);
+  const double rcond = M.reciprocal_condition_number();
+
+  deallog << rcond << 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 = {{1,3,11}};
+  for (const auto &s : sizes)
+    {
+      deallog << "size=" << s << std::endl;
+      // test<float>(s);
+      test<double>(s);
+    }
+}

diff --git a/tests/lapack/full_matrix_22.output b/tests/lapack/full_matrix_22.output
new file mode 100644 (file)
index 0000000..1ea812f
--- /dev/null
+++ b/tests/lapack/full_matrix_22.output
@@ -0,0 +1,7 @@
+
+DEAL::size=1
+DEAL::1.00000
+DEAL::size=3
+DEAL::0.0555556
+DEAL::size=11
+DEAL::0.00137741