From: Reza Rastak <rastak@stanford.edu>
Date: Mon, 4 May 2020 00:23:47 +0000 (-0700)
Subject: Use lapack gesvd() for orthogonal_project
X-Git-Tag: v9.3.0-rc1~1380^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10027%2Fhead;p=dealii.git

Use lapack gesvd() for orthogonal_project
---

diff --git a/include/deal.II/base/tensor.h b/include/deal.II/base/tensor.h
index dfaed33216..956bc7e674 100644
--- a/include/deal.II/base/tensor.h
+++ b/include/deal.II/base/tensor.h
@@ -25,8 +25,6 @@
 #include <deal.II/base/tensor_accessors.h>
 #include <deal.II/base/utilities.h>
 
-#include <deal.II/lac/lapack_full_matrix.h>
-
 #ifdef DEAL_II_WITH_ADOLC
 #  include <adolc/adouble.h> // Taped double
 #endif
@@ -2647,34 +2645,16 @@ cofactor(const Tensor<2, dim, Number> &t)
  * where $\mathbf U$ and $\mathbf V$ are computed from the SVD decomposition:
  * $\mathbf U  \mathbf S \mathbf V^T$,
  * effectively replacing $\mathbf S$ with the identity matrix.
- * @param tensor The tensor which to find the closest orthogonal
- * tensor to.
+ * @param A The tensor which to find the closest orthogonal tensor to.
+ * @pre @p Number must be either `float` or `double`.
+ * @pre In order to use this function, this program must be linked with the
+ *      LAPACK library.
+ * @pre @p A must not be singular.
  * @relatesalso Tensor
  */
 template <int dim, typename Number>
 Tensor<2, dim, Number>
-project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor)
-{
-  Tensor<2, dim, Number>   output_tensor;
-  FullMatrix<Number>       matrix(dim);
-  LAPACKFullMatrix<Number> lapack_matrix(dim);
-  LAPACKFullMatrix<Number> result(dim);
-
-  // todo: find or add dealii functionality to copy in one step.
-  matrix.copy_from(tensor);
-  lapack_matrix.copy_from(matrix);
-
-  // now compute the svd of the matrices
-  lapack_matrix.compute_svd();
-
-  // Use the SVD results to orthogonalize: $U V^T$
-  lapack_matrix.get_svd_u().mmult(result, lapack_matrix.get_svd_vt());
-
-  // todo: find or add dealii functionality to copy in one step.
-  matrix = result;
-  matrix.copy_to(output_tensor);
-  return output_tensor;
-}
+project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &A);
 
 
 /**
diff --git a/include/deal.II/lac/lapack_support.h b/include/deal.II/lac/lapack_support.h
index 75b001352a..6870e4b4c9 100644
--- a/include/deal.II/lac/lapack_support.h
+++ b/include/deal.II/lac/lapack_support.h
@@ -155,6 +155,10 @@ namespace LAPACKSupport
    * Character constant.
    */
   static const char N = 'N';
+  /**
+   * Character constant.
+   */
+  static const char O = 'O';
   /**
    * Character constant.
    */
