ARTS  2.4.0(git:4fb77825)
lin_alg.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Claudia Emde <claudia.emde@dlr.de>
2  This program is free software; you can redistribute it and/or modify it
3  under the terms of the GNU General Public License as published by the
4  Free Software Foundation; either version 2, or (at your option) any
5  later version.
6 
7  This program is distributed in the hope that it will be useful,
8  but WITHOUT ANY WARRANTY; without even the implied warranty of
9  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  GNU General Public License for more details.
11 
12  You should have received a copy of the GNU General Public License
13  along with this program; if not, write to the Free Software
14  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
15  USA.
16 */
17 
29 #include "lin_alg.h"
30 
31 /*===========================================================================
32  === External declarations
33  ===========================================================================*/
34 
35 #include <Eigen/Dense>
36 #include <Eigen/Eigenvalues>
37 #include <cmath>
38 #include <new>
39 #include <stdexcept>
40 #include "array.h"
41 #include "arts.h"
42 #include "arts_omp.h"
43 #include "lapack.h"
44 #include "logic.h"
45 #include "matpackIII.h"
46 
48 
56 void ludcmp(Matrix& LU, ArrayOfIndex& indx, ConstMatrixView A) {
57  // Assert that A is quadratic.
58  const Index n = A.nrows();
59  assert(is_size(LU, n, n));
60 
61  int n_int, info;
62  int* ipiv = new int[n];
63 
64  LU = transpose(A);
65 
66  // Standard case: The arts matrix is not transposed, the leading
67  // dimension is the row stride of the matrix.
68  n_int = (int)n;
69 
70  // Compute LU decomposition using LAPACK dgetrf_.
71  lapack::dgetrf_(&n_int, &n_int, LU.mdata, &n_int, ipiv, &info);
72 
73  // Copy pivot array to pivot vector.
74  for (Index i = 0; i < n; i++) {
75  indx[i] = ipiv[i];
76  }
77  delete[] ipiv;
78 }
79 
81 
91 void lubacksub(VectorView x,
92  ConstMatrixView LU,
93  ConstVectorView b,
94  const ArrayOfIndex& indx) {
95  Index n = LU.nrows();
96 
97  /* Check if the dimensions of the input matrix and vectors agree and if LU
98  is a quadratic matrix.*/
99  DEBUG_ONLY(Index column_stride = LU.mcr.get_stride());
100  DEBUG_ONLY(Index vec_stride = b.mrange.get_stride());
101 
102  assert(is_size(LU, n, n));
103  assert(is_size(b, n));
104  assert(is_size(indx, n));
105  assert(column_stride == 1);
106  assert(vec_stride == 1);
107 
108  char trans = 'N';
109  int n_int = (int)n;
110  int one = (int)1;
111  int ipiv[n];
112  double rhs[n];
113  int info;
114 
115  for (Index i = 0; i < n; i++) {
116  ipiv[i] = (int)indx[i];
117  rhs[i] = b[i];
118  }
119 
120  lapack::dgetrs_(
121  &trans, &n_int, &one, LU.mdata, &n_int, ipiv, rhs, &n_int, &info);
122 
123  for (Index i = 0; i < n; i++) {
124  x[i] = rhs[i];
125  }
126 }
127 
129 
138 void solve(VectorView x, ConstMatrixView A, ConstVectorView b) {
139  Index n = A.ncols();
140 
141  // Check dimensions of the system.
142  assert(n == A.nrows());
143  assert(n == x.nelem());
144  assert(n == b.nelem());
145 
146  // Allocate matrix and index vector for the LU decomposition.
147  Matrix LU = Matrix(n, n);
148  ArrayOfIndex indx(n);
149 
150  // Perform LU decomposition.
151  ludcmp(LU, indx, A);
152 
153  // Solve the system using backsubstitution.
154  lubacksub(x, LU, b, indx);
155 }
156 
158 
167 void inv(MatrixView Ainv, ConstMatrixView A) {
168  // A must be a square matrix.
169  assert(A.ncols() == A.nrows());
170 
171  Index n = A.ncols();
172  Matrix LU(A);
173 
174  int n_int, info;
175  int* ipiv = new int[n];
176 
177  n_int = (int)n;
178 
179  // Compute LU decomposition using LAPACK dgetrf_.
180  lapack::dgetrf_(&n_int, &n_int, LU.mdata, &n_int, ipiv, &info);
181 
182  // Invert matrix.
183  int lwork = n_int;
184  double* work;
185 
186  try {
187  work = new double[lwork];
188  } catch (std::bad_alloc& ba) {
189  throw runtime_error(
190  "Error inverting matrix: Could not allocate workspace memory.");
191  }
192 
193  lapack::dgetri_(&n_int, LU.mdata, &n_int, ipiv, work, &lwork, &info);
194  delete[] work;
195  delete[] ipiv;
196 
197  // Check for success.
198  if (info == 0) {
199  Ainv = LU;
200  } else {
201  throw runtime_error("Error inverting matrix: Matrix not of full rank.");
202  }
203 }
204 
205 void inv(ComplexMatrixView Ainv, const ConstComplexMatrixView A) {
206  // A must be a square matrix.
207  assert(A.ncols() == A.nrows());
208 
209  Index n = A.ncols();
210 
211  // Workdata
212  Ainv = A;
213  int n_int = int(n), lwork = n_int, info;
214  std::vector<int> ipiv(n);
215  ComplexVector work(lwork);
216 
217  // Compute LU decomposition using LAPACK dgetrf_.
218  lapack::zgetrf_(&n_int, &n_int, Ainv.get_c_array(), &n_int, ipiv.data(), &info);
219  lapack::zgetri_(&n_int, Ainv.get_c_array(), &n_int, ipiv.data(), work.get_c_array(), &lwork, &info);
220 
221  // Check for success.
222  if (info not_eq 0) {
223  throw runtime_error("Error inverting matrix: Matrix not of full rank.");
224  }
225 }
226 
228 
245 void diagonalize(MatrixView P,
246  VectorView WR,
247  VectorView WI,
248  ConstMatrixView A) {
249  Index n = A.ncols();
250 
251  // A must be a square matrix.
252  assert(n == A.nrows());
253  assert(n == WR.nelem());
254  assert(n == WI.nelem());
255  assert(n == P.nrows());
256  assert(n == P.ncols());
257 
258  Matrix A_tmp = A;
259  Matrix P2 = P;
260  Vector WR2 = WR;
261  Vector WI2 = WI;
262 
263  // Integers
264  int LDA, LDA_L, LDA_R, n_int, info;
265  n_int = (int)n;
266  LDA = LDA_L = LDA_R = (int)A.mcr.get_extent();
267 
268  // We want to calculate RP not LP
269  char l_eig = 'N', r_eig = 'V';
270 
271  // Work matrix
272  int lwork = 2 * n_int + n_int * n_int;
273  double *work, *rwork;
274  try {
275  rwork = new double[2 * n_int];
276  work = new double[lwork];
277  } catch (std::bad_alloc& ba) {
278  throw std::runtime_error(
279  "Error diagonalizing: Could not allocate workspace memory.");
280  }
281 
282  // Memory references
283  double* adata = A_tmp.mdata;
284  double* rpdata = P2.mdata;
285  double* lpdata = new double[0]; //To not confuse the compiler
286  double* wrdata = WR2.mdata;
287  double* widata = WI2.mdata;
288 
289  // Main calculations. Note that errors in the output is ignored
290  lapack::dgeev_(&l_eig,
291  &r_eig,
292  &n_int,
293  adata,
294  &LDA,
295  wrdata,
296  widata,
297  lpdata,
298  &LDA_L,
299  rpdata,
300  &LDA_R,
301  work,
302  &lwork,
303  rwork,
304  &info);
305 
306  // Free memory. Can these be sent in to speed things up?
307  delete[] work;
308  delete[] rwork;
309  delete[] lpdata;
310 
311  // Re-order. This can be done better?
312  for (Index i = 0; i < n; i++)
313  for (Index j = 0; j < n; j++) P(j, i) = P2(i, j);
314  WI = WI2;
315  WR = WR2;
316 }
317 
319 
330 void diagonalize(ComplexMatrixView P,
331  ComplexVectorView W,
332  const ConstComplexMatrixView A) {
333  Index n = A.ncols();
334 
335  // A must be a square matrix.
336  assert(n == A.nrows());
337  assert(n == W.nelem());
338  assert(n == P.nrows());
339  assert(n == P.ncols());
340 
341  ComplexMatrix A_tmp = A;
342 
343  // Integers
344  int LDA=int(A.get_column_extent()), LDA_L=int(A.get_column_extent()), LDA_R=int(A.get_column_extent()), n_int=int(n), info;
345 
346  // We want to calculate RP not LP
347  char l_eig = 'N', r_eig = 'V';
348 
349  // Work matrix
350  int lwork = 2 * n_int + n_int * n_int;
351  ComplexVector work(lwork);
352  ComplexVector lpdata(0);
353  Vector rwork(2 * n_int);
354 
355  // Main calculations. Note that errors in the output is ignored
356  lapack::zgeev_(&l_eig,
357  &r_eig,
358  &n_int,
359  A_tmp.get_c_array(),
360  &LDA,
361  W.get_c_array(),
362  lpdata.get_c_array(),
363  &LDA_L,
364  P.get_c_array(),
365  &LDA_R,
366  work.get_c_array(),
367  &lwork,
368  rwork.get_c_array(),
369  &info);
370 
371  for (Index i = 0; i < n; i++)
372  for (Index j = 0; j <= i; j++)
373  std::swap(P(j, i), P(i, j));
374 }
375 
377 
391 void matrix_exp(MatrixView F, ConstMatrixView A, const Index& q) {
392  const Index n = A.ncols();
393 
394  /* Check if A and F are a quadratic and of the same dimension. */
395  assert(is_size(A, n, n));
396  assert(is_size(F, n, n));
397 
398  Numeric A_norm_inf, c;
399  Numeric j;
400  Matrix D(n, n), N(n, n), X(n, n), cX(n, n, 0.0), B(n, n, 0.0);
401  Vector N_col_vec(n, 0.), F_col_vec(n, 0.);
402 
403  A_norm_inf = norm_inf(A);
404 
405  // This formular is derived in the book by Golub and Van Loan.
406  j = 1 + floor(1. / log(2.) * log(A_norm_inf));
407 
408  if (j < 0) j = 0.;
409  Index j_index = (Index)(j);
410 
411  // Scale matrix
412  F = A;
413  F /= pow(2, j);
414 
415  /* The higher q the more accurate is the computation,
416  see user guide for accuracy */
417  // Index q = 8;
418  Numeric q_n = (Numeric)(q);
419  id_mat(D);
420  id_mat(N);
421  id_mat(X);
422  c = 1.;
423 
424  for (Index k = 0; k < q; k++) {
425  Numeric k_n = (Numeric)(k + 1);
426  c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
427  mult(B, F, X); // X = F * X
428  X = B;
429  cX = X;
430  cX *= c; // cX = X*c
431  N += cX; // N = N + X*c
432  cX *= pow(-1, k_n); // cX = (-1)^k*c*X
433  D += cX; // D = D + (-1)^k*c*X
434  }
435 
436  /*Solve the equation system DF=N for F using LU decomposition,
437  use the backsubstitution routine for columns of N*/
438 
439  /* Now use X for the LU decomposition matrix of D.*/
440  ArrayOfIndex indx(n);
441 
442  ludcmp(X, indx, D);
443 
444  for (Index i = 0; i < n; i++) {
445  N_col_vec = N(joker, i); // extract column vectors of N
446  lubacksub(F_col_vec, X, N_col_vec, indx);
447  F(joker, i) = F_col_vec; // construct F matrix from column vectors
448  }
449 
450  /* The square of F gives the result. */
451  for (Index k = 0; k < j_index; k++) {
452  mult(B, F, F); // F = F^2
453  F = B;
454  }
455 }
456 
457 //Matrix exponent with decomposition
458 void matrix_exp2(MatrixView F, ConstMatrixView A) {
459  const Index n = F.nrows();
460 
461  assert(is_size(F, n, n));
462  assert(is_size(A, n, n));
463 
464  Matrix P(n, n), invP(n, n);
465  Vector WR(n), WI(n);
466 
467  diagonalize(P, WR, WI, A);
468  inv(invP, P);
469 
470  for (Index k = 0; k < n; k++) {
471  if (WI[k] != 0)
472  throw std::runtime_error(
473  "We require real eigenvalues for your chosen method.");
474  P(joker, k) *= exp(WR[k]);
475  }
476 
477  //mult(tmp,Wdiag,invP);
478  mult(F, P, invP);
479 }
480 
481 void matrix_exp_4x4(MatrixView arts_f, ConstMatrixView A, const Index& q) {
482  /* Check if A and F are a quadratic and of the same dimension. */
483  assert(is_size(A, 4, 4));
484  assert(is_size(arts_f, 4, 4));
485 
486  // Set constants and Numerics
487  const Numeric A_norm_inf = norm_inf(A);
488  const Numeric e = 1.0 + floor(1.0 / log(2.0) * log(A_norm_inf));
489  const Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
490  const Numeric inv_pow2 = 1. / pow(2, r);
491  Numeric c = 0.5;
492 
493  const Eigen::Matrix4d M = MapToEigen4x4(A) * inv_pow2;
494  Eigen::Matrix4d X = M;
495  Eigen::Matrix4d cX = c * X;
496  Matrix4x4ViewMap F = MapToEigen4x4(arts_f);
497  Eigen::Matrix4d D = Eigen::Matrix4d::Identity();
498  F.setIdentity();
499  F.noalias() += cX;
500  D.noalias() -= cX;
501 
502  bool p = true;
503  for (Index k = 2; k <= q; k++) {
504  c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
505  X = M * X;
506  cX.noalias() = c * X;
507  F.noalias() += cX;
508  if (p)
509  D.noalias() += cX;
510  else
511  D.noalias() -= cX;
512  p = not p;
513  }
514 
515  F = D.inverse() * F;
516 
517  for (Index k = 1; k <= r; k++) F *= F;
518 }
519 
520 //Matrix exponent with decomposition
521 void matrix_exp2_4x4(MatrixView arts_f, ConstMatrixView A) {
522  assert(is_size(arts_f, 4, 4));
523  assert(is_size(A, 4, 4));
524 
525  Matrix4x4ViewMap F = MapToEigen4x4(arts_f);
526  Eigen::EigenSolver<Eigen::Matrix4d> es;
527  es.compute(MapToEigen4x4(A));
528 
529  const Eigen::Matrix4cd U = es.eigenvectors();
530  const Eigen::Vector4cd q = es.eigenvalues();
531  const Eigen::Vector4cd eq = q.array().exp();
532  Eigen::Matrix4cd res = U * eq.asDiagonal();
533  res *= U.inverse();
534  F = res.real();
535 }
536 
538 
559 void special_matrix_exp_and_dmatrix_exp_dx_for_rt(MatrixView F,
560  Tensor3View dF_upp,
561  Tensor3View dF_low,
562  ConstMatrixView A,
563  ConstTensor3View dA_upp,
564  ConstTensor3View dA_low,
565  const Index& q) {
566  const Index n_partials = dA_upp.npages();
567 
568  const Index n = A.ncols();
569 
570  /* Check if A and F are a quadratic and of the same dimension. */
571  assert(is_size(A, n, n));
572  assert(is_size(F, n, n));
573  assert(n_partials == dF_upp.npages());
574  assert(n_partials == dF_low.npages());
575  assert(n_partials == dA_low.npages());
576  for (Index ii = 0; ii < n_partials; ii++) {
577  assert(is_size(dA_upp(ii, joker, joker), n, n));
578  assert(is_size(dA_low(ii, joker, joker), n, n));
579  assert(is_size(dF_upp(ii, joker, joker), n, n));
580  assert(is_size(dF_low(ii, joker, joker), n, n));
581  }
582 
583  // This sets up some constants
584  Numeric A_norm_inf, e;
585  A_norm_inf = norm_inf(A);
586  e = 1. + floor(1. / log(2.) * log(A_norm_inf));
587  Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
588  Numeric pow2rm1 = 1. / pow(2, r);
589  Numeric c = 0.5;
590 
591  // For non-partials
592  Matrix M = A, X(n, n), cX(n, n), D(n, n);
593  M *= pow2rm1;
594  X = M;
595  cX = X;
596  cX *= c;
597  id_mat(F);
598  F += cX;
599  id_mat(D);
600  D -= cX;
601 
602  // For partials in parallel
603  Tensor3 dM_upp(n_partials, n, n), dM_low(n_partials, n, n),
604  dD_upp(n_partials, n, n), dD_low(n_partials, n, n),
605  Y_upp(n_partials, n, n), Y_low(n_partials, n, n),
606  cY_upp(n_partials, n, n), cY_low(n_partials, n, n);
607  for (Index ii = 0; ii < n_partials; ii++) {
608  for (Index jj = 0; jj < n; jj++) {
609  for (Index kk = 0; kk < n; kk++) {
610  dM_upp(ii, jj, kk) = dA_upp(ii, jj, kk) * pow2rm1;
611  dM_low(ii, jj, kk) = dA_low(ii, jj, kk) * pow2rm1;
612 
613  Y_upp(ii, jj, kk) = dM_upp(ii, jj, kk);
614  Y_low(ii, jj, kk) = dM_low(ii, jj, kk);
615 
616  cY_upp(ii, jj, kk) = c * Y_upp(ii, jj, kk);
617  cY_low(ii, jj, kk) = c * Y_low(ii, jj, kk);
618 
619  dF_upp(ii, jj, kk) = cY_upp(ii, jj, kk);
620  dF_low(ii, jj, kk) = cY_low(ii, jj, kk);
621 
622  dD_upp(ii, jj, kk) = -cY_upp(ii, jj, kk);
623  dD_low(ii, jj, kk) = -cY_low(ii, jj, kk);
624  }
625  }
626  }
627 
628  // NOTE: MATLAB paper sets q = 6 but we allow other numbers
629  Matrix tmp1(n, n), tmp2(n, n);
630 
631  for (Index k = 2; k <= q; k++) {
632  c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
633 
634  // For partials in parallel
635  for (Index ii = 0; ii < n_partials; ii++) {
636  Matrix tmp_low1(n, n), tmp_upp1(n, n), tmp_low2(n, n), tmp_upp2(n, n);
637 
638  // Y = dM*X + M*Y
639  mult(tmp_upp1, dM_upp(ii, joker, joker), X);
640  mult(tmp_upp2, M, Y_upp(ii, joker, joker));
641  mult(tmp_low1, dM_low(ii, joker, joker), X);
642  mult(tmp_low2, M, Y_low(ii, joker, joker));
643 
644  for (Index jj = 0; jj < n; jj++) {
645  for (Index kk = 0; kk < n; kk++) {
646  Y_upp(ii, jj, kk) = tmp_upp1(jj, kk) + tmp_upp2(jj, kk);
647  Y_low(ii, jj, kk) = tmp_low1(jj, kk) + tmp_low2(jj, kk);
648 
649  // cY = c*Y
650  cY_upp(ii, jj, kk) = c * Y_upp(ii, jj, kk);
651  cY_low(ii, jj, kk) = c * Y_low(ii, jj, kk);
652 
653  // dF = dF + cY
654  dF_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
655  dF_low(ii, jj, kk) += cY_low(ii, jj, kk);
656 
657  if (k % 2 == 0) //For even numbers, add.
658  {
659  // dD = dD + cY
660  dD_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
661  dD_low(ii, jj, kk) += cY_low(ii, jj, kk);
662  } else {
663  // dD = dD - cY
664  dD_upp(ii, jj, kk) -= cY_upp(ii, jj, kk);
665  dD_low(ii, jj, kk) -= cY_low(ii, jj, kk);
666  }
667  }
668  }
669  }
670 
671  //For full derivative (Note X in partials above)
672  // X=M*X
673  mult(tmp1, M, X);
674  for (Index jj = 0; jj < n; jj++) {
675  for (Index kk = 0; kk < n; kk++) {
676  X(jj, kk) = tmp1(jj, kk);
677  cX(jj, kk) = tmp1(jj, kk) * c;
678  F(jj, kk) += cX(jj, kk);
679 
680  if (k % 2 == 0) //For even numbers, D = D + cX
681  D(jj, kk) += cX(jj, kk);
682  else //For odd numbers, D = D - cX
683  D(jj, kk) -= cX(jj, kk);
684  }
685  }
686  }
687 
688  // D^-1
689  inv(tmp1, D);
690 
691  // F = D\F, or D^-1*F
692  mult(tmp2, tmp1, F);
693  F = tmp2;
694 
695  // For partials in parallel
696  for (Index ii = 0; ii < n_partials; ii++) {
697  Matrix tmp_low(n, n), tmp_upp(n, n);
698  //dF = D \ (dF - dD*F), or D^-1 * (dF - dD*F)
699  mult(tmp_upp, dD_upp(ii, joker, joker), F);
700  mult(tmp_low, dD_low(ii, joker, joker), F);
701 
702  for (Index jj = 0; jj < n; jj++) {
703  for (Index kk = 0; kk < n; kk++) {
704  dF_upp(ii, jj, kk) -= tmp_upp(jj, kk); // dF - dD * F
705  dF_low(ii, jj, kk) -= tmp_low(jj, kk); // dF - dD * F
706  }
707  }
708 
709  mult(tmp_upp, tmp1, dF_upp(ii, joker, joker));
710  mult(tmp_low, tmp1, dF_low(ii, joker, joker));
711 
712  for (Index jj = 0; jj < n; jj++) {
713  for (Index kk = 0; kk < n; kk++) {
714  dF_upp(ii, jj, kk) = tmp_upp(jj, kk);
715  dF_low(ii, jj, kk) = tmp_low(jj, kk);
716  }
717  }
718  }
719 
720  for (Index k = 1; k <= r; k++) {
721  // For partials in parallel
722  for (Index ii = 0; ii < n_partials; ii++) {
723  Matrix tmp_low1(n, n), tmp_upp1(n, n), tmp_low2(n, n), tmp_upp2(n, n);
724 
725  // dF=F*dF+dF*F
726  mult(tmp_upp1, F, dF_upp(ii, joker, joker)); //F*dF
727  mult(tmp_upp2, dF_upp(ii, joker, joker), F); //dF*F
728  mult(tmp_low1, F, dF_low(ii, joker, joker)); //F*dF
729  mult(tmp_low2, dF_low(ii, joker, joker), F); //dF*F
730 
731  for (Index jj = 0; jj < n; jj++) {
732  for (Index kk = 0; kk < n; kk++) {
733  dF_upp(ii, jj, kk) = tmp_upp1(jj, kk) + tmp_upp2(jj, kk);
734  dF_low(ii, jj, kk) = tmp_low1(jj, kk) + tmp_low2(jj, kk);
735  }
736  }
737  }
738 
739  // F=F*F
740  mult(tmp1, F, F);
741  F = tmp1;
742  }
743 }
744 
745 Matrix4x4ViewMap MapToEigen4x4(Tensor3View& A, const Index& i) {
746  MatrixView B = A(i, joker, joker);
747  return MapToEigen4x4(B);
748 }
749 
750 MatrixViewMap MapToEigen(Tensor3View& A, const Index& i) {
751  MatrixView B = A(i, joker, joker);
752  return MapToEigen(B);
753 }
754 
789 void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(
790  MatrixView F, ConstMatrixView A) {
791  static const Numeric sqrt_05 = sqrt(0.5);
792 
793  const Numeric a = A(0, 0), b = A(0, 1), c = A(0, 2), d = A(0, 3), u = A(1, 2),
794  v = A(1, 3), w = A(2, 3);
795  const Numeric exp_a = exp(a);
796  const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
797  w2 = w * w;
798 
799  const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
800 
801  Numeric Const1;
802  Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
803  Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
804  Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
805  Const1 += u2 * (u2 * 0.5 + v2 + w2);
806  Const1 += v2 * (v2 * 0.5 + w2);
807  Const1 *= 2;
808  Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
809  Const1 += w2 * w2;
810 
811  if (Const1 > 0.0)
812  Const1 = sqrt(Const1);
813  else
814  Const1 = 0.0;
815 
816  if (Const1 == 0 and
817  Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
818  {
819  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
820  return;
821  }
822 
823  const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
824  const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
825  const Numeric x = sqrt_BpA.real() * sqrt_05;
826  const Numeric y = sqrt_BmA.imag() * sqrt_05;
827  const Numeric x2 = x * x;
828  const Numeric y2 = y * y;
829  const Numeric cos_y = cos(y);
830  const Numeric sin_y = sin(y);
831  const Numeric cosh_x = cosh(x);
832  const Numeric sinh_x = sinh(x);
833  const Numeric x2y2 = x2 + y2;
834  const Numeric inv_x2y2 = 1.0 / x2y2;
835 
836  Numeric C0, C1, C2, C3;
837  Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
838 
839  // X and Y cannot both be zero
840  if (x == 0.0) {
841  inv_y = 1.0 / y;
842  C0 = 1.0;
843  C1 = 1.0;
844  C2 = (1.0 - cos_y) * inv_x2y2;
845  C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
846  } else if (y == 0.0) {
847  inv_x = 1.0 / x;
848  C0 = 1.0;
849  C1 = 1.0;
850  C2 = (cosh_x - 1.0) * inv_x2y2;
851  C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
852  } else {
853  inv_x = 1.0 / x;
854  inv_y = 1.0 / y;
855 
856  C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
857  C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
858  C2 = (cosh_x - cos_y) * inv_x2y2;
859  C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
860  }
861 
862  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
863 
864  F(0, 0) += C2 * (b2 + c2 + d2);
865 
866  F(1, 1) += C2 * (b2 - u2 - v2);
867 
868  F(2, 2) += C2 * (c2 - u2 - w2);
869 
870  F(3, 3) += C2 * (d2 - v2 - w2);
871 
872  F(0, 1) = F(1, 0) = C1 * b;
873 
874  F(0, 1) +=
875  C2 * (-c * u - d * v) +
876  C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
877 
878  F(1, 0) +=
879  C2 * (c * u + d * v) +
880  C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) + d * (b * d + u * w));
881 
882  F(0, 2) = F(2, 0) = C1 * c;
883 
884  F(0, 2) +=
885  C2 * (b * u - d * w) +
886  C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
887 
888  F(2, 0) +=
889  C2 * (-b * u + d * w) +
890  C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) + d * (c * d - u * v));
891 
892  F(0, 3) = F(3, 0) = C1 * d;
893 
894  F(0, 3) +=
895  C2 * (b * v + c * w) +
896  C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
897 
898  F(3, 0) +=
899  C2 * (-b * v - c * w) +
900  C3 * (b * (b * d + u * w) + c * (c * d - u * v) - d * (-d2 + v2 + w2));
901 
902  F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
903 
904  F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
905  w * (b * d + u * w));
906 
907  F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
908  v * (c * d - u * v));
909 
910  F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
911 
912  F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
913  w * (b * c - v * w));
914 
915  F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
916  v * (-d2 + v2 + w2));
917 
918  F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
919 
920  F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
921  w * (-c2 + u2 + w2));
922 
923  F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
924  w * (-d2 + v2 + w2));
925 
926  F *= exp_a;
927 }
928 
963 void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(
964  MatrixView F, ConstMatrixView A) {
965  Matrix4x4ViewMap eigF = MapToEigen4x4(F);
966  Eigen::Matrix4d eigA = MapToEigen4x4(A);
967  // eigA << 0, A(0, 1), A(0, 2), A(0, 3),
968  // A(0, 1), 0, A(1, 2), A(1, 3),
969  // A(0, 2), -A(1, 2), 0, A(2, 3),
970  // A(0, 3), -A(1, 3), -A(2, 3), 0;
971  eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
972 
973  Eigen::Matrix4d eigA2 = eigA;
974  eigA2 *= eigA;
975  Eigen::Matrix4d eigA3 = eigA2;
976  eigA3 *= eigA;
977 
978  static const Numeric sqrt_05 = sqrt(0.5);
979 
980  const Numeric a = A(0, 0), b = A(0, 1), c = A(0, 2), d = A(0, 3), u = A(1, 2),
981  v = A(1, 3), w = A(2, 3);
982  const Numeric exp_a = exp(a);
983  const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
984  w2 = w * w;
985 
986  const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
987 
988  Numeric Const1;
989  Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
990  Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
991  Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
992  Const1 += u2 * (u2 * 0.5 + v2 + w2);
993  Const1 += v2 * (v2 * 0.5 + w2);
994  Const1 *= 2;
995  Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
996  Const1 += w2 * w2;
997 
998  if (Const1 > 0.0)
999  Const1 = sqrt(Const1);
1000  else
1001  Const1 = 0.0;
1002 
1003  if (Const1 == 0 and
1004  Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1005  {
1006  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1007  return;
1008  }
1009 
1010  const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1011  const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1012  const Numeric x = sqrt_BpA.real() * sqrt_05;
1013  const Numeric y = sqrt_BmA.imag() * sqrt_05;
1014  const Numeric x2 = x * x;
1015  const Numeric y2 = y * y;
1016  const Numeric cos_y = cos(y);
1017  const Numeric sin_y = sin(y);
1018  const Numeric cosh_x = cosh(x);
1019  const Numeric sinh_x = sinh(x);
1020  const Numeric x2y2 = x2 + y2;
1021  const Numeric inv_x2y2 = 1.0 / x2y2;
1022 
1023  Numeric C0, C1, C2, C3;
1024  Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1025 
1026  // X and Y cannot both be zero
1027  if (x == 0.0) {
1028  inv_y = 1.0 / y;
1029  C0 = 1.0;
1030  C1 = 1.0;
1031  C2 = (1.0 - cos_y) * inv_x2y2;
1032  C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1033  } else if (y == 0.0) {
1034  inv_x = 1.0 / x;
1035  C0 = 1.0;
1036  C1 = 1.0;
1037  C2 = (cosh_x - 1.0) * inv_x2y2;
1038  C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1039  } else {
1040  inv_x = 1.0 / x;
1041  inv_y = 1.0 / y;
1042 
1043  C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1044  C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1045  C2 = (cosh_x - cos_y) * inv_x2y2;
1046  C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1047  }
1048 
1049  eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1050  eigF(0, 0) += C0;
1051  eigF(1, 1) += C0;
1052  eigF(2, 2) += C0;
1053  eigF(3, 3) += C0;
1054  eigF *= exp_a;
1055 }
1056 
1111 void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(
1112  MatrixView F,
1113  Tensor3View dF_upp,
1114  Tensor3View dF_low,
1115  ConstMatrixView A,
1116  ConstTensor3View dA_upp,
1117  ConstTensor3View dA_low) {
1118  static const Numeric sqrt_05 = sqrt(0.5);
1119  const Index npd = dF_upp.npages();
1120 
1121  const Numeric a = A(0, 0), b = A(0, 1), c = A(0, 2), d = A(0, 3), u = A(1, 2),
1122  v = A(1, 3), w = A(2, 3);
1123  const Numeric exp_a = exp(a);
1124  const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
1125  w2 = w * w;
1126 
1127  const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1128 
1129  Numeric Const1;
1130  Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1131  Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1132  Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
1133  Const1 += u2 * (u2 * 0.5 + v2 + w2);
1134  Const1 += v2 * (v2 * 0.5 + w2);
1135  Const1 *= 2;
1136  Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
1137  Const1 += w2 * w2;
1138 
1139  if (Const1 > 0.0)
1140  Const1 = sqrt(Const1);
1141  else
1142  Const1 = 0.0;
1143 
1144  if (Const1 == 0 and
1145  Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1146  {
1147  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1148  for (Index i = 0; i < npd; i++) {
1149  // mult incase the derivative is still non-zero, e.g., as when mag-field is considered 0 but derivative is computed
1150  mult(dF_upp(i, joker, joker), F, dA_upp(i, joker, joker));
1151  mult(dF_low(i, joker, joker), F, dA_low(i, joker, joker));
1152  }
1153  return;
1154  }
1155 
1156  const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1157  const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1158  const Numeric x = sqrt_BpA.real() * sqrt_05;
1159  const Numeric y = sqrt_BmA.imag() * sqrt_05;
1160  const Numeric x2 = x * x;
1161  const Numeric y2 = y * y;
1162  const Numeric cos_y = cos(y);
1163  const Numeric sin_y = sin(y);
1164  const Numeric cosh_x = cosh(x);
1165  const Numeric sinh_x = sinh(x);
1166  const Numeric x2y2 = x2 + y2;
1167  const Numeric inv_x2y2 = 1.0 / x2y2;
1168 
1169  Numeric C0, C1, C2, C3;
1170  Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1171 
1172  // X and Y cannot both be zero
1173  if (x == 0.0) {
1174  inv_y = 1.0 / y;
1175  C0 = 1.0;
1176  C1 = 1.0;
1177  C2 = (1.0 - cos_y) * inv_x2y2;
1178  C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1179  } else if (y == 0.0) {
1180  inv_x = 1.0 / x;
1181  C0 = 1.0;
1182  C1 = 1.0;
1183  C2 = (cosh_x - 1.0) * inv_x2y2;
1184  C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1185  } else {
1186  inv_x = 1.0 / x;
1187  inv_y = 1.0 / y;
1188 
1189  C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1190  C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1191  C2 = (cosh_x - cos_y) * inv_x2y2;
1192  C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1193  }
1194 
1195  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
1196 
1197  F(0, 0) += C2 * (b2 + c2 + d2);
1198 
1199  F(1, 1) += C2 * (b2 - u2 - v2);
1200 
1201  F(2, 2) += C2 * (c2 - u2 - w2);
1202 
1203  F(3, 3) += C2 * (d2 - v2 - w2);
1204 
1205  F(0, 1) = F(1, 0) = C1 * b;
1206 
1207  F(0, 1) +=
1208  C2 * (-c * u - d * v) +
1209  C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
1210 
1211  F(1, 0) +=
1212  C2 * (c * u + d * v) +
1213  C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) + d * (b * d + u * w));
1214 
1215  F(0, 2) = F(2, 0) = C1 * c;
1216 
1217  F(0, 2) +=
1218  C2 * (b * u - d * w) +
1219  C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
1220 
1221  F(2, 0) +=
1222  C2 * (-b * u + d * w) +
1223  C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) + d * (c * d - u * v));
1224 
1225  F(0, 3) = F(3, 0) = C1 * d;
1226 
1227  F(0, 3) +=
1228  C2 * (b * v + c * w) +
1229  C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
1230 
1231  F(3, 0) +=
1232  C2 * (-b * v - c * w) +
1233  C3 * (b * (b * d + u * w) + c * (c * d - u * v) - d * (-d2 + v2 + w2));
1234 
1235  F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
1236 
1237  F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1238  w * (b * d + u * w));
1239 
1240  F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1241  v * (c * d - u * v));
1242 
1243  F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
1244 
1245  F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1246  w * (b * c - v * w));
1247 
1248  F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1249  v * (-d2 + v2 + w2));
1250 
1251  F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
1252 
1253  F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1254  w * (-c2 + u2 + w2));
1255 
1256  F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1257  w * (-d2 + v2 + w2));
1258 
1259  F *= exp_a;
1260 
1261  if (npd) {
1262  const Numeric inv_x2 = inv_x * inv_x;
1263  const Numeric inv_y2 = inv_y * inv_y;
1264 
1265  for (Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1266  for (Index i = 0; i < npd; i++) {
1267  MatrixView dF =
1268  upp_or_low ? dF_low(i, joker, joker) : dF_upp(i, joker, joker);
1269  ConstMatrixView dA =
1270  upp_or_low ? dA_low(i, joker, joker) : dA_upp(i, joker, joker);
1271 
1272  const Numeric da = dA(0, 0), db = dA(0, 1), dc = dA(0, 2),
1273  dd = dA(0, 3), du = dA(1, 2), dv = dA(1, 3),
1274  dw = dA(2, 3);
1275 
1276  const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1277  du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1278 
1279  const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1280 
1281  Numeric dConst1;
1282  if (Const1 > 0.) {
1283  dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1284  dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1285 
1286  dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1287  dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1288 
1289  dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1290  dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1291 
1292  dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1293  dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1294 
1295  dConst1 += dv2 * (v2 * 0.5 + w2);
1296  dConst1 += v2 * (dv2 * 0.5 + dw2);
1297 
1298  dConst1 += 4 * ((db * d * u * w - db * c * v * w - dc * d * u * v +
1299  b * dd * u * w - b * dc * v * w - c * dd * u * v +
1300  b * d * du * w - b * c * dv * w - c * d * du * v +
1301  b * d * u * dw - b * c * v * dw - c * d * u * dv));
1302  dConst1 += dw2 * w2;
1303  dConst1 /= Const1;
1304  } else
1305  dConst1 = 0.0;
1306 
1307  Numeric dC0, dC1, dC2, dC3;
1308  if (x == 0.0) {
1309  const Numeric dy =
1310  (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1311 
1312  dC0 = 0.0;
1313  dC1 = 0.0;
1314  dC2 = -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1315  dC3 = -2 * y * dy * C3 * inv_x2y2 +
1316  (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1317  ;
1318  } else if (y == 0.0) {
1319  const Numeric dx =
1320  (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1321 
1322  dC0 = 0.0;
1323  dC1 = 0.0;
1324  dC2 = -2 * x * dx * C2 * inv_x2y2 + dx * sinh_x * inv_x2y2;
1325  dC3 = -2 * x * dx * C3 * inv_x2y2 +
1326  (cosh_x * dx * inv_x - dx * sinh_x * inv_x2) * inv_x2y2;
1327  } else {
1328  const Numeric dx =
1329  (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1330  const Numeric dy =
1331  (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1332  const Numeric dy2 = 2 * y * dy;
1333  const Numeric dx2 = 2 * x * dx;
1334  const Numeric dx2dy2 = dx2 + dy2;
1335 
1336  dC0 = -dx2dy2 * C0 * inv_x2y2 +
1337  (2 * cos_y * dx * x + 2 * cosh_x * dy * y + dx * sinh_x * y2 -
1338  dy * sin_y * x2) *
1339  inv_x2y2;
1340 
1341  dC1 = -dx2dy2 * C1 * inv_x2y2 +
1342  (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1343  dy * sin_y * x2 * inv_y2 - dx * sinh_x * y2 * inv_x2 +
1344  cosh_x * dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1345  inv_x2y2;
1346 
1347  dC2 = -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1348 
1349  dC3 = -dx2dy2 * C3 * inv_x2y2 +
1350  (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1351  cosh_x * dx * inv_x - dx * sinh_x * inv_x2) *
1352  inv_x2y2;
1353  }
1354 
1355  dF(0, 0) = dF(1, 1) = dF(2, 2) = dF(3, 3) = dC0;
1356 
1357  dF(0, 0) += dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2);
1358 
1359  dF(1, 1) += dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2);
1360 
1361  dF(2, 2) += dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2);
1362 
1363  dF(3, 3) += dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2);
1364 
1365  dF(0, 1) = dF(1, 0) = db * C1 + b * dC1;
1366 
1367  dF(0, 1) += dC2 * (-c * u - d * v) +
1368  C2 * (-dc * u - dd * v - c * du - d * dv) +
1369  dC3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) -
1370  v * (b * v + c * w)) +
1371  C3 * (db * (b2 + c2 + d2) - du * (b * u - d * w) -
1372  dv * (b * v + c * w) + b * (db2 + dc2 + dd2) -
1373  u * (db * u - dd * w) - v * (db * v + dc * w) -
1374  u * (b * du - d * dw) - v * (b * dv + c * dw));
1375  dF(1, 0) += dC2 * (c * u + d * v) +
1376  C2 * (dc * u + dd * v + c * du + d * dv) +
1377  dC3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) +
1378  d * (b * d + u * w)) +
1379  C3 * (-db * (-b2 + u2 + v2) + dc * (b * c - v * w) +
1380  dd * (b * d + u * w) - b * (-db2 + du2 + dv2) +
1381  c * (db * c - dv * w) + d * (db * d + du * w) +
1382  c * (b * dc - v * dw) + d * (b * dd + u * dw));
1383 
1384  dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
1385 
1386  dF(0, 2) += dC2 * (b * u - d * w) +
1387  C2 * (db * u - dd * w + b * du - d * dw) +
1388  dC3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) -
1389  w * (b * v + c * w)) +
1390  C3 * (dc * (b2 + c2 + d2) - du * (c * u + d * v) -
1391  dw * (b * v + c * w) + c * (db2 + dc2 + dd2) -
1392  u * (dc * u + dd * v) - w * (db * v + dc * w) -
1393  u * (c * du + d * dv) - w * (b * dv + c * dw));
1394  dF(2, 0) += dC2 * (-b * u + d * w) +
1395  C2 * (-db * u + dd * w - b * du + d * dw) +
1396  dC3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) +
1397  d * (c * d - u * v)) +
1398  C3 * (db * (b * c - v * w) - dc * (-c2 + u2 + w2) +
1399  dd * (c * d - u * v) + b * (db * c - dv * w) -
1400  c * (-dc2 + du2 + dw2) + d * (dc * d - du * v) +
1401  b * (b * dc - v * dw) + d * (c * dd - u * dv));
1402 
1403  dF(0, 3) = dF(3, 0) = dC1 * d + C1 * dd;
1404 
1405  dF(0, 3) += dC2 * (b * v + c * w) +
1406  C2 * (db * v + dc * w + b * dv + c * dw) +
1407  dC3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) +
1408  w * (b * u - d * w)) +
1409  C3 * (dd * (b2 + c2 + d2) - dv * (c * u + d * v) +
1410  dw * (b * u - d * w) + d * (db2 + dc2 + dd2) -
1411  v * (dc * u + dd * v) + w * (db * u - dd * w) -
1412  v * (c * du + d * dv) + w * (b * du - d * dw));
1413  dF(3, 0) += dC2 * (-b * v - c * w) +
1414  C2 * (-db * v - dc * w - b * dv - c * dw) +
1415  dC3 * (b * (b * d + u * w) + c * (c * d - u * v) -
1416  d * (-d2 + v2 + w2)) +
1417  C3 * (db * (b * d + u * w) + dc * (c * d - u * v) -
1418  dd * (-d2 + v2 + w2) + b * (db * d + du * w) +
1419  c * (dc * d - du * v) - d * (-dd2 + dv2 + dw2) +
1420  b * (b * dd + u * dw) + c * (c * dd - u * dv));
1421 
1422  dF(1, 2) = dF(2, 1) =
1423  dC2 * (b * c - v * w) + C2 * (db * c + b * dc - dv * w - v * dw);
1424 
1425  dF(1, 2) += dC1 * u + C1 * du +
1426  dC3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1427  w * (b * d + u * w)) +
1428  C3 * (dc * (c * u + d * v) - du * (-b2 + u2 + v2) -
1429  dw * (b * d + u * w) + c * (dc * u + dd * v) -
1430  u * (-db2 + du2 + dv2) - w * (db * d + du * w) +
1431  c * (c * du + d * dv) - w * (b * dd + u * dw));
1432  dF(2, 1) += -dC1 * u - C1 * du +
1433  dC3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1434  v * (c * d - u * v)) +
1435  C3 * (-db * (b * u - d * w) + du * (-c2 + u2 + w2) -
1436  dv * (c * d - u * v) - b * (db * u - dd * w) +
1437  u * (-dc2 + du2 + dw2) - v * (dc * d - du * v) -
1438  b * (b * du - d * dw) - v * (c * dd - u * dv));
1439 
1440  dF(1, 3) = dF(3, 1) =
1441  dC2 * (b * d + u * w) + C2 * (db * d + b * dd + du * w + u * dw);
1442 
1443  dF(1, 3) += dC1 * v + C1 * dv +
1444  dC3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1445  w * (b * c - v * w)) +
1446  C3 * (dd * (c * u + d * v) - dv * (-b2 + u2 + v2) +
1447  dw * (b * c - v * w) + d * (dc * u + dd * v) -
1448  v * (-db2 + du2 + dv2) + w * (db * c - dv * w) +
1449  d * (c * du + d * dv) + w * (b * dc - v * dw));
1450  dF(3, 1) += -dC1 * v - C1 * dv +
1451  dC3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1452  v * (-d2 + v2 + w2)) +
1453  C3 * (-db * (b * v + c * w) - du * (c * d - u * v) +
1454  dv * (-d2 + v2 + w2) - b * (db * v + dc * w) -
1455  u * (dc * d - du * v) + v * (-dd2 + dv2 + dw2) -
1456  b * (b * dv + c * dw) - u * (c * dd - u * dv));
1457 
1458  dF(2, 3) = dF(3, 2) =
1459  dC2 * (c * d - u * v) + C2 * (dc * d + c * dd - du * v - u * dv);
1460 
1461  dF(2, 3) += dC1 * w + C1 * dw +
1462  dC3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1463  w * (-c2 + u2 + w2)) +
1464  C3 * (-dd * (b * u - d * w) + dv * (b * c - v * w) -
1465  dw * (-c2 + u2 + w2) - d * (db * u - dd * w) +
1466  v * (db * c - dv * w) - w * (-dc2 + du2 + dw2) -
1467  d * (b * du - d * dw) + v * (b * dc - v * dw));
1468  dF(3, 2) += -dC1 * w - C1 * dw +
1469  dC3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1470  w * (-d2 + v2 + w2)) +
1471  C3 * (-dc * (b * v + c * w) + du * (b * d + u * w) +
1472  dw * (-d2 + v2 + w2) - c * (db * v + dc * w) +
1473  u * (db * d + du * w) + w * (-dd2 + dv2 + dw2) -
1474  c * (b * dv + c * dw) + u * (b * dd + u * dw));
1475 
1476  dF *= exp_a;
1477 
1478  // Finalize derivation by the chian rule
1479  dF(0, 0) += F(0, 0) * da;
1480  dF(0, 1) += F(0, 1) * da;
1481  dF(0, 2) += F(0, 2) * da;
1482  dF(0, 3) += F(0, 3) * da;
1483  dF(1, 0) += F(1, 0) * da;
1484  dF(1, 1) += F(1, 1) * da;
1485  dF(1, 2) += F(1, 2) * da;
1486  dF(1, 3) += F(1, 3) * da;
1487  dF(2, 0) += F(2, 0) * da;
1488  dF(2, 1) += F(2, 1) * da;
1489  dF(2, 2) += F(2, 2) * da;
1490  dF(2, 3) += F(2, 3) * da;
1491  dF(3, 0) += F(3, 0) * da;
1492  dF(3, 1) += F(3, 1) * da;
1493  dF(3, 2) += F(3, 2) * da;
1494  dF(3, 3) += F(3, 3) * da;
1495  }
1496  }
1497  }
1498 }
1499 
1550 void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(
1551  MatrixView F,
1552  Tensor3View dF_upp,
1553  Tensor3View dF_low,
1554  ConstMatrixView A,
1555  ConstTensor3View dA_upp,
1556  ConstTensor3View dA_low) {
1557  Matrix4x4ViewMap eigF = MapToEigen4x4(F);
1558  Eigen::Matrix4d eigA = MapToEigen4x4(A);
1559  // eigA << 0, A(0, 1), A(0, 2), A(0, 3),
1560  // A(0, 1), 0, A(1, 2), A(1, 3),
1561  // A(0, 2), -A(1, 2), 0, A(2, 3),
1562  // A(0, 3), -A(1, 3), -A(2, 3), 0;
1563  eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
1564 
1565  Eigen::Matrix4d eigA2 = eigA;
1566  eigA2 *= eigA;
1567  Eigen::Matrix4d eigA3 = eigA2;
1568  eigA3 *= eigA;
1569 
1570  static const Numeric sqrt_05 = sqrt(0.5);
1571  const Index npd = dF_low.npages();
1572 
1573  const Numeric a = A(0, 0), b = A(0, 1), c = A(0, 2), d = A(0, 3), u = A(1, 2),
1574  v = A(1, 3), w = A(2, 3);
1575  const Numeric exp_a = exp(a);
1576  const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
1577  w2 = w * w;
1578 
1579  const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1580 
1581  Numeric Const1;
1582  Const1 = b2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1583  Const1 += c2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1584  Const1 += d2 * (d2 * 0.5 + u2 - v2 - w2);
1585  Const1 += u2 * (u2 * 0.5 + v2 + w2);
1586  Const1 += v2 * (v2 * 0.5 + w2);
1587  Const1 *= 2;
1588  Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
1589  Const1 += w2 * w2;
1590 
1591  if (Const1 > 0.0)
1592  Const1 = sqrt(Const1);
1593  else
1594  Const1 = 0.0;
1595 
1596  if (Const1 == 0 and
1597  Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1598  {
1599  F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1600  for (Index i = 0; i < npd; i++) {
1601  // mult incase the derivative is still non-zero, e.g., as when mag-field is considered 0 but derivative is computed
1602  mult(dF_upp(i, joker, joker), F, dA_upp(i, joker, joker));
1603  mult(dF_low(i, joker, joker), F, dA_low(i, joker, joker));
1604  }
1605  return;
1606  }
1607 
1608  const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1609  const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1610  const Numeric x = sqrt_BpA.real() * sqrt_05;
1611  const Numeric y = sqrt_BmA.imag() * sqrt_05;
1612  const Numeric x2 = x * x;
1613  const Numeric y2 = y * y;
1614  const Numeric cos_y = cos(y);
1615  const Numeric sin_y = sin(y);
1616  const Numeric cosh_x = cosh(x);
1617  const Numeric sinh_x = sinh(x);
1618  const Numeric x2y2 = x2 + y2;
1619  const Numeric inv_x2y2 = 1.0 / x2y2;
1620 
1621  Numeric C0, C1, C2, C3;
1622  Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1623 
1624  // X and Y cannot both be zero
1625  if (x == 0.0) {
1626  inv_y = 1.0 / y;
1627  C0 = 1.0;
1628  C1 = 1.0;
1629  C2 = (1.0 - cos_y) * inv_x2y2;
1630  C3 = (1.0 - sin_y * inv_y) * inv_x2y2;
1631  } else if (y == 0.0) {
1632  inv_x = 1.0 / x;
1633  C0 = 1.0;
1634  C1 = 1.0;
1635  C2 = (cosh_x - 1.0) * inv_x2y2;
1636  C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1637  } else {
1638  inv_x = 1.0 / x;
1639  inv_y = 1.0 / y;
1640 
1641  C0 = (cos_y * x2 + cosh_x * y2) * inv_x2y2;
1642  C1 = (sin_y * x2 * inv_y + sinh_x * y2 * inv_x) * inv_x2y2;
1643  C2 = (cosh_x - cos_y) * inv_x2y2;
1644  C3 = (sinh_x * inv_x - sin_y * inv_y) * inv_x2y2;
1645  }
1646 
1647  eigF(0, 0) = eigF(1, 1) = eigF(2, 2) = eigF(3, 3) = C0;
1648  eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1649  eigF(0, 0) += C0;
1650  eigF(1, 1) += C0;
1651  eigF(2, 2) += C0;
1652  eigF(3, 3) += C0;
1653  eigF *= exp_a;
1654 
1655  if (npd) {
1656  const Numeric inv_x2 = inv_x * inv_x;
1657  const Numeric inv_y2 = inv_y * inv_y;
1658 
1659  for (Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1660  for (Index i = 0; i < npd; i++) {
1661  MatrixView tmp =
1662  upp_or_low ? dF_low(i, joker, joker) : dF_upp(i, joker, joker);
1663  Matrix4x4ViewMap dF = MapToEigen4x4(tmp);
1664  Eigen::Matrix4d dA = MapToEigen4x4(
1665  upp_or_low ? dA_low(i, joker, joker) : dA_upp(i, joker, joker));
1666 
1667  const Numeric da = dA(0, 0), db = dA(0, 1), dc = dA(0, 2),
1668  dd = dA(0, 3), du = dA(1, 2), dv = dA(1, 3),
1669  dw = dA(2, 3);
1670 
1671  dA(0, 0) = dA(1, 1) = dA(2, 2) = dA(3, 3) = 0.0;
1672 
1673  Eigen::Matrix4d dA2;
1674  ;
1675  dA2.noalias() = eigA * dA;
1676  dA2.noalias() += dA * eigA;
1677 
1678  Eigen::Matrix4d dA3;
1679  dA3.noalias() = dA2 * eigA;
1680  dA3.noalias() += eigA2 * dA;
1681 
1682  const Numeric db2 = 2 * db * b, dc2 = 2 * dc * c, dd2 = 2 * dd * d,
1683  du2 = 2 * du * u, dv2 = 2 * dv * v, dw2 = 2 * dw * w;
1684 
1685  const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1686 
1687  Numeric dConst1;
1688  if (Const1 > 0.) {
1689  dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1690  dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1691 
1692  dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1693  dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1694 
1695  dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1696  dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1697 
1698  dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1699  dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1700 
1701  dConst1 += dv2 * (v2 * 0.5 + w2);
1702  dConst1 += v2 * (dv2 * 0.5 + dw2);
1703 
1704  dConst1 += 4 * ((db * d * u * w - db * c * v * w - dc * d * u * v +
1705  b * dd * u * w - b * dc * v * w - c * dd * u * v +
1706  b * d * du * w - b * c * dv * w - c * d * du * v +
1707  b * d * u * dw - b * c * v * dw - c * d * u * dv));
1708  dConst1 += dw2 * w2;
1709  dConst1 /= Const1;
1710  } else
1711  dConst1 = 0.0;
1712 
1713  if (x == 0.0) {
1714  const Numeric dy =
1715  (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1716 
1717  const Numeric dC2 =
1718  -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1719  const Numeric dC3 =
1720  -2 * y * dy * C3 * inv_x2y2 +
1721  (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1722 
1723  dF.noalias() = dC2 * eigA2;
1724  +C2* dA2 + dC3* eigA3 + C3* dA3;
1725  dF *= exp_a;
1726  dF.noalias() += eigF * da;
1727  } else if (y == 0.0) {
1728  const Numeric dx =
1729  (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1730 
1731  const Numeric dC2 =
1732  -2 * x * dx * C2 * inv_x2y2 + dx * sinh_x * inv_x2y2;
1733  const Numeric dC3 =
1734  -2 * x * dx * C3 * inv_x2y2 +
1735  (cosh_x * dx * inv_x - dx * sinh_x * inv_x2) * inv_x2y2;
1736 
1737  dF.noalias() = dC2 * eigA2 + C2 * dA2 + dC3 * eigA3 + C3 * dA3;
1738  dF *= exp_a;
1739  dF.noalias() += eigF * da;
1740  } else {
1741  const Numeric dx =
1742  (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1743  const Numeric dy =
1744  (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1745  const Numeric dy2 = 2 * y * dy;
1746  const Numeric dx2 = 2 * x * dx;
1747  const Numeric dx2dy2 = dx2 + dy2;
1748 
1749  const Numeric dC0 = -dx2dy2 * C0 * inv_x2y2 +
1750  (2 * cos_y * dx * x + 2 * cosh_x * dy * y +
1751  dx * sinh_x * y2 - dy * sin_y * x2) *
1752  inv_x2y2;
1753 
1754  const Numeric dC1 =
1755  -dx2dy2 * C1 * inv_x2y2 +
1756  (cos_y * dy * x2 * inv_y + dx2 * sin_y * inv_y -
1757  dy * sin_y * x2 * inv_y2 - dx * sinh_x * y2 * inv_x2 +
1758  cosh_x * dx * y2 * inv_x + dy2 * sinh_x * inv_x) *
1759  inv_x2y2;
1760 
1761  const Numeric dC2 =
1762  -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1763 
1764  const Numeric dC3 = -dx2dy2 * C3 * inv_x2y2 +
1765  (dy * sin_y * inv_y2 - cos_y * dy * inv_y +
1766  cosh_x * dx * inv_x - dx * sinh_x * inv_x2) *
1767  inv_x2y2;
1768 
1769  dF.noalias() = dC1 * eigA + C1 * dA + dC2 * eigA2 + C2 * dA2 +
1770  dC3 * eigA3 + C3 * dA3;
1771  dF(0, 0) += dC0;
1772  dF(1, 1) += dC0;
1773  dF(2, 2) += dC0;
1774  dF(3, 3) += dC0;
1775  dF *= exp_a;
1776  dF.noalias() += eigF * da;
1777  }
1778  }
1779  }
1780  }
1781 }
1782 
1783 void propmat4x4_to_transmat4x4(MatrixView F,
1784  Tensor3View dF_upp,
1785  Tensor3View dF_low,
1786  ConstMatrixView A,
1787  ConstTensor3View dA_upp,
1788  ConstTensor3View dA_low,
1789  const Index& q) {
1790  const Index n_partials = dA_upp.npages();
1791 
1792  /* Check if A and F are a quadratic and of the same dimension. */
1793  assert(is_size(A, 4, 4));
1794  assert(is_size(F, 4, 4));
1795  assert(n_partials == dF_upp.npages());
1796  assert(n_partials == dF_low.npages());
1797  assert(n_partials == dA_low.npages());
1798  for (Index ii = 0; ii < n_partials; ii++) {
1799  assert(is_size(dA_upp(ii, joker, joker), 4, 4));
1800  assert(is_size(dA_low(ii, joker, joker), 4, 4));
1801  assert(is_size(dF_upp(ii, joker, joker), 4, 4));
1802  assert(is_size(dF_low(ii, joker, joker), 4, 4));
1803  }
1804 
1805  // Set constants and Numerics
1806  const Numeric A_norm_inf = norm_inf(A);
1807  const Numeric e = 1.0 + floor(1.0 / log(2.0) * log(A_norm_inf));
1808  const Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
1809  const Numeric inv_pow2 = 1. / pow(2, r);
1810  Numeric c = 0.5;
1811 
1812  // Create M, dM, X and Y
1813  Eigen::Matrix4d M = MapToEigen4x4(A);
1814  M *= inv_pow2;
1815  Eigen::Matrix4d X = M;
1816  Eigen::Matrix4d cX = c * X;
1817  Eigen::Matrix4d D = Eigen::Matrix4d::Identity();
1818  Matrix4x4ViewMap eigF = MapToEigen4x4(F);
1819  eigF.setIdentity();
1820  eigF.noalias() += cX;
1821  D.noalias() -= cX;
1822 
1823  Array<Eigen::Matrix4d> dMu(n_partials), dMl(n_partials), Yu(n_partials),
1824  Yl(n_partials);
1825  Array<Eigen::Matrix4d> cYu(n_partials), cYl(n_partials), dDu(n_partials),
1826  dDl(n_partials);
1827 
1828  Array<Matrix4x4ViewMap> dFu, dFl;
1829  dFu.reserve(n_partials);
1830  dFl.reserve(n_partials);
1831 
1832  for (Index i = 0; i < n_partials; i++) {
1833  dMu[i].noalias() = MapToEigen4x4(dA_upp(i, joker, joker)) * inv_pow2;
1834  dMl[i].noalias() = MapToEigen4x4(dA_low(i, joker, joker)) * inv_pow2;
1835 
1836  Yu[i] = dMu[i];
1837  Yl[i] = dMl[i];
1838 
1839  cYu[i].noalias() = c * dMu[i];
1840  cYl[i].noalias() = c * dMl[i];
1841 
1842  dFu.push_back(MapToEigen4x4(dF_upp, i));
1843  dFl.push_back(MapToEigen4x4(dF_low, i));
1844  dFu[i] = cYu[i];
1845  dFl[i] = cYl[i];
1846 
1847  dDu[i] = -cYu[i];
1848  dDl[i] = -cYl[i];
1849  }
1850 
1851  bool p = true;
1852  for (Index k = 2; k <= q; k++) {
1853  c *= (Numeric)(q - k + 1) / (Numeric)(k * (2 * q - k + 1));
1854  for (Index i = 0; i < n_partials; i++) {
1855  Yu[i] = dMu[i] * X + M * Yu[i];
1856  Yl[i] = dMl[i] * X + M * Yl[i];
1857 
1858  cYu[i].noalias() = c * Yu[i];
1859  cYl[i].noalias() = c * Yl[i];
1860 
1861  dFu[i].noalias() += cYu[i];
1862  dFl[i].noalias() += cYl[i];
1863  }
1864 
1865  X = M * X;
1866  cX.noalias() = c * X;
1867  eigF.noalias() += cX;
1868 
1869  if (p) {
1870  D.noalias() += cX;
1871 
1872  for (Index i = 0; i < n_partials; i++) {
1873  dDu[i].noalias() += cYu[i];
1874  dDl[i].noalias() += cYl[i];
1875  }
1876  } else {
1877  D.noalias() -= cX;
1878 
1879  for (Index i = 0; i < n_partials; i++) {
1880  dDu[i].noalias() -= cYu[i];
1881  dDl[i].noalias() -= cYl[i];
1882  }
1883  }
1884  p = not p;
1885  }
1886 
1887  const Eigen::Matrix4d invD = D.inverse();
1888 
1889  eigF = invD * eigF;
1890  for (Index i = 0; i < n_partials; i++) {
1891  dFu[i] = invD * (dFu[i] - dDu[i] * eigF);
1892  dFl[i] = invD * (dFl[i] - dDl[i] * eigF);
1893  }
1894 
1895  for (Index k = 1; k <= r; k++) {
1896  for (Index i = 0; i < n_partials; i++) {
1897  dFu[i] = dFu[i] * eigF + eigF * dFu[i];
1898  dFl[i] = dFl[i] * eigF + eigF * dFl[i];
1899  }
1900  eigF = eigF * eigF;
1901  }
1902 }
1903 
1905 
1926 void matrix_exp_dmatrix_exp(MatrixView F,
1927  Tensor3View dF,
1928  ConstMatrixView A,
1929  ConstTensor3View dA,
1930  const Index& q) {
1931  const Index n_partials = dA.npages();
1932 
1933  const Index n = A.ncols();
1934 
1935  /* Check if A and F are a quadratic and of the same dimension. */
1936  assert(is_size(A, n, n));
1937  assert(is_size(F, n, n));
1938  assert(n_partials == dF.npages());
1939  for (Index ii = 0; ii < n_partials; ii++) {
1940  assert(is_size(dA(ii, joker, joker), n, n));
1941  assert(is_size(dF(ii, joker, joker), n, n));
1942  }
1943 
1944  // This sets up some cnstants
1945  Numeric A_norm_inf, e;
1946  A_norm_inf = norm_inf(A);
1947  e = 1. + floor(1. / log(2.) * log(A_norm_inf));
1948  Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
1949  Numeric pow2rm1 = 1. / pow(2, r);
1950  Numeric c = 0.5;
1951 
1952  // For non-derivatives
1953  Matrix M = A, X(n, n), cX(n, n), D(n, n);
1954  M *= pow2rm1;
1955  X = M;
1956  cX = X;
1957  cX *= c;
1958  id_mat(F);
1959  F += cX;
1960  id_mat(D);
1961  D -= cX;
1962 
1963  // For derivatives
1964  Tensor3 dM(n_partials, n, n), Y(n_partials, n, n), cY(n_partials, n, n),
1965  dD(n_partials, n, n);
1966  for (Index ii = 0; ii < n_partials; ii++) {
1967  for (Index jj = 0; jj < n; jj++) {
1968  for (Index kk = 0; kk < n; kk++) {
1969  dM(ii, jj, kk) = dA(ii, jj, kk) * pow2rm1;
1970 
1971  Y(ii, jj, kk) = dM(ii, jj, kk);
1972 
1973  cY(ii, jj, kk) = c * Y(ii, jj, kk);
1974 
1975  dF(ii, jj, kk) = cY(ii, jj, kk);
1976 
1977  dD(ii, jj, kk) = -cY(ii, jj, kk);
1978  }
1979  }
1980  }
1981 
1982  // NOTE: MATLAB paper sets q = 6 but we allow other numbers
1983 
1984  Matrix tmp1(n, n), tmp2(n, n);
1985 
1986  for (Index k = 2; k <= q; k++) {
1987  c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
1988 
1989  // For partials
1990  for (Index ii = 0; ii < n_partials; ii++) {
1991  // Y = dM*X + M*Y
1992  mult(tmp1, dM(ii, joker, joker), X);
1993  mult(tmp2, M, Y(ii, joker, joker));
1994 
1995  for (Index jj = 0; jj < n; jj++)
1996  for (Index kk = 0; kk < n; kk++) {
1997  Y(ii, jj, kk) = tmp1(jj, kk) + tmp2(jj, kk);
1998 
1999  // cY = c*Y
2000  cY(ii, jj, kk) = c * Y(ii, jj, kk);
2001 
2002  // dF = dF + cY
2003  dF(ii, jj, kk) += cY(ii, jj, kk);
2004 
2005  if (k % 2 == 0) //For even numbers, add.
2006  {
2007  // dD = dD + cY
2008  dD(ii, jj, kk) += cY(ii, jj, kk);
2009  } else {
2010  // dD = dD - cY
2011  dD(ii, jj, kk) -= cY(ii, jj, kk);
2012  }
2013  }
2014  }
2015 
2016  //For full derivative (Note X in partials above)
2017  // X=M*X
2018  mult(tmp1, M, X);
2019  X = tmp1;
2020 
2021  //cX = c*X
2022  cX = X;
2023  cX *= c;
2024 
2025  // F = F + cX
2026  F += cX;
2027 
2028  if (k % 2 == 0) //For even numbers, D = D + cX
2029  D += cX;
2030  else //For odd numbers, D = D - cX
2031  D -= cX;
2032  }
2033 
2034  // D^-1
2035  inv(tmp1, D);
2036 
2037  // F = D\F, or D^-1*F
2038  mult(tmp2, tmp1, F);
2039  F = tmp2;
2040 
2041  // For partials
2042  for (Index ii = 0; ii < n_partials; ii++) {
2043  //dF = D \ (dF - dF*F), or D^-1 * (dF - dF*F)
2044  mult(tmp2, dD(ii, joker, joker), F); // dF * F
2045  dF(ii, joker, joker) -= tmp2; // dF - dF * F
2046  mult(tmp2, tmp1, dF(ii, joker, joker));
2047  dF(ii, joker, joker) = tmp2;
2048  }
2049 
2050  for (Index k = 1; k <= r; k++) {
2051  for (Index ii = 0; ii < n_partials; ii++) {
2052  // dF=F*dF+dF*F
2053  mult(tmp1, F, dF(ii, joker, joker)); //F*dF
2054  mult(tmp2, dF(ii, joker, joker), F); //dF*F
2055  dF(ii, joker, joker) = tmp1;
2056  dF(ii, joker, joker) += tmp2;
2057  }
2058 
2059  // F=F*F
2060  mult(tmp1, F, F);
2061  F = tmp1;
2062  }
2063 }
2064 
2066 
2087 void matrix_exp_dmatrix_exp(MatrixView F,
2088  MatrixView dF,
2089  ConstMatrixView A,
2090  ConstMatrixView dA,
2091  const Index& q) {
2092  const Index n = A.ncols();
2093 
2094  /* Check if A and F are a quadratic and of the same dimension. */
2095  assert(is_size(A, n, n));
2096  assert(is_size(F, n, n));
2097  assert(is_size(dA, n, n));
2098  assert(is_size(dF, n, n));
2099 
2100  // This is the definition of how to scale
2101  Numeric A_norm_inf, e;
2102  A_norm_inf = norm_inf(A);
2103  e = 1. + floor(1. / log(2.) * log(A_norm_inf));
2104  Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
2105  Numeric pow2rm1 = 1. / pow(2, r);
2106 
2107  // A and dA are scaled
2108  Matrix M = A, dM = dA;
2109  M *= pow2rm1;
2110  dM *= pow2rm1;
2111 
2112  // These variables will hold the multiplication calculations
2113  Matrix X(n, n), Y(n, n);
2114  X = M;
2115  Y = dM;
2116 
2117  // cX = c*M
2118  // cY = c*dM
2119  Matrix cX = X, cY = Y;
2120  Numeric c = 0.5;
2121  cX *= c;
2122  cY *= c;
2123 
2124  Matrix D(n, n), dD(n, n);
2125  // F = I + c*M
2126  id_mat(F);
2127  F += cX;
2128 
2129  // dF = c*dM;
2130  dF = cY;
2131 
2132  //D = I -c*M
2133  id_mat(D);
2134  D -= cX;
2135 
2136  // dD = -c*dM
2137  dD = cY;
2138  dD *= -1.;
2139 
2140  // NOTE: MATLAB paper sets q = 6 but we allow other numbers
2141 
2142  Matrix tmp1(n, n), tmp2(n, n);
2143 
2144  for (Index k = 2; k <= q; k++) {
2145  c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
2146 
2147  // Y = dM*X + M*Y
2148  mult(tmp1, dM, X);
2149  mult(tmp2, M, Y);
2150  Y = tmp1;
2151  Y += tmp2;
2152 
2153  // X=M*X
2154  mult(tmp1, M, X);
2155  X = tmp1;
2156 
2157  //cX = c*X
2158  cX = X;
2159  cX *= c;
2160 
2161  // cY = c*Y
2162  cY = Y;
2163  cY *= c;
2164 
2165  // F = F + cX
2166  F += cX;
2167 
2168  // dF = dF + cY
2169  dF += cY;
2170 
2171  if (k % 2 == 0) //For even numbers, add.
2172  {
2173  // D = D + cX
2174  D += cX;
2175 
2176  // dD = dD + cY
2177  dD += cY;
2178  } else //For odd numbers, subtract
2179  {
2180  // D = D - cX
2181  D -= cX;
2182 
2183  // dD = dD - cY
2184  dD -= cY;
2185  }
2186  }
2187 
2188  // D^-1
2189  inv(tmp1, D);
2190 
2191  // F = D\F, or D^-1*F
2192  mult(tmp2, tmp1, F);
2193  F = tmp2;
2194 
2195  //dF = D \ (dF - dF*F), or D^-1 * (dF - dF*F)
2196  mult(tmp2, dD, F); // dF * F
2197  dF -= tmp2; // dF - dF * F
2198  mult(tmp2, tmp1, dF);
2199  dF = tmp2;
2200 
2201  for (Index k = 1; k <= r; k++) {
2202  // dF=F*dF+dF*F
2203  mult(tmp1, F, dF); //F*dF
2204  mult(tmp2, dF, F); //dF*F
2205  dF = tmp1;
2206  dF += tmp2;
2207 
2208  // F=F*F
2209  mult(tmp1, F, F);
2210  F = tmp1;
2211  }
2212 }
2213 
2215 
2223 Numeric norm_inf(ConstMatrixView A) {
2224  Numeric norm_inf = 0;
2225 
2226  for (Index j = 0; j < A.nrows(); j++) {
2227  Numeric row_sum = 0;
2228  //Calculate the row sum for all rows
2229  for (Index i = 0; i < A.ncols(); i++) row_sum += abs(A(i, j));
2230  //Pick out the row with the highest row sum
2231  if (norm_inf < row_sum) norm_inf = row_sum;
2232  }
2233  return norm_inf;
2234 }
2235 
2237 
2240 void id_mat(MatrixView I) {
2241  const Index n = I.ncols();
2242  assert(n == I.nrows());
2243 
2244  I = 0;
2245  for (Index i = 0; i < n; i++) I(i, i) = 1.;
2246 }
2247 
2256 Numeric det(ConstMatrixView A) {
2257  const Index dim = A.nrows();
2258  assert(dim == A.ncols());
2259 
2260  if (dim == 3)
2261  return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +
2262  A(0, 2) * A(1, 0) * A(2, 1) - A(0, 2) * A(1, 1) * A(2, 0) -
2263  A(0, 1) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 2) * A(2, 1);
2264  else if (dim == 2)
2265  return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
2266  else if (dim == 1)
2267  return A(0, 0);
2268 
2269  Numeric ret_val = 0.;
2270 
2271  for (Index j = 0; j < dim; j++) {
2272  Matrix temp(dim - 1, dim - 1);
2273  for (Index I = 1; I < dim; I++)
2274  for (Index J = 0; J < dim; J++) {
2275  if (J < j)
2276  temp(I - 1, J) = A(I, J);
2277  else if (J > j)
2278  temp(I - 1, J - 1) = A(I, J);
2279  }
2280 
2281  Numeric tempNum = det(temp);
2282 
2283  ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
2284  }
2285  return ret_val;
2286 }
2287 
2302 void linreg(Vector& p, ConstVectorView x, ConstVectorView y) {
2303  const Index n = x.nelem();
2304 
2305  assert(y.nelem() == n);
2306 
2307  p.resize(2);
2308 
2309  // Basic algorithm found at e.g.
2310  // http://en.wikipedia.org/wiki/Simple_linear_regression
2311  // The basic algorithm is as follows:
2312  /*
2313  Numeric s1=0, s2=0, s3=0, s4=0;
2314  for( Index i=0; i<n; i++ )
2315  {
2316  s1 += x[i] * y[i];
2317  s2 += x[i];
2318  s3 += y[i];
2319  s4 += x[i] * x[i];
2320  }
2321 
2322  p[1] = ( s1 - (s2*s3)/n ) / ( s4 - s2*s2/n );
2323  p[0] = s3/n - p[1]*s2/n;
2324  */
2325 
2326  // A version abit more numerical stable:
2327  // Mean value of x is removed before the fit: x' = (x-mean(x))
2328  // This corresponds to that s2 in version above becomes 0
2329  // y = a + b*x'
2330  // p[1] = b
2331  // p[0] = a - p[1]*mean(x)
2332 
2333  Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
2334 
2335  for (Index i = 0; i < n; i++) {
2336  xm += x[i] / Numeric(n);
2337  }
2338 
2339  for (Index i = 0; i < n; i++) {
2340  const Numeric xv = x[i] - xm;
2341  s1 += xv * y[i];
2342  s3 += y[i];
2343  s4 += xv * xv;
2344  }
2345 
2346  p[1] = s1 / s4;
2347  p[0] = s3 / Numeric(n) - p[1] * xm;
2348 }
2349 
2350 Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept {
2351  // Size of the problem
2352  const Index n = x.nelem();
2353  Matrix AT, ATA(n, n);
2354  Vector ATy(n);
2355 
2356  // Solver
2357  AT = transpose(A);
2358  mult(ATA, AT, A);
2359  mult(ATy, AT, y);
2360  solve(x, ATA, ATy);
2361 
2362  // Residual
2363  if (residual) {
2364  Vector r(n);
2365  mult(r, ATA, x);
2366  r -= ATy;
2367  return r * r;
2368  } else {
2369  return 0;
2370  }
2371 }
2372 
2373 Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eig(const Eigen::Ref<Eigen::MatrixXcd> A)
2374 {
2375  Eigen::ComplexEigenSolver<Eigen::MatrixXcd> ces;
2376  return ces.compute(A);
2377 }
matrix_exp2
void matrix_exp2(MatrixView F, ConstMatrixView A)
Definition: lin_alg.cc:458
Matrix
The Matrix class.
Definition: matpackI.h:1193
special_matrix_exp_and_dmatrix_exp_dx_for_rt
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.
Definition: lin_alg.cc:559
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
ComplexVectorView::get_c_array
const Complex * get_c_array() const
Conversion to plain C-array.
Definition: complex.cc:421
lapack.h
Interface for the LAPACK library.
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
linreg
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:2302
joker
const Joker joker
temp
#define temp
Definition: legacy_continua.cc:20951
matrix_exp2_4x4
void matrix_exp2_4x4(MatrixView arts_f, ConstMatrixView A)
Definition: lin_alg.cc:521
ConstComplexVectorView::nelem
Index nelem() const
Definition: complex.h:280
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:34
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ComplexVector
The ComplexVector class.
Definition: complex.h:573
Matrix4x4ViewMap
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:113
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
swap
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Definition: complex.cc:731
ConstMatrixView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1081
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
DEBUG_ONLY
#define DEBUG_ONLY(...)
Definition: debug.h:36
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
ConstVectorView::mrange
Range mrange
The range of mdata that is actually used.
Definition: matpackI.h:594
MapToEigen
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:750
ComplexMatrixView::get_c_array
const Complex * get_c_array() const
Conversion to plain C-array.
Definition: complex.cc:1112
propmat4x4_to_transmat4x4
void propmat4x4_to_transmat4x4(MatrixView F, Tensor3View dF_upp, Tensor3View dF_low, ConstMatrixView A, ConstTensor3View dA_upp, ConstTensor3View dA_low, const Index &q)
Definition: lin_alg.cc:1783
b2
#define b2
Definition: complex.h:59
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
array.h
This file contains the definition of Array.
matrix_exp
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:391
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
ConstMatrixView::mcr
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1079
MapToEigen4x4
Matrix4x4ViewMap MapToEigen4x4(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:745
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
lapack::dgeev_
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)
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
matrix_exp_dmatrix_exp
void matrix_exp_dmatrix_exp(MatrixView F, Tensor3View dF, ConstMatrixView A, ConstTensor3View dA, const Index &q)
General exponential of a Matrix with their derivatives.
Definition: lin_alg.cc:1926
ConstComplexMatrixView::nrows
Index nrows() const
Definition: complex.h:630
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
solve
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Definition: lin_alg.cc:138
matrix_exp_4x4
void matrix_exp_4x4(MatrixView arts_f, ConstMatrixView A, const Index &q)
Definition: lin_alg.cc:481
dw
constexpr Complex dw(Complex z, Complex w) noexcept
The Faddeeva function partial derivative.
Definition: linefunctions.cc:45
ConstComplexMatrixView::get_column_extent
Index get_column_extent() const
Get the extent of the underlying data.
Definition: complex.h:665
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
lapack::zgetri_
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
ConstComplexMatrixView::ncols
Index ncols() const
Definition: complex.h:631
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
diagonalize
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:245
lapack::dgetri_
void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info)
Matrix inversion.
lapack::dgetrf_
void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info)
LU decomposition.
ComplexMatrix
The ComplexMatrix class.
Definition: complex.h:858
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
lin_alg.h
Linear algebra functions.
matpackIII.h
ConstVectorView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:596
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
Range::get_extent
constexpr Index get_extent() const
Returns the extent of the range.
Definition: matpackI.h:329
oem::Identity
invlib::MatrixIdentity< Matrix > Identity
Definition: oem.h:39
N
#define N
Definition: rng.cc:164
Range::get_stride
constexpr Index get_stride() const
Returns the stride of the range.
Definition: matpackI.h:331
lapack::dgetrs_
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.
ComplexMatrixView
The ComplexMatrixView class.
Definition: complex.h:731
q
#define q
Definition: legacy_continua.cc:21712
MatrixViewMap
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:110
logic.h
Header file for logic.cc.
lapack::zgeev_
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)
lapack::zgetrf_
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
dx
#define dx
Definition: legacy_continua.cc:22180
ConstComplexMatrixView
A constant view of a ComplexMatrix.
Definition: complex.h:618
det
Numeric det(ConstMatrixView A)
Definition: lin_alg.cc:2256
eig
Eigen::ComplexEigenSolver< Eigen::MatrixXcd > eig(const Eigen::Ref< Eigen::MatrixXcd > A)
Return the Eigen decomposition of the eigen matrix.
Definition: lin_alg.cc:2373
cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(MatrixView F, ConstMatrixView A)
Definition: lin_alg.cc:789
lsf
Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept
Least squares fitting by solving x for known A and y.
Definition: lin_alg.cc:2350
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: complex.cc:1509
ludcmp
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:56
M
#define M
Definition: rng.cc:165
cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen
void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(MatrixView F, ConstMatrixView A)
Definition: lin_alg.cc:963
Vector
The Vector class.
Definition: matpackI.h:860
arts_omp.h
Header file for helper functions for OpenMP.
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
Zeeman::dM
constexpr Index dM(Polarization type) noexcept
Gives the change of M given a polarization type.
Definition: zeemandata.h:51
x2
#define x2
Definition: lineshapemodel.h:415
ComplexVectorView
The ComplexVectorView class.
Definition: complex.h:367
norm_inf
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
Definition: lin_alg.cc:2223
VectorView::get_c_array
const Numeric * get_c_array() const
Conversion to plain C-array, const-version.
Definition: matpackI.cc:272
lubacksub
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:91
arts.h
The global header file for ARTS.