Reference documentation for deal.II version GIT f3d24b0aad 2022-12-01 05:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
tensor.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2005 - 2022 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>
20 
21 #include <array>
22 
24 
25 namespace
26 {
27  template <int dim, typename Number>
28  void
29  calculate_svd_in_place(Tensor<2, dim, Number> &A_in_VT_out,
31  {
32  // inputs: A
33  // outputs: V^T, U
34  // SVD: A = U S V^T
35  // Since Tensor stores data in row major order and lapack expects column
36  // major ordering, we take the SVD of A^T by running the gesvd command.
37  // The results (V^T)^T and U^T are provided in column major that we use
38  // as row major results V^T and U.
39  // It essentially computs A^T = (V^T)^T S U^T and gives us V^T and U.
40  // This trick gives what we originally wanted (A = U S V^T) but the order
41  // of U and V^T is reversed.
42  std::array<Number, dim> S;
43  const types::blas_int N = dim;
44  // lwork must be >= max(1, 3*min(m,n)+max(m,n), 5*min(m,n))
45  const types::blas_int lwork = 5 * dim;
46  std::array<Number, lwork> work;
47  types::blas_int info;
48  gesvd(&LAPACKSupport::O, // replace VT in place
50  &N,
51  &N,
52  A_in_VT_out.begin_raw(),
53  &N,
54  S.data(),
55  A_in_VT_out.begin_raw(),
56  &N,
57  U.begin_raw(),
58  &N,
59  work.data(),
60  &lwork,
61  &info);
62  Assert(info == 0, LAPACKSupport::ExcErrorCode("gesvd", info));
63  Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());
64  }
65 } // namespace
66 
67 
68 
69 template <int dim, typename Number>
72 {
74  calculate_svd_in_place(VT, U);
75  return U * VT;
76 }
77 
78 
79 
80 template Tensor<2, 1, float>
82 template Tensor<2, 2, float>
84 template Tensor<2, 3, float>
86 template Tensor<2, 1, double>
88 template Tensor<2, 2, double>
90 template Tensor<2, 3, double>
92 
Number * begin_raw()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:458
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:459
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
#define Assert(cond, exc)
Definition: exceptions.h:1501
static ::ExceptionBase & ExcSingular()
void gesvd(const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, number3 *, const ::types::blas_int *, number4 *, const ::types::blas_int *, number5 *, const ::types::blas_int *, ::types::blas_int *)
static const char U
static const char A
static const char N
static const char O
int blas_int
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition: tensor.cc:71