diff --git a/source/base/CMakeLists.txt b/source/base/CMakeLists.txt
index c16cdecb9e..a046beecc8 100644
--- a/source/base/CMakeLists.txt
+++ b/source/base/CMakeLists.txt
@@ -83,6 +83,7 @@ SET(_unity_include_src
   subscriptor.cc
   symbolic_function.cc
   table_handler.cc
+  tensor.cc
   tensor_function.cc
   tensor_function_parser.cc
   tensor_product_polynomials.cc
diff --git a/source/base/tensor.cc b/source/base/tensor.cc
new file mode 100644
index 0000000000..758f89a506
--- /dev/null
+++ b/source/base/tensor.cc
@@ -0,0 +1,90 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2005 - 2019 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/tensor.h>
+
+#include <deal.II/lac/exceptions.h>
+#include <deal.II/lac/lapack_templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace
+{
+  template <int dim, typename Number>
+  void calculate_svd_in_place(Tensor<2, dim, Number> &A_in_VT_out,
+                              Tensor<2, dim, Number> &U)
+  {
+    // inputs: A
+    // outputs: V^T, U
+    // SVD: A = U S V^T
+    // Since Tensor stores data in row major order and lapack expects column
+    // major ordering, we take the SVD of A^T by running the gesvd command.
+    // The results (V^T)^T and U^T are provided in column major that we use
+    // as row major results V^T and U.
+    // It essentially computs A^T = (V^T)^T S U^T and gives us V^T and U.
+    // This trick gives what we originally wanted (A = U S V^T) but the order
+    // of U and V^T is reversed.
+    std::array<Number, dim> S;
+    const types::blas_int   N = dim;
+    // lwork must be >= max(1, 3*min(m,n)+max(m,n), 5*min(m,n))
+    const types::blas_int     lwork = 5 * dim;
+    std::array<Number, lwork> work;
+    types::blas_int           info;
+    gesvd(&LAPACKSupport::O, // replace VT in place
+          &LAPACKSupport::A,
+          &N,
+          &N,
+          A_in_VT_out.begin_raw(),
+          &N,
+          S.data(),
+          A_in_VT_out.begin_raw(),
+          &N,
+          U.begin_raw(),
+          &N,
+          work.data(),
+          &lwork,
+          &info);
+    Assert(info == 0, LAPACKSupport::ExcErrorCode("gesvd", info));
+    Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());
+  }
+} // namespace
+
+
+
+template <int dim, typename Number>
+Tensor<2, dim, Number>
+project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &A)
+{
+  Tensor<2, dim, Number> VT(A), U;
+  calculate_svd_in_place(VT, U);
+  return U * VT;
+}
+
+
+
+template Tensor<2, 1, float>
+project_onto_orthogonal_tensors(const Tensor<2, 1, float> &);
+template Tensor<2, 2, float>
+project_onto_orthogonal_tensors(const Tensor<2, 2, float> &);
+template Tensor<2, 3, float>
+project_onto_orthogonal_tensors(const Tensor<2, 3, float> &);
+template Tensor<2, 1, double>
+project_onto_orthogonal_tensors(const Tensor<2, 1, double> &);
+template Tensor<2, 2, double>
+project_onto_orthogonal_tensors(const Tensor<2, 2, double> &);
+template Tensor<2, 3, double>
+project_onto_orthogonal_tensors(const Tensor<2, 3, double> &);
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/tests/tensors/project_orthogonal.cc b/tests/tensors/project_orthogonal.cc
index 5187b41892..a9cf7dfb50 100644
--- a/tests/tensors/project_orthogonal.cc
+++ b/tests/tensors/project_orthogonal.cc
@@ -17,8 +17,6 @@
 
 #include <deal.II/base/tensor.h>
 
-#include <deal.II/lac/full_matrix.h>
-
 #include "../tests.h"
 
 template <int dim, typename number>
