36#include <Eigen/Eigenvalues>
62 int* ipiv =
new int[n];
74 for (
Index i = 0; i < n; i++) {
115 for (
Index i = 0; i < n; i++) {
116 ipiv[i] = (int)indx[i];
121 &trans, &n_int, &one, LU.
mdata, &n_int, ipiv, rhs, &n_int, &info);
123 for (
Index i = 0; i < n; i++) {
175 int* ipiv =
new int[n];
187 work =
new double[lwork];
188 }
catch (std::bad_alloc& ba) {
190 "Error inverting matrix: Could not allocate workspace memory.");
198 ARTS_USER_ERROR_IF (info not_eq 0,
"Error inverting matrix: Matrix not of full rank.");
210 int n_int = int(n), lwork = n_int, info;
211 std::vector<int> ipiv(n);
220 "Error inverting matrix: Matrix not of full rank.");
260 int LDA, LDA_L, LDA_R, n_int, info;
265 char l_eig =
'N', r_eig =
'V';
268 int lwork = 2 * n_int + n_int * n_int;
269 double *work, *rwork;
271 rwork =
new double[2 * n_int];
272 work =
new double[lwork];
273 }
catch (std::bad_alloc& ba) {
275 "Error diagonalizing: Could not allocate workspace memory.");
279 double* adata = A_tmp.
mdata;
280 double* rpdata = P2.
mdata;
281 double* lpdata =
new double[0];
282 double* wrdata = WR2.
mdata;
283 double* widata = WI2.
mdata;
308 for (
Index i = 0; i < n; i++)
309 for (
Index j = 0; j < n; j++) P(j, i) = P2(i, j);
343 char l_eig =
'N', r_eig =
'V';
346 int lwork = 2 * n_int + n_int * n_int;
367 for (
Index i = 0; i < n; i++)
368 for (
Index j = 0; j <= i; j++)
396 Matrix D(n, n),
N(n, n), X(n, n), cX(n, n, 0.0), B(n, n, 0.0);
397 Vector N_col_vec(n, 0.), F_col_vec(n, 0.);
402 j = 1 + floor(1. / log(2.) * log(A_norm_inf));
420 for (
Index k = 0; k <
q; k++) {
422 c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
440 for (
Index i = 0; i < n; i++) {
442 lubacksub(F_col_vec, X, N_col_vec, indx);
443 F(
joker, i) = F_col_vec;
447 for (
Index k = 0; k < j_index; k++) {
460 Matrix P(n, n), invP(n, n);
466 for (
Index k = 0; k < n; k++) {
468 "We require real eigenvalues for your chosen method.");
469 P(
joker, k) *= exp(WR[k]);
483 const Numeric e = 1.0 + floor(1.0 / log(2.0) * log(A_norm_inf));
484 const Index r = (e + 1.) > 0. ? (
Index)(e + 1.) : 0;
489 Eigen::Matrix4d X =
M;
490 Eigen::Matrix4d cX =
c * X;
498 for (
Index k = 2; k <=
q; k++) {
501 cX.noalias() =
c * X;
512 for (
Index k = 1; k <= r; k++) F *= F;
521 Eigen::EigenSolver<Eigen::Matrix4d> es;
524 const Eigen::Matrix4cd U = es.eigenvectors();
525 const Eigen::Vector4cd
q = es.eigenvalues();
526 const Eigen::Vector4cd eq =
q.array().exp();
527 Eigen::Matrix4cd res = U * eq.asDiagonal();
571 for (
Index ii = 0; ii < n_partials; ii++) {
581 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
582 Index r = (e + 1.) > 0. ? (
Index)(e + 1.) : 0;
587 Matrix M = A, X(n, n), cX(n, n), D(n, n);
598 Tensor3 dM_upp(n_partials, n, n), dM_low(n_partials, n, n),
599 dD_upp(n_partials, n, n), dD_low(n_partials, n, n),
600 Y_upp(n_partials, n, n), Y_low(n_partials, n, n),
601 cY_upp(n_partials, n, n), cY_low(n_partials, n, n);
602 for (
Index ii = 0; ii < n_partials; ii++) {
603 for (
Index jj = 0; jj < n; jj++) {
605 dM_upp(ii, jj, kk) = dA_upp(ii, jj, kk) * pow2rm1;
606 dM_low(ii, jj, kk) = dA_low(ii, jj, kk) * pow2rm1;
608 Y_upp(ii, jj, kk) = dM_upp(ii, jj, kk);
609 Y_low(ii, jj, kk) = dM_low(ii, jj, kk);
611 cY_upp(ii, jj, kk) =
c * Y_upp(ii, jj, kk);
612 cY_low(ii, jj, kk) =
c * Y_low(ii, jj, kk);
614 dF_upp(ii, jj, kk) = cY_upp(ii, jj, kk);
615 dF_low(ii, jj, kk) = cY_low(ii, jj, kk);
617 dD_upp(ii, jj, kk) = -cY_upp(ii, jj, kk);
618 dD_low(ii, jj, kk) = -cY_low(ii, jj, kk);
624 Matrix tmp1(n, n), tmp2(n, n);
626 for (
Index k = 2; k <=
q; k++) {
630 for (
Index ii = 0; ii < n_partials; ii++) {
631 Matrix tmp_low1(n, n), tmp_upp1(n, n), tmp_low2(n, n), tmp_upp2(n, n);
639 for (
Index jj = 0; jj < n; jj++) {
640 for (
Index kk = 0; kk < n; kk++) {
641 Y_upp(ii, jj, kk) = tmp_upp1(jj, kk) + tmp_upp2(jj, kk);
642 Y_low(ii, jj, kk) = tmp_low1(jj, kk) + tmp_low2(jj, kk);
645 cY_upp(ii, jj, kk) =
c * Y_upp(ii, jj, kk);
646 cY_low(ii, jj, kk) =
c * Y_low(ii, jj, kk);
649 dF_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
650 dF_low(ii, jj, kk) += cY_low(ii, jj, kk);
655 dD_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
656 dD_low(ii, jj, kk) += cY_low(ii, jj, kk);
659 dD_upp(ii, jj, kk) -= cY_upp(ii, jj, kk);
660 dD_low(ii, jj, kk) -= cY_low(ii, jj, kk);
669 for (
Index jj = 0; jj < n; jj++) {
670 for (
Index kk = 0; kk < n; kk++) {
671 X(jj, kk) = tmp1(jj, kk);
672 cX(jj, kk) = tmp1(jj, kk) *
c;
673 F(jj, kk) += cX(jj, kk);
676 D(jj, kk) += cX(jj, kk);
678 D(jj, kk) -= cX(jj, kk);
691 for (
Index ii = 0; ii < n_partials; ii++) {
692 Matrix tmp_low(n, n), tmp_upp(n, n);
697 for (
Index jj = 0; jj < n; jj++) {
698 for (
Index kk = 0; kk < n; kk++) {
699 dF_upp(ii, jj, kk) -= tmp_upp(jj, kk);
700 dF_low(ii, jj, kk) -= tmp_low(jj, kk);
707 for (
Index jj = 0; jj < n; jj++) {
708 for (
Index kk = 0; kk < n; kk++) {
709 dF_upp(ii, jj, kk) = tmp_upp(jj, kk);
710 dF_low(ii, jj, kk) = tmp_low(jj, kk);
715 for (
Index k = 1; k <= r; k++) {
717 for (
Index ii = 0; ii < n_partials; ii++) {
718 Matrix tmp_low1(n, n), tmp_upp1(n, n), tmp_low2(n, n), tmp_upp2(n, n);
726 for (
Index jj = 0; jj < n; jj++) {
727 for (
Index kk = 0; kk < n; kk++) {
728 dF_upp(ii, jj, kk) = tmp_upp1(jj, kk) + tmp_upp2(jj, kk);
729 dF_low(ii, jj, kk) = tmp_low1(jj, kk) + tmp_low2(jj, kk);
788 const Numeric a = A(0, 0),
b = A(0, 1),
c = A(0, 2),
d = A(0, 3),
u = A(1, 2),
789 v = A(1, 3),
w = A(2, 3);
794 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
797 Const1 =
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
798 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
799 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
800 Const1 += u2 * (u2 * 0.5 + v2 + w2);
801 Const1 += v2 * (v2 * 0.5 + w2);
803 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
807 Const1 =
sqrt(Const1);
814 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
820 const Numeric x = sqrt_BpA.real() * sqrt_05;
821 const Numeric y = sqrt_BmA.imag() * sqrt_05;
826 const Numeric cosh_x = cosh(x);
827 const Numeric sinh_x = sinh(x);
829 const Numeric inv_x2y2 = 1.0 / x2y2;
832 Numeric inv_y = 0.0, inv_x = 0.0;
839 C2 = (1.0 - cos_y) * inv_x2y2;
840 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
841 }
else if (y == 0.0) {
845 C2 = (cosh_x - 1.0) * inv_x2y2;
846 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
851 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
852 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
853 C2 = (cosh_x - cos_y) * inv_x2y2;
854 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
857 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
859 F(0, 0) += C2 * (
b2 + c2 + d2);
861 F(1, 1) += C2 * (
b2 - u2 - v2);
863 F(2, 2) += C2 * (c2 - u2 - w2);
865 F(3, 3) += C2 * (d2 - v2 - w2);
867 F(0, 1) = F(1, 0) = C1 *
b;
870 C2 * (-
c *
u -
d *
v) +
871 C3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
v * (
b *
v +
c *
w));
874 C2 * (
c *
u +
d *
v) +
875 C3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
d * (
b *
d +
u *
w));
877 F(0, 2) = F(2, 0) = C1 *
c;
880 C2 * (
b *
u -
d *
w) +
881 C3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
w * (
b *
v +
c *
w));
884 C2 * (-
b *
u +
d *
w) +
885 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
d * (
c *
d -
u *
v));
887 F(0, 3) = F(3, 0) = C1 *
d;
890 C2 * (
b *
v +
c *
w) +
891 C3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
w * (
b *
u -
d *
w));
894 C2 * (-
b *
v -
c *
w) +
895 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
d * (-d2 + v2 + w2));
897 F(1, 2) = F(2, 1) = C2 * (
b *
c -
v *
w);
899 F(1, 2) += C1 *
u + C3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
900 w * (
b *
d +
u *
w));
902 F(2, 1) += -C1 *
u + C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
903 v * (
c *
d -
u *
v));
905 F(1, 3) = F(3, 1) = C2 * (
b *
d +
u *
w);
907 F(1, 3) += C1 *
v + C3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
908 w * (
b *
c -
v *
w));
910 F(3, 1) += -C1 *
v + C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
911 v * (-d2 + v2 + w2));
913 F(2, 3) = F(3, 2) = C2 * (
c *
d -
u *
v);
915 F(2, 3) += C1 *
w + C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
916 w * (-c2 + u2 + w2));
918 F(3, 2) += -C1 *
w + C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
919 w * (-d2 + v2 + w2));
966 eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
968 Eigen::Matrix4d eigA2 = eigA;
970 Eigen::Matrix4d eigA3 = eigA2;
975 const Numeric a = A(0, 0),
b = A(0, 1),
c = A(0, 2),
d = A(0, 3),
u = A(1, 2),
976 v = A(1, 3),
w = A(2, 3);
981 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
984 Const1 =
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
985 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
986 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
987 Const1 += u2 * (u2 * 0.5 + v2 + w2);
988 Const1 += v2 * (v2 * 0.5 + w2);
990 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
994 Const1 =
sqrt(Const1);
1001 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1007 const Numeric x = sqrt_BpA.real() * sqrt_05;
1008 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1013 const Numeric cosh_x = cosh(x);
1014 const Numeric sinh_x = sinh(x);
1016 const Numeric inv_x2y2 = 1.0 / x2y2;
1019 Numeric inv_y = 0.0, inv_x = 0.0;
1026 C2 = (1.0 - cos_y) * inv_x2y2;
1027 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1028 }
else if (y == 0.0) {
1032 C2 = (cosh_x - 1.0) * inv_x2y2;
1033 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1038 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1039 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1040 C2 = (cosh_x - cos_y) * inv_x2y2;
1041 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1044 eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1116 const Numeric a = A(0, 0),
b = A(0, 1),
c = A(0, 2),
d = A(0, 3),
u = A(1, 2),
1117 v = A(1, 3),
w = A(2, 3);
1122 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
1125 Const1 =
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1126 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1127 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
1128 Const1 += u2 * (u2 * 0.5 + v2 + w2);
1129 Const1 += v2 * (v2 * 0.5 + w2);
1135 Const1 =
sqrt(Const1);
1142 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1143 for (
Index i = 0; i < npd; i++) {
1153 const Numeric x = sqrt_BpA.real() * sqrt_05;
1154 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1159 const Numeric cosh_x = cosh(x);
1160 const Numeric sinh_x = sinh(x);
1162 const Numeric inv_x2y2 = 1.0 / x2y2;
1165 Numeric inv_y = 0.0, inv_x = 0.0;
1172 C2 = (1.0 - cos_y) * inv_x2y2;
1173 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1174 }
else if (y == 0.0) {
1178 C2 = (cosh_x - 1.0) * inv_x2y2;
1179 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1184 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1185 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1186 C2 = (cosh_x - cos_y) * inv_x2y2;
1187 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1190 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
1192 F(0, 0) += C2 * (
b2 + c2 + d2);
1194 F(1, 1) += C2 * (
b2 - u2 - v2);
1196 F(2, 2) += C2 * (c2 - u2 - w2);
1198 F(3, 3) += C2 * (d2 - v2 - w2);
1200 F(0, 1) = F(1, 0) = C1 *
b;
1203 C2 * (-
c *
u -
d *
v) +
1204 C3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
v * (
b *
v +
c *
w));
1207 C2 * (
c *
u +
d *
v) +
1208 C3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
d * (
b *
d +
u *
w));
1210 F(0, 2) = F(2, 0) = C1 *
c;
1213 C2 * (
b *
u -
d *
w) +
1214 C3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
w * (
b *
v +
c *
w));
1217 C2 * (-
b *
u +
d *
w) +
1218 C3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
d * (
c *
d -
u *
v));
1220 F(0, 3) = F(3, 0) = C1 *
d;
1223 C2 * (
b *
v +
c *
w) +
1224 C3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
w * (
b *
u -
d *
w));
1227 C2 * (-
b *
v -
c *
w) +
1228 C3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
d * (-d2 + v2 + w2));
1230 F(1, 2) = F(2, 1) = C2 * (
b *
c -
v *
w);
1232 F(1, 2) += C1 *
u + C3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
1233 w * (
b *
d +
u *
w));
1235 F(2, 1) += -C1 *
u + C3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1236 v * (
c *
d -
u *
v));
1238 F(1, 3) = F(3, 1) = C2 * (
b *
d +
u *
w);
1240 F(1, 3) += C1 *
v + C3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
1241 w * (
b *
c -
v *
w));
1243 F(3, 1) += -C1 *
v + C3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1244 v * (-d2 + v2 + w2));
1246 F(2, 3) = F(3, 2) = C2 * (
c *
d -
u *
v);
1248 F(2, 3) += C1 *
w + C3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1249 w * (-c2 + u2 + w2));
1251 F(3, 2) += -C1 *
w + C3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1252 w * (-d2 + v2 + w2));
1257 const Numeric inv_x2 = inv_x * inv_x;
1258 const Numeric inv_y2 = inv_y * inv_y;
1260 for (
Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1261 for (
Index i = 0; i < npd; i++) {
1267 const Numeric da = dA(0, 0), db = dA(0, 1), dc = dA(0, 2),
1268 dd = dA(0, 3), du = dA(1, 2), dv = dA(1, 3),
1271 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1272 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1274 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1278 dConst1 = db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1279 dConst1 +=
b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1281 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1282 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1284 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1285 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1287 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1288 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1290 dConst1 += dv2 * (v2 * 0.5 + w2);
1291 dConst1 += v2 * (dv2 * 0.5 + dw2);
1293 dConst1 += 4 * ((db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1294 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1295 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1296 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1297 dConst1 += dw2 * w2;
1305 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1309 dC2 = -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1310 dC3 = -2 * y * dy * C3 * inv_x2y2 +
1311 (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1313 }
else if (y == 0.0) {
1315 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1319 dC2 = -2 * x *
dx * C2 * inv_x2y2 +
dx * sinh_x * inv_x2y2;
1320 dC3 = -2 * x *
dx * C3 * inv_x2y2 +
1321 (cosh_x *
dx * inv_x -
dx * sinh_x * inv_x2) * inv_x2y2;
1324 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1326 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1327 const Numeric dy2 = 2 * y * dy;
1329 const Numeric dx2dy2 = dx2 + dy2;
1331 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1332 (2 * cos_y *
dx * x + 2 * cosh_x * dy * y +
dx * sinh_x * y2 -
1336 dC1 = -dx2dy2 * C1 * inv_x2y2 +
1337 (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1338 dy * sin_y * x2 * inv_y2 -
dx * sinh_x * y2 * inv_x2 +
1339 cosh_x *
dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1342 dC2 = -dx2dy2 * C2 * inv_x2y2 + (
dx * sinh_x + dy * sin_y) * inv_x2y2;
1344 dC3 = -dx2dy2 * C3 * inv_x2y2 +
1345 (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1346 cosh_x *
dx * inv_x -
dx * sinh_x * inv_x2) *
1350 dF(0, 0) = dF(1, 1) = dF(2, 2) = dF(3, 3) = dC0;
1352 dF(0, 0) += dC2 * (
b2 + c2 + d2) + C2 * (db2 + dc2 + dd2);
1354 dF(1, 1) += dC2 * (
b2 - u2 - v2) + C2 * (db2 - du2 - dv2);
1356 dF(2, 2) += dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2);
1358 dF(3, 3) += dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2);
1360 dF(0, 1) = dF(1, 0) = db * C1 +
b * dC1;
1362 dF(0, 1) += dC2 * (-
c *
u -
d *
v) +
1363 C2 * (-dc *
u - dd *
v -
c * du -
d * dv) +
1364 dC3 * (
b * (
b2 + c2 + d2) -
u * (
b *
u -
d *
w) -
1365 v * (
b *
v +
c *
w)) +
1366 C3 * (db * (
b2 + c2 + d2) - du * (
b *
u -
d *
w) -
1367 dv * (
b *
v +
c *
w) +
b * (db2 + dc2 + dd2) -
1368 u * (db *
u - dd *
w) -
v * (db *
v + dc *
w) -
1369 u * (
b * du -
d * dw) -
v * (
b * dv +
c * dw));
1370 dF(1, 0) += dC2 * (
c *
u +
d *
v) +
1371 C2 * (dc *
u + dd *
v +
c * du +
d * dv) +
1372 dC3 * (-
b * (-
b2 + u2 + v2) +
c * (
b *
c -
v *
w) +
1373 d * (
b *
d +
u *
w)) +
1374 C3 * (-db * (-
b2 + u2 + v2) + dc * (
b *
c -
v *
w) +
1375 dd * (
b *
d +
u *
w) -
b * (-db2 + du2 + dv2) +
1376 c * (db *
c - dv *
w) +
d * (db *
d + du *
w) +
1377 c * (
b * dc -
v * dw) +
d * (
b * dd +
u * dw));
1379 dF(0, 2) = dF(2, 0) = dC1 *
c + C1 * dc;
1381 dF(0, 2) += dC2 * (
b *
u -
d *
w) +
1382 C2 * (db *
u - dd *
w +
b * du -
d * dw) +
1383 dC3 * (
c * (
b2 + c2 + d2) -
u * (
c *
u +
d *
v) -
1384 w * (
b *
v +
c *
w)) +
1385 C3 * (dc * (
b2 + c2 + d2) - du * (
c *
u +
d *
v) -
1386 dw * (
b *
v +
c *
w) +
c * (db2 + dc2 + dd2) -
1387 u * (dc *
u + dd *
v) -
w * (db *
v + dc *
w) -
1388 u * (
c * du +
d * dv) -
w * (
b * dv +
c * dw));
1389 dF(2, 0) += dC2 * (-
b *
u +
d *
w) +
1390 C2 * (-db *
u + dd *
w -
b * du +
d * dw) +
1391 dC3 * (
b * (
b *
c -
v *
w) -
c * (-c2 + u2 + w2) +
1392 d * (
c *
d -
u *
v)) +
1393 C3 * (db * (
b *
c -
v *
w) - dc * (-c2 + u2 + w2) +
1394 dd * (
c *
d -
u *
v) +
b * (db *
c - dv *
w) -
1395 c * (-dc2 + du2 + dw2) +
d * (dc *
d - du *
v) +
1396 b * (
b * dc -
v * dw) +
d * (
c * dd -
u * dv));
1398 dF(0, 3) = dF(3, 0) = dC1 *
d + C1 * dd;
1400 dF(0, 3) += dC2 * (
b *
v +
c *
w) +
1401 C2 * (db *
v + dc *
w +
b * dv +
c * dw) +
1402 dC3 * (
d * (
b2 + c2 + d2) -
v * (
c *
u +
d *
v) +
1403 w * (
b *
u -
d *
w)) +
1404 C3 * (dd * (
b2 + c2 + d2) - dv * (
c *
u +
d *
v) +
1405 dw * (
b *
u -
d *
w) +
d * (db2 + dc2 + dd2) -
1406 v * (dc *
u + dd *
v) +
w * (db *
u - dd *
w) -
1407 v * (
c * du +
d * dv) +
w * (
b * du -
d * dw));
1408 dF(3, 0) += dC2 * (-
b *
v -
c *
w) +
1409 C2 * (-db *
v - dc *
w -
b * dv -
c * dw) +
1410 dC3 * (
b * (
b *
d +
u *
w) +
c * (
c *
d -
u *
v) -
1411 d * (-d2 + v2 + w2)) +
1412 C3 * (db * (
b *
d +
u *
w) + dc * (
c *
d -
u *
v) -
1413 dd * (-d2 + v2 + w2) +
b * (db *
d + du *
w) +
1414 c * (dc *
d - du *
v) -
d * (-dd2 + dv2 + dw2) +
1415 b * (
b * dd +
u * dw) +
c * (
c * dd -
u * dv));
1417 dF(1, 2) = dF(2, 1) =
1418 dC2 * (
b *
c -
v *
w) + C2 * (db *
c +
b * dc - dv *
w -
v * dw);
1420 dF(1, 2) += dC1 *
u + C1 * du +
1421 dC3 * (
c * (
c *
u +
d *
v) -
u * (-
b2 + u2 + v2) -
1422 w * (
b *
d +
u *
w)) +
1423 C3 * (dc * (
c *
u +
d *
v) - du * (-
b2 + u2 + v2) -
1424 dw * (
b *
d +
u *
w) +
c * (dc *
u + dd *
v) -
1425 u * (-db2 + du2 + dv2) -
w * (db *
d + du *
w) +
1426 c * (
c * du +
d * dv) -
w * (
b * dd +
u * dw));
1427 dF(2, 1) += -dC1 *
u - C1 * du +
1428 dC3 * (-
b * (
b *
u -
d *
w) +
u * (-c2 + u2 + w2) -
1429 v * (
c *
d -
u *
v)) +
1430 C3 * (-db * (
b *
u -
d *
w) + du * (-c2 + u2 + w2) -
1431 dv * (
c *
d -
u *
v) -
b * (db *
u - dd *
w) +
1432 u * (-dc2 + du2 + dw2) -
v * (dc *
d - du *
v) -
1433 b * (
b * du -
d * dw) -
v * (
c * dd -
u * dv));
1435 dF(1, 3) = dF(3, 1) =
1436 dC2 * (
b *
d +
u *
w) + C2 * (db *
d +
b * dd + du *
w +
u * dw);
1438 dF(1, 3) += dC1 *
v + C1 * dv +
1439 dC3 * (
d * (
c *
u +
d *
v) -
v * (-
b2 + u2 + v2) +
1440 w * (
b *
c -
v *
w)) +
1441 C3 * (dd * (
c *
u +
d *
v) - dv * (-
b2 + u2 + v2) +
1442 dw * (
b *
c -
v *
w) +
d * (dc *
u + dd *
v) -
1443 v * (-db2 + du2 + dv2) +
w * (db *
c - dv *
w) +
1444 d * (
c * du +
d * dv) +
w * (
b * dc -
v * dw));
1445 dF(3, 1) += -dC1 *
v - C1 * dv +
1446 dC3 * (-
b * (
b *
v +
c *
w) -
u * (
c *
d -
u *
v) +
1447 v * (-d2 + v2 + w2)) +
1448 C3 * (-db * (
b *
v +
c *
w) - du * (
c *
d -
u *
v) +
1449 dv * (-d2 + v2 + w2) -
b * (db *
v + dc *
w) -
1450 u * (dc *
d - du *
v) +
v * (-dd2 + dv2 + dw2) -
1451 b * (
b * dv +
c * dw) -
u * (
c * dd -
u * dv));
1453 dF(2, 3) = dF(3, 2) =
1454 dC2 * (
c *
d -
u *
v) + C2 * (dc *
d +
c * dd - du *
v -
u * dv);
1456 dF(2, 3) += dC1 *
w + C1 * dw +
1457 dC3 * (-
d * (
b *
u -
d *
w) +
v * (
b *
c -
v *
w) -
1458 w * (-c2 + u2 + w2)) +
1459 C3 * (-dd * (
b *
u -
d *
w) + dv * (
b *
c -
v *
w) -
1460 dw * (-c2 + u2 + w2) -
d * (db *
u - dd *
w) +
1461 v * (db *
c - dv *
w) -
w * (-dc2 + du2 + dw2) -
1462 d * (
b * du -
d * dw) +
v * (
b * dc -
v * dw));
1463 dF(3, 2) += -dC1 *
w - C1 * dw +
1464 dC3 * (-
c * (
b *
v +
c *
w) +
u * (
b *
d +
u *
w) +
1465 w * (-d2 + v2 + w2)) +
1466 C3 * (-dc * (
b *
v +
c *
w) + du * (
b *
d +
u *
w) +
1467 dw * (-d2 + v2 + w2) -
c * (db *
v + dc *
w) +
1468 u * (db *
d + du *
w) +
w * (-dd2 + dv2 + dw2) -
1469 c * (
b * dv +
c * dw) +
u * (
b * dd +
u * dw));
1474 dF(0, 0) += F(0, 0) * da;
1475 dF(0, 1) += F(0, 1) * da;
1476 dF(0, 2) += F(0, 2) * da;
1477 dF(0, 3) += F(0, 3) * da;
1478 dF(1, 0) += F(1, 0) * da;
1479 dF(1, 1) += F(1, 1) * da;
1480 dF(1, 2) += F(1, 2) * da;
1481 dF(1, 3) += F(1, 3) * da;
1482 dF(2, 0) += F(2, 0) * da;
1483 dF(2, 1) += F(2, 1) * da;
1484 dF(2, 2) += F(2, 2) * da;
1485 dF(2, 3) += F(2, 3) * da;
1486 dF(3, 0) += F(3, 0) * da;
1487 dF(3, 1) += F(3, 1) * da;
1488 dF(3, 2) += F(3, 2) * da;
1489 dF(3, 3) += F(3, 3) * da;
1558 eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
1560 Eigen::Matrix4d eigA2 = eigA;
1562 Eigen::Matrix4d eigA3 = eigA2;
1568 const Numeric a = A(0, 0),
b = A(0, 1),
c = A(0, 2),
d = A(0, 3),
u = A(1, 2),
1569 v = A(1, 3),
w = A(2, 3);
1574 const Numeric Const2 =
b2 + c2 + d2 - u2 - v2 - w2;
1577 Const1 =
b2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1578 Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1579 Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
1580 Const1 += u2 * (u2 * 0.5 + v2 + w2);
1581 Const1 += v2 * (v2 * 0.5 + w2);
1583 Const1 += 8 * (
b *
d *
u *
w -
b *
c *
v *
w -
c *
d *
u *
v);
1587 Const1 =
sqrt(Const1);
1594 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1595 for (
Index i = 0; i < npd; i++) {
1605 const Numeric x = sqrt_BpA.real() * sqrt_05;
1606 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1611 const Numeric cosh_x = cosh(x);
1612 const Numeric sinh_x = sinh(x);
1614 const Numeric inv_x2y2 = 1.0 / x2y2;
1617 Numeric inv_y = 0.0, inv_x = 0.0;
1624 C2 = (1.0 - cos_y) * inv_x2y2;
1625 C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1626 }
else if (y == 0.0) {
1630 C2 = (cosh_x - 1.0) * inv_x2y2;
1631 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1636 C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1637 C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1638 C2 = (cosh_x - cos_y) * inv_x2y2;
1639 C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1642 eigF(0, 0) = eigF(1, 1) = eigF(2, 2) = eigF(3, 3) = C0;
1643 eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1651 const Numeric inv_x2 = inv_x * inv_x;
1652 const Numeric inv_y2 = inv_y * inv_y;
1654 for (
Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1655 for (
Index i = 0; i < npd; i++) {
1662 const Numeric da = dA(0, 0), db = dA(0, 1), dc = dA(0, 2),
1663 dd = dA(0, 3), du = dA(1, 2), dv = dA(1, 3),
1666 dA(0, 0) = dA(1, 1) = dA(2, 2) = dA(3, 3) = 0.0;
1668 Eigen::Matrix4d dA2;
1670 dA2.noalias() = eigA * dA;
1671 dA2.noalias() += dA * eigA;
1673 Eigen::Matrix4d dA3;
1674 dA3.noalias() = dA2 * eigA;
1675 dA3.noalias() += eigA2 * dA;
1677 const Numeric db2 = 2 * db *
b, dc2 = 2 * dc *
c, dd2 = 2 * dd *
d,
1678 du2 = 2 * du *
u, dv2 = 2 * dv *
v, dw2 = 2 * dw *
w;
1680 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1684 dConst1 = db2 * (
b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1685 dConst1 +=
b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1687 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1688 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1690 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1691 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1693 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1694 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1696 dConst1 += dv2 * (v2 * 0.5 + w2);
1697 dConst1 += v2 * (dv2 * 0.5 + dw2);
1699 dConst1 += 4 * ((db *
d *
u *
w - db *
c *
v *
w - dc *
d *
u *
v +
1700 b * dd *
u *
w -
b * dc *
v *
w -
c * dd *
u *
v +
1701 b *
d * du *
w -
b *
c * dv *
w -
c *
d * du *
v +
1702 b *
d *
u * dw -
b *
c *
v * dw -
c *
d *
u * dv));
1703 dConst1 += dw2 * w2;
1710 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1713 -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1715 -2 * y * dy * C3 * inv_x2y2 +
1716 (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1718 dF.noalias() = dC2 * eigA2;
1719 +C2* dA2 + dC3* eigA3 + C3* dA3;
1721 dF.noalias() += eigF * da;
1722 }
else if (y == 0.0) {
1724 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1727 -2 * x *
dx * C2 * inv_x2y2 +
dx * sinh_x * inv_x2y2;
1729 -2 * x *
dx * C3 * inv_x2y2 +
1730 (cosh_x *
dx * inv_x -
dx * sinh_x * inv_x2) * inv_x2y2;
1732 dF.noalias() = dC2 * eigA2 + C2 * dA2 + dC3 * eigA3 + C3 * dA3;
1734 dF.noalias() += eigF * da;
1737 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1739 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1740 const Numeric dy2 = 2 * y * dy;
1742 const Numeric dx2dy2 = dx2 + dy2;
1744 const Numeric dC0 = -dx2dy2 * C0 * inv_x2y2 +
1745 (2 * cos_y *
dx * x + 2 * cosh_x * dy * y +
1746 dx * sinh_x * y2 - dy * sin_y * x2) *
1750 -dx2dy2 * C1 * inv_x2y2 +
1751 (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1752 dy * sin_y * x2 * inv_y2 -
dx * sinh_x * y2 * inv_x2 +
1753 cosh_x *
dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1757 -dx2dy2 * C2 * inv_x2y2 + (
dx * sinh_x + dy * sin_y) * inv_x2y2;
1759 const Numeric dC3 = -dx2dy2 * C3 * inv_x2y2 +
1760 (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1761 cosh_x *
dx * inv_x -
dx * sinh_x * inv_x2) *
1764 dF.noalias() = dC1 * eigA + C1 * dA + dC2 * eigA2 + C2 * dA2 +
1765 dC3 * eigA3 + C3 * dA3;
1771 dF.noalias() += eigF * da;
1793 for (
Index ii = 0; ii < n_partials; ii++) {
1802 const Numeric e = 1.0 + floor(1.0 / log(2.0) * log(A_norm_inf));
1803 const Index r = (e + 1.) > 0. ? (
Index)(e + 1.) : 0;
1810 Eigen::Matrix4d X =
M;
1811 Eigen::Matrix4d cX =
c * X;
1815 eigF.noalias() += cX;
1824 dFu.reserve(n_partials);
1825 dFl.reserve(n_partials);
1827 for (
Index i = 0; i < n_partials; i++) {
1834 cYu[i].noalias() =
c * dMu[i];
1835 cYl[i].noalias() =
c * dMl[i];
1847 for (
Index k = 2; k <=
q; k++) {
1849 for (
Index i = 0; i < n_partials; i++) {
1850 Yu[i] = dMu[i] * X +
M * Yu[i];
1851 Yl[i] = dMl[i] * X +
M * Yl[i];
1853 cYu[i].noalias() =
c * Yu[i];
1854 cYl[i].noalias() =
c * Yl[i];
1856 dFu[i].noalias() += cYu[i];
1857 dFl[i].noalias() += cYl[i];
1861 cX.noalias() =
c * X;
1862 eigF.noalias() += cX;
1867 for (
Index i = 0; i < n_partials; i++) {
1868 dDu[i].noalias() += cYu[i];
1869 dDl[i].noalias() += cYl[i];
1874 for (
Index i = 0; i < n_partials; i++) {
1875 dDu[i].noalias() -= cYu[i];
1876 dDl[i].noalias() -= cYl[i];
1882 const Eigen::Matrix4d invD = D.inverse();
1885 for (
Index i = 0; i < n_partials; i++) {
1886 dFu[i] = invD * (dFu[i] - dDu[i] * eigF);
1887 dFl[i] = invD * (dFl[i] - dDl[i] * eigF);
1890 for (
Index k = 1; k <= r; k++) {
1891 for (
Index i = 0; i < n_partials; i++) {
1892 dFu[i] = dFu[i] * eigF + eigF * dFu[i];
1893 dFl[i] = dFl[i] * eigF + eigF * dFl[i];
1934 for (
Index ii = 0; ii < n_partials; ii++) {
1942 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
1943 Index r = (e + 1.) > 0. ? (
Index)(e + 1.) : 0;
1948 Matrix M = A, X(n, n), cX(n, n), D(n, n);
1959 Tensor3 dM(n_partials, n, n),
Y(n_partials, n, n), cY(n_partials, n, n),
1960 dD(n_partials, n, n);
1961 for (
Index ii = 0; ii < n_partials; ii++) {
1962 for (
Index jj = 0; jj < n; jj++) {
1963 for (
Index kk = 0; kk < n; kk++) {
1964 dM(ii, jj, kk) = dA(ii, jj, kk) * pow2rm1;
1966 Y(ii, jj, kk) =
dM(ii, jj, kk);
1968 cY(ii, jj, kk) =
c *
Y(ii, jj, kk);
1970 dF(ii, jj, kk) = cY(ii, jj, kk);
1972 dD(ii, jj, kk) = -cY(ii, jj, kk);
1979 Matrix tmp1(n, n), tmp2(n, n);
1981 for (
Index k = 2; k <=
q; k++) {
1985 for (
Index ii = 0; ii < n_partials; ii++) {
1990 for (
Index jj = 0; jj < n; jj++)
1991 for (
Index kk = 0; kk < n; kk++) {
1992 Y(ii, jj, kk) = tmp1(jj, kk) + tmp2(jj, kk);
1995 cY(ii, jj, kk) =
c *
Y(ii, jj, kk);
1998 dF(ii, jj, kk) += cY(ii, jj, kk);
2003 dD(ii, jj, kk) += cY(ii, jj, kk);
2006 dD(ii, jj, kk) -= cY(ii, jj, kk);
2033 mult(tmp2, tmp1, F);
2037 for (
Index ii = 0; ii < n_partials; ii++) {
2045 for (
Index k = 1; k <= r; k++) {
2046 for (
Index ii = 0; ii < n_partials; ii++) {
2098 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
2099 Index r = (e + 1.) > 0. ? (
Index)(e + 1.) : 0;
2119 Matrix D(n, n), dD(n, n);
2137 Matrix tmp1(n, n), tmp2(n, n);
2139 for (
Index k = 2; k <=
q; k++) {
2187 mult(tmp2, tmp1, F);
2193 mult(tmp2, tmp1, dF);
2196 for (
Index k = 1; k <= r; k++) {
2224 for (
Index i = 0; i < A.
ncols(); i++) row_sum +=
abs(A(i, j));
2240 for (
Index i = 0; i < n; i++) I(i, i) = 1.;
2256 return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +
2257 A(0, 2) * A(1, 0) * A(2, 1) - A(0, 2) * A(1, 1) * A(2, 0) -
2258 A(0, 1) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 2) * A(2, 1);
2260 return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
2266 for (
Index j = 0; j < dim; j++) {
2268 for (
Index I = 1; I < dim; I++)
2269 for (
Index J = 0; J < dim; J++) {
2271 temp(I - 1, J) = A(I, J);
2273 temp(I - 1, J - 1) = A(I, J);
2278 ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
2328 Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
2330 for (
Index i = 0; i < n; i++) {
2334 for (
Index i = 0; i < n; i++) {
2342 p[0] = s3 /
Numeric(n) - p[1] * xm;
2347 const Index n = x.nelem();
2368Eigen::ComplexEigenSolver<Eigen::MatrixXcd>
eig(
const Eigen::Ref<Eigen::MatrixXcd> A)
2370 Eigen::ComplexEigenSolver<Eigen::MatrixXcd> ces;
2371 return ces.compute(A);
This file contains the definition of Array.
The global header file for ARTS.
Header file for helper functions for OpenMP.
The ComplexMatrixView class.
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
The ComplexVectorView class.
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
A constant view of a ComplexMatrix.
Index ncols() const noexcept
Index nrows() const noexcept
Index get_column_extent() const
Get the extent of the underlying data.
Index nelem() const noexcept
A constant view of a Matrix.
Index nrows() const noexcept
Index ncols() const noexcept
Numeric * mdata
Pointer to the plain C array that holds the data.
Range mcr
The column range of mdata that is actually used.
A constant view of a Tensor3.
Index npages() const
Returns the number of pages.
A constant view of a Vector.
Numeric * mdata
Pointer to the plain C array that holds the data.
Index nelem() const noexcept
Returns the number of elements.
constexpr Index get_extent() const noexcept
Returns the extent of the range.
constexpr Index get_stride() const noexcept
Returns the stride of the range.
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
void resize(Index n)
Resize function.
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
Interface for the LAPACK library.
void matrix_exp2(MatrixView F, ConstMatrixView A)
void matrix_exp_4x4(MatrixView arts_f, ConstMatrixView A, const Index &q)
void propmat4x4_to_transmat4x4(MatrixView F, Tensor3View dF_upp, Tensor3View dF_low, ConstMatrixView A, ConstTensor3View dA_upp, ConstTensor3View dA_low, const Index &q)
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
void id_mat(MatrixView I)
Identity Matrix.
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(MatrixView F, ConstMatrixView A)
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
void matrix_exp2_4x4(MatrixView arts_f, ConstMatrixView A)
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Numeric det(ConstMatrixView A)
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Matrix4x4ViewMap MapToEigen4x4(Tensor3View &A, const Index &i)
void special_matrix_exp_and_dmatrix_exp_dx_for_rt(MatrixView F, Tensor3View dF_upp, Tensor3View dF_low, ConstMatrixView A, ConstTensor3View dA_upp, ConstTensor3View dA_low, const Index &q)
Special exponential of a Matrix with their derivatives.
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept
Least squares fitting by solving x for known A and y.
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Eigen::ComplexEigenSolver< Eigen::MatrixXcd > eig(const Eigen::Ref< Eigen::MatrixXcd > A)
Return the Eigen decomposition of the eigen matrix.
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(MatrixView F, ConstMatrixView A)
void matrix_exp_dmatrix_exp(MatrixView F, Tensor3View dF, ConstMatrixView A, ConstTensor3View dA, const Index &q)
General exponential of a Matrix with their derivatives.
Linear algebra functions.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Header file for logic.cc.
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
std::complex< Numeric > Complex
constexpr Index dM(Polarization type) noexcept
Gives the change of M given a polarization type.
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
void zgeev_(char *jobvl, char *jobvr, int *n, std::complex< double > *A, int *lda, std::complex< double > *W, std::complex< double > *VL, int *ldvl, std::complex< double > *VR, int *ldvr, std::complex< double > *work, int *lwork, double *rwork, int *info)
void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info)
Matrix inversion.
void dgeev_(char *jobvl, char *jobvr, int *n, double *A, int *lda, double *WR, double *WI, double *VL, int *ldvl, double *VR, int *ldvr, double *work, int *lwork, double *rwork, int *info)
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
void dgetrs_(char *trans, int *n, int *nrhs, double *A, int *lda, int *ipiv, double *b, int *ldb, int *info)
Solve linear system of equations.
void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info)
LU decomposition.
invlib::MatrixIdentity< Matrix > Identity
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Numeric pow(const Rational base, Numeric exp)
Power of.
Numeric sqrt(const Rational r)
Square root.