17 #ifndef __deal2__full_matrix_templates_h
18 #define __deal2__full_matrix_templates_h
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>
36 DEAL_II_NAMESPACE_OPEN
39 template <
typename number>
46 template <
typename number>
50 Table<2,number> (m, n)
54 template <
typename number>
57 const number *entries)
59 Table<2,number> (m, n)
65 template <
typename number>
73 template <
typename number>
76 Table<2,number> (id.m(), id.n())
83 template <
typename number>
92 template <
typename number>
93 template <
typename number2>
103 template <
typename number>
107 this->reinit (
id.m(),
id.n());
116 template <
typename number>
117 template <
typename number2>
125 (*
this)(i,j) = M(i,j);
132 template <
typename number>
136 Assert (!this->empty(), ExcEmptyMatrix());
138 const number *p = &this->values[0];
139 const number *
const e = &this->values[0] + this->n_elements();
141 if (*p++ != number(0.0))
149 template <
typename number>
156 number *p = &(*this)(0,0);
157 const number *e = &(*this)(0,0) + n()*m();
166 template <
typename number>
173 number *p = &(*this)(0,0);
174 const number *e = &(*this)(0,0) + n()*m();
176 const number factor_inv = number(1.)/factor;
188 template <
typename number>
189 template <
typename number2>
193 const bool adding)
const
195 Assert (!this->empty(), ExcEmptyMatrix());
200 Assert (&src != &dst, ExcSourceEqualsDestination());
202 const number *e = &this->values[0];
207 const size_type size_m = m(), size_n = n();
210 number2 s = adding ? dst(i) : 0.;
212 s += src_ptr[j] * number2(*(e++));
219 template <
typename number>
220 template <
typename number2>
223 const bool adding)
const
225 Assert (!this->empty(), ExcEmptyMatrix());
230 Assert (&src != &dst, ExcSourceEqualsDestination());
232 const number *e = &this->values[0];
233 number2 *dst_ptr = &dst(0);
234 const size_type size_m = m(), size_n = n();
245 const number2 d = src(i);
247 dst_ptr[j] += d * number2(*(e++));
252 template <
typename number>
253 template <
typename number2,
typename number3>
258 Assert (!this->empty(), ExcEmptyMatrix());
264 Assert (&src != &dst, ExcSourceEqualsDestination());
271 number s = number(right(i));
273 s -= number(src(j)) * (*this)(i,j);
277 return std::sqrt(res);
282 template <
typename number>
283 template <
typename number2>
287 Assert (!this->empty(), ExcEmptyMatrix());
296 number s = number(src(i));
298 s -= number(dst(j)) * (*this)(i,j);
299 dst(i) = s/(*this)(i,i);
306 template <
typename number>
307 template <
typename number2>
311 Assert (!this->empty(), ExcEmptyMatrix());
315 for (
int i=nu-1; i>=0; --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));
327 template <
typename number>
328 template <
typename number2>
335 Assert (dst_offset_i < m(),
337 Assert (dst_offset_j < n(),
339 Assert (src_offset_i < src.
m(),
341 Assert (src_offset_j < src.
n(),
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);
352 (*
this)(dst_offset_i+i,dst_offset_j+j)
353 = src(src_offset_i+i,src_offset_j+j);
357 template <
typename number>
358 template <
typename number2>
360 const std::vector<size_type> &p_rows,
361 const std::vector<size_type> &p_cols)
363 Assert (p_rows.size() == this->n_rows(),
365 Assert (p_cols.size() == this->n_cols(),
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]);
385 template <
typename number>
390 Assert (!this->empty(), ExcEmptyMatrix());
393 (*
this)(i,k) += s*(*
this)(j,k);
397 template <
typename number>
404 Assert (!this->empty(), ExcEmptyMatrix());
408 (*
this)(i,l) += s*(*
this)(j,l) + t*(*
this)(k,l);
412 template <
typename number>
416 Assert (!this->empty(), ExcEmptyMatrix());
419 (*
this)(k,i) += s*(*
this)(k,j);
423 template <
typename number>
428 Assert (!this->empty(), ExcEmptyMatrix());
430 for (
size_t l=0; l<n(); ++l)
431 (*
this)(l,i) += s*(*
this)(l,j) + t*(*
this)(l,k);
436 template <
typename number>
440 Assert (!this->empty(), ExcEmptyMatrix());
443 std::swap ((*
this)(i,k),
448 template <
typename number>
452 Assert (!this->empty(), ExcEmptyMatrix());
455 std::swap ((*
this)(k,i),
460 template <
typename number>
463 Assert (!this->empty(), ExcEmptyMatrix());
471 template <
typename number>
472 template <
typename number2>
476 Assert (!this->empty(), ExcEmptyMatrix());
483 (*
this)(i,j) = a * number(A(i,j));
487 template <
typename number>
488 template <
typename number2>
495 Assert (!this->empty(), ExcEmptyMatrix());
504 (*
this)(i,j) = a * number(A(i,j)) + b * number(B(i,j));
508 template <
typename number>
509 template <
typename number2>
518 Assert (!this->empty(), ExcEmptyMatrix());
529 (*
this)(i,j) = a * number(A(i,j)) +
536 template <
typename number>
537 template <
typename number2>
540 const bool adding)
const
542 Assert (!this->empty(), ExcEmptyMatrix());
550 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
556 if (this->n()*this->m()*src.
n() > 300)
570 const int m = src.
n();
571 const int n = this->m();
572 const int k = this->n();
573 const char *notrans =
"n";
575 const number alpha = 1.;
576 const number beta = (adding ==
true) ? 1. : 0.;
581 gemm(notrans, notrans, &m, &n, &k, &alpha, &src(0,0), &m,
582 &this->values[0], &k, &beta, &dst(0,0), &m);
589 const size_type m = this->m(), n = src.
n(), l = this->n();
597 number2 add_value = adding ? dst(i,j) : 0.;
599 add_value += (number2)(*this)(i,k) * (number2)(src(k,j));
600 dst(i,j) = add_value;
606 template <
typename number>
607 template <
typename number2>
610 const bool adding)
const
612 Assert (!this->empty(), ExcEmptyMatrix());
621 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
627 if (this->n()*this->m()*src.
n() > 300)
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";
647 const number alpha = 1.;
648 const number beta = (adding ==
true) ? 1. : 0.;
652 gemm(notrans, trans, &m, &n, &k, &alpha, &src(0,0), &m,
653 &this->values[0], &n, &beta, &dst(0,0), &m);
660 const size_type m = n(), n = src.
n(), l = this->m();
667 number2 add_value = 0.;
669 add_value += (number2)(*this)(k,i) * (number2)(*this)(k,j);
672 dst(i,j) += add_value;
674 dst(j,i) += add_value;
677 dst(i,j) = dst(j,i) = add_value;
688 number2 add_value = adding ? dst(i,j) : 0.;
690 add_value += (number2)(*this)(k,i) * (number2)(src(k,j));
691 dst(i,j) = add_value;
697 template <
typename number>
698 template <
typename number2>
701 const bool adding)
const
703 Assert (!this->empty(), ExcEmptyMatrix());
711 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
717 if (this->n()*this->m()*src.
m() > 300)
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";
737 const number alpha = 1.;
738 const number beta = (adding ==
true) ? 1. : 0.;
742 gemm(trans, notrans, &m, &n, &k, &alpha, &src(0,0), &k,
743 &this->values[0], &k, &beta, &dst(0,0), &m);
750 const size_type m = this->m(), n = src.
m(), l = this->n();
757 number2 add_value = 0.;
759 add_value += (number2)(*this)(i,k) * (number2)(*this)(j,k);
762 dst(i,j) += add_value;
764 dst(j,i) += add_value;
767 dst(i,j) = dst(j,i) = add_value;
775 number2 add_value = adding ? dst(i,j) : 0.;
777 add_value += (number2)(*this)(i,k) * (number2)(src(j,k));
778 dst(i,j) = add_value;
784 template <
typename number>
785 template <
typename number2>
788 const bool adding)
const
790 Assert (!this->empty(), ExcEmptyMatrix());
799 #if defined(HAVE_DGEMM_) && defined (HAVE_SGEMM_)
805 if (this->n()*this->m()*src.
m() > 300)
819 const int m = src.
n();
820 const int n = this->n();
821 const int k = this->m();
822 const char *trans =
"t";
824 const number alpha = 1.;
825 const number beta = (adding ==
true) ? 1. : 0.;
829 gemm(trans, trans, &m, &n, &k, &alpha, &src(0,0), &k,
830 &this->values[0], &n, &beta, &dst(0,0), &m);
837 const size_type m = n(), n = src.
m(), l = this->m();
847 number2 add_value = adding ? dst(i,j) : 0.;
849 add_value += (number2)(*this)(k,i) * (number2)(src(j,k));
850 dst(i,j) = add_value;
855 template <
typename number>
861 const bool transpose_B,
862 const bool transpose_D,
863 const number scaling)
895 ADij += A(i,k)*D(j,k);
898 ADij += A(i,k)*D(k,j);
904 this->
operator()(k,j) += scaling * ADij * B(i,k);
907 this->
operator()(k,j) += scaling * ADij * B(k,i);
912 template <
typename number>
913 template <
typename number2>
917 Assert (!this->empty(), ExcEmptyMatrix());
924 const number *val_ptr = &this->values[0];
925 const number2 *v_ptr;
930 const number *
const val_end_of_row = val_ptr+n_rows;
932 while (val_ptr != val_end_of_row)
933 s += number(*val_ptr++) * number(*v_ptr++);
942 template <
typename number>
943 template <
typename number2>
948 Assert (!this->empty(), ExcEmptyMatrix());
956 const number *val_ptr = &this->values[0];
957 const number2 *v_ptr;
962 const number *
const val_end_of_row = val_ptr+n_cols;
964 while (val_ptr != val_end_of_row)
965 s += number(*val_ptr++) * number(*v_ptr++);
967 sum += s * number(u(row));
975 template <
typename number>
979 Assert (m() == n(), ExcNotQuadratic());
985 const number t = ((*this)(i,j) + (*
this)(j,i)) / number(2.);
986 (*this)(i,j) = (*
this)(j,i) = t;
991 template <
typename number>
995 Assert (!this->empty(), ExcEmptyMatrix());
998 const size_type n_rows = m(), n_cols = n();
1000 for (
size_type col=0; col<n_cols; ++col)
1003 for (
size_type row=0; row<n_rows; ++row)
1004 sum += std::abs((*
this)(row,col));
1013 template <
typename number>
1017 Assert (!this->empty(), ExcEmptyMatrix());
1020 const size_type n_rows = m(), n_cols = n();
1022 for (
size_type row=0; row<n_rows; ++row)
1025 for (
size_type col=0; col<n_cols; ++col)
1026 sum += std::abs((*
this)(row,col));
1035 template <
typename number>
1036 template <
typename number2>
1041 Assert (!this->empty(), ExcEmptyMatrix());
1048 (*
this)(i,j) += a * number(A(i,j));
1052 template <
typename number>
1053 template <
typename number2>
1060 Assert (!this->empty(), ExcEmptyMatrix());
1069 (*
this)(i,j) += a * number(A(i,j)) + b * number(B(i,j));
1074 template <
typename number>
1075 template <
typename number2>
1084 Assert (!this->empty(), ExcEmptyMatrix());
1096 (*
this)(i,j) += a * number(A(i,j)) +
1097 b * number(B(i,j)) +
1103 template <
typename number>
1104 template <
typename number2>
1106 const number factor,
1112 Assert (dst_offset_i < m(),
1114 Assert (dst_offset_j < n(),
1116 Assert (src_offset_i < src.
m(),
1118 Assert (src_offset_j < src.
n(),
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);
1127 (*
this)(dst_offset_i+i,dst_offset_j+j)
1128 += factor * number(src(src_offset_i+i,src_offset_j+j));
1133 template <
typename number>
1134 template <
typename number2>
1136 const number factor,
1142 Assert (dst_offset_i < m(),
1144 Assert (dst_offset_j < n(),
1146 Assert (src_offset_i < src.
n(),
1148 Assert (src_offset_j < src.
m(),
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);
1159 (*
this)(dst_offset_i+i,dst_offset_j+j)
1160 += factor * number(src(src_offset_i+j,src_offset_j+i));
1165 template <
typename number>
1166 template <
typename number2>
1171 Assert (!this->empty(), ExcEmptyMatrix());
1173 Assert (m() == n(), ExcNotQuadratic());
1179 (*
this)(i,j) += a * number(A(j,i));
1183 template <
typename number>
1192 template <
typename number>
1196 Assert (!this->empty(), ExcEmptyMatrix());
1198 Assert (this->n_cols() == this->n_rows(),
1201 switch (this->n_cols())
1204 return (*
this)(0,0);
1206 return (*
this)(0,0)*(*
this)(1,1) - (*
this)(1,0)*(*
this)(0,1);
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));
1222 template <
typename number>
1226 Assert (!this->empty(), ExcEmptyMatrix());
1228 Assert (this->n_cols() == this->n_rows(),
1232 for (
size_type i=0; i<this->n_rows(); ++i)
1240 template <
typename number>
1244 Assert (!this->empty(), ExcEmptyMatrix());
1247 for (
size_type i=0; i<this->n_rows()*this->n_cols(); ++i)
1249 return std::sqrt(s);
1254 template <
typename number>
1258 Assert (!this->empty(), ExcEmptyMatrix());
1262 for (
size_type i=0; i<this->n_rows(); ++i)
1263 for (
size_type j=0; j<this->n_cols(); ++j)
1265 const number x_ij = (*this)(i,j);
1266 const number x_ji = (*this)(j,i);
1273 return std::sqrt(a)/std::sqrt(s);
1279 template <
typename number>
1280 template <
typename number2>
1284 Assert (!this->empty(), ExcEmptyMatrix());
1286 Assert (this->n_cols() == this->n_rows(),
1288 Assert (this->n_cols() == M.n_cols(),
1290 Assert (this->n_rows() == M.n_rows(),
1301 switch (this->n_cols())
1304 (*this)(0,0) = number2(1.0)/M(0,0);
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;
1320 const number2 t4 = M(0,0)*M(1,1),
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;
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;
1453 template <
typename number>
1454 template <
typename number2>
1474 this->reinit(A.
m(), A.
n());
1476 double SLik2 = 0.0, SLikLjk = 0.0;
1477 for (
size_type i=0; i< this->n_cols(); i++)
1485 SLikLjk += (*this)(i,k)*(*
this)(j,k);
1487 (*this)(i,j) = (1./(*
this)(j,j))*(A(i,j) - SLikLjk);
1488 SLik2 += (*this)(i,j)*(*
this)(i,j);
1491 ExcMatrixNotPositiveDefinite());
1493 (*this)(i,i) = std::sqrt(A(i,i) - SLik2);
1499 template <
typename number>
1500 template <
typename number2>
1510 for (
size_type j = 0; j< this->n(); j++)
1512 (*this)(i,j) = V(i)*W(j);
1518 template <
typename number>
1519 template <
typename number2>
1534 A_t.mmult(A_t_times_A,A);
1535 if (number(A_t_times_A.determinant())==number(0))
1536 Assert(
false, ExcSingular())
1539 A_t_times_A_inv.invert(A_t_times_A);
1540 A_t_times_A_inv.mmult(left_inv,A_t);
1546 template <
typename number>
1547 template <
typename number2>
1562 A.
mmult(A_times_A_t,A_t);
1563 if (number(A_times_A_t.determinant())==number(0))
1564 Assert(
false, ExcSingular())
1567 A_times_A_t_inv.invert(A_times_A_t);
1568 A_t.mmult(right_inv,A_times_A_t_inv);
1575 template <
typename number>
1587 Assert (!this->empty(), ExcEmptyMatrix());
1588 Assert(this->m()-dst_r>src_r_j-src_r_i,
1590 Assert(this->n()-dst_c>src_c_j-src_c_i,
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]);
1604 template <
typename number>
1615 Assert (!this->empty(), ExcEmptyMatrix());
1616 Assert(dim-dst_r>src_r_j-src_r_i,
1618 Assert(dim-dst_c>src_c_j-src_c_i,
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));
1633 template <
typename number>
1634 template <
typename somenumber>
1638 const number om)
const
1640 Assert (m() == n(), ExcNotQuadratic());
1644 const size_t n = src.
size();
1645 somenumber *dst_ptr = dst.
begin();
1646 const somenumber *src_ptr = src.
begin();
1648 for (
size_type i=0; i<n; ++i, ++dst_ptr, ++src_ptr)
1649 *dst_ptr = somenumber(om) **src_ptr / somenumber((*
this)(i,i));
1654 template <
typename number>
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
1665 unsigned int width = width_;
1667 Assert ((!this->empty()) || (this->n_cols()+this->n_rows()==0),
1672 std::ios::fmtflags old_flags = out.flags();
1673 unsigned int old_precision = out.precision (precision);
1677 out.setf (std::ios::scientific, std::ios::floatfield);
1679 width = precision+7;
1683 out.setf (std::ios::fixed, std::ios::floatfield);
1685 width = precision+2;
1691 if (std::abs((*
this)(i,j)) > threshold)
1692 out << std::setw(width)
1693 << (*this)(i,j) * number(denominator) <<
' ';
1695 out << std::setw(width) << zero_string <<
' ';
1701 out.flags (old_flags);
1702 out.precision(old_precision);
1706 template <
typename number>
1710 Assert (!this->empty(), ExcEmptyMatrix());
1711 Assert (this->n_cols() == this->n_rows(), ExcNotQuadratic());
1718 #if defined(HAVE_DGETRF_) && defined (HAVE_SGETRF_) && \
1719 defined(HAVE_DGETRI_) && defined (HAVE_SGETRI_)
1723 if (this->n_cols() > 15)
1744 const int nn = this->n();
1747 std::vector<int> ipiv(nn);
1752 getrf(&nn, &nn, &this->values[0], &nn, &ipiv[0], &info);
1755 Assert(info == 0, LACExceptions::ExcSingular());
1758 std::vector<number> inv_work (nn);
1763 getri(&nn, &this->values[0], &nn, &ipiv[0], &inv_work[0], &nn, &info);
1766 Assert(info == 0, LACExceptions::ExcSingular());
1786 double diagonal_sum = 0;
1788 diagonal_sum += std::abs((*
this)(i,i));
1789 const double typical_diagonal_element = diagonal_sum/N;
1794 std::vector<size_type> p(N);
1808 if (std::abs((*
this)(i,j)) > max)
1810 max = std::abs((*
this)(i,j));
1816 Assert(max > 1.e-16*typical_diagonal_element,
1817 ExcNotRegular(max));
1823 std::swap ((*
this)(j,k), (*
this)(r,k));
1825 std::swap (p[j], p[r]);
1829 const number hr = number(1.)/(*this)(j,j);
1837 (*this)(i,k) -= (*
this)(i,j)*(*
this)(j,k)*hr;
1843 (*this)(j,i) *= -hr;
1848 std::vector<number> hv(N);
1852 hv[p[k]] = (*
this)(i,k);
1854 (*
this)(i,k) = hv[k];
1860 template <
typename number>
1870 DEAL_II_NAMESPACE_CLOSE
unsigned int n_rows() const
real_type l1_norm() const
void diagadd(const number s)
FullMatrix< number > & operator=(const FullMatrix< number > &)
#define AssertDimension(dim1, dim2)
void backward(Vector< number2 > &dst, const Vector< number2 > &src) const
FullMatrix & operator*=(const number factor)
bool operator==(const FullMatrix< number > &) const
void cholesky(const FullMatrix< number2 > &A)
void add_row(const size_type i, const number s, const size_type j)
real_type linfty_norm() const
::ExceptionBase & ExcMessage(std::string arg1)
number residual(Vector< number2 > &dst, const Vector< number2 > &x, const Vector< number3 > &b) const
void outer_product(const Vector< number2 > &V, const Vector< number2 > &W)
void right_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)
void invert(const FullMatrix< number2 > &M)
TableBase< N, T > & operator=(const TableBase< N, T > &src)
void left_invert(const FullMatrix< number2 > &M)
void Tadd(const number s, const FullMatrix< number2 > &B)
#define AssertThrow(cond, exc)
void vmult(Vector< number2 > &w, const Vector< number2 > &v, const bool adding=false) const
bool is_finite(const double x)
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 forward(Vector< number2 > &dst, const Vector< number2 > &src) const
real_type relative_symmetry_norm2() const
number2 matrix_scalar_product(const Vector< number2 > &u, const Vector< number2 > &v) const
::ExceptionBase & ExcIO()
std::size_t memory_consumption() const
numbers::NumberTraits< number >::real_type real_type
void swap_row(const size_type i, const size_type j)
number2 matrix_norm_square(const Vector< number2 > &v) const
static real_type abs_square(const number &x)
#define Assert(cond, exc)
real_type frobenius_norm() const
static bool equal(const T *p1, const T *p2)
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
void Tvmult(Vector< number2 > &w, const Vector< number2 > &v, 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
void add_col(const size_type i, const number s, const size_type j)
void mTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
FullMatrix & operator/=(const number factor)
void swap_col(const size_type i, const size_type j)
::ExceptionBase & ExcNumberNotFinite()
void add(const number a, const FullMatrix< number2 > &A)
void TmTmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void Tmmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
FullMatrix(const size_type n=0)
SymmetricTensor< 4, dim, Number > invert(const SymmetricTensor< 4, dim, Number > &t)
number determinant() const
void mmult(FullMatrix< number2 > &C, const FullMatrix< number2 > &B, const bool adding=false) const
void copy_from(const MATRIX &)
::ExceptionBase & ExcInternalError()
void equ(const number a, const FullMatrix< number2 > &A)
unsigned int n_cols() const
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)
VectorizedArray< Number > max(const ::VectorizedArray< Number > &x, const ::VectorizedArray< Number > &y)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
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
bool operator==(const TableBase< N, T > &T2) const