@@ -53,18 +51,24 @@ main()
 {
   initlog();
   // not orthogonal
-  test(Tensor<2, 3, double>{{{1, 2, 3}, {4, 5, 6}, {7, 8, 9}}});
+  test(Tensor<2, 3, float>{{{1, 2, 3}, {2, 1, 4}, {3, 4, 1}}});
+  test(Tensor<2, 3, double>{{{1, 2, 3}, {2, 1, 4}, {3, 4, 1}}});
   // already orthogonal
+  test(Tensor<2, 3, float>{{{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}});
   test(Tensor<2, 3, double>{{{1, 0, 0}, {0, 0, 1}, {0, 1, 0}}});
   // not orthogonal, but det = 1
   test(Tensor<2, 3, double>{{{1, 2, 0}, {0, 1, 0}, {0, 0, 1}}});
   // 2D not orthogonal
+  test(Tensor<2, 2, float>{{{1, 2}, {3, 4}}});
   test(Tensor<2, 2, double>{{{1, 2}, {3, 4}}});
   // 2D already orthogonal
+  test(Tensor<2, 2, float>{{{0, 1}, {1, 0}}});
   test(Tensor<2, 2, double>{{{0, 1}, {1, 0}}});
   // 1D not orthogonal
+  test(Tensor<2, 1, float>{{{2.5}}});
   test(Tensor<2, 1, double>{{{2.5}}});
   // 1D already orthogonal
+  test(Tensor<2, 1, float>{{{1.}}});
   test(Tensor<2, 1, double>{{{1.}}});
 
   deallog << "OK" << std::endl;
diff --git a/tests/tensors/project_orthogonal.with_lapack=true.output b/tests/tensors/project_orthogonal.with_lapack=true.output
index cec2818537..a7a807d4c0 100644
--- a/tests/tensors/project_orthogonal.with_lapack=true.output
+++ b/tests/tensors/project_orthogonal.with_lapack=true.output
@@ -1,11 +1,27 @@
 
+DEAL:: Tensor<2, 3, float>
+DEAL::original tensor A = 1.00000 2.00000 3.00000 2.00000 1.00000 4.00000 3.00000 4.00000 1.00000
+DEAL::det A = 20.0000
+DEAL::A A^T = 14.0000 16.0000 14.0000 16.0000 21.0000 14.0000 14.0000 14.0000 26.0000
+DEAL::projected to Q = -0.488363 0.591135 0.641920 0.591135 -0.317015 0.741661 0.641920 0.741661 -0.194623
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000 -1.78814e-07 1.49012e-07 -1.78814e-07 1.00000 -1.49012e-08 1.49012e-07 -1.49012e-08 1.00000
+DEAL::
 DEAL:: Tensor<2, 3, double>
-DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000 5.00000 6.00000 7.00000 8.00000 9.00000
-DEAL::det A = 0.00000
-DEAL::A A^T = 14.0000 32.0000 50.0000 32.0000 77.0000 122.000 50.0000 122.000 194.000
-DEAL::projected to Q = -0.752720 0.389148 0.531015 0.389148 -0.387594 0.835664 0.531015 0.835664 0.140313
+DEAL::original tensor A = 1.00000 2.00000 3.00000 2.00000 1.00000 4.00000 3.00000 4.00000 1.00000
+DEAL::det A = 20.0000
+DEAL::A A^T = 14.0000 16.0000 14.0000 16.0000 21.0000 14.0000 14.0000 14.0000 26.0000
+DEAL::projected to Q = -0.488363 0.591135 0.641920 0.591135 -0.317014 0.741661 0.641920 0.741661 -0.194623
 DEAL::det Q = 1.00000
-DEAL::Q Q^T = 1.00000 5.55112e-17 2.77556e-16 5.55112e-17 1.00000 5.55112e-17 2.77556e-16 5.55112e-17 1.00000
+DEAL::Q Q^T = 1.00000 7.21645e-16 2.08167e-16 7.21645e-16 1.00000 -2.77556e-16 2.08167e-16 -2.77556e-16 1.00000
+DEAL::
+DEAL:: Tensor<2, 3, float>
+DEAL::original tensor A = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.00000 0.00000
+DEAL::det A = -1.00000
+DEAL::A A^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.00000 0.00000
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
 DEAL::
 DEAL:: Tensor<2, 3, double>
 DEAL::original tensor A = 1.00000 0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 1.00000 0.00000
@@ -23,6 +39,14 @@ DEAL::projected to Q = 0.707107 0.707107 0.00000 -0.707107 0.707107 0.00000 0.00
 DEAL::det Q = 1.00000
 DEAL::Q Q^T = 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000
 DEAL::
+DEAL:: Tensor<2, 2, float>
+DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000
+DEAL::det A = -2.00000
+DEAL::A A^T = 5.00000 11.0000 11.0000 25.0000
+DEAL::projected to Q = -0.514496 0.857493 0.857493 0.514496
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 -5.96046e-08 -5.96046e-08 1.00000
+DEAL::
 DEAL:: Tensor<2, 2, double>
 DEAL::original tensor A = 1.00000 2.00000 3.00000 4.00000
 DEAL::det A = -2.00000
@@ -31,6 +55,14 @@ DEAL::projected to Q = -0.514496 0.857493 0.857493 0.514496
 DEAL::det Q = -1.00000
 DEAL::Q Q^T = 1.00000 2.77556e-16 2.77556e-16 1.00000
 DEAL::
+DEAL:: Tensor<2, 2, float>
+DEAL::original tensor A = 0.00000 1.00000 1.00000 0.00000
+DEAL::det A = -1.00000
+DEAL::A A^T = 1.00000 0.00000 0.00000 1.00000
+DEAL::projected to Q = 0.00000 1.00000 1.00000 0.00000
+DEAL::det Q = -1.00000
+DEAL::Q Q^T = 1.00000 0.00000 0.00000 1.00000
+DEAL::
 DEAL:: Tensor<2, 2, double>
 DEAL::original tensor A = 0.00000 1.00000 1.00000 0.00000
 DEAL::det A = -1.00000
@@ -39,6 +71,14 @@ DEAL::projected to Q = 0.00000 1.00000 1.00000 0.00000
 DEAL::det Q = -1.00000
 DEAL::Q Q^T = 1.00000 0.00000 0.00000 1.00000
 DEAL::
+DEAL:: Tensor<2, 1, float>
+DEAL::original tensor A = 2.50000
+DEAL::det A = 2.50000
+DEAL::A A^T = 6.25000
+DEAL::projected to Q = 1.00000
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000
+DEAL::
 DEAL:: Tensor<2, 1, double>
 DEAL::original tensor A = 2.50000
 DEAL::det A = 2.50000
@@ -47,6 +87,14 @@ DEAL::projected to Q = 1.00000
 DEAL::det Q = 1.00000
 DEAL::Q Q^T = 1.00000
 DEAL::
+DEAL:: Tensor<2, 1, float>
+DEAL::original tensor A = 1.00000
+DEAL::det A = 1.00000
+DEAL::A A^T = 1.00000
+DEAL::projected to Q = 1.00000
+DEAL::det Q = 1.00000
+DEAL::Q Q^T = 1.00000
+DEAL::
 DEAL:: Tensor<2, 1, double>
 DEAL::original tensor A = 1.00000
 DEAL::det A = 1.00000