ARTS 2.5.0 (git: 9ee3ac6c)
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
56void ludcmp(Matrix& LU, ArrayOfIndex& indx, ConstMatrixView A) {
57 // Assert that A is quadratic.
58 const Index n = A.nrows();
59 ARTS_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
91void 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 ARTS_ASSERT(is_size(LU, n, n));
103 ARTS_ASSERT(is_size(b, n));
104 ARTS_ASSERT(is_size(indx, n));
105 ARTS_ASSERT(column_stride == 1);
106 ARTS_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
138void solve(VectorView x, ConstMatrixView A, ConstVectorView b) {
139 Index n = A.ncols();
140
141 // Check dimensions of the system.
142 ARTS_ASSERT(n == A.nrows());
143 ARTS_ASSERT(n == x.nelem());
144 ARTS_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
167void inv(MatrixView Ainv, ConstMatrixView A) {
168 // A must be a square matrix.
169 ARTS_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 ARTS_USER_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 ARTS_USER_ERROR_IF (info not_eq 0, "Error inverting matrix: Matrix not of full rank.");
199 Ainv = LU;
200}
201
202void inv(ComplexMatrixView Ainv, const ConstComplexMatrixView A) {
203 // A must be a square matrix.
204 ARTS_ASSERT(A.ncols() == A.nrows());
205
206 Index n = A.ncols();
207
208 // Workdata
209 Ainv = A;
210 int n_int = int(n), lwork = n_int, info;
211 std::vector<int> ipiv(n);
212 ComplexVector work(lwork);
213
214 // Compute LU decomposition using LAPACK dgetrf_.
215 lapack::zgetrf_(&n_int, &n_int, Ainv.get_c_array(), &n_int, ipiv.data(), &info);
216 lapack::zgetri_(&n_int, Ainv.get_c_array(), &n_int, ipiv.data(), work.get_c_array(), &lwork, &info);
217
218 // Check for success.
219 ARTS_USER_ERROR_IF (info not_eq 0,
220 "Error inverting matrix: Matrix not of full rank.");
221}
222
224
241void diagonalize(MatrixView P,
242 VectorView WR,
243 VectorView WI,
244 ConstMatrixView A) {
245 Index n = A.ncols();
246
247 // A must be a square matrix.
248 ARTS_ASSERT(n == A.nrows());
249 ARTS_ASSERT(n == WR.nelem());
250 ARTS_ASSERT(n == WI.nelem());
251 ARTS_ASSERT(n == P.nrows());
252 ARTS_ASSERT(n == P.ncols());
253
254 Matrix A_tmp = A;
255 Matrix P2 = P;
256 Vector WR2 = WR;
257 Vector WI2 = WI;
258
259 // Integers
260 int LDA, LDA_L, LDA_R, n_int, info;
261 n_int = (int)n;
262 LDA = LDA_L = LDA_R = (int)A.mcr.get_extent();
263
264 // We want to calculate RP not LP
265 char l_eig = 'N', r_eig = 'V';
266
267 // Work matrix
268 int lwork = 2 * n_int + n_int * n_int;
269 double *work, *rwork;
270 try {
271 rwork = new double[2 * n_int];
272 work = new double[lwork];
273 } catch (std::bad_alloc& ba) {
274 ARTS_USER_ERROR (
275 "Error diagonalizing: Could not allocate workspace memory.");
276 }
277
278 // Memory references
279 double* adata = A_tmp.mdata;
280 double* rpdata = P2.mdata;
281 double* lpdata = new double[0]; //To not confuse the compiler
282 double* wrdata = WR2.mdata;
283 double* widata = WI2.mdata;
284
285 // Main calculations. Note that errors in the output is ignored
286 lapack::dgeev_(&l_eig,
287 &r_eig,
288 &n_int,
289 adata,
290 &LDA,
291 wrdata,
292 widata,
293 lpdata,
294 &LDA_L,
295 rpdata,
296 &LDA_R,
297 work,
298 &lwork,
299 rwork,
300 &info);
301
302 // Free memory. Can these be sent in to speed things up?
303 delete[] work;
304 delete[] rwork;
305 delete[] lpdata;
306
307 // Re-order. This can be done better?
308 for (Index i = 0; i < n; i++)
309 for (Index j = 0; j < n; j++) P(j, i) = P2(i, j);
310 WI = WI2;
311 WR = WR2;
312}
313
315
326void diagonalize(ComplexMatrixView P,
327 ComplexVectorView W,
328 const ConstComplexMatrixView A) {
329 Index n = A.ncols();
330
331 // A must be a square matrix.
332 ARTS_ASSERT(n == A.nrows());
333 ARTS_ASSERT(n == W.nelem());
334 ARTS_ASSERT(n == P.nrows());
335 ARTS_ASSERT(n == P.ncols());
336
337 ComplexMatrix A_tmp = A;
338
339 // Integers
340 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;
341
342 // We want to calculate RP not LP
343 char l_eig = 'N', r_eig = 'V';
344
345 // Work matrix
346 int lwork = 2 * n_int + n_int * n_int;
347 ComplexVector work(lwork);
348 ComplexVector lpdata(0);
349 Vector rwork(2 * n_int);
350
351 // Main calculations. Note that errors in the output is ignored
352 lapack::zgeev_(&l_eig,
353 &r_eig,
354 &n_int,
355 A_tmp.get_c_array(),
356 &LDA,
357 W.get_c_array(),
358 lpdata.get_c_array(),
359 &LDA_L,
360 P.get_c_array(),
361 &LDA_R,
362 work.get_c_array(),
363 &lwork,
364 rwork.get_c_array(),
365 &info);
366
367 for (Index i = 0; i < n; i++)
368 for (Index j = 0; j <= i; j++)
369 std::swap(P(j, i), P(i, j));
370}
371
373
387void matrix_exp(MatrixView F, ConstMatrixView A, const Index& q) {
388 const Index n = A.ncols();
389
390 /* Check if A and F are a quadratic and of the same dimension. */
391 ARTS_ASSERT(is_size(A, n, n));
392 ARTS_ASSERT(is_size(F, n, n));
393
394 Numeric A_norm_inf, c;
395 Numeric 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.);
398
399 A_norm_inf = norm_inf(A);
400
401 // This formular is derived in the book by Golub and Van Loan.
402 j = 1 + floor(1. / log(2.) * log(A_norm_inf));
403
404 if (j < 0) j = 0.;
405 Index j_index = (Index)(j);
406
407 // Scale matrix
408 F = A;
409 F /= pow(2, j);
410
411 /* The higher q the more accurate is the computation,
412 see user guide for accuracy */
413 // Index q = 8;
414 Numeric q_n = (Numeric)(q);
415 id_mat(D);
416 id_mat(N);
417 id_mat(X);
418 c = 1.;
419
420 for (Index k = 0; k < q; k++) {
421 Numeric k_n = (Numeric)(k + 1);
422 c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
423 mult(B, F, X); // X = F * X
424 X = B;
425 cX = X;
426 cX *= c; // cX = X*c
427 N += cX; // N = N + X*c
428 cX *= pow(-1, k_n); // cX = (-1)^k*c*X
429 D += cX; // D = D + (-1)^k*c*X
430 }
431
432 /*Solve the equation system DF=N for F using LU decomposition,
433 use the backsubstitution routine for columns of N*/
434
435 /* Now use X for the LU decomposition matrix of D.*/
436 ArrayOfIndex indx(n);
437
438 ludcmp(X, indx, D);
439
440 for (Index i = 0; i < n; i++) {
441 N_col_vec = N(joker, i); // extract column vectors of N
442 lubacksub(F_col_vec, X, N_col_vec, indx);
443 F(joker, i) = F_col_vec; // construct F matrix from column vectors
444 }
445
446 /* The square of F gives the result. */
447 for (Index k = 0; k < j_index; k++) {
448 mult(B, F, F); // F = F^2
449 F = B;
450 }
451}
452
453//Matrix exponent with decomposition
454void matrix_exp2(MatrixView F, ConstMatrixView A) {
455 const Index n = F.nrows();
456
457 ARTS_ASSERT(is_size(F, n, n));
458 ARTS_ASSERT(is_size(A, n, n));
459
460 Matrix P(n, n), invP(n, n);
461 Vector WR(n), WI(n);
462
463 diagonalize(P, WR, WI, A);
464 inv(invP, P);
465
466 for (Index k = 0; k < n; k++) {
467 ARTS_USER_ERROR_IF (WI[k] != 0,
468 "We require real eigenvalues for your chosen method.");
469 P(joker, k) *= exp(WR[k]);
470 }
471
472 //mult(tmp,Wdiag,invP);
473 mult(F, P, invP);
474}
475
476void matrix_exp_4x4(MatrixView arts_f, ConstMatrixView A, const Index& q) {
477 /* Check if A and F are a quadratic and of the same dimension. */
478 ARTS_ASSERT(is_size(A, 4, 4));
479 ARTS_ASSERT(is_size(arts_f, 4, 4));
480
481 // Set constants and Numerics
482 const Numeric A_norm_inf = norm_inf(A);
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;
485 const Numeric inv_pow2 = 1. / pow(2, r);
486 Numeric c = 0.5;
487
488 const Eigen::Matrix4d M = MapToEigen4x4(A) * inv_pow2;
489 Eigen::Matrix4d X = M;
490 Eigen::Matrix4d cX = c * X;
491 Matrix4x4ViewMap F = MapToEigen4x4(arts_f);
492 Eigen::Matrix4d D = Eigen::Matrix4d::Identity();
493 F.setIdentity();
494 F.noalias() += cX;
495 D.noalias() -= cX;
496
497 bool p = true;
498 for (Index k = 2; k <= q; k++) {
499 c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
500 X = M * X;
501 cX.noalias() = c * X;
502 F.noalias() += cX;
503 if (p)
504 D.noalias() += cX;
505 else
506 D.noalias() -= cX;
507 p = not p;
508 }
509
510 F = D.inverse() * F;
511
512 for (Index k = 1; k <= r; k++) F *= F;
513}
514
515//Matrix exponent with decomposition
516void matrix_exp2_4x4(MatrixView arts_f, ConstMatrixView A) {
517 ARTS_ASSERT(is_size(arts_f, 4, 4));
518 ARTS_ASSERT(is_size(A, 4, 4));
519
520 Matrix4x4ViewMap F = MapToEigen4x4(arts_f);
521 Eigen::EigenSolver<Eigen::Matrix4d> es;
522 es.compute(MapToEigen4x4(A));
523
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();
528 res *= U.inverse();
529 F = res.real();
530}
531
533
554void special_matrix_exp_and_dmatrix_exp_dx_for_rt(MatrixView F,
555 Tensor3View dF_upp,
556 Tensor3View dF_low,
557 ConstMatrixView A,
558 ConstTensor3View dA_upp,
559 ConstTensor3View dA_low,
560 const Index& q) {
561 const Index n_partials = dA_upp.npages();
562
563 const Index n = A.ncols();
564
565 /* Check if A and F are a quadratic and of the same dimension. */
566 ARTS_ASSERT(is_size(A, n, n));
567 ARTS_ASSERT(is_size(F, n, n));
568 ARTS_ASSERT(n_partials == dF_upp.npages());
569 ARTS_ASSERT(n_partials == dF_low.npages());
570 ARTS_ASSERT(n_partials == dA_low.npages());
571 for (Index ii = 0; ii < n_partials; ii++) {
572 ARTS_ASSERT(is_size(dA_upp(ii, joker, joker), n, n));
573 ARTS_ASSERT(is_size(dA_low(ii, joker, joker), n, n));
574 ARTS_ASSERT(is_size(dF_upp(ii, joker, joker), n, n));
575 ARTS_ASSERT(is_size(dF_low(ii, joker, joker), n, n));
576 }
577
578 // This sets up some constants
579 Numeric A_norm_inf, e;
580 A_norm_inf = norm_inf(A);
581 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
582 Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
583 Numeric pow2rm1 = 1. / pow(2, r);
584 Numeric c = 0.5;
585
586 // For non-partials
587 Matrix M = A, X(n, n), cX(n, n), D(n, n);
588 M *= pow2rm1;
589 X = M;
590 cX = X;
591 cX *= c;
592 id_mat(F);
593 F += cX;
594 id_mat(D);
595 D -= cX;
596
597 // For partials in parallel
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++) {
604 for (Index kk = 0; kk < n; kk++) {
605 dM_upp(ii, jj, kk) = dA_upp(ii, jj, kk) * pow2rm1;
606 dM_low(ii, jj, kk) = dA_low(ii, jj, kk) * pow2rm1;
607
608 Y_upp(ii, jj, kk) = dM_upp(ii, jj, kk);
609 Y_low(ii, jj, kk) = dM_low(ii, jj, kk);
610
611 cY_upp(ii, jj, kk) = c * Y_upp(ii, jj, kk);
612 cY_low(ii, jj, kk) = c * Y_low(ii, jj, kk);
613
614 dF_upp(ii, jj, kk) = cY_upp(ii, jj, kk);
615 dF_low(ii, jj, kk) = cY_low(ii, jj, kk);
616
617 dD_upp(ii, jj, kk) = -cY_upp(ii, jj, kk);
618 dD_low(ii, jj, kk) = -cY_low(ii, jj, kk);
619 }
620 }
621 }
622
623 // NOTE: MATLAB paper sets q = 6 but we allow other numbers
624 Matrix tmp1(n, n), tmp2(n, n);
625
626 for (Index k = 2; k <= q; k++) {
627 c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
628
629 // For partials in parallel
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);
632
633 // Y = dM*X + M*Y
634 mult(tmp_upp1, dM_upp(ii, joker, joker), X);
635 mult(tmp_upp2, M, Y_upp(ii, joker, joker));
636 mult(tmp_low1, dM_low(ii, joker, joker), X);
637 mult(tmp_low2, M, Y_low(ii, joker, joker));
638
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);
643
644 // cY = c*Y
645 cY_upp(ii, jj, kk) = c * Y_upp(ii, jj, kk);
646 cY_low(ii, jj, kk) = c * Y_low(ii, jj, kk);
647
648 // dF = dF + cY
649 dF_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
650 dF_low(ii, jj, kk) += cY_low(ii, jj, kk);
651
652 if (k % 2 == 0) //For even numbers, add.
653 {
654 // dD = dD + cY
655 dD_upp(ii, jj, kk) += cY_upp(ii, jj, kk);
656 dD_low(ii, jj, kk) += cY_low(ii, jj, kk);
657 } else {
658 // dD = dD - cY
659 dD_upp(ii, jj, kk) -= cY_upp(ii, jj, kk);
660 dD_low(ii, jj, kk) -= cY_low(ii, jj, kk);
661 }
662 }
663 }
664 }
665
666 //For full derivative (Note X in partials above)
667 // X=M*X
668 mult(tmp1, M, X);
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);
674
675 if (k % 2 == 0) //For even numbers, D = D + cX
676 D(jj, kk) += cX(jj, kk);
677 else //For odd numbers, D = D - cX
678 D(jj, kk) -= cX(jj, kk);
679 }
680 }
681 }
682
683 // D^-1
684 inv(tmp1, D);
685
686 // F = D\F, or D^-1*F
687 mult(tmp2, tmp1, F);
688 F = tmp2;
689
690 // For partials in parallel
691 for (Index ii = 0; ii < n_partials; ii++) {
692 Matrix tmp_low(n, n), tmp_upp(n, n);
693 //dF = D \ (dF - dD*F), or D^-1 * (dF - dD*F)
694 mult(tmp_upp, dD_upp(ii, joker, joker), F);
695 mult(tmp_low, dD_low(ii, joker, joker), F);
696
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); // dF - dD * F
700 dF_low(ii, jj, kk) -= tmp_low(jj, kk); // dF - dD * F
701 }
702 }
703
704 mult(tmp_upp, tmp1, dF_upp(ii, joker, joker));
705 mult(tmp_low, tmp1, dF_low(ii, joker, joker));
706
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);
711 }
712 }
713 }
714
715 for (Index k = 1; k <= r; k++) {
716 // For partials in parallel
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);
719
720 // dF=F*dF+dF*F
721 mult(tmp_upp1, F, dF_upp(ii, joker, joker)); //F*dF
722 mult(tmp_upp2, dF_upp(ii, joker, joker), F); //dF*F
723 mult(tmp_low1, F, dF_low(ii, joker, joker)); //F*dF
724 mult(tmp_low2, dF_low(ii, joker, joker), F); //dF*F
725
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);
730 }
731 }
732 }
733
734 // F=F*F
735 mult(tmp1, F, F);
736 F = tmp1;
737 }
738}
739
740Matrix4x4ViewMap MapToEigen4x4(Tensor3View& A, const Index& i) {
741 MatrixView B = A(i, joker, joker);
742 return MapToEigen4x4(B);
743}
744
745MatrixViewMap MapToEigen(Tensor3View& A, const Index& i) {
746 MatrixView B = A(i, joker, joker);
747 return MapToEigen(B);
748}
749
784void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(
785 MatrixView F, ConstMatrixView A) {
786 static const Numeric sqrt_05 = sqrt(0.5);
787
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);
790 const Numeric exp_a = exp(a);
791 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
792 w2 = w * w;
793
794 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
795
796 Numeric Const1;
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);
802 Const1 *= 2;
803 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
804 Const1 += w2 * w2;
805
806 if (Const1 > 0.0)
807 Const1 = sqrt(Const1);
808 else
809 Const1 = 0.0;
810
811 if (Const1 == 0 and
812 Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
813 {
814 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
815 return;
816 }
817
818 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
819 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
820 const Numeric x = sqrt_BpA.real() * sqrt_05;
821 const Numeric y = sqrt_BmA.imag() * sqrt_05;
822 const Numeric x2 = x * x;
823 const Numeric y2 = y * y;
824 const Numeric cos_y = cos(y);
825 const Numeric sin_y = sin(y);
826 const Numeric cosh_x = cosh(x);
827 const Numeric sinh_x = sinh(x);
828 const Numeric x2y2 = x2 + y2;
829 const Numeric inv_x2y2 = 1.0 / x2y2;
830
831 Numeric C0, C1, C2, C3;
832 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
833
834 // X and Y cannot both be zero
835 if (x == 0.0) {
836 inv_y = 1.0 / y;
837 C0 = 1.0;
838 C1 = 1.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) {
842 inv_x = 1.0 / x;
843 C0 = 1.0;
844 C1 = 1.0;
845 C2 = (cosh_x - 1.0) * inv_x2y2;
846 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
847 } else {
848 inv_x = 1.0 / x;
849 inv_y = 1.0 / y;
850
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;
855 }
856
857 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
858
859 F(0, 0) += C2 * (b2 + c2 + d2);
860
861 F(1, 1) += C2 * (b2 - u2 - v2);
862
863 F(2, 2) += C2 * (c2 - u2 - w2);
864
865 F(3, 3) += C2 * (d2 - v2 - w2);
866
867 F(0, 1) = F(1, 0) = C1 * b;
868
869 F(0, 1) +=
870 C2 * (-c * u - d * v) +
871 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
872
873 F(1, 0) +=
874 C2 * (c * u + d * v) +
875 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) + d * (b * d + u * w));
876
877 F(0, 2) = F(2, 0) = C1 * c;
878
879 F(0, 2) +=
880 C2 * (b * u - d * w) +
881 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
882
883 F(2, 0) +=
884 C2 * (-b * u + d * w) +
885 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) + d * (c * d - u * v));
886
887 F(0, 3) = F(3, 0) = C1 * d;
888
889 F(0, 3) +=
890 C2 * (b * v + c * w) +
891 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
892
893 F(3, 0) +=
894 C2 * (-b * v - c * w) +
895 C3 * (b * (b * d + u * w) + c * (c * d - u * v) - d * (-d2 + v2 + w2));
896
897 F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
898
899 F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
900 w * (b * d + u * w));
901
902 F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
903 v * (c * d - u * v));
904
905 F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
906
907 F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
908 w * (b * c - v * w));
909
910 F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
911 v * (-d2 + v2 + w2));
912
913 F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
914
915 F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
916 w * (-c2 + u2 + w2));
917
918 F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
919 w * (-d2 + v2 + w2));
920
921 F *= exp_a;
922}
923
958void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(
959 MatrixView F, ConstMatrixView A) {
960 Matrix4x4ViewMap eigF = MapToEigen4x4(F);
961 Eigen::Matrix4d eigA = MapToEigen4x4(A);
962 // eigA << 0, A(0, 1), A(0, 2), A(0, 3),
963 // A(0, 1), 0, A(1, 2), A(1, 3),
964 // A(0, 2), -A(1, 2), 0, A(2, 3),
965 // A(0, 3), -A(1, 3), -A(2, 3), 0;
966 eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
967
968 Eigen::Matrix4d eigA2 = eigA;
969 eigA2 *= eigA;
970 Eigen::Matrix4d eigA3 = eigA2;
971 eigA3 *= eigA;
972
973 static const Numeric sqrt_05 = sqrt(0.5);
974
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);
977 const Numeric exp_a = exp(a);
978 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
979 w2 = w * w;
980
981 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
982
983 Numeric Const1;
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);
989 Const1 *= 2;
990 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
991 Const1 += w2 * w2;
992
993 if (Const1 > 0.0)
994 Const1 = sqrt(Const1);
995 else
996 Const1 = 0.0;
997
998 if (Const1 == 0 and
999 Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1000 {
1001 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1002 return;
1003 }
1004
1005 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1006 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1007 const Numeric x = sqrt_BpA.real() * sqrt_05;
1008 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1009 const Numeric x2 = x * x;
1010 const Numeric y2 = y * y;
1011 const Numeric cos_y = cos(y);
1012 const Numeric sin_y = sin(y);
1013 const Numeric cosh_x = cosh(x);
1014 const Numeric sinh_x = sinh(x);
1015 const Numeric x2y2 = x2 + y2;
1016 const Numeric inv_x2y2 = 1.0 / x2y2;
1017
1018 Numeric C0, C1, C2, C3;
1019 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1020
1021 // X and Y cannot both be zero
1022 if (x == 0.0) {
1023 inv_y = 1.0 / y;
1024 C0 = 1.0;
1025 C1 = 1.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) {
1029 inv_x = 1.0 / x;
1030 C0 = 1.0;
1031 C1 = 1.0;
1032 C2 = (cosh_x - 1.0) * inv_x2y2;
1033 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1034 } else {
1035 inv_x = 1.0 / x;
1036 inv_y = 1.0 / y;
1037
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;
1042 }
1043
1044 eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1045 eigF(0, 0) += C0;
1046 eigF(1, 1) += C0;
1047 eigF(2, 2) += C0;
1048 eigF(3, 3) += C0;
1049 eigF *= exp_a;
1050}
1051
1106void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__explicit(
1107 MatrixView F,
1108 Tensor3View dF_upp,
1109 Tensor3View dF_low,
1110 ConstMatrixView A,
1111 ConstTensor3View dA_upp,
1112 ConstTensor3View dA_low) {
1113 static const Numeric sqrt_05 = sqrt(0.5);
1114 const Index npd = dF_upp.npages();
1115
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);
1118 const Numeric exp_a = exp(a);
1119 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
1120 w2 = w * w;
1121
1122 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1123
1124 Numeric Const1;
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);
1130 Const1 *= 2;
1131 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
1132 Const1 += w2 * w2;
1133
1134 if (Const1 > 0.0)
1135 Const1 = sqrt(Const1);
1136 else
1137 Const1 = 0.0;
1138
1139 if (Const1 == 0 and
1140 Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1141 {
1142 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1143 for (Index i = 0; i < npd; i++) {
1144 // mult incase the derivative is still non-zero, e.g., as when mag-field is considered 0 but derivative is computed
1145 mult(dF_upp(i, joker, joker), F, dA_upp(i, joker, joker));
1146 mult(dF_low(i, joker, joker), F, dA_low(i, joker, joker));
1147 }
1148 return;
1149 }
1150
1151 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1152 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1153 const Numeric x = sqrt_BpA.real() * sqrt_05;
1154 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1155 const Numeric x2 = x * x;
1156 const Numeric y2 = y * y;
1157 const Numeric cos_y = cos(y);
1158 const Numeric sin_y = sin(y);
1159 const Numeric cosh_x = cosh(x);
1160 const Numeric sinh_x = sinh(x);
1161 const Numeric x2y2 = x2 + y2;
1162 const Numeric inv_x2y2 = 1.0 / x2y2;
1163
1164 Numeric C0, C1, C2, C3;
1165 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1166
1167 // X and Y cannot both be zero
1168 if (x == 0.0) {
1169 inv_y = 1.0 / y;
1170 C0 = 1.0;
1171 C1 = 1.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) {
1175 inv_x = 1.0 / x;
1176 C0 = 1.0;
1177 C1 = 1.0;
1178 C2 = (cosh_x - 1.0) * inv_x2y2;
1179 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1180 } else {
1181 inv_x = 1.0 / x;
1182 inv_y = 1.0 / y;
1183
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;
1188 }
1189
1190 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = C0;
1191
1192 F(0, 0) += C2 * (b2 + c2 + d2);
1193
1194 F(1, 1) += C2 * (b2 - u2 - v2);
1195
1196 F(2, 2) += C2 * (c2 - u2 - w2);
1197
1198 F(3, 3) += C2 * (d2 - v2 - w2);
1199
1200 F(0, 1) = F(1, 0) = C1 * b;
1201
1202 F(0, 1) +=
1203 C2 * (-c * u - d * v) +
1204 C3 * (b * (b2 + c2 + d2) - u * (b * u - d * w) - v * (b * v + c * w));
1205
1206 F(1, 0) +=
1207 C2 * (c * u + d * v) +
1208 C3 * (-b * (-b2 + u2 + v2) + c * (b * c - v * w) + d * (b * d + u * w));
1209
1210 F(0, 2) = F(2, 0) = C1 * c;
1211
1212 F(0, 2) +=
1213 C2 * (b * u - d * w) +
1214 C3 * (c * (b2 + c2 + d2) - u * (c * u + d * v) - w * (b * v + c * w));
1215
1216 F(2, 0) +=
1217 C2 * (-b * u + d * w) +
1218 C3 * (b * (b * c - v * w) - c * (-c2 + u2 + w2) + d * (c * d - u * v));
1219
1220 F(0, 3) = F(3, 0) = C1 * d;
1221
1222 F(0, 3) +=
1223 C2 * (b * v + c * w) +
1224 C3 * (d * (b2 + c2 + d2) - v * (c * u + d * v) + w * (b * u - d * w));
1225
1226 F(3, 0) +=
1227 C2 * (-b * v - c * w) +
1228 C3 * (b * (b * d + u * w) + c * (c * d - u * v) - d * (-d2 + v2 + w2));
1229
1230 F(1, 2) = F(2, 1) = C2 * (b * c - v * w);
1231
1232 F(1, 2) += C1 * u + C3 * (c * (c * u + d * v) - u * (-b2 + u2 + v2) -
1233 w * (b * d + u * w));
1234
1235 F(2, 1) += -C1 * u + C3 * (-b * (b * u - d * w) + u * (-c2 + u2 + w2) -
1236 v * (c * d - u * v));
1237
1238 F(1, 3) = F(3, 1) = C2 * (b * d + u * w);
1239
1240 F(1, 3) += C1 * v + C3 * (d * (c * u + d * v) - v * (-b2 + u2 + v2) +
1241 w * (b * c - v * w));
1242
1243 F(3, 1) += -C1 * v + C3 * (-b * (b * v + c * w) - u * (c * d - u * v) +
1244 v * (-d2 + v2 + w2));
1245
1246 F(2, 3) = F(3, 2) = C2 * (c * d - u * v);
1247
1248 F(2, 3) += C1 * w + C3 * (-d * (b * u - d * w) + v * (b * c - v * w) -
1249 w * (-c2 + u2 + w2));
1250
1251 F(3, 2) += -C1 * w + C3 * (-c * (b * v + c * w) + u * (b * d + u * w) +
1252 w * (-d2 + v2 + w2));
1253
1254 F *= exp_a;
1255
1256 if (npd) {
1257 const Numeric inv_x2 = inv_x * inv_x;
1258 const Numeric inv_y2 = inv_y * inv_y;
1259
1260 for (Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1261 for (Index i = 0; i < npd; i++) {
1262 MatrixView dF =
1263 upp_or_low ? dF_low(i, joker, joker) : dF_upp(i, joker, joker);
1264 ConstMatrixView dA =
1265 upp_or_low ? dA_low(i, joker, joker) : dA_upp(i, joker, joker);
1266
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),
1269 dw = dA(2, 3);
1270
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;
1273
1274 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1275
1276 Numeric dConst1;
1277 if (Const1 > 0.) {
1278 dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1279 dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1280
1281 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1282 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1283
1284 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1285 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1286
1287 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1288 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1289
1290 dConst1 += dv2 * (v2 * 0.5 + w2);
1291 dConst1 += v2 * (dv2 * 0.5 + dw2);
1292
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;
1298 dConst1 /= Const1;
1299 } else
1300 dConst1 = 0.0;
1301
1302 Numeric dC0, dC1, dC2, dC3;
1303 if (x == 0.0) {
1304 const Numeric dy =
1305 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1306
1307 dC0 = 0.0;
1308 dC1 = 0.0;
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;
1312 ;
1313 } else if (y == 0.0) {
1314 const Numeric dx =
1315 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1316
1317 dC0 = 0.0;
1318 dC1 = 0.0;
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;
1322 } else {
1323 const Numeric dx =
1324 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1325 const Numeric dy =
1326 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1327 const Numeric dy2 = 2 * y * dy;
1328 const Numeric dx2 = 2 * x * dx;
1329 const Numeric dx2dy2 = dx2 + dy2;
1330
1331 dC0 = -dx2dy2 * C0 * inv_x2y2 +
1332 (2 * cos_y * dx * x + 2 * cosh_x * dy * y + dx * sinh_x * y2 -
1333 dy * sin_y * x2) *
1334 inv_x2y2;
1335
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) *
1340 inv_x2y2;
1341
1342 dC2 = -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1343
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) *
1347 inv_x2y2;
1348 }
1349
1350 dF(0, 0) = dF(1, 1) = dF(2, 2) = dF(3, 3) = dC0;
1351
1352 dF(0, 0) += dC2 * (b2 + c2 + d2) + C2 * (db2 + dc2 + dd2);
1353
1354 dF(1, 1) += dC2 * (b2 - u2 - v2) + C2 * (db2 - du2 - dv2);
1355
1356 dF(2, 2) += dC2 * (c2 - u2 - w2) + C2 * (dc2 - du2 - dw2);
1357
1358 dF(3, 3) += dC2 * (d2 - v2 - w2) + C2 * (dd2 - dv2 - dw2);
1359
1360 dF(0, 1) = dF(1, 0) = db * C1 + b * dC1;
1361
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));
1378
1379 dF(0, 2) = dF(2, 0) = dC1 * c + C1 * dc;
1380
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));
1397
1398 dF(0, 3) = dF(3, 0) = dC1 * d + C1 * dd;
1399
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));
1416
1417 dF(1, 2) = dF(2, 1) =
1418 dC2 * (b * c - v * w) + C2 * (db * c + b * dc - dv * w - v * dw);
1419
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));
1434
1435 dF(1, 3) = dF(3, 1) =
1436 dC2 * (b * d + u * w) + C2 * (db * d + b * dd + du * w + u * dw);
1437
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));
1452
1453 dF(2, 3) = dF(3, 2) =
1454 dC2 * (c * d - u * v) + C2 * (dc * d + c * dd - du * v - u * dv);
1455
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));
1470
1471 dF *= exp_a;
1472
1473 // Finalize derivation by the chian rule
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;
1490 }
1491 }
1492 }
1493}
1494
1545void cayley_hamilton_fitted_method_4x4_propmat_to_transmat__eigen(
1546 MatrixView F,
1547 Tensor3View dF_upp,
1548 Tensor3View dF_low,
1549 ConstMatrixView A,
1550 ConstTensor3View dA_upp,
1551 ConstTensor3View dA_low) {
1552 Matrix4x4ViewMap eigF = MapToEigen4x4(F);
1553 Eigen::Matrix4d eigA = MapToEigen4x4(A);
1554 // eigA << 0, A(0, 1), A(0, 2), A(0, 3),
1555 // A(0, 1), 0, A(1, 2), A(1, 3),
1556 // A(0, 2), -A(1, 2), 0, A(2, 3),
1557 // A(0, 3), -A(1, 3), -A(2, 3), 0;
1558 eigA(0, 0) = eigA(1, 1) = eigA(2, 2) = eigA(3, 3) = 0.0;
1559
1560 Eigen::Matrix4d eigA2 = eigA;
1561 eigA2 *= eigA;
1562 Eigen::Matrix4d eigA3 = eigA2;
1563 eigA3 *= eigA;
1564
1565 static const Numeric sqrt_05 = sqrt(0.5);
1566 const Index npd = dF_low.npages();
1567
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);
1570 const Numeric exp_a = exp(a);
1571 const Numeric b2 = b * b, c2 = c * c, d2 = d * d, u2 = u * u, v2 = v * v,
1572 w2 = w * w;
1573
1574 const Numeric Const2 = b2 + c2 + d2 - u2 - v2 - w2;
1575
1576 Numeric Const1;
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);
1582 Const1 *= 2;
1583 Const1 += 8 * (b * d * u * w - b * c * v * w - c * d * u * v);
1584 Const1 += w2 * w2;
1585
1586 if (Const1 > 0.0)
1587 Const1 = sqrt(Const1);
1588 else
1589 Const1 = 0.0;
1590
1591 if (Const1 == 0 and
1592 Const2 == 0) // Diagonal matrix... ---> x and y will be zero...
1593 {
1594 F(0, 0) = F(1, 1) = F(2, 2) = F(3, 3) = exp_a;
1595 for (Index i = 0; i < npd; i++) {
1596 // mult incase the derivative is still non-zero, e.g., as when mag-field is considered 0 but derivative is computed
1597 mult(dF_upp(i, joker, joker), F, dA_upp(i, joker, joker));
1598 mult(dF_low(i, joker, joker), F, dA_low(i, joker, joker));
1599 }
1600 return;
1601 }
1602
1603 const Complex sqrt_BpA = sqrt(Complex(Const2 + Const1, 0.0));
1604 const Complex sqrt_BmA = sqrt(Complex(Const2 - Const1, 0.0));
1605 const Numeric x = sqrt_BpA.real() * sqrt_05;
1606 const Numeric y = sqrt_BmA.imag() * sqrt_05;
1607 const Numeric x2 = x * x;
1608 const Numeric y2 = y * y;
1609 const Numeric cos_y = cos(y);
1610 const Numeric sin_y = sin(y);
1611 const Numeric cosh_x = cosh(x);
1612 const Numeric sinh_x = sinh(x);
1613 const Numeric x2y2 = x2 + y2;
1614 const Numeric inv_x2y2 = 1.0 / x2y2;
1615
1616 Numeric C0, C1, C2, C3;
1617 Numeric inv_y = 0.0, inv_x = 0.0; // Init'd to remove warnings
1618
1619 // X and Y cannot both be zero
1620 if (x == 0.0) {
1621 inv_y = 1.0 / y;
1622 C0 = 1.0;
1623 C1 = 1.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) {
1627 inv_x = 1.0 / x;
1628 C0 = 1.0;
1629 C1 = 1.0;
1630 C2 = (cosh_x - 1.0) * inv_x2y2;
1631 C3 = (sinh_x * inv_x - 1.0) * inv_x2y2;
1632 } else {
1633 inv_x = 1.0 / x;
1634 inv_y = 1.0 / y;
1635
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;
1640 }
1641
1642 eigF(0, 0) = eigF(1, 1) = eigF(2, 2) = eigF(3, 3) = C0;
1643 eigF.noalias() = C1 * eigA + C2 * eigA2 + C3 * eigA3;
1644 eigF(0, 0) += C0;
1645 eigF(1, 1) += C0;
1646 eigF(2, 2) += C0;
1647 eigF(3, 3) += C0;
1648 eigF *= exp_a;
1649
1650 if (npd) {
1651 const Numeric inv_x2 = inv_x * inv_x;
1652 const Numeric inv_y2 = inv_y * inv_y;
1653
1654 for (Index upp_or_low = 0; upp_or_low < 2; upp_or_low++) {
1655 for (Index i = 0; i < npd; i++) {
1656 MatrixView tmp =
1657 upp_or_low ? dF_low(i, joker, joker) : dF_upp(i, joker, joker);
1658 Matrix4x4ViewMap dF = MapToEigen4x4(tmp);
1659 Eigen::Matrix4d dA = MapToEigen4x4(
1660 upp_or_low ? dA_low(i, joker, joker) : dA_upp(i, joker, joker));
1661
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),
1664 dw = dA(2, 3);
1665
1666 dA(0, 0) = dA(1, 1) = dA(2, 2) = dA(3, 3) = 0.0;
1667
1668 Eigen::Matrix4d dA2;
1669 ;
1670 dA2.noalias() = eigA * dA;
1671 dA2.noalias() += dA * eigA;
1672
1673 Eigen::Matrix4d dA3;
1674 dA3.noalias() = dA2 * eigA;
1675 dA3.noalias() += eigA2 * dA;
1676
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;
1679
1680 const Numeric dConst2 = db2 + dc2 + dd2 - du2 - dv2 - dw2;
1681
1682 Numeric dConst1;
1683 if (Const1 > 0.) {
1684 dConst1 = db2 * (b2 * 0.5 + c2 + d2 - u2 - v2 + w2);
1685 dConst1 += b2 * (db2 * 0.5 + dc2 + dd2 - du2 - dv2 + dw2);
1686
1687 dConst1 += dc2 * (c2 * 0.5 + d2 - u2 + v2 - w2);
1688 dConst1 += c2 * (dc2 * 0.5 + dd2 - du2 + dv2 - dw2);
1689
1690 dConst1 += dd2 * (d2 * 0.5 + u2 - v2 - w2);
1691 dConst1 += d2 * (dd2 * 0.5 + du2 - dv2 - dw2);
1692
1693 dConst1 += du2 * (u2 * 0.5 + v2 + w2);
1694 dConst1 += u2 * (du2 * 0.5 + dv2 + dw2);
1695
1696 dConst1 += dv2 * (v2 * 0.5 + w2);
1697 dConst1 += v2 * (dv2 * 0.5 + dw2);
1698
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;
1704 dConst1 /= Const1;
1705 } else
1706 dConst1 = 0.0;
1707
1708 if (x == 0.0) {
1709 const Numeric dy =
1710 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1711
1712 const Numeric dC2 =
1713 -2 * y * dy * C2 * inv_x2y2 + dy * sin_y * inv_x2y2;
1714 const Numeric dC3 =
1715 -2 * y * dy * C3 * inv_x2y2 +
1716 (dy * sin_y * inv_y2 - cos_y * dy * inv_y) * inv_x2y2;
1717
1718 dF.noalias() = dC2 * eigA2;
1719 +C2* dA2 + dC3* eigA3 + C3* dA3;
1720 dF *= exp_a;
1721 dF.noalias() += eigF * da;
1722 } else if (y == 0.0) {
1723 const Numeric dx =
1724 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1725
1726 const Numeric dC2 =
1727 -2 * x * dx * C2 * inv_x2y2 + dx * sinh_x * inv_x2y2;
1728 const Numeric dC3 =
1729 -2 * x * dx * C3 * inv_x2y2 +
1730 (cosh_x * dx * inv_x - dx * sinh_x * inv_x2) * inv_x2y2;
1731
1732 dF.noalias() = dC2 * eigA2 + C2 * dA2 + dC3 * eigA3 + C3 * dA3;
1733 dF *= exp_a;
1734 dF.noalias() += eigF * da;
1735 } else {
1736 const Numeric dx =
1737 (0.5 * (dConst2 + dConst1) / sqrt_BpA).real() * sqrt_05;
1738 const Numeric dy =
1739 (0.5 * (dConst2 - dConst1) / sqrt_BmA).imag() * sqrt_05;
1740 const Numeric dy2 = 2 * y * dy;
1741 const Numeric dx2 = 2 * x * dx;
1742 const Numeric dx2dy2 = dx2 + dy2;
1743
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) *
1747 inv_x2y2;
1748
1749 const Numeric dC1 =
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) *
1754 inv_x2y2;
1755
1756 const Numeric dC2 =
1757 -dx2dy2 * C2 * inv_x2y2 + (dx * sinh_x + dy * sin_y) * inv_x2y2;
1758
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) *
1762 inv_x2y2;
1763
1764 dF.noalias() = dC1 * eigA + C1 * dA + dC2 * eigA2 + C2 * dA2 +
1765 dC3 * eigA3 + C3 * dA3;
1766 dF(0, 0) += dC0;
1767 dF(1, 1) += dC0;
1768 dF(2, 2) += dC0;
1769 dF(3, 3) += dC0;
1770 dF *= exp_a;
1771 dF.noalias() += eigF * da;
1772 }
1773 }
1774 }
1775 }
1776}
1777
1778void propmat4x4_to_transmat4x4(MatrixView F,
1779 Tensor3View dF_upp,
1780 Tensor3View dF_low,
1781 ConstMatrixView A,
1782 ConstTensor3View dA_upp,
1783 ConstTensor3View dA_low,
1784 const Index& q) {
1785 const Index n_partials = dA_upp.npages();
1786
1787 /* Check if A and F are a quadratic and of the same dimension. */
1788 ARTS_ASSERT(is_size(A, 4, 4));
1789 ARTS_ASSERT(is_size(F, 4, 4));
1790 ARTS_ASSERT(n_partials == dF_upp.npages());
1791 ARTS_ASSERT(n_partials == dF_low.npages());
1792 ARTS_ASSERT(n_partials == dA_low.npages());
1793 for (Index ii = 0; ii < n_partials; ii++) {
1794 ARTS_ASSERT(is_size(dA_upp(ii, joker, joker), 4, 4));
1795 ARTS_ASSERT(is_size(dA_low(ii, joker, joker), 4, 4));
1796 ARTS_ASSERT(is_size(dF_upp(ii, joker, joker), 4, 4));
1797 ARTS_ASSERT(is_size(dF_low(ii, joker, joker), 4, 4));
1798 }
1799
1800 // Set constants and Numerics
1801 const Numeric A_norm_inf = norm_inf(A);
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;
1804 const Numeric inv_pow2 = 1. / pow(2, r);
1805 Numeric c = 0.5;
1806
1807 // Create M, dM, X and Y
1808 Eigen::Matrix4d M = MapToEigen4x4(A);
1809 M *= inv_pow2;
1810 Eigen::Matrix4d X = M;
1811 Eigen::Matrix4d cX = c * X;
1812 Eigen::Matrix4d D = Eigen::Matrix4d::Identity();
1813 Matrix4x4ViewMap eigF = MapToEigen4x4(F);
1814 eigF.setIdentity();
1815 eigF.noalias() += cX;
1816 D.noalias() -= cX;
1817
1818 Array<Eigen::Matrix4d> dMu(n_partials), dMl(n_partials), Yu(n_partials),
1819 Yl(n_partials);
1820 Array<Eigen::Matrix4d> cYu(n_partials), cYl(n_partials), dDu(n_partials),
1821 dDl(n_partials);
1822
1823 Array<Matrix4x4ViewMap> dFu, dFl;
1824 dFu.reserve(n_partials);
1825 dFl.reserve(n_partials);
1826
1827 for (Index i = 0; i < n_partials; i++) {
1828 dMu[i].noalias() = MapToEigen4x4(dA_upp(i, joker, joker)) * inv_pow2;
1829 dMl[i].noalias() = MapToEigen4x4(dA_low(i, joker, joker)) * inv_pow2;
1830
1831 Yu[i] = dMu[i];
1832 Yl[i] = dMl[i];
1833
1834 cYu[i].noalias() = c * dMu[i];
1835 cYl[i].noalias() = c * dMl[i];
1836
1837 dFu.push_back(MapToEigen4x4(dF_upp, i));
1838 dFl.push_back(MapToEigen4x4(dF_low, i));
1839 dFu[i] = cYu[i];
1840 dFl[i] = cYl[i];
1841
1842 dDu[i] = -cYu[i];
1843 dDl[i] = -cYl[i];
1844 }
1845
1846 bool p = true;
1847 for (Index k = 2; k <= q; k++) {
1848 c *= (Numeric)(q - k + 1) / (Numeric)(k * (2 * q - k + 1));
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];
1852
1853 cYu[i].noalias() = c * Yu[i];
1854 cYl[i].noalias() = c * Yl[i];
1855
1856 dFu[i].noalias() += cYu[i];
1857 dFl[i].noalias() += cYl[i];
1858 }
1859
1860 X = M * X;
1861 cX.noalias() = c * X;
1862 eigF.noalias() += cX;
1863
1864 if (p) {
1865 D.noalias() += cX;
1866
1867 for (Index i = 0; i < n_partials; i++) {
1868 dDu[i].noalias() += cYu[i];
1869 dDl[i].noalias() += cYl[i];
1870 }
1871 } else {
1872 D.noalias() -= cX;
1873
1874 for (Index i = 0; i < n_partials; i++) {
1875 dDu[i].noalias() -= cYu[i];
1876 dDl[i].noalias() -= cYl[i];
1877 }
1878 }
1879 p = not p;
1880 }
1881
1882 const Eigen::Matrix4d invD = D.inverse();
1883
1884 eigF = invD * eigF;
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);
1888 }
1889
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];
1894 }
1895 eigF = eigF * eigF;
1896 }
1897}
1898
1900
1921void matrix_exp_dmatrix_exp(MatrixView F,
1922 Tensor3View dF,
1923 ConstMatrixView A,
1924 ConstTensor3View dA,
1925 const Index& q) {
1926 const Index n_partials = dA.npages();
1927
1928 const Index n = A.ncols();
1929
1930 /* Check if A and F are a quadratic and of the same dimension. */
1931 ARTS_ASSERT(is_size(A, n, n));
1932 ARTS_ASSERT(is_size(F, n, n));
1933 ARTS_ASSERT(n_partials == dF.npages());
1934 for (Index ii = 0; ii < n_partials; ii++) {
1935 ARTS_ASSERT(is_size(dA(ii, joker, joker), n, n));
1936 ARTS_ASSERT(is_size(dF(ii, joker, joker), n, n));
1937 }
1938
1939 // This sets up some cnstants
1940 Numeric A_norm_inf, e;
1941 A_norm_inf = norm_inf(A);
1942 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
1943 Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
1944 Numeric pow2rm1 = 1. / pow(2, r);
1945 Numeric c = 0.5;
1946
1947 // For non-derivatives
1948 Matrix M = A, X(n, n), cX(n, n), D(n, n);
1949 M *= pow2rm1;
1950 X = M;
1951 cX = X;
1952 cX *= c;
1953 id_mat(F);
1954 F += cX;
1955 id_mat(D);
1956 D -= cX;
1957
1958 // For derivatives
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;
1965
1966 Y(ii, jj, kk) = dM(ii, jj, kk);
1967
1968 cY(ii, jj, kk) = c * Y(ii, jj, kk);
1969
1970 dF(ii, jj, kk) = cY(ii, jj, kk);
1971
1972 dD(ii, jj, kk) = -cY(ii, jj, kk);
1973 }
1974 }
1975 }
1976
1977 // NOTE: MATLAB paper sets q = 6 but we allow other numbers
1978
1979 Matrix tmp1(n, n), tmp2(n, n);
1980
1981 for (Index k = 2; k <= q; k++) {
1982 c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
1983
1984 // For partials
1985 for (Index ii = 0; ii < n_partials; ii++) {
1986 // Y = dM*X + M*Y
1987 mult(tmp1, dM(ii, joker, joker), X);
1988 mult(tmp2, M, Y(ii, joker, joker));
1989
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);
1993
1994 // cY = c*Y
1995 cY(ii, jj, kk) = c * Y(ii, jj, kk);
1996
1997 // dF = dF + cY
1998 dF(ii, jj, kk) += cY(ii, jj, kk);
1999
2000 if (k % 2 == 0) //For even numbers, add.
2001 {
2002 // dD = dD + cY
2003 dD(ii, jj, kk) += cY(ii, jj, kk);
2004 } else {
2005 // dD = dD - cY
2006 dD(ii, jj, kk) -= cY(ii, jj, kk);
2007 }
2008 }
2009 }
2010
2011 //For full derivative (Note X in partials above)
2012 // X=M*X
2013 mult(tmp1, M, X);
2014 X = tmp1;
2015
2016 //cX = c*X
2017 cX = X;
2018 cX *= c;
2019
2020 // F = F + cX
2021 F += cX;
2022
2023 if (k % 2 == 0) //For even numbers, D = D + cX
2024 D += cX;
2025 else //For odd numbers, D = D - cX
2026 D -= cX;
2027 }
2028
2029 // D^-1
2030 inv(tmp1, D);
2031
2032 // F = D\F, or D^-1*F
2033 mult(tmp2, tmp1, F);
2034 F = tmp2;
2035
2036 // For partials
2037 for (Index ii = 0; ii < n_partials; ii++) {
2038 //dF = D \ (dF - dF*F), or D^-1 * (dF - dF*F)
2039 mult(tmp2, dD(ii, joker, joker), F); // dF * F
2040 dF(ii, joker, joker) -= tmp2; // dF - dF * F
2041 mult(tmp2, tmp1, dF(ii, joker, joker));
2042 dF(ii, joker, joker) = tmp2;
2043 }
2044
2045 for (Index k = 1; k <= r; k++) {
2046 for (Index ii = 0; ii < n_partials; ii++) {
2047 // dF=F*dF+dF*F
2048 mult(tmp1, F, dF(ii, joker, joker)); //F*dF
2049 mult(tmp2, dF(ii, joker, joker), F); //dF*F
2050 dF(ii, joker, joker) = tmp1;
2051 dF(ii, joker, joker) += tmp2;
2052 }
2053
2054 // F=F*F
2055 mult(tmp1, F, F);
2056 F = tmp1;
2057 }
2058}
2059
2061
2082void matrix_exp_dmatrix_exp(MatrixView F,
2083 MatrixView dF,
2084 ConstMatrixView A,
2085 ConstMatrixView dA,
2086 const Index& q) {
2087 const Index n = A.ncols();
2088
2089 /* Check if A and F are a quadratic and of the same dimension. */
2090 ARTS_ASSERT(is_size(A, n, n));
2091 ARTS_ASSERT(is_size(F, n, n));
2092 ARTS_ASSERT(is_size(dA, n, n));
2093 ARTS_ASSERT(is_size(dF, n, n));
2094
2095 // This is the definition of how to scale
2096 Numeric A_norm_inf, e;
2097 A_norm_inf = norm_inf(A);
2098 e = 1. + floor(1. / log(2.) * log(A_norm_inf));
2099 Index r = (e + 1.) > 0. ? (Index)(e + 1.) : 0;
2100 Numeric pow2rm1 = 1. / pow(2, r);
2101
2102 // A and dA are scaled
2103 Matrix M = A, dM = dA;
2104 M *= pow2rm1;
2105 dM *= pow2rm1;
2106
2107 // These variables will hold the multiplication calculations
2108 Matrix X(n, n), Y(n, n);
2109 X = M;
2110 Y = dM;
2111
2112 // cX = c*M
2113 // cY = c*dM
2114 Matrix cX = X, cY = Y;
2115 Numeric c = 0.5;
2116 cX *= c;
2117 cY *= c;
2118
2119 Matrix D(n, n), dD(n, n);
2120 // F = I + c*M
2121 id_mat(F);
2122 F += cX;
2123
2124 // dF = c*dM;
2125 dF = cY;
2126
2127 //D = I -c*M
2128 id_mat(D);
2129 D -= cX;
2130
2131 // dD = -c*dM
2132 dD = cY;
2133 dD *= -1.;
2134
2135 // NOTE: MATLAB paper sets q = 6 but we allow other numbers
2136
2137 Matrix tmp1(n, n), tmp2(n, n);
2138
2139 for (Index k = 2; k <= q; k++) {
2140 c *= (Numeric)(q - k + 1) / (Numeric)((k) * (2 * q - k + 1));
2141
2142 // Y = dM*X + M*Y
2143 mult(tmp1, dM, X);
2144 mult(tmp2, M, Y);
2145 Y = tmp1;
2146 Y += tmp2;
2147
2148 // X=M*X
2149 mult(tmp1, M, X);
2150 X = tmp1;
2151
2152 //cX = c*X
2153 cX = X;
2154 cX *= c;
2155
2156 // cY = c*Y
2157 cY = Y;
2158 cY *= c;
2159
2160 // F = F + cX
2161 F += cX;
2162
2163 // dF = dF + cY
2164 dF += cY;
2165
2166 if (k % 2 == 0) //For even numbers, add.
2167 {
2168 // D = D + cX
2169 D += cX;
2170
2171 // dD = dD + cY
2172 dD += cY;
2173 } else //For odd numbers, subtract
2174 {
2175 // D = D - cX
2176 D -= cX;
2177
2178 // dD = dD - cY
2179 dD -= cY;
2180 }
2181 }
2182
2183 // D^-1
2184 inv(tmp1, D);
2185
2186 // F = D\F, or D^-1*F
2187 mult(tmp2, tmp1, F);
2188 F = tmp2;
2189
2190 //dF = D \ (dF - dF*F), or D^-1 * (dF - dF*F)
2191 mult(tmp2, dD, F); // dF * F
2192 dF -= tmp2; // dF - dF * F
2193 mult(tmp2, tmp1, dF);
2194 dF = tmp2;
2195
2196 for (Index k = 1; k <= r; k++) {
2197 // dF=F*dF+dF*F
2198 mult(tmp1, F, dF); //F*dF
2199 mult(tmp2, dF, F); //dF*F
2200 dF = tmp1;
2201 dF += tmp2;
2202
2203 // F=F*F
2204 mult(tmp1, F, F);
2205 F = tmp1;
2206 }
2207}
2208
2210
2218Numeric norm_inf(ConstMatrixView A) {
2219 Numeric norm_inf = 0;
2220
2221 for (Index j = 0; j < A.nrows(); j++) {
2222 Numeric row_sum = 0;
2223 //Calculate the row sum for all rows
2224 for (Index i = 0; i < A.ncols(); i++) row_sum += abs(A(i, j));
2225 //Pick out the row with the highest row sum
2226 if (norm_inf < row_sum) norm_inf = row_sum;
2227 }
2228 return norm_inf;
2229}
2230
2232
2235void id_mat(MatrixView I) {
2236 const Index n = I.ncols();
2237 ARTS_ASSERT(n == I.nrows());
2238
2239 I = 0;
2240 for (Index i = 0; i < n; i++) I(i, i) = 1.;
2241}
2242
2251Numeric det(ConstMatrixView A) {
2252 const Index dim = A.nrows();
2253 ARTS_ASSERT(dim == A.ncols());
2254
2255 if (dim == 3)
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);
2259 else if (dim == 2)
2260 return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
2261 else if (dim == 1)
2262 return A(0, 0);
2263
2264 Numeric ret_val = 0.;
2265
2266 for (Index j = 0; j < dim; j++) {
2267 Matrix temp(dim - 1, dim - 1);
2268 for (Index I = 1; I < dim; I++)
2269 for (Index J = 0; J < dim; J++) {
2270 if (J < j)
2271 temp(I - 1, J) = A(I, J);
2272 else if (J > j)
2273 temp(I - 1, J - 1) = A(I, J);
2274 }
2275
2276 Numeric tempNum = det(temp);
2277
2278 ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
2279 }
2280 return ret_val;
2281}
2282
2297void linreg(Vector& p, ConstVectorView x, ConstVectorView y) {
2298 const Index n = x.nelem();
2299
2300 ARTS_ASSERT(y.nelem() == n);
2301
2302 p.resize(2);
2303
2304 // Basic algorithm found at e.g.
2305 // http://en.wikipedia.org/wiki/Simple_linear_regression
2306 // The basic algorithm is as follows:
2307 /*
2308 Numeric s1=0, s2=0, s3=0, s4=0;
2309 for( Index i=0; i<n; i++ )
2310 {
2311 s1 += x[i] * y[i];
2312 s2 += x[i];
2313 s3 += y[i];
2314 s4 += x[i] * x[i];
2315 }
2316
2317 p[1] = ( s1 - (s2*s3)/n ) / ( s4 - s2*s2/n );
2318 p[0] = s3/n - p[1]*s2/n;
2319 */
2320
2321 // A version abit more numerical stable:
2322 // Mean value of x is removed before the fit: x' = (x-mean(x))
2323 // This corresponds to that s2 in version above becomes 0
2324 // y = a + b*x'
2325 // p[1] = b
2326 // p[0] = a - p[1]*mean(x)
2327
2328 Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
2329
2330 for (Index i = 0; i < n; i++) {
2331 xm += x[i] / Numeric(n);
2332 }
2333
2334 for (Index i = 0; i < n; i++) {
2335 const Numeric xv = x[i] - xm;
2336 s1 += xv * y[i];
2337 s3 += y[i];
2338 s4 += xv * xv;
2339 }
2340
2341 p[1] = s1 / s4;
2342 p[0] = s3 / Numeric(n) - p[1] * xm;
2343}
2344
2345Numeric lsf(VectorView x, ConstMatrixView A, ConstVectorView y, bool residual) noexcept {
2346 // Size of the problem
2347 const Index n = x.nelem();
2348 Matrix AT, ATA(n, n);
2349 Vector ATy(n);
2350
2351 // Solver
2352 AT = transpose(A);
2353 mult(ATA, AT, A);
2354 mult(ATy, AT, y);
2355 solve(x, ATA, ATy);
2356
2357 // Residual
2358 if (residual) {
2359 Vector r(n);
2360 mult(r, ATA, x);
2361 r -= ATy;
2362 return r * r;
2363 } else {
2364 return 0;
2365 }
2366}
2367
2368Eigen::ComplexEigenSolver<Eigen::MatrixXcd> eig(const Eigen::Ref<Eigen::MatrixXcd> A)
2369{
2370 Eigen::ComplexEigenSolver<Eigen::MatrixXcd> ces;
2371 return ces.compute(A);
2372}
array.h
This file contains the definition of Array.
arts.h
The global header file for ARTS.
arts_omp.h
Header file for helper functions for OpenMP.
ComplexMatrixView
The ComplexMatrixView class.
Definition: matpack_complex.h:757
ComplexMatrixView::get_c_array
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpack_complex.cc:1107
ComplexMatrix
The ComplexMatrix class.
Definition: matpack_complex.h:884
ComplexVectorView
The ComplexVectorView class.
Definition: matpack_complex.h:388
ComplexVectorView::get_c_array
const Complex * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array.
Definition: matpack_complex.cc:418
ComplexVector
The ComplexVector class.
Definition: matpack_complex.h:592
ConstComplexMatrixView
A constant view of a ComplexMatrix.
Definition: matpack_complex.h:644
ConstComplexMatrixView::ncols
Index ncols() const noexcept
Definition: matpack_complex.h:656
ConstComplexMatrixView::nrows
Index nrows() const noexcept
Definition: matpack_complex.h:655
ConstComplexMatrixView::get_column_extent
Index get_column_extent() const
Get the extent of the underlying data.
Definition: matpack_complex.h:691
ConstComplexVectorView::nelem
Index nelem() const noexcept
Definition: matpack_complex.h:301
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:1014
ConstMatrixView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1113
ConstMatrixView::mcr
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1111
ConstMatrixView::nrows
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
ConstMatrixView::ncols
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:489
ConstVectorView::mdata
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:612
ConstVectorView::nelem
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
MatrixView
The MatrixView class.
Definition: matpackI.h:1125
Matrix
The Matrix class.
Definition: matpackI.h:1225
Range::get_stride
constexpr Index get_stride() const ARTS_NOEXCEPT
Returns the stride of the range.
Definition: matpackI.h:336
Range::get_extent
constexpr Index get_extent() const ARTS_NOEXCEPT
Returns the extent of the range.
Definition: matpackI.h:334
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
VectorView
The VectorView class.
Definition: matpackI.h:626
VectorView::get_c_array
const Numeric * get_c_array() const ARTS_NOEXCEPT
Conversion to plain C-array, const-version.
Definition: matpackI.cc:271
Vector
The Vector class.
Definition: matpackI.h:876
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
mult
void mult(MatrixView C, ConstMatrixView A, const Block &B)
Definition: covariance_matrix.cc:38
DEBUG_ONLY
#define DEBUG_ONLY(...)
Definition: debug.h:52
ARTS_ASSERT
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
ARTS_USER_ERROR
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
ARTS_USER_ERROR_IF
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
lapack.h
Interface for the LAPACK library.
abs
#define abs(x)
Definition: legacy_continua.cc:20567
q
#define q
Definition: legacy_continua.cc:21653
temp
#define temp
Definition: legacy_continua.cc:20892
dx
#define dx
Definition: legacy_continua.cc:22121
matrix_exp2
void matrix_exp2(MatrixView F, ConstMatrixView A)
Definition: lin_alg.cc:454
matrix_exp_4x4
void matrix_exp_4x4(MatrixView arts_f, ConstMatrixView A, const Index &q)
Definition: lin_alg.cc:476
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:1778
linreg
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:2297
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2235
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:958
solve
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Definition: lin_alg.cc:138
norm_inf
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
Definition: lin_alg.cc:2218
diagonalize
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:241
lubacksub
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:91
matrix_exp2_4x4
void matrix_exp2_4x4(MatrixView arts_f, ConstMatrixView A)
Definition: lin_alg.cc:516
matrix_exp
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:387
ludcmp
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:56
det
Numeric det(ConstMatrixView A)
Definition: lin_alg.cc:2251
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
MapToEigen4x4
Matrix4x4ViewMap MapToEigen4x4(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:740
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:554
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:2345
MapToEigen
MatrixViewMap MapToEigen(Tensor3View &A, const Index &i)
Definition: lin_alg.cc:745
eig
Eigen::ComplexEigenSolver< Eigen::MatrixXcd > eig(const Eigen::Ref< Eigen::MatrixXcd > A)
Return the Eigen decomposition of the eigen matrix.
Definition: lin_alg.cc:2368
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:784
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:1921
lin_alg.h
Linear algebra functions.
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
logic.h
Header file for logic.cc.
MatrixViewMap
Eigen::Map< MatrixType, 0, StrideType > MatrixViewMap
Definition: matpackI.h:115
Matrix4x4ViewMap
Eigen::Map< Matrix4x4Type, 0, StrideType > Matrix4x4ViewMap
Definition: matpackI.h:118
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
transpose
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
Definition: matpack_complex.cc:1499
swap
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Definition: matpack_complex.cc:726
joker
const Joker joker
Complex
std::complex< Numeric > Complex
Definition: matpack_complex.h:35
b2
#define b2
Definition: matpack_complex.h:74
Zeeman::dM
constexpr Index dM(Polarization type) noexcept
Gives the change of M given a polarization type.
Definition: zeemandata.h:52
lapack::zgetri_
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
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::dgetri_
void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info)
Matrix inversion.
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)
lapack::zgetrf_
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
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.
lapack::dgetrf_
void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info)
LU decomposition.
oem::Identity
invlib::MatrixIdentity< Matrix > Identity
Definition: oem.h:38
oem::Matrix
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
u
#define u
d
#define d
v
#define v
w
#define w
a
#define a
c
#define c
b
#define b
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
N
#define N
Definition: rng.cc:164
M
#define M
Definition: rng.cc:165