From 67df80eec292cfd295f75b65611dba7a775ee7c2 Mon Sep 17 00:00:00 2001 From: Reza Rastak Date: Sun, 3 May 2020 17:23:47 -0700 Subject: [PATCH] Use lapack gesvd() for orthogonal_project --- include/deal.II/base/tensor.h | 32 ++----- include/deal.II/lac/lapack_support.h | 4 + source/base/CMakeLists.txt | 1 + source/base/tensor.cc | 90 +++++++++++++++++++ tests/tensors/project_orthogonal.cc | 10 ++- ...project_orthogonal.with_lapack=true.output | 58 ++++++++++-- 6 files changed, 161 insertions(+), 34 deletions(-) create mode 100644 source/base/tensor.cc 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 #include -#include - #ifdef DEAL_II_WITH_ADOLC # include // 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 Tensor<2, dim, Number> -project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &tensor) -{ - Tensor<2, dim, Number> output_tensor; - FullMatrix matrix(dim); - LAPACKFullMatrix lapack_matrix(dim); - LAPACKFullMatrix 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 + +#include +#include + +DEAL_II_NAMESPACE_OPEN + +namespace +{ + template + 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 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 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 +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 -#include - #include "../tests.h" template @@ -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 -- 2.39.5