1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2019 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #include <deal.II/base/tensor.h>
17
18 #include <deal.II/lac/exceptions.h>
19 #include <deal.II/lac/lapack_templates.h>
20
21 DEAL_II_NAMESPACE_OPEN
22
23 namespace
24 {
25 template <int dim, typename Number>
calculate_svd_in_place(Tensor<2,dim,Number> & A_in_VT_out,Tensor<2,dim,Number> & U)26 void calculate_svd_in_place(Tensor<2, dim, Number> &A_in_VT_out,
27 Tensor<2, dim, Number> &U)
28 {
29 // inputs: A
30 // outputs: V^T, U
31 // SVD: A = U S V^T
32 // Since Tensor stores data in row major order and lapack expects column
33 // major ordering, we take the SVD of A^T by running the gesvd command.
34 // The results (V^T)^T and U^T are provided in column major that we use
35 // as row major results V^T and U.
36 // It essentially computs A^T = (V^T)^T S U^T and gives us V^T and U.
37 // This trick gives what we originally wanted (A = U S V^T) but the order
38 // of U and V^T is reversed.
39 std::array<Number, dim> S;
40 const types::blas_int N = dim;
41 // lwork must be >= max(1, 3*min(m,n)+max(m,n), 5*min(m,n))
42 const types::blas_int lwork = 5 * dim;
43 std::array<Number, lwork> work;
44 types::blas_int info;
45 gesvd(&LAPACKSupport::O, // replace VT in place
46 &LAPACKSupport::A,
47 &N,
48 &N,
49 A_in_VT_out.begin_raw(),
50 &N,
51 S.data(),
52 A_in_VT_out.begin_raw(),
53 &N,
54 U.begin_raw(),
55 &N,
56 work.data(),
57 &lwork,
58 &info);
59 Assert(info == 0, LAPACKSupport::ExcErrorCode("gesvd", info));
60 Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());
61 }
62 } // namespace
63
64
65
66 template <int dim, typename Number>
67 Tensor<2, dim, Number>
project_onto_orthogonal_tensors(const Tensor<2,dim,Number> & A)68 project_onto_orthogonal_tensors(const Tensor<2, dim, Number> &A)
69 {
70 Tensor<2, dim, Number> VT(A), U;
71 calculate_svd_in_place(VT, U);
72 return U * VT;
73 }
74
75
76
77 template Tensor<2, 1, float>
78 project_onto_orthogonal_tensors(const Tensor<2, 1, float> &);
79 template Tensor<2, 2, float>
80 project_onto_orthogonal_tensors(const Tensor<2, 2, float> &);
81 template Tensor<2, 3, float>
82 project_onto_orthogonal_tensors(const Tensor<2, 3, float> &);
83 template Tensor<2, 1, double>
84 project_onto_orthogonal_tensors(const Tensor<2, 1, double> &);
85 template Tensor<2, 2, double>
86 project_onto_orthogonal_tensors(const Tensor<2, 2, double> &);
87 template Tensor<2, 3, double>
88 project_onto_orthogonal_tensors(const Tensor<2, 3, double> &);
89
90 DEAL_II_NAMESPACE_CLOSE
91