Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
full_matrix.templates.h
1 // ---------------------------------------------------------------------
2 // @f$Id: full_matrix.templates.h 30040 2013-07-18 17:06:48Z maier @f$
3 //
4 // Copyright (C) 1999 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__full_matrix_templates_h
18 #define __deal2__full_matrix_templates_h
19 
20 
21 //TODO: this file has a lot of operations between matrices and matrices or matrices and vectors of different precision. we should go through the file and in each case pick the more accurate data type for intermediate results. currently, the choice is pretty much random. this may also allow us some operations where one operand is complex and the other is not
22 
23 #include <deal.II/base/config.h>
24 #include <deal.II/base/template_constraints.h>
25 #include <deal.II/lac/vector.h>
26 #include <deal.II/lac/full_matrix.h>
27 #include <deal.II/lac/lapack_full_matrix.h>
28 #include <deal.II/lac/lapack_templates.h>
29 
30 #include <vector>
31 #include <cmath>
32 #include <cstdlib>
33 #include <cstdio>
34 #include <algorithm>
35 
37 
38 
39 template <typename number>
41  :
42  Table<2,number> (n,n)
43 {}
44 
45 
46 template <typename number>
48  const size_type n)
49  :
50  Table<2,number> (m, n)
51 {}
52 
53 
54 template <typename number>
56  const size_type n,
57  const number *entries)
58  :
59  Table<2,number> (m, n)
60 {
61  this->fill (entries);
62 }
63 
64 
65 template <typename number>
67  :
68  Table<2,number> (m)
69 {}
70 
71 
72 
73 template <typename number>
75  :
76  Table<2,number> (id.m(), id.n())
77 {
78  for (size_type i=0; i<id.m(); ++i)
79  (*this)(i,i) = 1;
80 }
81 
82 
83 template <typename number>
86 {
88  return *this;
89 }
90 
91 
92 template <typename number>
93 template <typename number2>
96 {
98  return *this;
99 }
100 
101 
102 
103 template <typename number>
106 {
107  this->reinit (id.m(), id.n());
108  for (size_type i=0; i<id.m(); ++i)
109  (*this)(i,i) = 1.;
110 
111  return *this;
112 }
113 
114 
115 
116 template <typename number>
117 template <typename number2>
120 {
121  Assert (this->m() == M.n_rows(), ExcDimensionMismatch(this->m(), M.n_rows()));
122  Assert (this->n() == M.n_cols(), ExcDimensionMismatch(this->n(), M.n_rows()));
123  for (size_type i=0; i<this->m(); ++i)
124  for (size_type j=0; j<this->n(); ++j)
125  (*this)(i,j) = M(i,j);
126 
127  return *this;
128 }
129 
130 
131 
132 template <typename number>
133 bool
135 {
136  Assert (!this->empty(), ExcEmptyMatrix());
137 
138  const number *p = &this->values[0];
139  const number *const e = &this->values[0] + this->n_elements();
140  while (p!=e)
141  if (*p++ != number(0.0))
142  return false;
143 
144  return true;
145 }
146 
147 
148 
149 template <typename number>
152 {
153 
155 
156  number *p = &(*this)(0,0);
157  const number *e = &(*this)(0,0) + n()*m();
158  while (p != e)
159  *p++ *= factor;
160 
161  return *this;
162 }
163 
164 
165 
166 template <typename number>
169 {
170 
172 
173  number *p = &(*this)(0,0);
174  const number *e = &(*this)(0,0) + n()*m();
175 
176  const number factor_inv = number(1.)/factor;
177 
179 
180  while (p != e)
181  *p++ *= factor_inv;
182 
183  return *this;
184 }
185 
186 
187 
188 template <typename number>
189 template <typename number2>
190 void
192  const Vector<number2> &src,
193  const bool adding) const
194 {
195  Assert (!this->empty(), ExcEmptyMatrix());
196 
197  Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
198  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
199 
200  Assert (&src != &dst, ExcSourceEqualsDestination());
201 
202  const number *e = &this->values[0];
203  // get access to the data in order to
204  // avoid copying it when using the ()
205  // operator
206  const number2 *src_ptr = &(*const_cast<Vector<number2>*>(&src))(0);
207  const size_type size_m = m(), size_n = n();
208  for (size_type i=0; i<size_m; ++i)
209  {
210  number2 s = adding ? dst(i) : 0.;
211  for (size_type j=0; j<size_n; ++j)
212  s += src_ptr[j] * number2(*(e++));
213  dst(i) = s;
214  }
215 }
216 
217 
218 
219 template <typename number>
220 template <typename number2>
222  const Vector<number2> &src,
223  const bool adding) const
224 {
225  Assert (!this->empty(), ExcEmptyMatrix());
226 
227  Assert(dst.size() == n(), ExcDimensionMismatch(dst.size(), n()));
228  Assert(src.size() == m(), ExcDimensionMismatch(src.size(), m()));
229 
230  Assert (&src != &dst, ExcSourceEqualsDestination());
231 
232  const number *e = &this->values[0];
233  number2 *dst_ptr = &dst(0);
234  const size_type size_m = m(), size_n = n();
235 
236  // zero out data if we are not adding
237  if (!adding)
238  for (size_type j=0; j<size_n; ++j)
239  dst_ptr[j] = 0.;
240 
241  // write the loop in a way that we can
242  // access the data contiguously
243  for (size_type i=0; i<size_m; ++i)
244  {
245  const number2 d = src(i);
246  for (size_type j=0; j<size_n; ++j)
247  dst_ptr[j] += d * number2(*(e++));
248  };
249 }
250 
251 
252 template <typename number>
253 template <typename number2, typename number3>
255  const Vector<number2> &src,
256  const Vector<number3> &right) const
257 {
258  Assert (!this->empty(), ExcEmptyMatrix());
259 
260  Assert(dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
261  Assert(src.size() == n(), ExcDimensionMismatch(src.size(), n()));
262  Assert(right.size() == m(), ExcDimensionMismatch(right.size(), m()));
263 
264  Assert (&src != &dst, ExcSourceEqualsDestination());
265 
266  number res = 0.;
267  const size_type size_m = m(),
268  size_n = n();
269  for (size_type i=0; i<size_n; ++i)
270  {
271  number s = number(right(i));
272  for (size_type j=0; j<size_m; ++j)
273  s -= number(src(j)) * (*this)(i,j);
274  dst(i) = s;
275  res += s*s;
276  }
277  return std::sqrt(res);
278 }
279 
280 
281 
282 template <typename number>
283 template <typename number2>
285  const Vector<number2> &src) const
286 {
287  Assert (!this->empty(), ExcEmptyMatrix());
288 
289  Assert (dst.size() == m(), ExcDimensionMismatch(dst.size(), m()));
290  Assert (src.size() == n(), ExcDimensionMismatch(src.size(), n()));
291 
292  size_type i,j;
293  size_type nu = ( (m()<n()) ? m() : n());
294  for (i=0; i<nu; ++i)
295  {
296  number s = number(src(i));
297  for (j=0; j<i; ++j)
298  s -= number(dst(j)) * (*this)(i,j);
299  dst(i) = s/(*this)(i,i);
301  }
302 }
303 
304 
305 
306 template <typename number>
307 template <typename number2>
309  const Vector<number2> &src) const
310 {
311  Assert (!this->empty(), ExcEmptyMatrix());
312 
313  size_type j;
314  size_type nu = (m()<n() ? m() : n());
315  for (int i=nu-1; i>=0; --i)
316  {
317  number2 s = src(i);
318  for (j=i+1; j<nu; ++j)
319  s -= dst(j) * number2((*this)(i,j));
320  dst(i) = s/number2((*this)(i,i));
322  }
323 }
324 
325 
326 
327 template <typename number>
328 template <typename number2>
330  const size_type dst_offset_i,
331  const size_type dst_offset_j,
332  const size_type src_offset_i,
333  const size_type src_offset_j)
334 {
335  Assert (dst_offset_i < m(),
336  ExcIndexRange (dst_offset_i, 0, m()));
337  Assert (dst_offset_j < n(),
338  ExcIndexRange (dst_offset_j, 0, n()));
339  Assert (src_offset_i < src.m(),
340  ExcIndexRange (src_offset_i, 0, src.m()));
341  Assert (src_offset_j < src.n(),
342  ExcIndexRange (src_offset_j, 0, src.n()));
343 
344  // Compute maximal size of copied block
345  const size_type rows = std::min (m() - dst_offset_i,
346  src.m() - src_offset_i);
347  const size_type cols = std::min (n() - dst_offset_j,
348  src.n() - src_offset_j);
349 
350  for (size_type i=0; i<rows ; ++i)
351  for (size_type j=0; j<cols ; ++j)
352  (*this)(dst_offset_i+i,dst_offset_j+j)
353  = src(src_offset_i+i,src_offset_j+j);
354 }
355 
356 
357 template <typename number>
358 template <typename number2>
360  const std::vector<size_type> &p_rows,
361  const std::vector<size_type> &p_cols)
362 {
363  Assert (p_rows.size() == this->n_rows(),
364  ExcDimensionMismatch (p_rows.size(), this->n_rows()));
365  Assert (p_cols.size() == this->n_cols(),
366  ExcDimensionMismatch (p_cols.size(), this->n_cols()));
367 
368  for (size_type i=0; i<this->n_rows(); ++i)
369  for (size_type j=0; j<this->n_cols(); ++j)
370  (*this)(i,j) = src(p_rows[i], p_cols[j]);
371 }
372 
373 
374 
375 /* template <typename number> */
376 /* template <typename number2> */
377 /* void FullMatrix<number>::fill (const number2* entries) */
378 /* { */
379 /* if (n_cols()*n_rows() != 0) */
380 /* std::copy (entries, entries+n_rows()*n_cols(), &this->values[0]); */
381 /* } */
382 
383 
384 
385 template <typename number>
387  const number s,
388  const size_type j)
389 {
390  Assert (!this->empty(), ExcEmptyMatrix());
391 
392  for (size_type k=0; k<m(); ++k)
393  (*this)(i,k) += s*(*this)(j,k);
394 }
395 
396 
397 template <typename number>
399  const number s,
400  const size_type j,
401  const number t,
402  const size_type k)
403 {
404  Assert (!this->empty(), ExcEmptyMatrix());
405 
406  const size_type size_m = m();
407  for (size_type l=0; l<size_m; ++l)
408  (*this)(i,l) += s*(*this)(j,l) + t*(*this)(k,l);
409 }
410 
411 
412 template <typename number>
413 void FullMatrix<number>::add_col (const size_type i, const number s,
414  const size_type j)
415 {
416  Assert (!this->empty(), ExcEmptyMatrix());
417 
418  for (size_type k=0; k<n(); ++k)
419  (*this)(k,i) += s*(*this)(k,j);
420 }
421 
422 
423 template <typename number>
424 void FullMatrix<number>::add_col (const size_type i, const number s,
425  const size_type j, const number t,
426  const size_type k)
427 {
428  Assert (!this->empty(), ExcEmptyMatrix());
429 
430  for (size_t l=0; l<n(); ++l)
431  (*this)(l,i) += s*(*this)(l,j) + t*(*this)(l,k);
432 }
433 
434 
435 
436 template <typename number>
438  const size_type j)
439 {
440  Assert (!this->empty(), ExcEmptyMatrix());
441 
442  for (size_type k=0; k<n(); ++k)
443  std::swap ((*this)(i,k),
444  (*this)(j,k));
445 }
446 
447 
448 template <typename number>
450  const size_type j)
451 {
452  Assert (!this->empty(), ExcEmptyMatrix());
453 
454  for (size_type k=0; k<m(); ++k)
455  std::swap ((*this)(k,i),
456  (*this)(k,j));
457 }
458 
459 
460 template <typename number>
461 void FullMatrix<number>::diagadd (const number src)
462 {
463  Assert (!this->empty(), ExcEmptyMatrix());
464  Assert (m() == n(), ExcDimensionMismatch(m(),n()));
465 
466  for (size_type i=0; i<n(); ++i)
467  (*this)(i,i) += src;
468 }
469 
470 
471 template <typename number>
472 template <typename number2>
473 void FullMatrix<number>::equ (const number a,
474  const FullMatrix<number2> &A)
475 {
476  Assert (!this->empty(), ExcEmptyMatrix());
477 
478  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
479  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
480 
481  for (size_type i=0; i<m(); ++i)
482  for (size_type j=0; j<n(); ++j)
483  (*this)(i,j) = a * number(A(i,j));
484 }
485 
486 
487 template <typename number>
488 template <typename number2>
489 void
490 FullMatrix<number>::equ (const number a,
491  const FullMatrix<number2> &A,
492  const number b,
493  const FullMatrix<number2> &B)
494 {
495  Assert (!this->empty(), ExcEmptyMatrix());
496 
497  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
498  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
499  Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
500  Assert (n() == B.n(), ExcDimensionMismatch(n(), B.n()));
501 
502  for (size_type i=0; i<m(); ++i)
503  for (size_type j=0; j<n(); ++j)
504  (*this)(i,j) = a * number(A(i,j)) + b * number(B(i,j));
505 }
506 
507 
508 template <typename number>
509 template <typename number2>
510 void
511 FullMatrix<number>::equ (const number a,
512  const FullMatrix<number2> &A,
513  const number b,
514  const FullMatrix<number2> &B,
515  const number c,
516  const FullMatrix<number2> &C)
517 {
518  Assert (!this->empty(), ExcEmptyMatrix());
519 
520  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
521  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
522  Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
523  Assert (n() == B.n(), ExcDimensionMismatch(n(), B.n()));
524  Assert (m() == C.m(), ExcDimensionMismatch(m(), C.m()));
525  Assert (n() == C.n(), ExcDimensionMismatch(n(), C.n()));
526 
527  for (size_type i=0; i<m(); ++i)
528  for (size_type j=0; j<n(); ++j)
529  (*this)(i,j) = a * number(A(i,j)) +
530  b * number(B(i,j)) +
531  c * number(C(i,j));
532 }
533 
534 
535 
536 template <typename number>
537 template <typename number2>
539  const FullMatrix<number2> &src,
540  const bool adding) const
541 {
542  Assert (!this->empty(), ExcEmptyMatrix());
543  Assert (n() == src.m(), ExcDimensionMismatch(n(), src.m()));
544  Assert (dst.n() == src.n(), ExcDimensionMismatch(dst.n(), src.n()));
545  Assert (dst.m() == m(), ExcDimensionMismatch(m(), dst.m()));
546 
547  // see if we can use BLAS algorithms for this and if the type for 'number'
548  // works for us (it is usually not efficient to use BLAS for very small
549  // matrices):
550 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
552  ||
554  &&
556  if (this->n()*this->m()*src.n() > 300)
557  {
558  // In case we have the BLAS function gemm detected at configure, we
559  // use that algorithm for matrix-matrix multiplication since it
560  // provides better performance than the deal.II native function (it
561  // uses cache and register blocking in order to access local data).
562  //
563  // Note that BLAS/LAPACK stores matrix elements column-wise (i.e., all
564  // values in one column, then all in the next, etc.), whereas the
565  // FullMatrix stores them row-wise. We ignore that difference, and
566  // give our row-wise data to BLAS, let BLAS build the product of
567  // transpose matrices, and read the result as if it were row-wise
568  // again. In other words, we calculate (B^T A^T)^T, which is AB.
569 
570  const int m = src.n();
571  const int n = this->m();
572  const int k = this->n();
573  const char *notrans = "n";
574 
575  const number alpha = 1.;
576  const number beta = (adding == true) ? 1. : 0.;
577 
578  // Use the BLAS function gemm for
579  // calculating the matrix-matrix
580  // product.
581  gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
582  &this->values[0], &k, &beta, &dst(0,0), &m);
583 
584  return;
585  }
586 
587 #endif
588 
589  const size_type m = this->m(), n = src.n(), l = this->n();
590 
591  // arrange the loops in a way that we keep write operations low, (writing is
592  // usually more costly than reading), even though we need to access the data
593  // in src not in a contiguous way.
594  for (size_type i=0; i<m; i++)
595  for (size_type j=0; j<n; j++)
596  {
597  number2 add_value = adding ? dst(i,j) : 0.;
598  for (size_type k=0; k<l; k++)
599  add_value += (number2)(*this)(i,k) * (number2)(src(k,j));
600  dst(i,j) = add_value;
601  }
602 }
603 
604 
605 
606 template <typename number>
607 template <typename number2>
609  const FullMatrix<number2> &src,
610  const bool adding) const
611 {
612  Assert (!this->empty(), ExcEmptyMatrix());
613  Assert (m() == src.m(), ExcDimensionMismatch(m(), src.m()));
614  Assert (n() == dst.m(), ExcDimensionMismatch(n(), dst.m()));
615  Assert (src.n() == dst.n(), ExcDimensionMismatch(src.n(), dst.n()));
616 
617 
618  // see if we can use BLAS algorithms for this and if the type for 'number'
619  // works for us (it is usually not efficient to use BLAS for very small
620  // matrices):
621 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
623  ||
625  &&
627  if (this->n()*this->m()*src.n() > 300)
628  {
629  // In case we have the BLAS function gemm detected at configure, we
630  // use that algorithm for matrix-matrix multiplication since it
631  // provides better performance than the deal.II native function (it
632  // uses cache and register blocking in order to access local data).
633  //
634  // Note that BLAS/LAPACK stores matrix elements column-wise (i.e., all
635  // values in one column, then all in the next, etc.), whereas the
636  // FullMatrix stores them row-wise. We ignore that difference, and
637  // give our row-wise data to BLAS, let BLAS build the product of
638  // transpose matrices, and read the result as if it were row-wise
639  // again. In other words, we calculate (B^T A)^T, which is A^T B.
640 
641  const int m = src.n();
642  const int n = this->n();
643  const int k = this->m();
644  const char *trans = "t";
645  const char *notrans = "n";
646 
647  const number alpha = 1.;
648  const number beta = (adding == true) ? 1. : 0.;
649 
650  // Use the BLAS function gemm for calculating the matrix-matrix
651  // product.
652  gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
653  &this->values[0], &n, &beta, &dst(0,0), &m);
654 
655  return;
656  }
657 
658 #endif
659 
660  const size_type m = n(), n = src.n(), l = this->m();
661 
662  // symmetric matrix if the two matrices are the same
663  if (PointerComparison::equal(this, &src))
664  for (size_type i=0; i<m; ++i)
665  for (size_type j=i; j<m; ++j)
666  {
667  number2 add_value = 0.;
668  for (size_type k=0; k<l; ++k)
669  add_value += (number2)(*this)(k,i) * (number2)(*this)(k,j);
670  if (adding)
671  {
672  dst(i,j) += add_value;
673  if (i<j)
674  dst(j,i) += add_value;
675  }
676  else
677  dst(i,j) = dst(j,i) = add_value;
678  }
679  // arrange the loops in a way that we keep write operations low, (writing is
680  // usually more costly than reading), even though we need to access the data
681  // in src not in a contiguous way. However, we should usually end up in the
682  // optimized gemm operation in case the matrix is big, so this shouldn't be
683  // too bad.
684  else
685  for (size_type i=0; i<m; i++)
686  for (size_type j=0; j<n; j++)
687  {
688  number2 add_value = adding ? dst(i,j) : 0.;
689  for (size_type k=0; k<l; k++)
690  add_value += (number2)(*this)(k,i) * (number2)(src(k,j));
691  dst(i,j) = add_value;
692  }
693 }
694 
695 
696 
697 template <typename number>
698 template <typename number2>
700  const FullMatrix<number2> &src,
701  const bool adding) const
702 {
703  Assert (!this->empty(), ExcEmptyMatrix());
704  Assert (n() == src.n(), ExcDimensionMismatch(n(), src.n()));
705  Assert (dst.n() == src.m(), ExcDimensionMismatch(dst.n(), src.m()));
706  Assert (dst.m() == m(), ExcDimensionMismatch(m(), dst.m()));
707 
708  // see if we can use BLAS algorithms for this and if the type for 'number'
709  // works for us (it is usually not efficient to use BLAS for very small
710  // matrices):
711 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
713  ||
715  &&
717  if (this->n()*this->m()*src.m() > 300)
718  {
719  // In case we have the BLAS function gemm detected at configure, we
720  // use that algorithm for matrix-matrix multiplication since it
721  // provides better performance than the deal.II native function (it
722  // uses cache and register blocking in order to access local data).
723  //
724  // Note that BLAS/LAPACK stores matrix elements column-wise (i.e., all
725  // values in one column, then all in the next, etc.), whereas the
726  // FullMatrix stores them row-wise. We ignore that difference, and
727  // give our row-wise data to BLAS, let BLAS build the product of
728  // transpose matrices, and read the result as if it were row-wise
729  // again. In other words, we calculate (B A^T)^T, which is AB^T.
730 
731  const int m = src.m();
732  const int n = this->m();
733  const int k = this->n();
734  const char *notrans = "n";
735  const char *trans = "t";
736 
737  const number alpha = 1.;
738  const number beta = (adding == true) ? 1. : 0.;
739 
740  // Use the BLAS function gemm for calculating the matrix-matrix
741  // product.
742  gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
743  &this->values[0], &k, &beta, &dst(0,0), &m);
744 
745  return;
746  }
747 
748 #endif
749 
750  const size_type m = this->m(), n = src.m(), l = this->n();
751 
752  // symmetric matrix if the two matrices are the same
753  if (PointerComparison::equal(this, &src))
754  for (size_type i=0; i<m; ++i)
755  for (size_type j=i; j<m; ++j)
756  {
757  number2 add_value = 0.;
758  for (size_type k=0; k<l; ++k)
759  add_value += (number2)(*this)(i,k) * (number2)(*this)(j,k);
760  if (adding)
761  {
762  dst(i,j) += add_value;
763  if (i<j)
764  dst(j,i) += add_value;
765  }
766  else
767  dst(i,j) = dst(j,i) = add_value;
768  }
769  else
770  // arrange the loops in a way that we keep write operations low, (writing is
771  // usually more costly than reading).
772  for (size_type i=0; i<m; i++)
773  for (size_type j=0; j<n; j++)
774  {
775  number2 add_value = adding ? dst(i,j) : 0.;
776  for (size_type k=0; k<l; k++)
777  add_value += (number2)(*this)(i,k) * (number2)(src(j,k));
778  dst(i,j) = add_value;
779  }
780 }
781 
782 
783 
784 template <typename number>
785 template <typename number2>
787  const FullMatrix<number2> &src,
788  const bool adding) const
789 {
790  Assert (!this->empty(), ExcEmptyMatrix());
791  Assert (m() == src.n(), ExcDimensionMismatch(m(), src.n()));
792  Assert (n() == dst.m(), ExcDimensionMismatch(n(), dst.m()));
793  Assert (src.m() == dst.n(), ExcDimensionMismatch(src.m(), dst.n()));
794 
795 
796  // see if we can use BLAS algorithms for this and if the type for 'number'
797  // works for us (it is usually not efficient to use BLAS for very small
798  // matrices):
799 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
801  ||
803  &&
805  if (this->n()*this->m()*src.m() > 300)
806  {
807  // In case we have the BLAS function gemm detected at configure, we
808  // use that algorithm for matrix-matrix multiplication since it
809  // provides better performance than the deal.II native function (it
810  // uses cache and register blocking in order to access local data).
811  //
812  // Note that BLAS/LAPACK stores matrix elements column-wise (i.e., all
813  // values in one column, then all in the next, etc.), whereas the
814  // FullMatrix stores them row-wise. We ignore that difference, and
815  // give our row-wise data to BLAS, let BLAS build the product of
816  // transpose matrices, and read the result as if it were row-wise
817  // again. In other words, we calculate (B A)^T, which is A^T B^T.
818 
819  const int m = src.n();
820  const int n = this->n();
821  const int k = this->m();
822  const char *trans = "t";
823 
824  const number alpha = 1.;
825  const number beta = (adding == true) ? 1. : 0.;
826 
827  // Use the BLAS function gemm for calculating the matrix-matrix
828  // product.
829  gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
830  &this->values[0], &n, &beta, &dst(0,0), &m);
831 
832  return;
833  }
834 
835 #endif
836 
837  const size_type m = n(), n = src.m(), l = this->m();
838 
839  // arrange the loops in a way that we keep write operations low, (writing is
840  // usually more costly than reading), even though we need to access the data
841  // in the calling matrix in a non-contiguous way, possibly leading to cache
842  // misses. However, we should usually end up in the optimized gemm operation
843  // in case the matrix is big, so this shouldn't be too bad.
844  for (size_type i=0; i<m; i++)
845  for (size_type j=0; j<n; j++)
846  {
847  number2 add_value = adding ? dst(i,j) : 0.;
848  for (size_type k=0; k<l; k++)
849  add_value += (number2)(*this)(k,i) * (number2)(src(j,k));
850  dst(i,j) = add_value;
851  }
852 }
853 
854 
855 template <typename number>
856 void
858  const FullMatrix<number> &A,
859  const FullMatrix<number> &B,
860  const FullMatrix<number> &D,
861  const bool transpose_B,
862  const bool transpose_D,
863  const number scaling)
864 {
865  if (transpose_B)
866  {
867  AssertDimension(B.m(), A.m());
868  AssertDimension(B.n(), m());
869  }
870  else
871  {
872  AssertDimension(B.n(), A.m());
873  AssertDimension(B.m(), m());
874  }
875  if (transpose_D)
876  {
877  AssertDimension(D.n(), A.n());
878  AssertDimension(D.m(), n());
879  }
880  else
881  {
882  AssertDimension(D.m(), A.n());
883  AssertDimension(D.n(), n());
884  }
885 
886  // For all entries of the product
887  // AD
888  for (size_type i=0; i<A.m(); ++i)
889  for (size_type j=0; j<n(); ++j)
890  {
891  // Compute the entry
892  number ADij = 0.;
893  if (transpose_D)
894  for (size_type k=0; k<A.n(); ++k)
895  ADij += A(i,k)*D(j,k);
896  else
897  for (size_type k=0; k<A.n(); ++k)
898  ADij += A(i,k)*D(k,j);
899  // And add it to this after
900  // multiplying with the right
901  // factor from B
902  if (transpose_B)
903  for (size_type k=0; k<m(); ++k)
904  this->operator()(k,j) += scaling * ADij * B(i,k);
905  else
906  for (size_type k=0; k<m(); ++k)
907  this->operator()(k,j) += scaling * ADij * B(k,i);
908  }
909 }
910 
911 
912 template <typename number>
913 template <typename number2>
914 number2
916 {
917  Assert (!this->empty(), ExcEmptyMatrix());
918 
919  Assert(m() == v.size(), ExcDimensionMismatch(m(),v.size()));
920  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
921 
922  number2 sum = 0.;
923  const size_type n_rows = m();
924  const number *val_ptr = &this->values[0];
925  const number2 *v_ptr;
926 
927  for (size_type row=0; row<n_rows; ++row)
928  {
929  number s = 0.;
930  const number *const val_end_of_row = val_ptr+n_rows;
931  v_ptr = v.begin();
932  while (val_ptr != val_end_of_row)
933  s += number(*val_ptr++) * number(*v_ptr++);
934 
935  sum += s * number(numbers::NumberTraits<number2>::conjugate(v(row)));
936  }
937 
938  return sum;
939 }
940 
941 
942 template <typename number>
943 template <typename number2>
944 number2
946  const Vector<number2> &v) const
947 {
948  Assert (!this->empty(), ExcEmptyMatrix());
949 
950  Assert(m() == u.size(), ExcDimensionMismatch(m(),v.size()));
951  Assert(n() == v.size(), ExcDimensionMismatch(n(),v.size()));
952 
953  number2 sum = 0.;
954  const size_type n_rows = m();
955  const size_type n_cols = n();
956  const number *val_ptr = &this->values[0];
957  const number2 *v_ptr;
958 
959  for (size_type row=0; row<n_rows; ++row)
960  {
961  number s = 0.;
962  const number *const val_end_of_row = val_ptr+n_cols;
963  v_ptr = v.begin();
964  while (val_ptr != val_end_of_row)
965  s += number(*val_ptr++) * number(*v_ptr++);
966 
967  sum += s * number(u(row));
968  }
969 
970  return sum;
971 }
972 
973 
974 
975 template <typename number>
976 void
978 {
979  Assert (m() == n(), ExcNotQuadratic());
980 
981  const size_type N = m();
982  for (size_type i=0; i<N; ++i)
983  for (size_type j=i+1; j<N; ++j)
984  {
985  const number t = ((*this)(i,j) + (*this)(j,i)) / number(2.);
986  (*this)(i,j) = (*this)(j,i) = t;
987  };
988 }
989 
990 
991 template <typename number>
994 {
995  Assert (!this->empty(), ExcEmptyMatrix());
996 
997  real_type sum=0, max=0;
998  const size_type n_rows = m(), n_cols = n();
999 
1000  for (size_type col=0; col<n_cols; ++col)
1001  {
1002  sum=0;
1003  for (size_type row=0; row<n_rows; ++row)
1004  sum += std::abs((*this)(row,col));
1005  if (sum > max)
1006  max = sum;
1007  }
1008  return max;
1009 }
1010 
1011 
1012 
1013 template <typename number>
1016 {
1017  Assert (!this->empty(), ExcEmptyMatrix());
1018 
1019  real_type sum=0, max=0;
1020  const size_type n_rows = m(), n_cols = n();
1021 
1022  for (size_type row=0; row<n_rows; ++row)
1023  {
1024  sum=0;
1025  for (size_type col=0; col<n_cols; ++col)
1026  sum += std::abs((*this)(row,col));
1027  if (sum > max)
1028  max = sum;
1029  }
1030  return max;
1031 }
1032 
1033 
1034 
1035 template <typename number>
1036 template <typename number2>
1037 void
1038 FullMatrix<number>::add (const number a,
1039  const FullMatrix<number2> &A)
1040 {
1041  Assert (!this->empty(), ExcEmptyMatrix());
1042 
1043  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
1044  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
1045 
1046  for (size_type i=0; i<m(); ++i)
1047  for (size_type j=0; j<n(); ++j)
1048  (*this)(i,j) += a * number(A(i,j));
1049 }
1050 
1051 
1052 template <typename number>
1053 template <typename number2>
1054 void
1055 FullMatrix<number>::add (const number a,
1056  const FullMatrix<number2> &A,
1057  const number b,
1058  const FullMatrix<number2> &B)
1059 {
1060  Assert (!this->empty(), ExcEmptyMatrix());
1061 
1062  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
1063  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
1064  Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1065  Assert (n() == B.n(), ExcDimensionMismatch(n(), B.n()));
1066 
1067  for (size_type i=0; i<m(); ++i)
1068  for (size_type j=0; j<n(); ++j)
1069  (*this)(i,j) += a * number(A(i,j)) + b * number(B(i,j));
1070 }
1071 
1072 
1073 
1074 template <typename number>
1075 template <typename number2>
1076 void
1077 FullMatrix<number>::add (const number a,
1078  const FullMatrix<number2> &A,
1079  const number b,
1080  const FullMatrix<number2> &B,
1081  const number c,
1082  const FullMatrix<number2> &C)
1083 {
1084  Assert (!this->empty(), ExcEmptyMatrix());
1085 
1086  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
1087  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
1088  Assert (m() == B.m(), ExcDimensionMismatch(m(), B.m()));
1089  Assert (n() == B.n(), ExcDimensionMismatch(n(), B.n()));
1090  Assert (m() == C.m(), ExcDimensionMismatch(m(), C.m()));
1091  Assert (n() == C.n(), ExcDimensionMismatch(n(), C.n()));
1092 
1093 
1094  for (size_type i=0; i<m(); ++i)
1095  for (size_type j=0; j<n(); ++j)
1096  (*this)(i,j) += a * number(A(i,j)) +
1097  b * number(B(i,j)) +
1098  c * number(C(i,j));
1099 }
1100 
1101 
1102 
1103 template <typename number>
1104 template <typename number2>
1106  const number factor,
1107  const size_type dst_offset_i,
1108  const size_type dst_offset_j,
1109  const size_type src_offset_i,
1110  const size_type src_offset_j)
1111 {
1112  Assert (dst_offset_i < m(),
1113  ExcIndexRange (dst_offset_i, 0, m()));
1114  Assert (dst_offset_j < n(),
1115  ExcIndexRange (dst_offset_j, 0, n()));
1116  Assert (src_offset_i < src.m(),
1117  ExcIndexRange (src_offset_i, 0, src.m()));
1118  Assert (src_offset_j < src.n(),
1119  ExcIndexRange (src_offset_j, 0, src.n()));
1120 
1121  // Compute maximal size of copied block
1122  const size_type rows = std::min (m() - dst_offset_i, src.m() - src_offset_i);
1123  const size_type cols = std::min (n() - dst_offset_j, src.n() - src_offset_j);
1124 
1125  for (size_type i=0; i<rows ; ++i)
1126  for (size_type j=0; j<cols ; ++j)
1127  (*this)(dst_offset_i+i,dst_offset_j+j)
1128  += factor * number(src(src_offset_i+i,src_offset_j+j));
1129 }
1130 
1131 
1132 
1133 template <typename number>
1134 template <typename number2>
1136  const number factor,
1137  const size_type dst_offset_i,
1138  const size_type dst_offset_j,
1139  const size_type src_offset_i,
1140  const size_type src_offset_j)
1141 {
1142  Assert (dst_offset_i < m(),
1143  ExcIndexRange (dst_offset_i, 0, m()));
1144  Assert (dst_offset_j < n(),
1145  ExcIndexRange (dst_offset_j, 0, n()));
1146  Assert (src_offset_i < src.n(),
1147  ExcIndexRange (src_offset_i, 0, src.n()));
1148  Assert (src_offset_j < src.m(),
1149  ExcIndexRange (src_offset_j, 0, src.m()));
1150 
1151  // Compute maximal size of copied block
1152  const size_type rows = std::min (m() - dst_offset_i, src.n() - src_offset_j);
1153  const size_type cols = std::min (n() - dst_offset_j,
1154  src.m() - src_offset_i);
1155 
1156 
1157  for (size_type i=0; i<rows ; ++i)
1158  for (size_type j=0; j<cols ; ++j)
1159  (*this)(dst_offset_i+i,dst_offset_j+j)
1160  += factor * number(src(src_offset_i+j,src_offset_j+i));
1161 }
1162 
1163 
1164 
1165 template <typename number>
1166 template <typename number2>
1167 void
1169  const FullMatrix<number2> &A)
1170 {
1171  Assert (!this->empty(), ExcEmptyMatrix());
1172 
1173  Assert (m() == n(), ExcNotQuadratic());
1174  Assert (m() == A.m(), ExcDimensionMismatch(m(), A.m()));
1175  Assert (n() == A.n(), ExcDimensionMismatch(n(), A.n()));
1176 
1177  for (size_type i=0; i<n(); ++i)
1178  for (size_type j=0; j<m(); ++j)
1179  (*this)(i,j) += a * number(A(j,i));
1180 }
1181 
1182 
1183 template <typename number>
1184 bool
1186 {
1187  // simply pass down to the base class
1188  return Table<2,number>::operator==(M);
1189 }
1190 
1191 
1192 template <typename number>
1193 number
1195 {
1196  Assert (!this->empty(), ExcEmptyMatrix());
1197 
1198  Assert (this->n_cols() == this->n_rows(),
1199  ExcDimensionMismatch(this->n_cols(), this->n_rows()));
1200 
1201  switch (this->n_cols())
1202  {
1203  case 1:
1204  return (*this)(0,0);
1205  case 2:
1206  return (*this)(0,0)*(*this)(1,1) - (*this)(1,0)*(*this)(0,1);
1207  case 3:
1208  return ((*this)(0,0)*(*this)(1,1)*(*this)(2,2)
1209  -(*this)(0,0)*(*this)(1,2)*(*this)(2,1)
1210  -(*this)(1,0)*(*this)(0,1)*(*this)(2,2)
1211  +(*this)(1,0)*(*this)(0,2)*(*this)(2,1)
1212  +(*this)(2,0)*(*this)(0,1)*(*this)(1,2)
1213  -(*this)(2,0)*(*this)(0,2)*(*this)(1,1));
1214  default:
1215  Assert (false, ExcNotImplemented());
1216  return 0;
1217  };
1218 }
1219 
1220 
1221 
1222 template <typename number>
1223 number
1225 {
1226  Assert (!this->empty(), ExcEmptyMatrix());
1227 
1228  Assert (this->n_cols() == this->n_rows(),
1229  ExcDimensionMismatch(this->n_cols(), this->n_rows()));
1230 
1231  number tr = 0;
1232  for (size_type i=0; i<this->n_rows(); ++i)
1233  tr += (*this)(i,i);
1234 
1235  return tr;
1236 }
1237 
1238 
1239 
1240 template <typename number>
1243 {
1244  Assert (!this->empty(), ExcEmptyMatrix());
1245 
1246  real_type s = 0.;
1247  for (size_type i=0; i<this->n_rows()*this->n_cols(); ++i)
1248  s += numbers::NumberTraits<number>::abs_square(this->values[i]);
1249  return std::sqrt(s);
1250 }
1251 
1252 
1253 
1254 template <typename number>
1257 {
1258  Assert (!this->empty(), ExcEmptyMatrix());
1259 
1260  real_type s = 0.;
1261  real_type a = 0.;
1262  for (size_type i=0; i<this->n_rows(); ++i)
1263  for (size_type j=0; j<this->n_cols(); ++j)
1264  {
1265  const number x_ij = (*this)(i,j);
1266  const number x_ji = (*this)(j,i);
1267 
1270  }
1271 
1272  if (s!=0.)
1273  return std::sqrt(a)/std::sqrt(s);
1274  return 0;
1275 }
1276 
1277 
1278 
1279 template <typename number>
1280 template <typename number2>
1281 void
1283 {
1284  Assert (!this->empty(), ExcEmptyMatrix());
1285 
1286  Assert (this->n_cols() == this->n_rows(),
1287  ExcNotQuadratic());
1288  Assert (this->n_cols() == M.n_cols(),
1289  ExcDimensionMismatch(this->n_cols(), M.n_cols()));
1290  Assert (this->n_rows() == M.n_rows(),
1291  ExcDimensionMismatch(this->n_rows(), M.n_rows()));
1292 
1293  if (PointerComparison::equal(&M, this))
1294  {
1295  // avoid overwriting source
1296  // by destination matrix:
1297  FullMatrix<number2> M2 = M;
1298  invert(M2);
1299  }
1300  else
1301  switch (this->n_cols())
1302  {
1303  case 1:
1304  (*this)(0,0) = number2(1.0)/M(0,0);
1305  return;
1306  case 2:
1307  // this is Maple output,
1308  // thus a bit unstructured
1309  {
1310  const number2 t4 = number2(1.0)/(M(0,0)*M(1,1)-M(0,1)*M(1,0));
1311  (*this)(0,0) = M(1,1)*t4;
1312  (*this)(0,1) = -M(0,1)*t4;
1313  (*this)(1,0) = -M(1,0)*t4;
1314  (*this)(1,1) = M(0,0)*t4;
1315  return;
1316  };
1317 
1318  case 3:
1319  {
1320  const number2 t4 = M(0,0)*M(1,1),
1321  t6 = M(0,0)*M(1,2),
1322  t8 = M(0,1)*M(1,0),
1323  t00 = M(0,2)*M(1,0),
1324  t01 = M(0,1)*M(2,0),
1325  t04 = M(0,2)*M(2,0),
1326  t07 = number2(1.0)/(t4*M(2,2)-t6*M(2,1)-t8*M(2,2)+
1327  t00*M(2,1)+t01*M(1,2)-t04*M(1,1));
1328  (*this)(0,0) = (M(1,1)*M(2,2)-M(1,2)*M(2,1))*t07;
1329  (*this)(0,1) = -(M(0,1)*M(2,2)-M(0,2)*M(2,1))*t07;
1330  (*this)(0,2) = -(-M(0,1)*M(1,2)+M(0,2)*M(1,1))*t07;
1331  (*this)(1,0) = -(M(1,0)*M(2,2)-M(1,2)*M(2,0))*t07;
1332  (*this)(1,1) = (M(0,0)*M(2,2)-t04)*t07;
1333  (*this)(1,2) = -(t6-t00)*t07;
1334  (*this)(2,0) = -(-M(1,0)*M(2,1)+M(1,1)*M(2,0))*t07;
1335  (*this)(2,1) = -(M(0,0)*M(2,1)-t01)*t07;
1336  (*this)(2,2) = (t4-t8)*t07;
1337  return;
1338  };
1339 
1340  case 4:
1341  {
1342  // with (linalg);
1343  // a:=matrix(4,4);
1344  // evalm(a);
1345  // ai:=inverse(a);
1346  // readlib(C);
1347  // C(ai,optimized,filename=x4);
1348 
1349  const number2 t14 = M(0,0)*M(1,1);
1350  const number2 t15 = M(2,2)*M(3,3);
1351  const number2 t17 = M(2,3)*M(3,2);
1352  const number2 t19 = M(0,0)*M(2,1);
1353  const number2 t20 = M(1,2)*M(3,3);
1354  const number2 t22 = M(1,3)*M(3,2);
1355  const number2 t24 = M(0,0)*M(3,1);
1356  const number2 t25 = M(1,2)*M(2,3);
1357  const number2 t27 = M(1,3)*M(2,2);
1358  const number2 t29 = M(1,0)*M(0,1);
1359  const number2 t32 = M(1,0)*M(2,1);
1360  const number2 t33 = M(0,2)*M(3,3);
1361  const number2 t35 = M(0,3)*M(3,2);
1362  const number2 t37 = M(1,0)*M(3,1);
1363  const number2 t38 = M(0,2)*M(2,3);
1364  const number2 t40 = M(0,3)*M(2,2);
1365  const number2 t42 = t14*t15-t14*t17-t19*t20+t19*t22+
1366  t24*t25-t24*t27-t29*t15+t29*t17+
1367  t32*t33-t32*t35-t37*t38+t37*t40;
1368  const number2 t43 = M(2,0)*M(0,1);
1369  const number2 t46 = M(2,0)*M(1,1);
1370  const number2 t49 = M(2,0)*M(3,1);
1371  const number2 t50 = M(0,2)*M(1,3);
1372  const number2 t52 = M(0,3)*M(1,2);
1373  const number2 t54 = M(3,0)*M(0,1);
1374  const number2 t57 = M(3,0)*M(1,1);
1375  const number2 t60 = M(3,0)*M(2,1);
1376  const number2 t63 = t43*t20-t43*t22-t46*t33+t46*t35+
1377  t49*t50-t49*t52-t54*t25+t54*t27+
1378  t57*t38-t57*t40-t60*t50+t60*t52;
1379  const number2 t65 = number2(1.)/(t42+t63);
1380  const number2 t71 = M(0,2)*M(2,1);
1381  const number2 t73 = M(0,3)*M(2,1);
1382  const number2 t75 = M(0,2)*M(3,1);
1383  const number2 t77 = M(0,3)*M(3,1);
1384  const number2 t81 = M(0,1)*M(1,2);
1385  const number2 t83 = M(0,1)*M(1,3);
1386  const number2 t85 = M(0,2)*M(1,1);
1387  const number2 t87 = M(0,3)*M(1,1);
1388  const number2 t101 = M(1,0)*M(2,2);
1389  const number2 t103 = M(1,0)*M(2,3);
1390  const number2 t105 = M(2,0)*M(1,2);
1391  const number2 t107 = M(2,0)*M(1,3);
1392  const number2 t109 = M(3,0)*M(1,2);
1393  const number2 t111 = M(3,0)*M(1,3);
1394  const number2 t115 = M(0,0)*M(2,2);
1395  const number2 t117 = M(0,0)*M(2,3);
1396  const number2 t119 = M(2,0)*M(0,2);
1397  const number2 t121 = M(2,0)*M(0,3);
1398  const number2 t123 = M(3,0)*M(0,2);
1399  const number2 t125 = M(3,0)*M(0,3);
1400  const number2 t129 = M(0,0)*M(1,2);
1401  const number2 t131 = M(0,0)*M(1,3);
1402  const number2 t133 = M(1,0)*M(0,2);
1403  const number2 t135 = M(1,0)*M(0,3);
1404  (*this)(0,0) = (M(1,1)*M(2,2)*M(3,3)-M(1,1)*M(2,3)*M(3,2)-
1405  M(2,1)*M(1,2)*M(3,3)+M(2,1)*M(1,3)*M(3,2)+
1406  M(3,1)*M(1,2)*M(2,3)-M(3,1)*M(1,3)*M(2,2))*t65;
1407  (*this)(0,1) = -(M(0,1)*M(2,2)*M(3,3)-M(0,1)*M(2,3)*M(3,2)-
1408  t71*M(3,3)+t73*M(3,2)+t75*M(2,3)-t77*M(2,2))*t65;
1409  (*this)(0,2) = (t81*M(3,3)-t83*M(3,2)-t85*M(3,3)+t87*M(3,2)+
1410  t75*M(1,3)-t77*M(1,2))*t65;
1411  (*this)(0,3) = -(t81*M(2,3)-t83*M(2,2)-t85*M(2,3)+t87*M(2,2)+
1412  t71*M(1,3)-t73*M(1,2))*t65;
1413  (*this)(1,0) = -(t101*M(3,3)-t103*M(3,2)-t105*M(3,3)+t107*M(3,2)+
1414  t109*M(2,3)-t111*M(2,2))*t65;
1415  (*this)(1,1) = (t115*M(3,3)-t117*M(3,2)-t119*M(3,3)+t121*M(3,2)+
1416  t123*M(2,3)-t125*M(2,2))*t65;
1417  (*this)(1,2) = -(t129*M(3,3)-t131*M(3,2)-t133*M(3,3)+t135*M(3,2)+
1418  t123*M(1,3)-t125*M(1,2))*t65;
1419  (*this)(1,3) = (t129*M(2,3)-t131*M(2,2)-t133*M(2,3)+t135*M(2,2)+
1420  t119*M(1,3)-t121*M(1,2))*t65;
1421  (*this)(2,0) = (t32*M(3,3)-t103*M(3,1)-t46*M(3,3)+t107*M(3,1)+
1422  t57*M(2,3)-t111*M(2,1))*t65;
1423  (*this)(2,1) = -(t19*M(3,3)-t117*M(3,1)-t43*M(3,3)+t121*M(3,1)+
1424  t54*M(2,3)-t125*M(2,1))*t65;
1425  (*this)(2,2) = (t14*M(3,3)-t131*M(3,1)-t29*M(3,3)+t135*M(3,1)+
1426  t54*M(1,3)-t125*M(1,1))*t65;
1427  (*this)(2,3) = -(t14*M(2,3)-t131*M(2,1)-t29*M(2,3)+t135*M(2,1)+
1428  t43*M(1,3)-t121*M(1,1))*t65;
1429  (*this)(3,0) = -(t32*M(3,2)-t101*M(3,1)-t46*M(3,2)+t105*M(3,1)+
1430  t57*M(2,2)-t109*M(2,1))*t65;
1431  (*this)(3,1) = (t19*M(3,2)-t115*M(3,1)-t43*M(3,2)+t119*M(3,1)+
1432  t54*M(2,2)-t123*M(2,1))*t65;
1433  (*this)(3,2) = -(t14*M(3,2)-t129*M(3,1)-t29*M(3,2)+t133*M(3,1)+
1434  t54*M(1,2)-t123*M(1,1))*t65;
1435  (*this)(3,3) = (t14*M(2,2)-t129*M(2,1)-t29*M(2,2)+t133*M(2,1)+
1436  t43*M(1,2)-t119*M(1,1))*t65;
1437 
1438  break;
1439  }
1440 
1441 
1442  default:
1443  // if no inversion is
1444  // hardcoded, fall back
1445  // to use the
1446  // Gauss-Jordan algorithm
1447  *this = M;
1448  gauss_jordan();
1449  };
1450 }
1451 
1452 
1453 template <typename number>
1454 template <typename number2>
1455 void
1457 {
1458  Assert (!A.empty(), ExcEmptyMatrix());
1459  Assert (A.n() == A.m(),
1460  ExcNotQuadratic());
1461  // Matrix must be symmetric.
1462  Assert(A.relative_symmetry_norm2() < 1.0e-10, ExcMessage("A must be symmetric."));
1463 
1464  if (PointerComparison::equal(&A, this))
1465  {
1466  // avoid overwriting source
1467  // by destination matrix:
1468  FullMatrix<number2> A2 = A;
1469  cholesky(A2);
1470  }
1471  else
1472  {
1473  /* reinit *this to 0 */
1474  this->reinit(A.m(), A.n());
1475 
1476  double SLik2 = 0.0, SLikLjk = 0.0;
1477  for (size_type i=0; i< this->n_cols(); i++)
1478  {
1479  SLik2 = 0.0;
1480  for (size_type j = 0; j < i; j++)
1481  {
1482  SLikLjk = 0.0;
1483  for (size_type k =0; k<j; k++)
1484  {
1485  SLikLjk += (*this)(i,k)*(*this)(j,k);
1486  };
1487  (*this)(i,j) = (1./(*this)(j,j))*(A(i,j) - SLikLjk);
1488  SLik2 += (*this)(i,j)*(*this)(i,j);
1489  }
1490  AssertThrow (A(i,i) - SLik2 >= 0,
1491  ExcMatrixNotPositiveDefinite());
1492 
1493  (*this)(i,i) = std::sqrt(A(i,i) - SLik2);
1494  }
1495  }
1496 }
1497 
1498 
1499 template <typename number>
1500 template <typename number2>
1501 void
1503  const Vector<number2> &W)
1504 {
1505  Assert (V.size() == W.size(), ExcMessage("Vectors V, W must be the same size."));
1506  this->reinit(V.size(), V.size());
1507 
1508  for (size_type i = 0; i<this->n(); i++)
1509  {
1510  for (size_type j = 0; j< this->n(); j++)
1511  {
1512  (*this)(i,j) = V(i)*W(j);
1513  }
1514  }
1515 }
1516 
1517 
1518 template <typename number>
1519 template <typename number2>
1520 void
1522 {
1523  Assert (!A.empty(), ExcEmptyMatrix());
1524  Assert(A.m()>A.n(), ExcDimensionMismatch(A.m(), A.n()));
1525  Assert(this->m()==A.n(), ExcDimensionMismatch(this->m(), A.n()));
1526  Assert(this->n()==A.m(), ExcDimensionMismatch(this->n(), A.m()));
1527 
1528  FullMatrix<number2> A_t(A.n(),A.m());
1529  FullMatrix<number2> A_t_times_A(A.n(),A.n());
1530  FullMatrix<number2> A_t_times_A_inv(A.n(),A.n());
1531  FullMatrix<number2> left_inv(A.n(),A.m());
1532 
1533  A_t.Tadd(A,1);
1534  A_t.mmult(A_t_times_A,A);
1535  if (number(A_t_times_A.determinant())==number(0))
1536  Assert(false, ExcSingular())
1537  else
1538  {
1539  A_t_times_A_inv.invert(A_t_times_A);
1540  A_t_times_A_inv.mmult(left_inv,A_t);
1541 
1542  *this=left_inv;
1543  }
1544 }
1545 
1546 template <typename number>
1547 template <typename number2>
1548 void
1550 {
1551  Assert (!A.empty(), ExcEmptyMatrix());
1552  Assert(A.n()>A.m(), ExcDimensionMismatch(A.n(), A.m()));
1553  Assert(this->m()==A.n(), ExcDimensionMismatch(this->m(), A.n()));
1554  Assert(this->n()==A.m(), ExcDimensionMismatch(this->n(), A.m()));
1555 
1556  FullMatrix<number> A_t(A.n(),A.m());
1557  FullMatrix<number> A_times_A_t(A.m(),A.m());
1558  FullMatrix<number> A_times_A_t_inv(A.m(),A.m());
1559  FullMatrix<number> right_inv(A.n(),A.m());
1560 
1561  A_t.Tadd(A,1);
1562  A.mmult(A_times_A_t,A_t);
1563  if (number(A_times_A_t.determinant())==number(0))
1564  Assert(false, ExcSingular())
1565  else
1566  {
1567  A_times_A_t_inv.invert(A_times_A_t);
1568  A_t.mmult(right_inv,A_times_A_t_inv);
1569 
1570  *this=right_inv;
1571  }
1572 }
1573 
1574 
1575 template <typename number>
1576 template <int dim>
1577 void
1579  const size_type src_r_i,
1580  const size_type src_r_j,
1581  const size_type src_c_i,
1582  const size_type src_c_j,
1583  const size_type dst_r,
1584  const size_type dst_c)
1585 {
1586 
1587  Assert (!this->empty(), ExcEmptyMatrix());
1588  Assert(this->m()-dst_r>src_r_j-src_r_i,
1589  ExcIndexRange(this->m()-dst_r,0,src_r_j-src_r_i));
1590  Assert(this->n()-dst_c>src_c_j-src_c_i,
1591  ExcIndexRange(this->n()-dst_c,0,src_c_j-src_c_i));
1592  Assert(dim>src_r_j, ExcIndexRange(dim,0,src_r_j));
1593  Assert(dim>src_c_j, ExcIndexRange(dim,0,src_r_j));
1594  Assert(src_r_j>=src_r_i, ExcIndexRange(src_r_j,0,src_r_i));
1595  Assert(src_c_j>=src_c_i, ExcIndexRange(src_r_j,0,src_r_i));
1596 
1597  for (size_type i=0; i<src_r_j-src_r_i+1; i++)
1598  for (size_type j=0; j<src_c_j-src_c_i+1; j++)
1599  (*this)(i+dst_r,j+dst_c) = number(T[i+src_r_i][j+src_c_i]);
1600 
1601 }
1602 
1603 
1604 template <typename number>
1605 template <int dim>
1606 void
1608  const size_type src_r_i,
1609  const size_type src_r_j,
1610  const size_type src_c_i,
1611  const size_type src_c_j,
1612  const size_type dst_r,
1613  const size_type dst_c) const
1614 {
1615  Assert (!this->empty(), ExcEmptyMatrix());
1616  Assert(dim-dst_r>src_r_j-src_r_i,
1617  ExcIndexRange(dim-dst_r,0,src_r_j-src_r_i));
1618  Assert(dim-dst_c>src_c_j-src_c_i,
1619  ExcIndexRange(dim-dst_c,0,src_c_j-src_c_i));
1620  Assert(this->m()>src_r_j, ExcIndexRange(dim,0,src_r_j));
1621  Assert(this->n()>src_c_j, ExcIndexRange(dim,0,src_r_j));
1622  Assert(src_r_j>=src_r_i, ExcIndexRange(src_r_j,0,src_r_i));
1623  Assert(src_c_j>=src_c_i, ExcIndexRange(src_r_j,0,src_r_i));
1624 
1625 
1626  for (size_type i=0; i<src_r_j-src_r_i+1; i++)
1627  for (size_type j=0; j<src_c_j-src_c_i+1; j++)
1628  T[i+dst_r][j+dst_c] = double ((*this)(i+src_r_i,j+src_c_i));
1629 }
1630 
1631 
1632 
1633 template <typename number>
1634 template <typename somenumber>
1635 void
1637  const Vector<somenumber> &src,
1638  const number om) const
1639 {
1640  Assert (m() == n(), ExcNotQuadratic());
1641  Assert (dst.size() == n(), ExcDimensionMismatch (dst.size(), n()));
1642  Assert (src.size() == n(), ExcDimensionMismatch (src.size(), n()));
1643 
1644  const size_t n = src.size();
1645  somenumber *dst_ptr = dst.begin();
1646  const somenumber *src_ptr = src.begin();
1647 
1648  for (size_type i=0; i<n; ++i, ++dst_ptr, ++src_ptr)
1649  *dst_ptr = somenumber(om) **src_ptr / somenumber((*this)(i,i));
1650 }
1651 
1652 
1653 
1654 template <typename number>
1655 void
1657  std::ostream &out,
1658  const unsigned int precision,
1659  const bool scientific,
1660  const unsigned int width_,
1661  const char *zero_string,
1662  const double denominator,
1663  const double threshold) const
1664 {
1665  unsigned int width = width_;
1666 
1667  Assert ((!this->empty()) || (this->n_cols()+this->n_rows()==0),
1668  ExcInternalError());
1669 
1670  // set output format, but store old
1671  // state
1672  std::ios::fmtflags old_flags = out.flags();
1673  unsigned int old_precision = out.precision (precision);
1674 
1675  if (scientific)
1676  {
1677  out.setf (std::ios::scientific, std::ios::floatfield);
1678  if (!width)
1679  width = precision+7;
1680  }
1681  else
1682  {
1683  out.setf (std::ios::fixed, std::ios::floatfield);
1684  if (!width)
1685  width = precision+2;
1686  }
1687 
1688  for (size_type i=0; i<m(); ++i)
1689  {
1690  for (size_type j=0; j<n(); ++j)
1691  if (std::abs((*this)(i,j)) > threshold)
1692  out << std::setw(width)
1693  << (*this)(i,j) * number(denominator) << ' ';
1694  else
1695  out << std::setw(width) << zero_string << ' ';
1696  out << std::endl;
1697  };
1698 
1699  AssertThrow (out, ExcIO());
1700  // reset output format
1701  out.flags (old_flags);
1702  out.precision(old_precision);
1703 }
1704 
1705 
1706 template <typename number>
1707 void
1709 {
1710  Assert (!this->empty(), ExcEmptyMatrix());
1711  Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
1712 
1713  // see if we can use Lapack algorithms
1714  // for this and if the type for 'number'
1715  // works for us (it is usually not
1716  // efficient to use Lapack for very small
1717  // matrices):
1718 #if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && \
1719  defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
1721  ||
1723  if (this->n_cols() > 15)
1724  {
1725  // In case we have the LAPACK functions
1726  // getrf and getri detected at configure,
1727  // we use these algorithms for inversion
1728  // since they provide better performance
1729  // than the deal.II native functions.
1730  //
1731  // Note that BLAS/LAPACK stores matrix
1732  // elements column-wise (i.e., all values in
1733  // one column, then all in the next, etc.),
1734  // whereas the FullMatrix stores them
1735  // row-wise.
1736  // We ignore that difference, and give our
1737  // row-wise data to LAPACK,
1738  // let LAPACK build the inverse of the
1739  // transpose matrix, and read the result as
1740  // if it were row-wise again. In other words,
1741  // we just got ((A^T)^{-1})^T, which is
1742  // A^{-1}.
1743 
1744  const int nn = this->n();
1745 
1746  // workspace for permutations
1747  std::vector<int> ipiv(nn);
1748  int info;
1749 
1750  // Use the LAPACK function getrf for
1751  // calculating the LU factorization.
1752  getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
1753 
1754  Assert(info >= 0, ExcInternalError());
1755  Assert(info == 0, LACExceptions::ExcSingular());
1756 
1757  // scratch array
1758  std::vector<number> inv_work (nn);
1759 
1760  // Use the LAPACK function getri for
1761  // calculating the actual inverse using
1762  // the LU factorization.
1763  getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
1764 
1765  Assert(info >= 0, ExcInternalError());
1766  Assert(info == 0, LACExceptions::ExcSingular());
1767 
1768  return;
1769  }
1770 
1771 #endif
1772 
1773  // otherwise do it by hand. use the
1774  // Gauss-Jordan-Algorithmus from
1775  // Stoer & Bulirsch I (4th Edition)
1776  // p. 153
1777  const size_type N = n();
1778 
1779  // first get an estimate of the
1780  // size of the elements of this
1781  // matrix, for later checks whether
1782  // the pivot element is large
1783  // enough, or whether we have to
1784  // fear that the matrix is not
1785  // regular
1786  double diagonal_sum = 0;
1787  for (size_type i=0; i<N; ++i)
1788  diagonal_sum += std::abs((*this)(i,i));
1789  const double typical_diagonal_element = diagonal_sum/N;
1790 
1791  // initialize the array that holds
1792  // the permutations that we find
1793  // during pivot search
1794  std::vector<size_type> p(N);
1795  for (size_type i=0; i<N; ++i)
1796  p[i] = i;
1797 
1798  for (size_type j=0; j<N; ++j)
1799  {
1800  // pivot search: search that
1801  // part of the line on and
1802  // right of the diagonal for
1803  // the largest element
1804  real_type max = std::abs((*this)(j,j));
1805  size_type r = j;
1806  for (size_type i=j+1; i<N; ++i)
1807  {
1808  if (std::abs((*this)(i,j)) > max)
1809  {
1810  max = std::abs((*this)(i,j));
1811  r = i;
1812  }
1813  }
1814  // check whether the pivot is
1815  // too small
1816  Assert(max > 1.e-16*typical_diagonal_element,
1817  ExcNotRegular(max));
1818 
1819  // row interchange
1820  if (r>j)
1821  {
1822  for (size_type k=0; k<N; ++k)
1823  std::swap ((*this)(j,k), (*this)(r,k));
1824 
1825  std::swap (p[j], p[r]);
1826  }
1827 
1828  // transformation
1829  const number hr = number(1.)/(*this)(j,j);
1830  (*this)(j,j) = hr;
1831  for (size_type k=0; k<N; ++k)
1832  {
1833  if (k==j) continue;
1834  for (size_type i=0; i<N; ++i)
1835  {
1836  if (i==j) continue;
1837  (*this)(i,k) -= (*this)(i,j)*(*this)(j,k)*hr;
1838  }
1839  }
1840  for (size_type i=0; i<N; ++i)
1841  {
1842  (*this)(i,j) *= hr;
1843  (*this)(j,i) *= -hr;
1844  }
1845  (*this)(j,j) = hr;
1846  }
1847  // column interchange
1848  std::vector<number> hv(N);
1849  for (size_type i=0; i<N; ++i)
1850  {
1851  for (size_type k=0; k<N; ++k)
1852  hv[p[k]] = (*this)(i,k);
1853  for (size_type k=0; k<N; ++k)
1854  (*this)(i,k) = hv[k];
1855  }
1856 }
1857 
1858 
1859 
1860 template <typename number>
1861 std::size_t
1863 {
1864  return (sizeof(*this) - sizeof (Table<2,number>)
1865  +
1867 }
1868 
1869 
1870 DEAL_II_NAMESPACE_CLOSE
1871 
1872 #endif
unsigned int n_rows() const
real_type frobenius_norm() const
number residual(Vector< number2 > &dst, const Vector< number2 > &x, const Vector< number3 > &b) const
void add_row(const size_type i, const number s, const size_type j)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
void copy_from(const MATRIX &)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
bool empty() 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
void equ(const number a, const FullMatrix< number2 > &A)
size_type m() const
numbers::NumberTraits< number >::real_type real_type
Definition: full_matrix.h:102
unsigned int size_type
Definition: full_matrix.h:74
::ExceptionBase & ExcMessage(std::string arg1)
FullMatrix & operator/=(const number factor)
void copy_to(Tensor< 2, dim > &T, const size_type src_r_i=0, const size_type src_r_j=dim-1, const size_type src_c_i=0, const size_type src_c_j=dim-1, const size_type dst_r=0, const size_type dst_c=0) const
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
bool all_zero() const
size_type n() const
void left_invert(const FullMatrix< number2 > &M)
void fill_permutation(const FullMatrix< number2 > &src, const std::vector< size_type > &p_rows, const std::vector< size_type > &p_cols)
bool is_finite(const double x)
void add_col(const size_type i, const number s, const size_type j)
void triple_product(const FullMatrix< number > &A, const FullMatrix< number > &B, const FullMatrix< number > &D, const bool transpose_B=false, const bool transpose_D=false, const number scaling=number(1.))
void cholesky(const FullMatrix< number2 > &A)
void Tmmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix< number > & operator=(const FullMatrix< number > &)
::ExceptionBase & ExcIO()
void swap_row(const size_type i, const size_type j)
FullMatrix(const size_type n=0)
void fill(const FullMatrix< number2 > &src, const size_type dst_offset_i=0, const size_type dst_offset_j=0, const size_type src_offset_i=0, const size_type src_offset_j=0)
number determinant() const
Tensor< 2, dim, Number > invert(const Tensor< 2, dim, Number > &t)
Definition: tensor.h:1781
void forward(Vector< number2 > &dst, const Vector< number2 > &src) const
static real_type abs_square(const number &x)
Definition: numbers.h:307
#define Assert(cond, exc)
Definition: exceptions.h:299
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
static bool equal(const T *p1, const T *p2)
void diagadd(const number s)
void invert(const FullMatrix< number2 > &M)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void Tadd(const number s, const FullMatrix< number2 > &B)
TableBase< N, T > & operator=(const TableBase< N, T > &src)
iterator begin()
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
real_type linfty_norm() const
void outer_product(const Vector< number2 > &V, const Vector< number2 > &W)
real_type relative_symmetry_norm2() const
void right_invert(const FullMatrix< number2 > &M)
bool operator==(const TableBase< N, T > &T2) const
void add(const number a, const FullMatrix< number2 > &A)
std::size_t memory_consumption() const
::ExceptionBase & ExcNumberNotFinite()
bool operator==(const FullMatrix< number > &) const
real_type l1_norm() const
void swap_col(const size_type i, const size_type j)
number2 matrix_norm_square(const Vector< number2 > &v) const
void TmTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix & operator*=(const number factor)
::ExceptionBase & ExcNotImplemented()
number trace() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
Definition: table.h:33
::ExceptionBase & ExcInternalError()
number2 matrix_scalar_product(const Vector< number2 > &u, const Vector< number2 > &v) const
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
unsigned int n_cols() const
VectorizedArray< Number > max(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
std::size_t size() const