deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00: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
lapack_full_matrix.cc
Go to the documentation of this file.
1// -----------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
4// Copyright (C) 2005 - 2025 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Detailed license information governing the source code and contributions
9// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
10//
11// -----------------------------------------------------------------------------
12
13#include <deal.II/base/numbers.h>
14
15#include <deal.II/lac/blas_extension_templates.h>
16#include <deal.II/lac/block_vector.h>
17#include <deal.II/lac/full_matrix.h>
18#include <deal.II/lac/lapack_full_matrix.h>
19#include <deal.II/lac/lapack_support.h>
20#include <deal.II/lac/lapack_templates.h>
21#include <deal.II/lac/sparse_matrix.h>
22#include <deal.II/lac/utilities.h>
23#include <deal.II/lac/vector.h>
24
25#include <iomanip>
26#include <iostream>
27
28DEAL_II_NAMESPACE_OPEN
29
30using namespace LAPACKSupport;
31
32namespace internal
33{
34 namespace LAPACKFullMatrixImplementation
35 {
36 // ZGEEV/CGEEV and DGEEV/SGEEV need different work arrays and different
37 // output arrays for eigenvalues. This makes working with generic scalar
38 // types a bit difficult. To get around this, geev_helper has the same
39 // signature for real and complex arguments, but it ignores some
40 // parameters when called with a real type and ignores different
41 // parameters when called with a complex type.
42 template <typename T>
43 void
44 geev_helper(const char vl,
45 const char vr,
46 AlignedVector<T> &matrix,
47 const types::blas_int n_rows,
48 std::vector<T> &real_part_eigenvalues,
49 std::vector<T> &imag_part_eigenvalues,
50 std::vector<T> &left_eigenvectors,
51 std::vector<T> &right_eigenvectors,
52 std::vector<T> &real_work,
53 std::vector<T> & /*complex_work*/,
54 const types::blas_int work_flag,
55 types::blas_int &info)
56 {
57 static_assert(std::is_same_v<T, double> || std::is_same_v<T, float>,
58 "Only implemented for double and float");
59 Assert(matrix.size() == static_cast<std::size_t>(n_rows * n_rows),
60 ExcInternalError());
61 Assert(static_cast<std::size_t>(n_rows) <= real_part_eigenvalues.size(),
62 ExcInternalError());
63 Assert(static_cast<std::size_t>(n_rows) <= imag_part_eigenvalues.size(),
64 ExcInternalError());
65 if (vl == 'V')
66 Assert(static_cast<std::size_t>(n_rows * n_rows) <=
67 left_eigenvectors.size(),
68 ExcInternalError());
69 if (vr == 'V')
70 Assert(static_cast<std::size_t>(n_rows * n_rows) <=
71 right_eigenvectors.size(),
72 ExcInternalError());
73 Assert(work_flag == -1 ||
74 static_cast<std::size_t>(2 * n_rows) <= real_work.size(),
75 ExcInternalError());
76 Assert(work_flag == -1 || std::max<long int>(1, 3 * n_rows) <= work_flag,
77 ExcInternalError());
78 geev(&vl,
79 &vr,
80 &n_rows,
81 matrix.data(),
82 &n_rows,
83 real_part_eigenvalues.data(),
84 imag_part_eigenvalues.data(),
85 left_eigenvectors.data(),
86 &n_rows,
87 right_eigenvectors.data(),
88 &n_rows,
89 real_work.data(),
90 &work_flag,
91 &info);
92 }
93
94
95
96 template <typename T>
97 void
98 geev_helper(const char vl,
99 const char vr,
100 AlignedVector<std::complex<T>> &matrix,
101 const types::blas_int n_rows,
102 std::vector<T> & /*real_part_eigenvalues*/,
103 std::vector<std::complex<T>> &eigenvalues,
104 std::vector<std::complex<T>> &left_eigenvectors,
105 std::vector<std::complex<T>> &right_eigenvectors,
106 std::vector<std::complex<T>> &complex_work,
107 std::vector<T> &real_work,
108 const types::blas_int work_flag,
109 types::blas_int &info)
110 {
111 static_assert(
112 std::is_same_v<T, double> || std::is_same_v<T, float>,
113 "Only implemented for std::complex<double> and std::complex<float>");
114 Assert(matrix.size() == static_cast<std::size_t>(n_rows * n_rows),
115 ExcInternalError());
116 Assert(static_cast<std::size_t>(n_rows) <= eigenvalues.size(),
117 ExcInternalError());
118 if (vl == 'V')
119 Assert(static_cast<std::size_t>(n_rows * n_rows) <=
120 left_eigenvectors.size(),
121 ExcInternalError());
122 if (vr == 'V')
123 Assert(static_cast<std::size_t>(n_rows * n_rows) <=
124 right_eigenvectors.size(),
125 ExcInternalError());
126 Assert(work_flag == -1 ||
127 std::max<std::size_t>(1, work_flag) <= real_work.size(),
128 ExcInternalError());
129 Assert(work_flag == -1 || std::max<long int>(1, 2 * n_rows) <= work_flag,
130 ExcInternalError());
131
132 geev(&vl,
133 &vr,
134 &n_rows,
135 matrix.data(),
136 &n_rows,
137 eigenvalues.data(),
138 left_eigenvectors.data(),
139 &n_rows,
140 right_eigenvectors.data(),
141 &n_rows,
142 complex_work.data(),
143 &work_flag,
144 real_work.data(),
145 &info);
146 }
147
148
149
150 template <typename T>
151 void
152 gesdd_helper(const char job,
153 const types::blas_int n_rows,
154 const types::blas_int n_cols,
155 AlignedVector<T> &matrix,
156 std::vector<T> &singular_values,
157 AlignedVector<T> &left_vectors,
158 AlignedVector<T> &right_vectors,
159 std::vector<T> &real_work,
160 std::vector<T> & /*complex work*/,
161 std::vector<types::blas_int> &integer_work,
162 const types::blas_int work_flag,
163 types::blas_int &info)
164 {
165 Assert(job == 'A' || job == 'S' || job == 'O' || job == 'N',
166 ExcInternalError());
167 Assert(static_cast<std::size_t>(n_rows * n_cols) == matrix.size(),
168 ExcInternalError());
169 Assert(std::min<std::size_t>(n_rows, n_cols) <= singular_values.size(),
170 ExcInternalError());
171 Assert(8 * std::min<std::size_t>(n_rows, n_cols) <= integer_work.size(),
172 ExcInternalError());
173 Assert(work_flag == -1 ||
174 static_cast<std::size_t>(work_flag) <= real_work.size(),
175 ExcInternalError());
176 gesdd(&job,
177 &n_rows,
178 &n_cols,
179 matrix.data(),
180 &n_rows,
181 singular_values.data(),
182 left_vectors.data(),
183 &n_rows,
184 right_vectors.data(),
185 &n_cols,
186 real_work.data(),
187 &work_flag,
188 integer_work.data(),
189 &info);
190 }
191
192
193
194 template <typename T>
195 void
196 gesdd_helper(const char job,
197 const types::blas_int n_rows,
198 const types::blas_int n_cols,
199 AlignedVector<std::complex<T>> &matrix,
200 std::vector<T> &singular_values,
201 AlignedVector<std::complex<T>> &left_vectors,
202 AlignedVector<std::complex<T>> &right_vectors,
203 std::vector<std::complex<T>> &work,
204 std::vector<T> &real_work,
205 std::vector<types::blas_int> &integer_work,
206 const types::blas_int &work_flag,
207 types::blas_int &info)
208 {
209 Assert(job == 'A' || job == 'S' || job == 'O' || job == 'N',
210 ExcInternalError());
211 Assert(static_cast<std::size_t>(n_rows * n_cols) == matrix.size(),
212 ExcInternalError());
213 Assert(static_cast<std::size_t>(std::min(n_rows, n_cols)) <=
214 singular_values.size(),
215 ExcInternalError());
216 Assert(8 * std::min<std::size_t>(n_rows, n_cols) <= integer_work.size(),
217 ExcInternalError());
218 Assert(work_flag == -1 ||
219 static_cast<std::size_t>(work_flag) <= real_work.size(),
220 ExcInternalError());
221
222 gesdd(&job,
223 &n_rows,
224 &n_cols,
225 matrix.data(),
226 &n_rows,
227 singular_values.data(),
228 left_vectors.data(),
229 &n_rows,
230 right_vectors.data(),
231 &n_cols,
232 work.data(),
233 &work_flag,
234 real_work.data(),
235 integer_work.data(),
236 &info);
237 }
238 } // namespace LAPACKFullMatrixImplementation
239} // namespace internal
240
241
242
243template <typename number>
244LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type n)
245 : TransposeTable<number>(n, n)
246 , state(matrix)
247 , property(general)
248{}
249
250
251
252template <typename number>
253LAPACKFullMatrix<number>::LAPACKFullMatrix(const size_type m, const size_type n)
254 : TransposeTable<number>(m, n)
255 , state(matrix)
256 , property(general)
257{}
258
259
260
261template <typename number>
262LAPACKFullMatrix<number>::LAPACKFullMatrix(const LAPACKFullMatrix &M)
263 : TransposeTable<number>(M)
264 , state(matrix)
265 , property(general)
266{}
267
268
269
270template <typename number>
271LAPACKFullMatrix<number> &
272LAPACKFullMatrix<number>::operator=(const LAPACKFullMatrix<number> &M)
273{
274 TransposeTable<number>::operator=(M);
275 state = M.state;
276 property = M.property;
277 return *this;
278}
279
280
281
282template <typename number>
283void
289
290
291
292template <typename number>
293void
294LAPACKFullMatrix<number>::grow_or_shrink(const size_type n)
295{
296 const size_type s = std::min(std::min(this->m(), n), this->n());
297 TransposeTable<number> copy(std::move(*this));
298 this->TransposeTable<number>::reinit(n, n);
299 for (size_type i = 0; i < s; ++i)
300 for (size_type j = 0; j < s; ++j)
301 (*this)(i, j) = copy(i, j);
302}
303
304
305
306template <typename number>
307void
308LAPACKFullMatrix<number>::apply_givens_rotation(
309 const std::array<number, 3> &csr,
310 const size_type i,
311 const size_type k,
312 const bool left)
313{
314 auto &A = *this;
315 // see Golub 2013 "Matrix computations", p241 5.1.9 Applying Givens
316 // Rotations but note the difference in notation, namely the sign of s: we
317 // have G * A, where G[1,1] = s
318 if (left)
319 {
320 for (size_type j = 0; j < A.n(); ++j)
321 {
322 const number t = A(i, j);
323 A(i, j) = csr[0] * A(i, j) + csr[1] * A(k, j);
324 A(k, j) = -csr[1] * t + csr[0] * A(k, j);
325 }
326 }
327 else
328 {
329 for (size_type j = 0; j < A.m(); ++j)
330 {
331 const number t = A(j, i);
332 A(j, i) = csr[0] * A(j, i) + csr[1] * A(j, k);
333 A(j, k) = -csr[1] * t + csr[0] * A(j, k);
334 }
335 }
336}
337
338
339
340template <typename number>
341void
342LAPACKFullMatrix<number>::remove_row_and_column(const size_type row,
343 const size_type col)
344{
345 AssertIndexRange(row, this->m());
346 AssertIndexRange(col, this->n());
347
348 const size_type nrows = this->m() - 1;
349 const size_type ncols = this->n() - 1;
350
351 TransposeTable<number> copy(std::move(*this));
352 this->TransposeTable<number>::reinit(nrows, ncols);
353
354 for (size_type j = 0; j < ncols; ++j)
355 {
356 const size_type jj = (j < col ? j : j + 1);
357 for (size_type i = 0; i < nrows; ++i)
358 {
359 const size_type ii = (i < row ? i : i + 1);
360 (*this)(i, j) = copy(ii, jj);
361 }
362 }
363}
364
365
366
367template <typename number>
368void
374
375
376
377template <typename number>
378template <typename number2>
379LAPACKFullMatrix<number> &
380LAPACKFullMatrix<number>::operator=(const FullMatrix<number2> &M)
381{
382 Assert(this->m() == M.m(), ExcDimensionMismatch(this->m(), M.m()));
383 Assert(this->n() == M.n(), ExcDimensionMismatch(this->n(), M.n()));
384 for (size_type i = 0; i < this->m(); ++i)
385 for (size_type j = 0; j < this->n(); ++j)
386 (*this)(i, j) = M(i, j);
387
388 state = LAPACKSupport::matrix;
389 property = LAPACKSupport::general;
390 return *this;
391}
392
393
394
395template <typename number>
396template <typename number2>
397LAPACKFullMatrix<number> &
398LAPACKFullMatrix<number>::operator=(const SparseMatrix<number2> &M)
399{
400 Assert(this->m() == M.n(), ExcDimensionMismatch(this->m(), M.n()));
401 Assert(this->n() == M.m(), ExcDimensionMismatch(this->n(), M.m()));
402 for (size_type i = 0; i < this->m(); ++i)
403 for (size_type j = 0; j < this->n(); ++j)
404 (*this)(i, j) = M.el(i, j);
405
406 state = LAPACKSupport::matrix;
407 property = LAPACKSupport::general;
408 return *this;
409}
410
411
412
413template <typename number>
414LAPACKFullMatrix<number> &
415LAPACKFullMatrix<number>::operator=(const number d)
416{
417 Assert(d == number(0), ExcScalarAssignmentOnlyForZeroValue());
418
419 if (this->n_elements() != 0)
420 this->reset_values();
421
422 state = LAPACKSupport::matrix;
423 return *this;
424}
425
426
427
428template <typename number>
429LAPACKFullMatrix<number> &
430LAPACKFullMatrix<number>::operator*=(const number factor)
431{
432 Assert(state == LAPACKSupport::matrix ||
433 state == LAPACKSupport::inverse_matrix,
434 ExcState(state));
435
436 AssertIsFinite(factor);
437 const char type = 'G';
438 const number cfrom = 1.;
439 const types::blas_int m = this->m();
440 const types::blas_int n = this->n();
441 const types::blas_int lda = this->m();
442 types::blas_int info = 0;
443 // kl and ku will not be referenced for type = G (dense matrices).
444 const types::blas_int kl = 0;
445 number *values = this->values.data();
446
447 lascl(&type, &kl, &kl, &cfrom, &factor, &m, &n, values, &lda, &info);
448
449 // Negative return value implies a wrong argument. This should be internal.
450 Assert(info >= 0, ExcInternalError());
451
452 return *this;
453}
454
455
456
457template <typename number>
458LAPACKFullMatrix<number> &
459LAPACKFullMatrix<number>::operator/=(const number factor)
460{
461 Assert(state == LAPACKSupport::matrix ||
462 state == LAPACKSupport::inverse_matrix,
463 ExcState(state));
464
465 AssertIsFinite(factor);
466 Assert(factor != number(0.), ExcZero());
467
468 const char type = 'G';
469 const number cto = 1.;
470 const types::blas_int m = this->m();
471 const types::blas_int n = this->n();
472 const types::blas_int lda = this->m();
473 types::blas_int info = 0;
474 // kl and ku will not be referenced for type = G (dense matrices).
475 const types::blas_int kl = 0;
476 number *values = this->values.data();
477
478 lascl(&type, &kl, &kl, &factor, &cto, &m, &n, values, &lda, &info);
479
480 // Negative return value implies a wrong argument. This should be internal.
481 Assert(info >= 0, ExcInternalError());
482
483 return *this;
484}
485
486
487
488template <typename number>
489void
490LAPACKFullMatrix<number>::add(const number a, const LAPACKFullMatrix<number> &A)
491{
492 Assert(state == LAPACKSupport::matrix ||
493 state == LAPACKSupport::inverse_matrix,
494 ExcState(state));
495
496 Assert(m() == A.m(), ExcDimensionMismatch(m(), A.m()));
497 Assert(n() == A.n(), ExcDimensionMismatch(n(), A.n()));
498
499 AssertIsFinite(a);
500
501 // BLAS does not offer functions to add matrices.
502 // LapackFullMatrix is stored in contiguous array
503 // ==> use BLAS 1 for adding vectors
504 const types::blas_int n = this->m() * this->n();
505 const types::blas_int inc = 1;
506 number *values = this->values.data();
507 const number *values_A = A.values.data();
508
509 axpy(&n, &a, values_A, &inc, values, &inc);
510}
511
512
513
514namespace
515{
516 template <typename number>
517 void
518 cholesky_rank1(LAPACKFullMatrix<number> &A,
519 const number a,
520 const Vector<number> &v)
521 {
522 const typename LAPACKFullMatrix<number>::size_type N = A.n();
523 Vector<number> z(v);
524 // Cholesky update / downdate, see
525 // 6.5.4 Cholesky Updating and Downdating, Golub 2013 Matrix computations
526 // Note that potrf() is called with LAPACKSupport::L , so the
527 // factorization is stored in lower triangular part.
528 // Also see discussion here
529 // http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=2646
530 if (a > 0.)
531 {
532 // simple update via a sequence of Givens rotations.
533 // Observe that
534 //
535 // | L^T |T | L^T |
536 // A <-- | | | | = L L^T + z z^T
537 // | z^T | | z^T |
538 //
539 // so we can get updated factor by doing a sequence of Givens
540 // rotations to make the matrix lower-triangular
541 // Also see LINPACK's dchud http://www.netlib.org/linpack/dchud.f
542 z *= std::sqrt(a);
543 for (typename LAPACKFullMatrix<number>::size_type k = 0; k < N; ++k)
544 {
545 const std::array<number, 3> csr =
546 Utilities::LinearAlgebra::givens_rotation(A(k, k), z(k));
547 A(k, k) = csr[2];
548 for (typename LAPACKFullMatrix<number>::size_type i = k + 1; i < N;
549 ++i)
550 {
551 const number t = A(i, k);
552 A(i, k) = csr[0] * A(i, k) + csr[1] * z(i);
553 z(i) = -csr[1] * t + csr[0] * z(i);
554 }
555 }
556 }
557 else
558 {
559 // downdating is not always possible as we may end up with
560 // negative definite matrix. If it's possible, then it boils
561 // down to application of hyperbolic rotations.
562 // Observe that
563 //
564 // | L^T |T | L^T |
565 // A <-- | | S | | = L L^T - z z^T
566 // | z^T | | z^T |
567 //
568 // |In 0 |
569 // S := | |
570 // |0 -1 |
571 //
572 // We are looking for H which is S-orthogonal (HSH^T=S) and
573 // can restore upper-triangular factor of the factorization of A above.
574 // We will use Hyperbolic rotations to do the job
575 //
576 // | c -s | | x1 | | r |
577 // | | = | | = | |, c^2 - s^2 = 1
578 // |-s c | | x2 | | 0 |
579 //
580 // which have real solution only if x2 <= x1.
581 // See also Linpack's http://www.netlib.org/linpack/dchdd.f and
582 // https://infoscience.epfl.ch/record/161468/files/cholupdate.pdf and
583 // "Analysis of a recursive Least Squares Hyperbolic Rotation Algorithm
584 // for Signal Processing", Alexander, Pan, Plemmons, 1988.
585 z *= std::sqrt(-a);
586 for (typename LAPACKFullMatrix<number>::size_type k = 0; k < N; ++k)
587 {
588 const std::array<number, 3> csr =
589 Utilities::LinearAlgebra::hyperbolic_rotation(A(k, k), z(k));
590 A(k, k) = csr[2];
591 for (typename LAPACKFullMatrix<number>::size_type i = k + 1; i < N;
592 ++i)
593 {
594 const number t = A(i, k);
595 A(i, k) = csr[0] * A(i, k) - csr[1] * z(i);
596 z(i) = -csr[1] * t + csr[0] * z(i);
597 }
598 }
599 }
600 }
601
602
603 template <typename number>
604 void
605 cholesky_rank1(LAPACKFullMatrix<std::complex<number>> & /*A*/,
606 const std::complex<number> /*a*/,
607 const Vector<std::complex<number>> & /*v*/)
608 {
609 AssertThrow(false, ExcNotImplemented());
610 }
611} // namespace
612
613
614
615template <typename number>
616void
617LAPACKFullMatrix<number>::rank1_update(const number a, const Vector<number> &v)
618{
619 Assert(property == LAPACKSupport::symmetric, ExcProperty(property));
620
621 Assert(n() == m(), ExcInternalError());
622 Assert(m() == v.size(), ExcDimensionMismatch(m(), v.size()));
623
624 AssertIsFinite(a);
625
626 if (state == LAPACKSupport::matrix)
627 {
628 {
629 const types::blas_int N = this->m();
630 const char uplo = LAPACKSupport::U;
631 const types::blas_int lda = N;
632 const types::blas_int incx = 1;
633
634 syr(&uplo, &N, &a, v.begin(), &incx, this->values.begin(), &lda);
635 }
636
637 const size_type N = this->m();
638 // FIXME: we should really only update upper or lower triangular parts
639 // of symmetric matrices and make sure the interface is consistent,
640 // for example operator(i,j) gives correct results regardless of storage.
641 for (size_type i = 0; i < N; ++i)
642 for (size_type j = 0; j < i; ++j)
643 (*this)(i, j) = (*this)(j, i);
644 }
645 else if (state == LAPACKSupport::cholesky)
646 {
647 cholesky_rank1(*this, a, v);
648 }
649 else
650 AssertThrow(false, ExcState(state));
651}
652
653
654
655template <typename number>
656void
657LAPACKFullMatrix<number>::vmult(Vector<number> &w,
658 const Vector<number> &v,
659 const bool adding) const
660{
661 const types::blas_int mm = this->m();
662 const types::blas_int nn = this->n();
663 const number alpha = 1.;
664 const number beta = (adding ? 1. : 0.);
665 const number null = 0.;
666
667 // use trmv for triangular matrices
668 if ((property == upper_triangular || property == lower_triangular) &&
669 (mm == nn) && state == matrix)
670 {
671 Assert(adding == false, ExcNotImplemented());
672
673 AssertDimension(v.size(), this->n());
674 AssertDimension(w.size(), this->m());
675
676 const char diag = 'N';
677 const char trans = 'N';
678 const char uplo =
679 (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
680
681 w = v;
682
683 const types::blas_int N = mm;
684 const types::blas_int lda = N;
685 const types::blas_int incx = 1;
686
687 trmv(
688 &uplo, &trans, &diag, &N, this->values.data(), &lda, w.data(), &incx);
689
690 return;
691 }
692
693 switch (state)
694 {
695 case matrix:
696 case inverse_matrix:
697 {
698 AssertDimension(v.size(), this->n());
699 AssertDimension(w.size(), this->m());
700
701 gemv("N",
702 &mm,
703 &nn,
704 &alpha,
705 this->values.data(),
706 &mm,
707 v.data(),
708 &one,
709 &beta,
710 w.data(),
711 &one);
712 break;
713 }
714 case svd:
715 {
716 std::scoped_lock lock(mutex);
717 AssertDimension(v.size(), this->n());
718 AssertDimension(w.size(), this->m());
719 // Compute V^T v
720 work.resize(std::max(mm, nn));
721 gemv("N",
722 &nn,
723 &nn,
724 &alpha,
725 svd_vt->values.data(),
726 &nn,
727 v.data(),
728 &one,
729 &null,
730 work.data(),
731 &one);
732 // Multiply by singular values
733 for (size_type i = 0; i < wr.size(); ++i)
734 work[i] *= wr[i];
735 // Multiply with U
736 gemv("N",
737 &mm,
738 &mm,
739 &alpha,
740 svd_u->values.data(),
741 &mm,
742 work.data(),
743 &one,
744 &beta,
745 w.data(),
746 &one);
747 break;
748 }
749 case inverse_svd:
750 {
751 std::scoped_lock lock(mutex);
752 AssertDimension(w.size(), this->n());
753 AssertDimension(v.size(), this->m());
754 // Compute U^T v
755 work.resize(std::max(mm, nn));
756 gemv("T",
757 &mm,
758 &mm,
759 &alpha,
760 svd_u->values.data(),
761 &mm,
762 v.data(),
763 &one,
764 &null,
765 work.data(),
766 &one);
767 // Multiply by singular values
768 for (size_type i = 0; i < wr.size(); ++i)
769 work[i] *= wr[i];
770 // Multiply with V
771 gemv("T",
772 &nn,
773 &nn,
774 &alpha,
775 svd_vt->values.data(),
776 &nn,
777 work.data(),
778 &one,
779 &beta,
780 w.data(),
781 &one);
782 break;
783 }
784 default:
785 Assert(false, ExcState(state));
786 }
787}
788
789
790
791template <typename number>
792void
793LAPACKFullMatrix<number>::Tvmult(Vector<number> &w,
794 const Vector<number> &v,
795 const bool adding) const
796{
797 const types::blas_int mm = this->m();
798 const types::blas_int nn = this->n();
799 const number alpha = 1.;
800 const number beta = (adding ? 1. : 0.);
801 const number null = 0.;
802
803 // use trmv for triangular matrices
804 if ((property == upper_triangular || property == lower_triangular) &&
805 (mm == nn) && state == matrix)
806 {
807 Assert(adding == false, ExcNotImplemented());
808
809 AssertDimension(v.size(), this->n());
810 AssertDimension(w.size(), this->m());
811
812 const char diag = 'N';
813 const char trans = 'T';
814 const char uplo =
815 (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
816
817 w = v;
818
819 const types::blas_int N = mm;
820 const types::blas_int lda = N;
821 const types::blas_int incx = 1;
822
823 trmv(
824 &uplo, &trans, &diag, &N, this->values.data(), &lda, w.data(), &incx);
825
826 return;
827 }
828
829
830 switch (state)
831 {
832 case matrix:
833 case inverse_matrix:
834 {
835 AssertDimension(w.size(), this->n());
836 AssertDimension(v.size(), this->m());
837
838 gemv("T",
839 &mm,
840 &nn,
841 &alpha,
842 this->values.data(),
843 &mm,
844 v.data(),
845 &one,
846 &beta,
847 w.data(),
848 &one);
849 break;
850 }
851 case svd:
852 {
853 std::scoped_lock lock(mutex);
854 AssertDimension(w.size(), this->n());
855 AssertDimension(v.size(), this->m());
856
857 // Compute U^T v
858 work.resize(std::max(mm, nn));
859 gemv("T",
860 &mm,
861 &mm,
862 &alpha,
863 svd_u->values.data(),
864 &mm,
865 v.data(),
866 &one,
867 &null,
868 work.data(),
869 &one);
870 // Multiply by singular values
871 for (size_type i = 0; i < wr.size(); ++i)
872 work[i] *= wr[i];
873 // Multiply with V
874 gemv("T",
875 &nn,
876 &nn,
877 &alpha,
878 svd_vt->values.data(),
879 &nn,
880 work.data(),
881 &one,
882 &beta,
883 w.data(),
884 &one);
885 break;
886 }
887 case inverse_svd:
888 {
889 std::scoped_lock lock(mutex);
890 AssertDimension(v.size(), this->n());
891 AssertDimension(w.size(), this->m());
892
893 // Compute V^T v
894 work.resize(std::max(mm, nn));
895 gemv("N",
896 &nn,
897 &nn,
898 &alpha,
899 svd_vt->values.data(),
900 &nn,
901 v.data(),
902 &one,
903 &null,
904 work.data(),
905 &one);
906 // Multiply by singular values
907 for (size_type i = 0; i < wr.size(); ++i)
908 work[i] *= wr[i];
909 // Multiply with U
910 gemv("N",
911 &mm,
912 &mm,
913 &alpha,
914 svd_u->values.data(),
915 &mm,
916 work.data(),
917 &one,
918 &beta,
919 w.data(),
920 &one);
921 break;
922 }
923 default:
924 Assert(false, ExcState(state));
925 }
926}
927
928
929
930template <typename number>
931void
932LAPACKFullMatrix<number>::vmult_add(Vector<number> &w,
933 const Vector<number> &v) const
934{
935 vmult(w, v, true);
936}
937
938
939
940template <typename number>
941void
942LAPACKFullMatrix<number>::Tvmult_add(Vector<number> &w,
943 const Vector<number> &v) const
944{
945 Tvmult(w, v, true);
946}
947
948
949
950template <typename number>
951void
952LAPACKFullMatrix<number>::mmult(LAPACKFullMatrix<number> &C,
953 const LAPACKFullMatrix<number> &B,
954 const bool adding) const
955{
956 Assert(state == matrix || state == inverse_matrix, ExcState(state));
957 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
958 Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state));
959 Assert(this->n() == B.m(), ExcDimensionMismatch(this->n(), B.m()));
960 Assert(C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
961 Assert(C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m()));
962 const types::blas_int mm = this->m();
963 const types::blas_int nn = B.n();
964 const types::blas_int kk = this->n();
965 const number alpha = 1.;
966 const number beta = (adding ? 1. : 0.);
967
968 gemm("N",
969 "N",
970 &mm,
971 &nn,
972 &kk,
973 &alpha,
974 this->values.data(),
975 &mm,
976 B.values.data(),
977 &kk,
978 &beta,
979 C.values.data(),
980 &mm);
981}
982
983
984
985template <typename number>
986void
987LAPACKFullMatrix<number>::mmult(FullMatrix<number> &C,
988 const LAPACKFullMatrix<number> &B,
989 const bool adding) const
990{
991 Assert(state == matrix || state == inverse_matrix, ExcState(state));
992 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
993 Assert(this->n() == B.m(), ExcDimensionMismatch(this->n(), B.m()));
994 Assert(C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
995 Assert(C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m()));
996 const types::blas_int mm = this->m();
997 const types::blas_int nn = B.n();
998 const types::blas_int kk = this->n();
999 const number alpha = 1.;
1000 const number beta = (adding ? 1. : 0.);
1001
1002 // since FullMatrix stores the matrix in transposed order compared to this
1003 // matrix, compute B^T * A^T = (A * B)^T
1004 gemm("T",
1005 "T",
1006 &nn,
1007 &mm,
1008 &kk,
1009 &alpha,
1010 B.values.data(),
1011 &kk,
1012 this->values.data(),
1013 &mm,
1014 &beta,
1015 &C(0, 0),
1016 &nn);
1017}
1018
1019
1020
1021template <typename number>
1022void
1023LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &C,
1024 const LAPACKFullMatrix<number> &B,
1025 const Vector<number> &V,
1026 const bool adding) const
1027{
1028 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1029 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1030 Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state));
1031
1032 const LAPACKFullMatrix<number> &A = *this;
1033
1034 Assert(A.m() == B.m(), ExcDimensionMismatch(A.m(), B.m()));
1035 Assert(C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
1036 Assert(C.m() == A.n(), ExcDimensionMismatch(A.n(), C.m()));
1037 Assert(V.size() == A.m(), ExcDimensionMismatch(A.m(), V.size()));
1038
1039 const types::blas_int mm = A.n();
1040 const types::blas_int nn = B.n();
1041 const types::blas_int kk = B.m();
1042
1043 // lapack does not have any triple product routines, including the case of
1044 // diagonal matrix in the middle, see
1045 // https://stackoverflow.com/questions/3548069/multiplying-three-matrices-in-blas-with-the-middle-one-being-diagonal
1046 // http://mathforum.org/kb/message.jspa?messageID=3546564
1047
1048 std::scoped_lock lock(mutex);
1049 // First, get V*B into "work" array
1050 work.resize(kk * nn);
1051 // following http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=2&t=768#p2577
1052 // do left-multiplication manually. Note that Xscal would require to first
1053 // copy the input vector as multiplication is done inplace.
1054 for (types::blas_int j = 0; j < nn; ++j)
1055 for (types::blas_int i = 0; i < kk; ++i)
1056 {
1057 Assert(j * kk + i < static_cast<types::blas_int>(work.size()),
1058 ExcInternalError());
1059 work[j * kk + i] = V(i) * B(i, j);
1060 }
1061
1062 // Now do the standard Tmmult:
1063 const number alpha = 1.;
1064 const number beta = (adding ? 1. : 0.);
1065
1066 gemm("T",
1067 "N",
1068 &mm,
1069 &nn,
1070 &kk,
1071 &alpha,
1072 this->values.data(),
1073 &kk,
1074 work.data(),
1075 &kk,
1076 &beta,
1077 C.values.data(),
1078 &mm);
1079}
1080
1081
1082
1083template <typename number>
1084void
1085LAPACKFullMatrix<number>::transpose(LAPACKFullMatrix<number> &B) const
1086{
1087 const LAPACKFullMatrix<number> &A = *this;
1088 AssertDimension(A.m(), B.n());
1089 AssertDimension(A.n(), B.m());
1090 const types::blas_int m = B.m();
1091 const types::blas_int n = B.n();
1092#ifdef DEAL_II_LAPACK_WITH_MKL
1093 const number one = 1.;
1094 omatcopy('C', 'C', n, m, one, A.values.data(), n, B.values.data(), m);
1095#else
1096 for (types::blas_int i = 0; i < m; ++i)
1097 for (types::blas_int j = 0; j < n; ++j)
1098 B(i, j) = numbers::NumberTraits<number>::conjugate(A(j, i));
1099#endif
1100}
1101
1102
1103
1104template <typename number>
1105void
1106LAPACKFullMatrix<number>::scale_rows(const Vector<number> &V)
1107{
1108 LAPACKFullMatrix<number> &A = *this;
1109 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1110 Assert(V.size() == A.m(), ExcDimensionMismatch(A.m(), V.size()));
1111
1112 const types::blas_int nn = A.n();
1113 const types::blas_int kk = A.m();
1114 for (types::blas_int j = 0; j < nn; ++j)
1115 for (types::blas_int i = 0; i < kk; ++i)
1116 A(i, j) *= V(i);
1117}
1118
1119
1120
1121template <typename number>
1122void
1123LAPACKFullMatrix<number>::Tmmult(LAPACKFullMatrix<number> &C,
1124 const LAPACKFullMatrix<number> &B,
1125 const bool adding) const
1126{
1127 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1128 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1129 Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state));
1130 Assert(this->m() == B.m(), ExcDimensionMismatch(this->m(), B.m()));
1131 Assert(C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
1132 Assert(C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m()));
1133 const types::blas_int mm = this->n();
1134 const types::blas_int nn = B.n();
1135 const types::blas_int kk = B.m();
1136 const number alpha = 1.;
1137 const number beta = (adding ? 1. : 0.);
1138
1139 if (PointerComparison::equal(this, &B))
1140 {
1141 syrk(&LAPACKSupport::U,
1142 "T",
1143 &nn,
1144 &kk,
1145 &alpha,
1146 this->values.data(),
1147 &kk,
1148 &beta,
1149 C.values.data(),
1150 &nn);
1151
1152 // fill-in lower triangular part
1153 for (types::blas_int j = 0; j < nn; ++j)
1154 for (types::blas_int i = 0; i < j; ++i)
1155 C(j, i) = C(i, j);
1156
1157 C.property = symmetric;
1158 }
1159 else
1160 {
1161 gemm("T",
1162 "N",
1163 &mm,
1164 &nn,
1165 &kk,
1166 &alpha,
1167 this->values.data(),
1168 &kk,
1169 B.values.data(),
1170 &kk,
1171 &beta,
1172 C.values.data(),
1173 &mm);
1174 }
1175}
1176
1177
1178
1179template <typename number>
1180void
1181LAPACKFullMatrix<number>::Tmmult(FullMatrix<number> &C,
1182 const LAPACKFullMatrix<number> &B,
1183 const bool adding) const
1184{
1185 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1186 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1187 Assert(this->m() == B.m(), ExcDimensionMismatch(this->m(), B.m()));
1188 Assert(C.n() == B.n(), ExcDimensionMismatch(C.n(), B.n()));
1189 Assert(C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m()));
1190 const types::blas_int mm = this->n();
1191 const types::blas_int nn = B.n();
1192 const types::blas_int kk = B.m();
1193 const number alpha = 1.;
1194 const number beta = (adding ? 1. : 0.);
1195
1196 // since FullMatrix stores the matrix in transposed order compared to this
1197 // matrix, compute B^T * A = (A^T * B)^T
1198 gemm("T",
1199 "N",
1200 &nn,
1201 &mm,
1202 &kk,
1203 &alpha,
1204 B.values.data(),
1205 &kk,
1206 this->values.data(),
1207 &kk,
1208 &beta,
1209 &C(0, 0),
1210 &nn);
1211}
1212
1213
1214
1215template <typename number>
1216void
1217LAPACKFullMatrix<number>::mTmult(LAPACKFullMatrix<number> &C,
1218 const LAPACKFullMatrix<number> &B,
1219 const bool adding) const
1220{
1221 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1222 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1223 Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state));
1224 Assert(this->n() == B.n(), ExcDimensionMismatch(this->n(), B.n()));
1225 Assert(C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m()));
1226 Assert(C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m()));
1227 const types::blas_int mm = this->m();
1228 const types::blas_int nn = B.m();
1229 const types::blas_int kk = B.n();
1230 const number alpha = 1.;
1231 const number beta = (adding ? 1. : 0.);
1232
1233 if (PointerComparison::equal(this, &B))
1234 {
1235 syrk(&LAPACKSupport::U,
1236 "N",
1237 &nn,
1238 &kk,
1239 &alpha,
1240 this->values.data(),
1241 &nn,
1242 &beta,
1243 C.values.data(),
1244 &nn);
1245
1246 // fill-in lower triangular part
1247 for (types::blas_int j = 0; j < nn; ++j)
1248 for (types::blas_int i = 0; i < j; ++i)
1249 C(j, i) = C(i, j);
1250
1251 C.property = symmetric;
1252 }
1253 else
1254 {
1255 gemm("N",
1256 "T",
1257 &mm,
1258 &nn,
1259 &kk,
1260 &alpha,
1261 this->values.data(),
1262 &mm,
1263 B.values.data(),
1264 &nn,
1265 &beta,
1266 C.values.data(),
1267 &mm);
1268 }
1269}
1270
1271
1272
1273template <typename number>
1274void
1275LAPACKFullMatrix<number>::mTmult(FullMatrix<number> &C,
1276 const LAPACKFullMatrix<number> &B,
1277 const bool adding) const
1278{
1279 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1280 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1281 Assert(this->n() == B.n(), ExcDimensionMismatch(this->n(), B.n()));
1282 Assert(C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m()));
1283 Assert(C.m() == this->m(), ExcDimensionMismatch(this->m(), C.m()));
1284 const types::blas_int mm = this->m();
1285 const types::blas_int nn = B.m();
1286 const types::blas_int kk = B.n();
1287 const number alpha = 1.;
1288 const number beta = (adding ? 1. : 0.);
1289
1290 // since FullMatrix stores the matrix in transposed order compared to this
1291 // matrix, compute B * A^T = (A * B^T)^T
1292 gemm("N",
1293 "T",
1294 &nn,
1295 &mm,
1296 &kk,
1297 &alpha,
1298 B.values.data(),
1299 &nn,
1300 this->values.data(),
1301 &mm,
1302 &beta,
1303 &C(0, 0),
1304 &nn);
1305}
1306
1307
1308
1309template <typename number>
1310void
1311LAPACKFullMatrix<number>::TmTmult(LAPACKFullMatrix<number> &C,
1312 const LAPACKFullMatrix<number> &B,
1313 const bool adding) const
1314{
1315 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1316 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1317 Assert(C.state == matrix || C.state == inverse_matrix, ExcState(C.state));
1318 Assert(this->m() == B.n(), ExcDimensionMismatch(this->m(), B.n()));
1319 Assert(C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m()));
1320 Assert(C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m()));
1321 const types::blas_int mm = this->n();
1322 const types::blas_int nn = B.m();
1323 const types::blas_int kk = B.n();
1324 const number alpha = 1.;
1325 const number beta = (adding ? 1. : 0.);
1326
1327 gemm("T",
1328 "T",
1329 &mm,
1330 &nn,
1331 &kk,
1332 &alpha,
1333 this->values.data(),
1334 &kk,
1335 B.values.data(),
1336 &nn,
1337 &beta,
1338 C.values.data(),
1339 &mm);
1340}
1341
1342
1343
1344template <typename number>
1345void
1346LAPACKFullMatrix<number>::TmTmult(FullMatrix<number> &C,
1347 const LAPACKFullMatrix<number> &B,
1348 const bool adding) const
1349{
1350 Assert(state == matrix || state == inverse_matrix, ExcState(state));
1351 Assert(B.state == matrix || B.state == inverse_matrix, ExcState(B.state));
1352 Assert(this->m() == B.n(), ExcDimensionMismatch(this->m(), B.n()));
1353 Assert(C.n() == B.m(), ExcDimensionMismatch(C.n(), B.m()));
1354 Assert(C.m() == this->n(), ExcDimensionMismatch(this->n(), C.m()));
1355 const types::blas_int mm = this->n();
1356 const types::blas_int nn = B.m();
1357 const types::blas_int kk = B.n();
1358 const number alpha = 1.;
1359 const number beta = (adding ? 1. : 0.);
1360
1361 // since FullMatrix stores the matrix in transposed order compared to this
1362 // matrix, compute B * A = (A^T * B^T)^T
1363 gemm("N",
1364 "N",
1365 &nn,
1366 &mm,
1367 &kk,
1368 &alpha,
1369 B.values.data(),
1370 &nn,
1371 this->values.data(),
1372 &kk,
1373 &beta,
1374 &C(0, 0),
1375 &nn);
1376}
1377
1378
1379
1380template <typename number>
1381void
1382LAPACKFullMatrix<number>::compute_lu_factorization()
1383{
1384 Assert(state == matrix, ExcState(state));
1385 state = LAPACKSupport::unusable;
1386
1387 const types::blas_int mm = this->m();
1388 const types::blas_int nn = this->n();
1389 number *const values = this->values.data();
1390 ipiv.resize(mm);
1391 types::blas_int info = 0;
1392 getrf(&mm, &nn, values, &mm, ipiv.data(), &info);
1393
1394 Assert(info >= 0, ExcInternalError());
1395
1396 // if info >= 0, the factorization has been completed
1397 state = lu;
1398
1399 AssertThrow(info == 0, LACExceptions::ExcSingular());
1400}
1401
1402
1403
1404template <typename number>
1405void
1406LAPACKFullMatrix<number>::set_property(const Property p)
1407{
1408 property = p;
1409}
1410
1411
1412
1413template <typename number>
1414number
1415LAPACKFullMatrix<number>::l1_norm() const
1416{
1417 const char type('O');
1418 return norm(type);
1419}
1420
1421
1422
1423template <typename number>
1424number
1425LAPACKFullMatrix<number>::linfty_norm() const
1426{
1427 const char type('I');
1428 return norm(type);
1429}
1430
1431
1432
1433template <typename number>
1434number
1435LAPACKFullMatrix<number>::frobenius_norm() const
1436{
1437 const char type('F');
1438 return norm(type);
1439}
1440
1441
1442
1443template <typename number>
1444number
1445LAPACKFullMatrix<number>::norm(const char type) const
1446{
1447 std::scoped_lock lock(mutex);
1448
1449 Assert(state == LAPACKSupport::matrix ||
1450 state == LAPACKSupport::inverse_matrix,
1451 ExcMessage("norms can be called in matrix state only."));
1452
1453 const types::blas_int N = this->n();
1454 const types::blas_int M = this->m();
1455 const number *const values = this->values.data();
1456 if (property == symmetric)
1457 {
1458 const types::blas_int lda = std::max<types::blas_int>(1, N);
1459 const types::blas_int lwork =
1460 (type == 'I' || type == 'O') ? std::max<types::blas_int>(1, N) : 0;
1461 work.resize(lwork);
1462 return lansy(&type, &LAPACKSupport::L, &N, values, &lda, work.data());
1463 }
1464 else
1465 {
1466 const types::blas_int lda = std::max<types::blas_int>(1, M);
1467 const types::blas_int lwork =
1468 (type == 'I') ? std::max<types::blas_int>(1, M) : 0;
1469 work.resize(lwork);
1470 return lange(&type, &M, &N, values, &lda, work.data());
1471 }
1472}
1473
1474
1475
1476template <typename number>
1477number
1478LAPACKFullMatrix<number>::trace() const
1479{
1480 Assert(state == LAPACKSupport::matrix ||
1481 state == LAPACKSupport::inverse_matrix,
1482 ExcMessage("Trace can be called in matrix state only."));
1483 Assert(this->n() == this->m(), ExcDimensionMismatch(this->n(), this->m()));
1484
1485 number tr = 0;
1486 for (size_type i = 0; i < this->m(); ++i)
1487 tr += (*this)(i, i);
1488
1489 return tr;
1490}
1491
1492
1493
1494template <typename number>
1495void
1496LAPACKFullMatrix<number>::compute_cholesky_factorization()
1497{
1498 Assert(state == matrix, ExcState(state));
1499 Assert(property == symmetric, ExcProperty(property));
1500 state = LAPACKSupport::unusable;
1501
1502 const types::blas_int mm = this->m();
1503 const types::blas_int nn = this->n();
1504 Assert(mm == nn, ExcDimensionMismatch(mm, nn));
1505
1506 number *const values = this->values.data();
1507 types::blas_int info = 0;
1508 const types::blas_int lda = std::max<types::blas_int>(1, nn);
1509 potrf(&LAPACKSupport::L, &nn, values, &lda, &info);
1510
1511 // info < 0 : the info-th argument had an illegal value
1512 Assert(info >= 0, ExcInternalError());
1513
1514 state = cholesky;
1515 AssertThrow(info == 0, LACExceptions::ExcSingular());
1516}
1517
1518
1519
1520template <typename number>
1521number
1522LAPACKFullMatrix<number>::reciprocal_condition_number(const number a_norm) const
1523{
1524 std::scoped_lock lock(mutex);
1525 Assert(state == cholesky, ExcState(state));
1526 number rcond = 0.;
1527
1528 const types::blas_int N = this->m();
1529 const number *values = this->values.data();
1530 types::blas_int info = 0;
1531 const types::blas_int lda = std::max<types::blas_int>(1, N);
1532 work.resize(3 * N);
1533 iwork.resize(N);
1534
1535 // use the same uplo as in Cholesky
1536 pocon(&LAPACKSupport::L,
1537 &N,
1538 values,
1539 &lda,
1540 &a_norm,
1541 &rcond,
1542 work.data(),
1543 iwork.data(),
1544 &info);
1545
1546 Assert(info >= 0, ExcInternalError());
1547
1548 return rcond;
1549}
1550
1551
1552
1553template <typename number>
1554number
1555LAPACKFullMatrix<number>::reciprocal_condition_number() const
1556{
1557 std::scoped_lock lock(mutex);
1558 Assert(property == upper_triangular || property == lower_triangular,
1559 ExcProperty(property));
1560 number rcond = 0.;
1561
1562 const types::blas_int N = this->m();
1563 const number *const values = this->values.data();
1564 types::blas_int info = 0;
1565 const types::blas_int lda = std::max<types::blas_int>(1, N);
1566 work.resize(3 * N);
1567 iwork.resize(N);
1568
1569 const char norm = '1';
1570 const char diag = 'N';
1571 const char uplo =
1572 (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
1573 trcon(&norm,
1574 &uplo,
1575 &diag,
1576 &N,
1577 values,
1578 &lda,
1579 &rcond,
1580 work.data(),
1581 iwork.data(),
1582 &info);
1583
1584 Assert(info >= 0, ExcInternalError());
1585
1586 return rcond;
1587}
1588
1589
1590
1591template <typename number>
1592void
1593LAPACKFullMatrix<number>::compute_svd()
1594{
1595 Assert(state == matrix, ExcState(state));
1596 state = LAPACKSupport::unusable;
1597
1598 const types::blas_int mm = this->m();
1599 const types::blas_int nn = this->n();
1600 wr.resize(std::max(mm, nn));
1601 std::fill(wr.begin(), wr.end(), 0.);
1602 ipiv.resize(8 * mm);
1603
1604 svd_u = std::make_unique<LAPACKFullMatrix<number>>(mm, mm);
1605 svd_vt = std::make_unique<LAPACKFullMatrix<number>>(nn, nn);
1606 types::blas_int info = 0;
1607
1608 // First determine optimal workspace size
1609 work.resize(1);
1610 types::blas_int lwork = -1;
1611
1612 // TODO double check size
1613 std::vector<typename numbers::NumberTraits<number>::real_type> real_work;
1614 if (numbers::NumberTraits<number>::is_complex)
1615 {
1616 // This array is only used by the complex versions.
1617 std::size_t min = std::min(this->m(), this->n());
1618 std::size_t max = std::max(this->m(), this->n());
1619 real_work.resize(
1620 std::max(5 * min * min + 5 * min, 2 * max * min + 2 * min * min + min));
1621 }
1622
1623 // make sure that the first entry in the work array is clear, in case the
1624 // routine does not completely overwrite the memory:
1625 work[0] = number();
1626 internal::LAPACKFullMatrixImplementation::gesdd_helper(LAPACKSupport::A,
1627 mm,
1628 nn,
1629 this->values,
1630 wr,
1631 svd_u->values,
1632 svd_vt->values,
1633 work,
1634 real_work,
1635 ipiv,
1636 lwork,
1637 info);
1638
1639 AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("gesdd", info));
1640 // Resize the work array. Add one to the size computed by LAPACK to be on
1641 // the safe side.
1642 lwork = static_cast<types::blas_int>(std::abs(work[0]) + 1);
1643
1644 work.resize(lwork);
1645 // Do the actual SVD.
1646 internal::LAPACKFullMatrixImplementation::gesdd_helper(LAPACKSupport::A,
1647 mm,
1648 nn,
1649 this->values,
1650 wr,
1651 svd_u->values,
1652 svd_vt->values,
1653 work,
1654 real_work,
1655 ipiv,
1656 lwork,
1657 info);
1658 AssertThrow(info == 0, LAPACKSupport::ExcErrorCode("gesdd", info));
1659
1660 work.resize(0);
1661 ipiv.resize(0);
1662
1663 state = LAPACKSupport::svd;
1664}
1665
1666
1667
1668template <typename number>
1669void
1670LAPACKFullMatrix<number>::compute_inverse_svd(const double threshold)
1671{
1672 if (state == LAPACKSupport::matrix)
1673 compute_svd();
1674
1675 Assert(state == LAPACKSupport::svd, ExcState(state));
1676
1677 const typename numbers::NumberTraits<number>::real_type one(1.0);
1678 const double lim = std::abs(wr[0]) * threshold;
1679 for (size_type i = 0; i < wr.size(); ++i)
1680 {
1681 if (std::abs(wr[i]) > lim)
1682 wr[i] = one / wr[i];
1683 else
1684 wr[i] = 0.;
1685 }
1686 state = LAPACKSupport::inverse_svd;
1687}
1688
1689
1690
1691template <typename number>
1692void
1693LAPACKFullMatrix<number>::compute_inverse_svd_with_kernel(
1694 const unsigned int kernel_size)
1695{
1696 if (state == LAPACKSupport::matrix)
1697 compute_svd();
1698
1699 Assert(state == LAPACKSupport::svd, ExcState(state));
1700
1701 const typename numbers::NumberTraits<number>::real_type one(1.0);
1702 const unsigned int n_wr = wr.size();
1703 for (size_type i = 0; i < n_wr - kernel_size; ++i)
1704 wr[i] = one / wr[i];
1705 for (size_type i = n_wr - kernel_size; i < n_wr; ++i)
1706 wr[i] = 0.;
1707 state = LAPACKSupport::inverse_svd;
1708}
1709
1710
1711
1712template <typename number>
1713void
1714LAPACKFullMatrix<number>::invert()
1715{
1716 Assert(state == matrix || state == lu || state == cholesky, ExcState(state));
1717 const types::blas_int mm = this->m();
1718 const types::blas_int nn = this->n();
1719 Assert(nn == mm, ExcNotQuadratic());
1720
1721 number *const values = this->values.data();
1722 types::blas_int info = 0;
1723
1724 if (property != symmetric)
1725 {
1726 if (state == matrix)
1727 compute_lu_factorization();
1728
1729 ipiv.resize(mm);
1730 inv_work.resize(mm);
1731 getri(&mm, values, &mm, ipiv.data(), inv_work.data(), &mm, &info);
1732 }
1733 else
1734 {
1735 if (state == matrix)
1736 compute_cholesky_factorization();
1737
1738 const types::blas_int lda = std::max<types::blas_int>(1, nn);
1739 potri(&LAPACKSupport::L, &nn, values, &lda, &info);
1740 // inverse is stored in lower diagonal, set the upper diagonal
1741 // appropriately:
1742 for (types::blas_int i = 0; i < nn; ++i)
1743 for (types::blas_int j = i + 1; j < nn; ++j)
1744 this->el(i, j) = this->el(j, i);
1745 }
1746
1747 Assert(info >= 0, ExcInternalError());
1748 AssertThrow(info == 0, LACExceptions::ExcSingular());
1749
1750 state = inverse_matrix;
1751}
1752
1753
1754
1755template <typename number>
1756void
1757LAPACKFullMatrix<number>::solve(Vector<number> &v, const bool transposed) const
1758{
1759 Assert(this->m() == this->n(), LACExceptions::ExcNotQuadratic());
1760 AssertDimension(this->m(), v.size());
1761 const char *trans = transposed ? &T : &N;
1762 const types::blas_int nn = this->n();
1763 const number *const values = this->values.data();
1764 const types::blas_int n_rhs = 1;
1765 types::blas_int info = 0;
1766
1767 if (state == lu)
1768 {
1769 getrs(
1770 trans, &nn, &n_rhs, values, &nn, ipiv.data(), v.begin(), &nn, &info);
1771 }
1772 else if (state == cholesky)
1773 {
1774 potrs(&LAPACKSupport::L, &nn, &n_rhs, values, &nn, v.begin(), &nn, &info);
1775 }
1776 else if (property == upper_triangular || property == lower_triangular)
1777 {
1778 const char uplo =
1779 (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
1780
1781 const types::blas_int lda = nn;
1782 const types::blas_int ldb = nn;
1783 trtrs(
1784 &uplo, trans, "N", &nn, &n_rhs, values, &lda, v.begin(), &ldb, &info);
1785 }
1786 else
1787 {
1788 Assert(false,
1789 ExcMessage(
1790 "The matrix has to be either factorized or triangular."));
1791 }
1792
1793 Assert(info == 0, ExcInternalError());
1794}
1795
1796
1797
1798template <typename number>
1799void
1800LAPACKFullMatrix<number>::solve(LAPACKFullMatrix<number> &B,
1801 const bool transposed) const
1802{
1803 Assert(B.state == matrix, ExcState(B.state));
1804
1805 Assert(this->m() == this->n(), LACExceptions::ExcNotQuadratic());
1806 AssertDimension(this->m(), B.m());
1807 const char *trans = transposed ? &T : &N;
1808 const types::blas_int nn = this->n();
1809 const number *const values = this->values.data();
1810 const types::blas_int n_rhs = B.n();
1811 types::blas_int info = 0;
1812
1813 if (state == lu)
1814 {
1815 getrs(trans,
1816 &nn,
1817 &n_rhs,
1818 values,
1819 &nn,
1820 ipiv.data(),
1821 B.values.data(),
1822 &nn,
1823 &info);
1824 }
1825 else if (state == cholesky)
1826 {
1827 potrs(&LAPACKSupport::L,
1828 &nn,
1829 &n_rhs,
1830 values,
1831 &nn,
1832 B.values.data(),
1833 &nn,
1834 &info);
1835 }
1836 else if (property == upper_triangular || property == lower_triangular)
1837 {
1838 const char uplo =
1839 (property == upper_triangular ? LAPACKSupport::U : LAPACKSupport::L);
1840
1841 const types::blas_int lda = nn;
1842 const types::blas_int ldb = nn;
1843 trtrs(&uplo,
1844 trans,
1845 "N",
1846 &nn,
1847 &n_rhs,
1848 values,
1849 &lda,
1850 B.values.data(),
1851 &ldb,
1852 &info);
1853 }
1854 else
1855 {
1856 Assert(false,
1857 ExcMessage(
1858 "The matrix has to be either factorized or triangular."));
1859 }
1860
1861 Assert(info == 0, ExcInternalError());
1862}
1863
1864
1865
1866template <typename number>
1867number
1868LAPACKFullMatrix<number>::determinant() const
1869{
1870 Assert(this->m() == this->n(), LACExceptions::ExcNotQuadratic());
1871
1872 // LAPACK doesn't offer a function to compute a matrix determinant.
1873 // This is due to the difficulty in maintaining numerical accuracy, as the
1874 // calculations are likely to overflow or underflow. See
1875 // http://www.netlib.org/lapack/faq.html#_are_there_routines_in_lapack_to_compute_determinants
1876 //
1877 // However, after a PLU decomposition one can compute this by multiplication
1878 // of the diagonal entries with one another. One must take into consideration
1879 // the number of permutations (row swaps) stored in the P matrix.
1880 //
1881 // See the implementations in the blaze library (detNxN)
1882 // https://bitbucket.org/blaze-lib/blaze
1883 // and also
1884 // https://dualm.wordpress.com/2012/01/06/computing-determinant-in-fortran/
1885 // http://icl.cs.utk.edu/lapack-forum/viewtopic.php?p=341&#p336
1886 // for further information.
1887 Assert(state == lu, ExcState(state));
1888 Assert(ipiv.size() == this->m(), ExcInternalError());
1889 number det = 1.0;
1890 for (size_type i = 0; i < this->m(); ++i)
1891 det *=
1892 (ipiv[i] == types::blas_int(i + 1) ? this->el(i, i) : -this->el(i, i));
1893 return det;
1894}
1895
1896
1897
1898template <typename number>
1899void
1900LAPACKFullMatrix<number>::compute_eigenvalues(const bool right, const bool left)
1901{
1902 Assert(state == matrix, ExcState(state));
1903 const types::blas_int nn = this->n();
1904 wr.resize(nn);
1905 wi.resize(nn);
1906 if (right)
1907 vr.resize(nn * nn);
1908 if (left)
1909 vl.resize(nn * nn);
1910
1911 types::blas_int info = 0;
1912 types::blas_int lwork = 1;
1913 const char jobvr = (right) ? V : N;
1914 const char jobvl = (left) ? V : N;
1915
1916 /*
1917 * The LAPACK routine xGEEV requires a sufficiently large work array; the
1918 * minimum requirement is
1919 *
1920 * work.size >= 4*nn.
1921 *
1922 * However, for better performance, a larger work array may be needed. The
1923 * first call determines the optimal work size and the second does the work.
1924 */
1925 lwork = -1;
1926 work.resize(1);
1927
1928 std::vector<typename numbers::NumberTraits<number>::real_type> real_work;
1929 if (numbers::NumberTraits<number>::is_complex)
1930 // This array is only used by the complex versions.
1931 real_work.resize(2 * this->m());
1932 internal::LAPACKFullMatrixImplementation::geev_helper(jobvl,
1933 jobvr,
1934 this->values,
1935 this->m(),
1936 wr,
1937 wi,
1938 vl,
1939 vr,
1940 work,
1941 real_work,
1942 lwork,
1943 info);
1944
1945 // geev returns info=0 on success. Since we only queried the optimal size
1946 // for work, everything else would not be acceptable.
1947 Assert(info == 0, ExcInternalError());
1948 // Allocate working array according to suggestion (same strategy as was
1949 // noted in compute_svd).
1950 lwork = static_cast<types::blas_int>(std::abs(work[0]) + 1);
1951
1952 // resize workspace array
1953 work.resize(lwork);
1954 real_work.resize(lwork);
1955
1956 // Finally compute the eigenvalues.
1957 internal::LAPACKFullMatrixImplementation::geev_helper(jobvl,
1958 jobvr,
1959 this->values,
1960 this->m(),
1961 wr,
1962 wi,
1963 vl,
1964 vr,
1965 work,
1966 real_work,
1967 lwork,
1968 info);
1969
1970 Assert(info >= 0, ExcInternalError());
1971 if (info < 0)
1972 {
1973 AssertThrow(info == 0,
1974 ExcMessage("Lapack error in geev: the " +
1975 std::to_string(-info) +
1976 "-th"
1977 " parameter had an illegal value."));
1978 }
1979 else
1980 {
1981 AssertThrow(
1982 info == 0,
1983 ExcMessage(
1984 "Lapack error in geev: the QR algorithm failed to compute "
1985 "all the eigenvalues, and no eigenvectors have been computed."));
1986 }
1987
1988 state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
1989}
1990
1991
1992
1993namespace
1994{
1995 // This function extracts complex eigenvectors from the underlying 'number'
1996 // array 'vr' of the LAPACK eigenvalue routine. For real-valued matrices
1997 // addressed by this function specialization, we might get complex
1998 // eigenvalues, which come in complex-conjugate pairs. In LAPACK, a compact
1999 // storage scheme is applied that stores the real and imaginary part of
2000 // eigenvectors only once, putting the real parts in one column and the
2001 // imaginary part in the next of a real-valued array. Here, we do the
2002 // unpacking into the usual complex values.
2003 template <typename RealNumber>
2004 void
2005 unpack_lapack_eigenvector_and_increment_index(
2006 const std::vector<RealNumber> &vr,
2007 const std::complex<RealNumber> &eigenvalue,
2008 FullMatrix<std::complex<RealNumber>> &result,
2009 unsigned int &index)
2010 {
2011 const std::size_t n = result.n();
2012 if (eigenvalue.imag() != 0.)
2013 {
2014 for (std::size_t j = 0; j < n; ++j)
2015 {
2016 result(j, index).real(vr[index * n + j]);
2017 result(j, index + 1).real(vr[index * n + j]);
2018 result(j, index).imag(vr[(index + 1) * n + j]);
2019 result(j, index + 1).imag(-vr[(index + 1) * n + j]);
2020 }
2021
2022 // we filled two columns with the complex-conjugate pair, so increment
2023 // returned index by 2
2024 index += 2;
2025 }
2026 else
2027 {
2028 for (unsigned int j = 0; j < n; ++j)
2029 result(j, index).real(vr[index * n + j]);
2030
2031 // real-valued case, we only filled one column
2032 ++index;
2033 }
2034 }
2035
2036 // This specialization fills the eigenvectors for complex-valued matrices,
2037 // in which case we simply read off the entry in the 'vr' array.
2038 template <typename ComplexNumber>
2039 void
2040 unpack_lapack_eigenvector_and_increment_index(
2041 const std::vector<ComplexNumber> &vr,
2042 const ComplexNumber &,
2043 FullMatrix<ComplexNumber> &result,
2044 unsigned int &index)
2045 {
2046 const std::size_t n = result.n();
2047 for (unsigned int j = 0; j < n; ++j)
2048 result(j, index) = vr[index * n + j];
2049
2050 // complex-valued case always only fills a single column
2051 ++index;
2052 }
2053} // namespace
2054
2055
2056
2057template <typename number>
2058FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
2059LAPACKFullMatrix<number>::get_right_eigenvectors() const
2060{
2061 Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
2062 Assert(vr.size() == this->n_rows() * this->n_cols(),
2063 ExcMessage("Right eigenvectors are not available! Did you "
2064 "set the associated flag in compute_eigenvalues()?"));
2065
2066 FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
2067 result(n(), n());
2068
2069 for (unsigned int i = 0; i < n();)
2070 unpack_lapack_eigenvector_and_increment_index(vr, eigenvalue(i), result, i);
2071
2072 return result;
2073}
2074
2075
2076
2077template <typename number>
2078FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
2079LAPACKFullMatrix<number>::get_left_eigenvectors() const
2080{
2081 Assert(state & LAPACKSupport::eigenvalues, ExcInvalidState());
2082 Assert(vl.size() == this->n_rows() * this->n_cols(),
2083 ExcMessage("Left eigenvectors are not available! Did you "
2084 "set the associated flag in compute_eigenvalues()?"));
2085
2086 FullMatrix<std::complex<typename numbers::NumberTraits<number>::real_type>>
2087 result(n(), n());
2088
2089 for (unsigned int i = 0; i < n();)
2090 unpack_lapack_eigenvector_and_increment_index(vl, eigenvalue(i), result, i);
2091
2092 return result;
2093}
2094
2095
2096
2097template <typename number>
2098void
2099LAPACKFullMatrix<number>::compute_eigenvalues_symmetric(
2100 const number lower_bound,
2101 const number upper_bound,
2102 const number abs_accuracy,
2103 Vector<number> &eigenvalues,
2104 FullMatrix<number> &eigenvectors)
2105{
2106 Assert(state == matrix, ExcState(state));
2107 const types::blas_int nn = (this->n() > 0 ? this->n() : 1);
2108 Assert(static_cast<size_type>(nn) == this->m(), ExcNotQuadratic());
2109
2110 wr.resize(nn);
2111 LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
2112
2113 number *const values_A = this->values.data();
2114 number *const values_eigenvectors = matrix_eigenvectors.values.data();
2115
2116 types::blas_int info(0), lwork(-1), n_eigenpairs(0);
2117 const char *const jobz(&V);
2118 const char *const uplo(&U);
2119 const char *const range(&V);
2120 const types::blas_int *const dummy(&one);
2121 std::vector<types::blas_int> iwork(static_cast<size_type>(5 * nn));
2122 std::vector<types::blas_int> ifail(static_cast<size_type>(nn));
2123
2124
2125 /*
2126 * The LAPACK routine xSYEVX requires a sufficiently large work array; the
2127 * minimum requirement is
2128 *
2129 * work.size >= 8*nn.
2130 *
2131 * However, for better performance, a larger work array may be needed. The
2132 * first call determines the optimal work size and the second does the work.
2133 */
2134 work.resize(1);
2135
2136 syevx(jobz,
2137 range,
2138 uplo,
2139 &nn,
2140 values_A,
2141 &nn,
2142 &lower_bound,
2143 &upper_bound,
2144 dummy,
2145 dummy,
2146 &abs_accuracy,
2147 &n_eigenpairs,
2148 wr.data(),
2149 values_eigenvectors,
2150 &nn,
2151 work.data(),
2152 &lwork,
2153 iwork.data(),
2154 ifail.data(),
2155 &info);
2156 // syevx returns info=0 on success. Since we only queried the optimal size
2157 // for work, everything else would not be acceptable.
2158 Assert(info == 0, ExcInternalError());
2159 // Allocate working array according to suggestion (same strategy as was noted
2160 // in compute_svd).
2161 lwork = static_cast<types::blas_int>(std::abs(work[0]) + 1);
2162 work.resize(static_cast<size_type>(lwork));
2163
2164 // Finally compute the eigenvalues.
2165 syevx(jobz,
2166 range,
2167 uplo,
2168 &nn,
2169 values_A,
2170 &nn,
2171 &lower_bound,
2172 &upper_bound,
2173 dummy,
2174 dummy,
2175 &abs_accuracy,
2176 &n_eigenpairs,
2177 wr.data(),
2178 values_eigenvectors,
2179 &nn,
2180 work.data(),
2181 &lwork,
2182 iwork.data(),
2183 ifail.data(),
2184 &info);
2185
2186 // Negative return value implies a wrong argument. This should be internal.
2187 Assert(info >= 0, ExcInternalError());
2188 if (info < 0)
2189 {
2190 AssertThrow(info == 0,
2191 ExcMessage("Lapack error in syevx: the " +
2192 std::to_string(-info) +
2193 "-th"
2194 " parameter had an illegal value."));
2195 }
2196 else if ((info > 0) && (info <= nn))
2197 {
2198 AssertThrow(info == 0,
2199 ExcMessage(
2200 "Lapack error in syevx: " + std::to_string(info) +
2201 " eigenvectors failed to converge."
2202 " (You may need to scale the abs_accuracy according"
2203 " to your matrix norm.)"));
2204 }
2205 else
2206 {
2207 AssertThrow(info == 0,
2208 ExcMessage("Lapack error in syevx: unknown error."));
2209 }
2210
2211 eigenvalues.reinit(n_eigenpairs);
2212 eigenvectors.reinit(nn, n_eigenpairs, true);
2213
2214 for (size_type i = 0; i < static_cast<size_type>(n_eigenpairs); ++i)
2215 {
2216 eigenvalues(i) = wr[i];
2217 size_type col_begin(i * nn);
2218 for (size_type j = 0; j < static_cast<size_type>(nn); ++j)
2219 {
2220 eigenvectors(j, i) = values_eigenvectors[col_begin + j];
2221 }
2222 }
2223
2224 state = LAPACKSupport::State(unusable);
2225}
2226
2227
2228
2229template <typename number>
2230void
2231LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
2232 LAPACKFullMatrix<number> &B,
2233 const number lower_bound,
2234 const number upper_bound,
2235 const number abs_accuracy,
2236 Vector<number> &eigenvalues,
2237 std::vector<Vector<number>> &eigenvectors,
2238 const types::blas_int itype)
2239{
2240 Assert(state == matrix, ExcState(state));
2241 const types::blas_int nn = (this->n() > 0 ? this->n() : 1);
2242 Assert(static_cast<size_type>(nn) == this->m(), ExcNotQuadratic());
2243 Assert(B.m() == B.n(), ExcNotQuadratic());
2244 Assert(static_cast<size_type>(nn) == B.n(), ExcDimensionMismatch(nn, B.n()));
2245
2246 wr.resize(nn);
2247 LAPACKFullMatrix<number> matrix_eigenvectors(nn, nn);
2248
2249 number *const values_A = this->values.data();
2250 number *const values_B = B.values.data();
2251 number *const values_eigenvectors = matrix_eigenvectors.values.data();
2252
2253 types::blas_int info(0), lwork(-1), n_eigenpairs(0);
2254 const char *const jobz(&V);
2255 const char *const uplo(&U);
2256 const char *const range(&V);
2257 const types::blas_int *const dummy(&one);
2258 iwork.resize(static_cast<size_type>(5 * nn));
2259 std::vector<types::blas_int> ifail(static_cast<size_type>(nn));
2260
2261
2262 /*
2263 * The LAPACK routine xSYGVX requires a sufficiently large work array; the
2264 * minimum requirement is
2265 *
2266 * work.size >= 8*nn.
2267 *
2268 * However, for better performance, a larger work array may be needed. The
2269 * first call determines the optimal work size and the second does the work.
2270 */
2271 work.resize(1);
2272
2273 sygvx(&itype,
2274 jobz,
2275 range,
2276 uplo,
2277 &nn,
2278 values_A,
2279 &nn,
2280 values_B,
2281 &nn,
2282 &lower_bound,
2283 &upper_bound,
2284 dummy,
2285 dummy,
2286 &abs_accuracy,
2287 &n_eigenpairs,
2288 wr.data(),
2289 values_eigenvectors,
2290 &nn,
2291 work.data(),
2292 &lwork,
2293 iwork.data(),
2294 ifail.data(),
2295 &info);
2296 // sygvx returns info=0 on success. Since we only queried the optimal size
2297 // for work, everything else would not be acceptable.
2298 Assert(info == 0, ExcInternalError());
2299 // Allocate working array according to suggestion (same strategy as was
2300 // noted in compute_svd).
2301 lwork = static_cast<types::blas_int>(std::abs(work[0]) + 1);
2302
2303 // resize workspace arrays
2304 work.resize(static_cast<size_type>(lwork));
2305
2306 // Finally compute the generalized eigenvalues.
2307 sygvx(&itype,
2308 jobz,
2309 range,
2310 uplo,
2311 &nn,
2312 values_A,
2313 &nn,
2314 values_B,
2315 &nn,
2316 &lower_bound,
2317 &upper_bound,
2318 dummy,
2319 dummy,
2320 &abs_accuracy,
2321 &n_eigenpairs,
2322 wr.data(),
2323 values_eigenvectors,
2324 &nn,
2325 work.data(),
2326 &lwork,
2327 iwork.data(),
2328 ifail.data(),
2329 &info);
2330
2331 // Negative return value implies a wrong argument. This should be internal.
2332 Assert(info >= 0, ExcInternalError());
2333 if (info < 0)
2334 {
2335 AssertThrow(info == 0,
2336 ExcMessage("Lapack error in sygvx: the " +
2337 std::to_string(-info) +
2338 "-th"
2339 " parameter had an illegal value."));
2340 }
2341 else if ((info > 0) && (info <= nn))
2342 {
2343 AssertThrow(
2344 info == 0,
2345 ExcMessage(
2346 "Lapack error in sygvx: ssyevx/dsyevx failed to converge, and " +
2347 std::to_string(info) +
2348 " eigenvectors failed to converge."
2349 " (You may need to scale the abs_accuracy"
2350 " according to the norms of matrices A and B.)"));
2351 }
2352 else if ((info > nn) && (info <= 2 * nn))
2353 {
2354 AssertThrow(info == 0,
2355 ExcMessage(
2356 "Lapack error in sygvx: the leading minor of order " +
2357 std::to_string(info - nn) +
2358 " of matrix B is not positive-definite."
2359 " The factorization of B could not be completed and"
2360 " no eigenvalues or eigenvectors were computed."));
2361 }
2362 else
2363 {
2364 AssertThrow(info == 0,
2365 ExcMessage("Lapack error in sygvx: unknown error."));
2366 }
2367
2368 eigenvalues.reinit(n_eigenpairs);
2369 eigenvectors.resize(n_eigenpairs);
2370
2371 for (size_type i = 0; i < static_cast<size_type>(n_eigenpairs); ++i)
2372 {
2373 eigenvalues(i) = wr[i];
2374 size_type col_begin(i * nn);
2375 eigenvectors[i].reinit(nn, true);
2376 for (size_type j = 0; j < static_cast<size_type>(nn); ++j)
2377 {
2378 eigenvectors[i](j) = values_eigenvectors[col_begin + j];
2379 }
2380 }
2381
2382 state = LAPACKSupport::State(unusable);
2383}
2384
2385
2386
2387template <typename number>
2388void
2389LAPACKFullMatrix<number>::compute_generalized_eigenvalues_symmetric(
2390 LAPACKFullMatrix<number> &B,
2391 std::vector<Vector<number>> &eigenvectors,
2392 const types::blas_int itype)
2393{
2394 Assert(state == matrix, ExcState(state));
2395 const types::blas_int nn = this->n();
2396 Assert(static_cast<size_type>(nn) == this->m(), ExcNotQuadratic());
2397 Assert(B.m() == B.n(), ExcNotQuadratic());
2398 Assert(static_cast<size_type>(nn) == B.n(), ExcDimensionMismatch(nn, B.n()));
2399 Assert(eigenvectors.size() <= static_cast<size_type>(nn),
2400 ExcMessage("eigenvectors.size() > matrix.n()"));
2401
2402 wr.resize(nn);
2403 wi.resize(nn); // This is set purely for consistency reasons with the
2404 // eigenvalues() function.
2405
2406 number *const values_A = this->values.data();
2407 number *const values_B = B.values.data();
2408
2409 types::blas_int info = 0;
2410 types::blas_int lwork = -1;
2411 const char *const jobz = (eigenvectors.size() > 0) ? (&V) : (&N);
2412 const char *const uplo = (&U);
2413
2414 /*
2415 * The LAPACK routine xSYGV requires a sufficiently large work array; the
2416 * minimum requirement is
2417 *
2418 * work.size >= 3*nn - 1.
2419 *
2420 * However, for better performance, a larger work array may be needed. The
2421 * first call determines the optimal work size and the second does the work.
2422 */
2423 work.resize(1);
2424
2425 sygv(&itype,
2426 jobz,
2427 uplo,
2428 &nn,
2429 values_A,
2430 &nn,
2431 values_B,
2432 &nn,
2433 wr.data(),
2434 work.data(),
2435 &lwork,
2436 &info);
2437 // sygv returns info=0 on success. Since we only queried the optimal size
2438 // for work, everything else would not be acceptable.
2439 Assert(info == 0, ExcInternalError());
2440 // Allocate working array according to suggestion (same strategy as was
2441 // noted in compute_svd).
2442 lwork = static_cast<types::blas_int>(std::abs(work[0]) + 1);
2443
2444 // resize workspace array
2445 work.resize(static_cast<size_type>(lwork));
2446
2447 // Finally compute the generalized eigenvalues.
2448 sygv(&itype,
2449 jobz,
2450 uplo,
2451 &nn,
2452 values_A,
2453 &nn,
2454 values_B,
2455 &nn,
2456 wr.data(),
2457 work.data(),
2458 &lwork,
2459 &info);
2460 // Negative return value implies a wrong argument. This should be internal.
2461
2462 Assert(info >= 0, ExcInternalError());
2463 if (info < 0)
2464 {
2465 AssertThrow(info == 0,
2466 ExcMessage("Lapack error in sygv: the " +
2467 std::to_string(-info) +
2468 "-th"
2469 " parameter had an illegal value."));
2470 }
2471 else if ((info > 0) && (info <= nn))
2472 {
2473 AssertThrow(
2474 info == 0,
2475 ExcMessage(
2476 "Lapack error in sygv: ssyev/dsyev failed to converge, and " +
2477 std::to_string(info) +
2478 " off-diagonal elements of an intermediate "
2479 " tridiagonal did not converge to zero."
2480 " (You may need to scale the abs_accuracy"
2481 " according to the norms of matrices A and B.)"));
2482 }
2483 else if ((info > nn) && (info <= 2 * nn))
2484 {
2485 AssertThrow(info == 0,
2486 ExcMessage(
2487 "Lapack error in sygv: the leading minor of order " +
2488 std::to_string(info - nn) +
2489 " of matrix B is not positive-definite."
2490 " The factorization of B could not be completed and"
2491 " no eigenvalues or eigenvectors were computed."));
2492 }
2493 else
2494 {
2495 AssertThrow(info == 0,
2496 ExcMessage("Lapack error in sygv: unknown error."));
2497 }
2498
2499 for (size_type i = 0; i < eigenvectors.size(); ++i)
2500 {
2501 size_type col_begin(i * nn);
2502 eigenvectors[i].reinit(nn, true);
2503 for (size_type j = 0; j < static_cast<size_type>(nn); ++j)
2504 {
2505 eigenvectors[i](j) = values_A[col_begin + j];
2506 }
2507 }
2508 state = LAPACKSupport::State(LAPACKSupport::eigenvalues | unusable);
2509}
2510
2511
2512
2513template <typename number>
2514void
2515LAPACKFullMatrix<number>::print_formatted(std::ostream &out,
2516 const unsigned int precision,
2517 const bool scientific,
2518 const unsigned int width_,
2519 const char *zero_string,
2520 const double denominator,
2521 const double threshold,
2522 const char *separator) const
2523{
2524 unsigned int width = width_;
2525
2526 Assert((!this->empty()) || (this->n() + this->m() == 0), ExcInternalError());
2527 Assert(state == LAPACKSupport::matrix ||
2528 state == LAPACKSupport::inverse_matrix ||
2529 state == LAPACKSupport::cholesky,
2530 ExcState(state));
2531
2532 // set output format, but store old
2533 // state
2534 std::ios::fmtflags old_flags = out.flags();
2535 std::streamsize old_precision = out.precision(precision);
2536
2537 if (scientific)
2538 {
2539 out.setf(std::ios::scientific, std::ios::floatfield);
2540 if (width == 0u)
2541 width = precision + 7;
2542 }
2543 else
2544 {
2545 out.setf(std::ios::fixed, std::ios::floatfield);
2546 if (width == 0u)
2547 width = precision + 2;
2548 }
2549
2550 for (size_type i = 0; i < this->m(); ++i)
2551 {
2552 // Cholesky is stored in lower triangular, so just output this part:
2553 const size_type nc = state == LAPACKSupport::cholesky ? i + 1 : this->n();
2554 for (size_type j = 0; j < nc; ++j)
2555 // we might have complex numbers, so use abs also to check for nan
2556 // since there is no isnan on complex numbers
2557 if (numbers::is_nan(std::abs((*this)(i, j))))
2558 out << std::setw(width) << (*this)(i, j) << separator;
2559 else if (std::abs(this->el(i, j)) > threshold)
2560 out << std::setw(width) << this->el(i, j) * denominator << separator;
2561 else
2562 out << std::setw(width) << zero_string << separator;
2563 out << std::endl;
2564 }
2565
2566 AssertThrow(out.fail() == false, ExcIO());
2567 // reset output format
2568 out.flags(old_flags);
2569 out.precision(old_precision);
2570}
2571
2572
2573
2574template <typename number>
2575LAPACKSupport::State
2576LAPACKFullMatrix<number>::get_state() const
2577{
2578 return this->state;
2579}
2580
2581
2582//----------------------------------------------------------------------//
2583
2584template <typename number>
2585void
2586PreconditionLU<number>::initialize(const LAPACKFullMatrix<number> &M)
2587{
2588 matrix = &M;
2589 mem = nullptr;
2590}
2591
2592
2593template <typename number>
2594void
2601
2602
2603template <typename number>
2604void
2605PreconditionLU<number>::vmult(Vector<number> &dst,
2606 const Vector<number> &src) const
2607{
2608 dst = src;
2609 matrix->solve(dst, false);
2610}
2611
2612
2613template <typename number>
2614void
2615PreconditionLU<number>::Tvmult(Vector<number> &dst,
2616 const Vector<number> &src) const
2617{
2618 dst = src;
2619 matrix->solve(dst, true);
2620}
2621
2622
2623template <typename number>
2624void
2625PreconditionLU<number>::vmult(BlockVector<number> &dst,
2626 const BlockVector<number> &src) const
2627{
2628 Assert(mem != nullptr, ExcNotInitialized());
2629 Vector<number> *aux = mem->alloc();
2630 *aux = src;
2631 matrix->solve(*aux, false);
2632 dst = *aux;
2633}
2634
2635
2636template <typename number>
2637void
2638PreconditionLU<number>::Tvmult(BlockVector<number> &dst,
2639 const BlockVector<number> &src) const
2640{
2641 Assert(mem != nullptr, ExcNotInitialized());
2642 Vector<number> *aux = mem->alloc();
2643 *aux = src;
2644 matrix->solve(*aux, true);
2645 dst = *aux;
2646}
2647
2648
2649
2650#include "lac/lapack_full_matrix.inst"
2651
2652
2653DEAL_II_NAMESPACE_CLOSE
blas_extension_templates.h
omatcopy
void omatcopy(char, char, ::types::blas_int, ::types::blas_int, const number1, const number2 *, ::types::blas_int, number3 *, ::types::blas_int)
Definition blas_extension_templates.h:41
block_vector.h
AlignedVector::data
pointer data()
EnableObserverPointer::operator=
EnableObserverPointer & operator=(const EnableObserverPointer &)
Definition enable_observer_pointer.h:309
FullMatrix::n
size_type n() const
FullMatrix::m
size_type m() const
LAPACKFullMatrix::operator*=
LAPACKFullMatrix< number > & operator*=(const number factor)
Definition lapack_full_matrix.cc:430
LAPACKFullMatrix::reciprocal_condition_number
number reciprocal_condition_number() const
Definition lapack_full_matrix.cc:1555
LAPACKFullMatrix::Tmmult
void Tmmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
Definition lapack_full_matrix.cc:1123
LAPACKFullMatrix::scale_rows
void scale_rows(const Vector< number > &V)
Definition lapack_full_matrix.cc:1106
LAPACKFullMatrix::get_right_eigenvectors
FullMatrix< std::complex< typename numbers::NumberTraits< number >::real_type > > get_right_eigenvectors() const
Definition lapack_full_matrix.cc:2059
LAPACKFullMatrix::add
void add(const number a, const LAPACKFullMatrix< number > &B)
Definition lapack_full_matrix.cc:490
LAPACKFullMatrix::Tvmult
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
Definition lapack_full_matrix.h:1156
LAPACKFullMatrix::transpose
void transpose(LAPACKFullMatrix< number > &B) const
Definition lapack_full_matrix.cc:1085
LAPACKFullMatrix::mTmult
void mTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
Definition lapack_full_matrix.cc:1217
LAPACKFullMatrix::get_state
LAPACKSupport::State get_state() const
Definition lapack_full_matrix.cc:2576
LAPACKFullMatrix::compute_eigenvalues_symmetric
void compute_eigenvalues_symmetric(const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, FullMatrix< number > &eigenvectors)
Definition lapack_full_matrix.cc:2099
LAPACKFullMatrix::vmult
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
Definition lapack_full_matrix.h:1129
LAPACKFullMatrix::reinit
void reinit(const size_type size)
Definition lapack_full_matrix.cc:284
LAPACKFullMatrix::operator=
LAPACKFullMatrix< number > & operator=(const LAPACKFullMatrix< number > &)
Definition lapack_full_matrix.cc:272
LAPACKFullMatrix::get_left_eigenvectors
FullMatrix< std::complex< typename numbers::NumberTraits< number >::real_type > > get_left_eigenvectors() const
Definition lapack_full_matrix.cc:2079
LAPACKFullMatrix::grow_or_shrink
void grow_or_shrink(const size_type size)
Definition lapack_full_matrix.cc:294
LAPACKFullMatrix::apply_givens_rotation
void apply_givens_rotation(const std::array< number, 3 > &csr, const size_type i, const size_type k, const bool left=true)
Definition lapack_full_matrix.cc:308
LAPACKFullMatrix::set_property
void set_property(const LAPACKSupport::Property property)
Definition lapack_full_matrix.cc:1406
LAPACKFullMatrix::norm
number norm(const char type) const
Definition lapack_full_matrix.cc:1445
LAPACKFullMatrix::solve
void solve(Vector< number > &v, const bool transposed=false) const
Definition lapack_full_matrix.cc:1757
LAPACKFullMatrix::compute_eigenvalues
void compute_eigenvalues(const bool right_eigenvectors=false, const bool left_eigenvectors=false)
Definition lapack_full_matrix.cc:1900
LAPACKFullMatrix::state
LAPACKSupport::State state
Definition lapack_full_matrix.h:935
LAPACKFullMatrix::size_type
std::make_unsigned_t< types::blas_int > size_type
Definition lapack_full_matrix.h:62
LAPACKFullMatrix::frobenius_norm
number frobenius_norm() const
Definition lapack_full_matrix.cc:1435
LAPACKFullMatrix::LAPACKFullMatrix
LAPACKFullMatrix(const size_type size=0)
Definition lapack_full_matrix.cc:244
LAPACKFullMatrix::property
LAPACKSupport::Property property
Definition lapack_full_matrix.h:941
LAPACKFullMatrix::m
size_type m() const
Definition lapack_full_matrix.h:1051
LAPACKFullMatrix::compute_inverse_svd
void compute_inverse_svd(const double threshold=0.)
Definition lapack_full_matrix.cc:1670
LAPACKFullMatrix::compute_generalized_eigenvalues_symmetric
void compute_generalized_eigenvalues_symmetric(LAPACKFullMatrix< number > &B, const number lower_bound, const number upper_bound, const number abs_accuracy, Vector< number > &eigenvalues, std::vector< Vector< number > > &eigenvectors, const types::blas_int itype=1)
Definition lapack_full_matrix.cc:2231
LAPACKFullMatrix::vmult_add
void vmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
Definition lapack_full_matrix.h:1143
LAPACKFullMatrix::n
size_type n() const
Definition lapack_full_matrix.h:1060
LAPACKFullMatrix::linfty_norm
number linfty_norm() const
Definition lapack_full_matrix.cc:1425
LAPACKFullMatrix::TmTmult
void TmTmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
Definition lapack_full_matrix.cc:1311
LAPACKFullMatrix::compute_inverse_svd_with_kernel
void compute_inverse_svd_with_kernel(const unsigned int kernel_size)
Definition lapack_full_matrix.cc:1693
LAPACKFullMatrix::Tvmult_add
void Tvmult_add(Vector< number2 > &w, const Vector< number2 > &v) const
Definition lapack_full_matrix.h:1170
LAPACKFullMatrix::rank1_update
void rank1_update(const number a, const Vector< number > &v)
Definition lapack_full_matrix.cc:617
LAPACKFullMatrix::remove_row_and_column
void remove_row_and_column(const size_type row, const size_type col)
Definition lapack_full_matrix.cc:342
LAPACKFullMatrix::operator/=
LAPACKFullMatrix< number > & operator/=(const number factor)
Definition lapack_full_matrix.cc:459
LAPACKFullMatrix::determinant
number determinant() const
Definition lapack_full_matrix.cc:1868
LAPACKFullMatrix::mmult
void mmult(LAPACKFullMatrix< number > &C, const LAPACKFullMatrix< number > &B, const bool adding=false) const
Definition lapack_full_matrix.cc:952
LAPACKFullMatrix::print_formatted
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1., const double threshold=0., const char *separator=" ") const
Definition lapack_full_matrix.cc:2515
PreconditionLU::vmult
void vmult(Vector< number > &, const Vector< number > &) const
Definition lapack_full_matrix.cc:2605
PreconditionLU::initialize
void initialize(const LAPACKFullMatrix< number > &)
Definition lapack_full_matrix.cc:2586
PreconditionLU::Tvmult
void Tvmult(Vector< number > &, const Vector< number > &) const
Definition lapack_full_matrix.cc:2615
SparseMatrix::el
number el(const size_type i, const size_type j) const
SparseMatrix::n
size_type n() const
SparseMatrix::m
size_type m() const
TableBase::values
AlignedVector< T > values
Definition table.h:793
TransposeTable::reinit
void reinit(const size_type size1, const size_type size2, const bool omit_default_initialization=false)
Vector::data
pointer data()
Vector::size
virtual size_type size() const override
Vector::begin
iterator begin()
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:38
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:39
full_matrix.h
StandardExceptions::ExcIO
static ::ExceptionBase & ExcIO()
StandardExceptions::ExcZero
static ::ExceptionBase & ExcZero()
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
LAPACKSupport::ExcErrorCode
static ::ExceptionBase & ExcErrorCode(std::string arg1, types::blas_int arg2)
StandardExceptions::ExcScalarAssignmentOnlyForZeroValue
static ::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
LAPACKSupport::ExcProperty
static ::ExceptionBase & ExcProperty(Property arg1)
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543
LACExceptions::ExcSingular
static ::ExceptionBase & ExcSingular()
AssertIsFinite
#define AssertIsFinite(number)
Definition exception_macros.h:932
AssertDimension
#define AssertDimension(dim1, dim2)
Definition exception_macros.h:829
AssertIndexRange
#define AssertIndexRange(index, range)
Definition exception_macros.h:901
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcDimensionMismatch
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
StandardExceptions::ExcNotInitialized
static ::ExceptionBase & ExcNotInitialized()
LACExceptions::ExcNotQuadratic
static ::ExceptionBase & ExcNotQuadratic()
StandardExceptions::ExcInvalidState
static ::ExceptionBase & ExcInvalidState()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
LAPACKSupport::ExcState
static ::ExceptionBase & ExcState(State arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
lapack_full_matrix.h
lapack_support.h
lapack_templates.h
getrs
void getrs(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:3308
syrk
void syrk(const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, number4 *, const ::types::blas_int *)
Definition lapack_templates.h:5380
geev
void geev(const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, number3 *, number4 *, const ::types::blas_int *, number5 *, const ::types::blas_int *, number6 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:1521
gemm
void gemm(const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, const ::types::blas_int *, const number4 *, number5 *, const ::types::blas_int *)
Definition lapack_templates.h:2251
pocon
void pocon(const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, const number2 *, number3 *, number4 *, ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:4033
lascl
void lascl(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const ::types::blas_int *, number3 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:3707
syr
void syr(const char *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, number3 *, const ::types::blas_int *)
Definition lapack_templates.h:5312
trtrs
void trtrs(const char *, const char *, const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:5796
syevx
void syevx(const char *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, const number2 *, const number3 *, const ::types::blas_int *, const ::types::blas_int *, const number4 *, ::types::blas_int *, number5 *, number6 *, const ::types::blas_int *, number7 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:4849
trcon
void trcon(const char *, const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, number3 *, ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:5526
gesdd
void gesdd(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 *, ::types::blas_int *)
Definition lapack_templates.h:2711
axpy
void axpy(const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, number3 *, const ::types::blas_int *)
Definition lapack_templates.h:1410
trmv
void trmv(const char *, const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *)
Definition lapack_templates.h:5672
sygvx
void sygvx(const ::types::blas_int *, const char *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, const number3 *, const number4 *, const ::types::blas_int *, const ::types::blas_int *, const number5 *, ::types::blas_int *, number6 *, number7 *, const ::types::blas_int *, number8 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:5121
gemv
void gemv(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const number2 *, const ::types::blas_int *, const number3 *, const ::types::blas_int *, const number4 *, number5 *, const ::types::blas_int *)
Definition lapack_templates.h:2428
lange
number1 lange(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *)
Definition lapack_templates.h:3485
potrs
void potrs(const char *, const ::types::blas_int *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:4364
sygv
void sygv(const ::types::blas_int *, const char *, const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, number2 *, const ::types::blas_int *, number3 *, number4 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:5022
potri
void potri(const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:4267
lansy
number1 lansy(const char *, const char *, const ::types::blas_int *, const number1 *, const ::types::blas_int *, number2 *)
Definition lapack_templates.h:3596
potrf
void potrf(const char *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:4170
getrf
void getrf(const ::types::blas_int *, const ::types::blas_int *, number1 *, const ::types::blas_int *, ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:3087
getri
void getri(const ::types::blas_int *, number1 *, const ::types::blas_int *, const ::types::blas_int *, number2 *, const ::types::blas_int *, ::types::blas_int *)
Definition lapack_templates.h:3193
LAPACKSupport::cholesky
@ cholesky
Contents is a Cholesky decomposition.
Definition lapack_support.h:62
LAPACKSupport::lu
@ lu
Contents is an LU decomposition.
Definition lapack_support.h:60
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition lapack_support.h:56
LAPACKSupport::unusable
@ unusable
Contents is something useless.
Definition lapack_support.h:70
LAPACKSupport::inverse_matrix
@ inverse_matrix
Contents is the inverse of a matrix.
Definition lapack_support.h:58
LAPACKSupport::svd
@ svd
Matrix contains singular value decomposition,.
Definition lapack_support.h:66
LAPACKSupport::inverse_svd
@ inverse_svd
Matrix is the inverse of a singular value decomposition.
Definition lapack_support.h:68
LAPACKSupport::eigenvalues
@ eigenvalues
Eigenvalue vector is filled.
Definition lapack_support.h:64
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition lapack_support.h:111
LAPACKSupport::upper_triangular
@ upper_triangular
Matrix is upper triangular.
Definition lapack_support.h:113
LAPACKSupport::lower_triangular
@ lower_triangular
Matrix is lower triangular.
Definition lapack_support.h:115
LAPACKSupport::general
@ general
No special properties.
Definition lapack_support.h:109
LAPACKSupport::L
constexpr char L
Definition lapack_support.h:171
LAPACKSupport::N
constexpr char N
Definition lapack_support.h:155
LAPACKSupport::U
constexpr char U
Definition lapack_support.h:167
LAPACKSupport::T
constexpr char T
Definition lapack_support.h:163
LAPACKSupport::V
constexpr char V
Definition lapack_support.h:175
LAPACKSupport::A
constexpr char A
Definition lapack_support.h:151
LAPACKSupport::one
constexpr types::blas_int one
Definition lapack_support.h:183
Utilities::LinearAlgebra::givens_rotation
std::array< NumberType, 3 > givens_rotation(const NumberType &x, const NumberType &y)
Utilities::LinearAlgebra::hyperbolic_rotation
std::array< NumberType, 3 > hyperbolic_rotation(const NumberType &x, const NumberType &y)
internal::LAPACKFullMatrixImplementation::gesdd_helper
void gesdd_helper(const char job, const types::blas_int n_rows, const types::blas_int n_cols, AlignedVector< T > &matrix, std::vector< T > &singular_values, AlignedVector< T > &left_vectors, AlignedVector< T > &right_vectors, std::vector< T > &real_work, std::vector< T > &, std::vector< types::blas_int > &integer_work, const types::blas_int work_flag, types::blas_int &info)
Definition lapack_full_matrix.cc:152
internal::LAPACKFullMatrixImplementation::geev_helper
void geev_helper(const char vl, const char vr, AlignedVector< T > &matrix, const types::blas_int n_rows, std::vector< T > &real_part_eigenvalues, std::vector< T > &imag_part_eigenvalues, std::vector< T > &left_eigenvectors, std::vector< T > &right_eigenvectors, std::vector< T > &real_work, std::vector< T > &, const types::blas_int work_flag, types::blas_int &info)
Definition lapack_full_matrix.cc:44
numbers::is_nan
bool is_nan(const double x)
Definition numbers.h:501
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596
std::max
::VectorizedArray< Number, width > max(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6579
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6505
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563
sparse_matrix.h
PointerComparison::equal
static bool equal(const T *p1, const T *p2)
Definition template_constraints.h:357
numbers::NumberTraits::conjugate
static constexpr const number & conjugate(const number &x)
Definition numbers.h:541
eigenvectors
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
Generated by doxygen 1.9.8