deal.II version GIT relicensing-2659-g040196caa3 2025-02-18 14:20:01+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\}}\)
Loading...
Searching...
No Matches
tensor.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 2020 - 2023 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
16#include <deal.II/base/tensor.h>
17
20
21#include <array>
22
24
25namespace
26{
27 template <int dim, typename Number>
28 void
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 computes 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;
48 constexpr std::size_t size =
50 std::array<Number, size> A_array;
51 A_in_VT_out.unroll(A_array.begin(), A_array.end());
52 std::array<Number, size> U_array;
53 U.unroll(U_array.begin(), U_array.end());
54 gesvd(&LAPACKSupport::O, // replace VT in place
56 &N,
57 &N,
58 A_array.data(),
59 &N,
60 S.data(),
61 A_array.data(),
62 &N,
63 U_array.data(),
64 &N,
65 work.data(),
66 &lwork,
67 &info);
69 Assert(S.back() / S.front() > 1.e-10, LACExceptions::ExcSingular());
70
74 }
75} // namespace
76
77
78
79template <int dim, typename Number>
87
88
89
102
ArrayView< std::remove_reference_t< typename std::iterator_traits< Iterator >::reference >, MemorySpaceType > make_array_view(const Iterator begin, const Iterator end)
Definition array_view.h:949
void unroll(const Iterator begin, const Iterator end) const
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:518
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:519
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 *)
std::size_t size
Definition mpi.cc:739
constexpr char O
constexpr char N
constexpr char U
constexpr char A
Tensor< 2, dim, Number > project_onto_orthogonal_tensors(const Tensor< 2, dim, Number > &A)
Definition tensor.cc:81