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