Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
tensor.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2005 - 2023 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
17#include <deal.II/base/tensor.h>
18
21
22#include <array>
23
25
26namespace
27{
28 template <int dim, typename Number>
29 void
30 calculate_svd_in_place(Tensor<2, dim, Number> &A_in_VT_out,
32 {
33 // inputs: A
34 // outputs: V^T, U
35 // SVD: A = U S V^T
36 // Since Tensor stores data in row major order and lapack expects column
37 // major ordering, we take the SVD of A^T by running the gesvd command.
38 // The results (V^T)^T and U^T are provided in column major that we use
39 // as row major results V^T and U.
40 // It essentially computs A^T = (V^T)^T S U^T and gives us V^T and U.
41 // This trick gives what we originally wanted (A = U S V^T) but the order
42 // of U and V^T is reversed.
43 std::array<Number, dim> S;
44 const types::blas_int N = dim;
45 // lwork must be >= max(1, 3*min(m,n)+max(m,n), 5*min(m,n))
46 const types::blas_int lwork = 5 * dim;
47 std::array<Number, lwork> work;
48 types::blas_int info;
49 constexpr std::size_t size =
51 std::array<Number, size> A_array;
52 A_in_VT_out.unroll(A_array.begin(), A_array.end());
53 std::array<Number, size> U_array;
54 U.unroll(U_array.begin(), U_array.end());
55 gesvd(&LAPACKSupport::O, // replace VT in place
57 &N,
58 &N,
59 A_array.data(),
60 &N,
61 S.data(),
62 A_array.data(),
63 &N,
64 U_array.data(),
65 &N,
66 work.data(),
67 &lwork,
68 &info);
69 Assert(info == 0, LAPACKSupport::ExcErrorCode("gesvd", info));
70 Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());
71
72 A_in_VT_out =
73 Tensor<2, dim, Number>(make_array_view(A_array.begin(), A_array.end()));
74 U = Tensor<2, dim, Number>(make_array_view(U_array.begin(), U_array.end()));
75 }
76} // namespace
77
78
79
80template <int dim, typename Number>
83{
85 calculate_svd_in_place(VT, U);
86 return U * VT;
87}
88
89
90
103
ArrayView< typename std::remove_reference< typename std::iterator_traits< Iterator >::reference >::type, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:704
void unroll(Vector< OtherNumber > &result) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
#define Assert(cond, exc)
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 A
static const char O
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition tensor.cc:82