ARTS  2.4.0(git:4fb77825)
test_linalg.cc
Go to the documentation of this file.
1 /* Copyright (C) 2002-2012 Claudia Emde <claudia.emde@dlr.de>
2 
3  This program is free software; you can redistribute it and/or modify it
4  under the terms of the GNU General Public License as published by the
5  Free Software Foundation; either version 2, or (at your option) any
6  later version.
7 
8  This program is distributed in the hope that it will be useful,
9  but WITHOUT ANY WARRANTY; without even the implied warranty of
10  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11  GNU General Public License for more details.
12 
13  You should have received a copy of the GNU General Public License
14  along with this program; if not, write to the Free Software
15  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
16  USA.
17 
18 */
19 
30 #include <stdlib.h>
31 #include <time.h>
32 #include <cmath>
33 #include <iostream>
34 #include <random>
35 #include "array.h"
36 #include "lin_alg.h"
37 #include "test_utils.h"
38 
39 using std::abs;
40 using std::cout;
41 using std::endl;
42 using std::setw;
43 
45 
49 void test_lusolve1D(void) {
50  Matrix a(1, 1);
51  ArrayOfIndex indx(1);
52  Matrix orig(1, 1);
53  Matrix b(1, 1);
54 
55  /* Assign test-matrix element. */
56  a(0, 0) = 3;
57 
58  /* ------------------------------------------------------------------------
59  Test the function ludcmp.
60  ----------------------------------------------------------------------- */
61  cout << "\n LU decomposition test \n";
62  cout << "initial matrix: \n";
63 
64  cout << " " << a(0, 0) << endl;
65 
66  /* input: Test-matrix a,
67  output: Decomposed matrix b (includes upper and lower triangle, cp.
68  Numerical Recipies)
69  and index which gives information about pivoting. */
70  ludcmp(b, indx, a);
71 
72  cout << "\n after decomposition: ";
73  cout << b(0, 0) << endl;
74 
75  /* Seperate b into the two triangular matrices. */
76  Matrix l(1, 1, 0.0);
77  Matrix u(1, 1, 0.0);
78  Matrix lu(1, 1, 0.0);
79 
80  l(0, 0) = 1.0;
81  u(0, 0) = b(0, 0);
82 
83  /*-------------------------------------------------------------------
84  end of ludcmp test
85  ------------------------------------------------------------------*/
86  /*--------------------------------------------------------------------
87  test backsubstitution routine lubacksub
88  -------------------------------------------------------------------*/
89 
90  Vector c(1);
91  c[0] = 6;
92  cout << indx[0] << " " << c[0] << endl;
93 
94  Vector x(1);
95  lubacksub(x, b, c, indx);
96 
97  cout << "\n solution vector x: ";
98  cout << x[0] << endl;
99 }
100 
102 
106 void test_lusolve4D(void) {
107  Matrix a(4, 4);
108  ArrayOfIndex indx(4);
109  Matrix orig(4, 4);
110  Matrix b(4, 4);
111 
112  /* Assign test-matrix elements. */
113 
114  a(0, 0) = 1;
115  a(0, 1) = 3;
116  a(0, 2) = 5;
117  a(0, 3) = 6;
118  a(1, 0) = 2;
119  a(1, 1) = 3;
120  a(1, 2) = 4;
121  a(1, 3) = 4;
122  a(2, 0) = 1;
123  a(2, 1) = 2;
124  a(2, 2) = 5;
125  a(2, 3) = 1;
126  a(3, 0) = 7;
127  a(3, 1) = 2;
128  a(3, 2) = 4;
129  a(3, 3) = 3;
130 
131  /* ------------------------------------------------------------------------
132  Test the function ludcmp.
133  ----------------------------------------------------------------------- */
134 
135  cout << "\n LU decomposition test \n";
136  cout << "initial matrix: \n";
137  for (Index i = 0; i < 4; i++) {
138  cout << "\n";
139  for (Index j = 0; j < 4; j++) cout << " " << a(i, j);
140  }
141  cout << "\n";
142 
143  /* input: Test-matrix a,
144  output: Decomposed matrix b (includes upper and lower triangle, cp.
145  Numerical Recipies)
146  and index which gives information about pivoting. */
147  ludcmp(b, indx, a);
148 
149  cout << "\n after decomposition";
150  for (Index i = 0; i < 4; i++) {
151  cout << "\n";
152  for (Index j = 0; j < 4; j++) cout << " " << b(i, j);
153  }
154  cout << "\n";
155 
156  /* Seperate b into the two triangular matrices. */
157  Matrix l(4, 4, 0.0);
158  Matrix u(4, 4, 0.0);
159  Matrix lu(4, 4, 0.0);
160 
161  for (Index i = 0; i < 4; i++) l(i, i) = 1.0;
162  l(1, 0) = b(1, 0);
163  l(2, Range(0, 2)) = b(2, Range(0, 2));
164  l(3, Range(0, 3)) = b(3, Range(0, 3));
165 
166  cout << "\n Matrix L";
167  for (Index i = 0; i < 4; i++) {
168  cout << "\n";
169  for (Index j = 0; j < 4; j++) cout << " " << l(i, j);
170  }
171  cout << "\n";
172 
173  u(0, Range(0, 4)) = b(0, Range(0, 4));
174  u(1, Range(1, 3)) = b(1, Range(1, 3));
175  u(2, Range(2, 2)) = b(2, Range(2, 2));
176  u(3, Range(3, 1)) = b(3, Range(3, 1));
177 
178  cout << "\n Matrix U";
179  for (Index i = 0; i < 4; i++) {
180  cout << "\n";
181  for (Index j = 0; j < 4; j++) cout << " " << u(i, j);
182  }
183  cout << "\n";
184 
185  /* Test, if LU = a. */
186  mult(lu, l, u);
187 
188  cout << "\n product L*U";
189  for (Index i = 0; i < 4; i++) {
190  cout << "\n";
191  for (Index j = 0; j < 4; j++) cout << " " << lu(i, j);
192  }
193  cout << "\n";
194 
195  /*-------------------------------------------------------------------
196  end of ludcmp test
197  ------------------------------------------------------------------*/
198 
199  /*--------------------------------------------------------------------
200  test backsubstitution routine lubacksub
201  -------------------------------------------------------------------*/
202 
203  Vector c(4);
204  c[0] = 2;
205  c[1] = 5;
206  c[2] = 6;
207  c[3] = 7;
208 
209  cout << "\n vector indx";
210  for (Index i = 0; i < 4; i++) {
211  cout << "\n";
212  cout << indx[i] << " " << c[i];
213  }
214  cout << endl;
215 
216  Vector x(4);
217  lubacksub(x, b, c, indx);
218 
219  cout << "\n solution vector x" << endl;
220  for (Index i = 0; i < 4; i++) {
221  cout << "\n";
222  cout << x[i];
223  }
224  cout << endl;
225 
226  cout << "\n test solution LU*x";
227  Vector y(4);
228  mult(y, a, x);
229  for (Index i = 0; i < 4; i++) {
230  cout << "\n";
231  cout << y[i];
232  }
233  cout << "\n";
234 }
235 
237 
252 void test_solve_linear_system(Index ntests, Index dim, bool verbose) {
253  Matrix A(dim, dim);
254  Matrix LU(dim, dim);
255  Vector x0(dim);
256  Vector x(dim);
257  Vector b(dim);
258  ArrayOfIndex indx(dim);
259 
260  // initialize random seed
261  srand((unsigned int)time(0));
262 
263  cout << endl << endl << "Testing linear system solution: n = " << dim;
264  cout << ", ntests = " << ntests << endl;
265  cout << endl << setw(10) << "Test no." << setw(20) << "lubacksub(...)";
266  cout << setw(20) << "solve(...)" << endl << endl;
267 
268  for (Index i = 0; i < ntests; i++) {
269  // Generate linear system, make sure the determinant
270  // is non-zero.
271  random_fill_matrix_pos_def(A, 10, false);
272  random_fill_vector(x0, 10, false);
273  mult(b, A, x0);
274 
275  // Test ludcmp/lubacksub.
276  ludcmp(LU, indx, A);
277  lubacksub(x, LU, b, indx);
278 
279  Numeric err = 0.0;
280  err = get_maximum_error(x, x0, true);
281 
282  cout << setw(10) << i << setw(20) << err;
283 
284  // Test solve.
285  solve(x, A, b);
286 
287  err = get_maximum_error(x, x0, true);
288 
289  cout << setw(20) << err << endl;
290 
291  if (verbose) {
292  cout << endl;
293  cout << "A:" << endl << A << endl << endl;
294  cout << "x0:" << endl << x0 << endl << endl;
295  cout << "x:" << endl << x << endl << endl;
296  cout << "Permutation Vector:" << endl << indx << endl;
297  cout << endl;
298  }
299  }
300 }
301 
303 
319 void test_inv(Index ntests, Index dim, bool verbose = false) {
320  Matrix A(dim, dim);
321  Matrix Ainv(dim, dim);
322  Matrix I0(dim, dim);
323  id_mat(I0);
324  Matrix I(dim, dim);
325 
326  // initialize random seed
327  srand((unsigned int)time(0));
328 
329  cout << endl << endl << "Testing matrix inversion: n = " << dim;
330  cout << ", ntests = " << ntests << endl << endl;
331  cout << setw(10) << "Test no." << setw(20) << "Max. rel. error" << endl
332  << endl;
333 
334  for (Index i = 0; i < ntests; i++) {
335  // Generate random matrix, make sure the determinant
336  // is non-zero.
337  random_fill_matrix_pos_def(A, 10, false);
338 
339  inv(Ainv, A);
340  mult(I, Ainv, A);
341 
342  Numeric err = get_maximum_error(I, I0, false);
343  // Print results.
344  cout << setw(10) << i << setw(20) << err << endl;
345 
346  if (verbose) {
347  cout << endl;
348  cout << "A:" << endl << A << endl << endl;
349  cout << "Ainv:" << endl << Ainv << endl << endl;
350  cout << "A*Ainv:" << endl << I << endl << endl;
351  }
352  }
353 }
354 
356 
359 void test_matrix_exp4D(void) {
360  Matrix A(4, 4);
361  Matrix F(4, 4);
362  A(0, 0) = 1;
363  A(0, 1) = 3;
364  A(0, 2) = 5;
365  A(0, 3) = 6;
366  A(1, 0) = 2;
367  A(1, 1) = 3;
368  A(1, 2) = 4;
369  A(1, 3) = 4;
370  A(2, 0) = 1;
371  A(2, 1) = 2;
372  A(2, 2) = 5;
373  A(2, 3) = 1;
374  A(3, 0) = 7;
375  A(3, 1) = 2;
376  A(3, 2) = 4;
377  A(3, 3) = 3;
378 
379  /* set parameter for accuracy */
380  Index q = 8;
381 
382  /*execute matrix exponential function*/
383  matrix_exp(F, A, q);
384 
385  cout << "\n Exponential of Matrix K";
386  for (Index i = 0; i < 4; i++) {
387  cout << "\n";
388  for (Index j = 0; j < 4; j++) cout << " " << F(i, j);
389  }
390  cout << "\n";
391 }
392 
394 
397 void test_matrix_exp1D(void) {
398  Matrix A(1, 1);
399  Matrix F(1, 1);
400  A(0, 0) = 5;
401 
402  /* set parameter for accuracy */
403  Index q = 8;
404 
405  /*execute matrix exponential function*/
406  matrix_exp(F, A, q);
407 
408  cout << "\n Exponential of Matrix A:\n";
409  cout << F(0, 0);
410  cout << "\n";
411 }
412 
414 
417 void test_matrix_exp3D(void) {
418  Matrix A(3, 3);
419  Matrix F(3, 3);
420  A(0, 0) = 1;
421  A(0, 1) = 3;
422  A(0, 2) = 5;
423  A(1, 0) = 2;
424  A(1, 1) = 3;
425  A(1, 2) = 4;
426  A(2, 0) = 1;
427  A(2, 1) = 2;
428  A(2, 2) = 5;
429 
430  /* set parameter for accuracy */
431  Index q = 8;
432 
433  /*execute matrix exponential function*/
434  matrix_exp(F, A, q);
435 
436  cout << "\n Exponential of Matrix A";
437  for (Index i = 0; i < 3; i++) {
438  cout << "\n";
439  for (Index j = 0; j < 3; j++) cout << " " << F(i, j);
440  }
441  cout << "\n";
442 }
443 
444 void test_real_diagonalize(Index ntests, Index dim) {
445  Matrix A(dim, dim), F1(dim, dim), F2(dim, dim), tmp1(dim, dim),
446  tmp2(dim, dim), P(dim, dim);
447  Vector Wr(dim), Wi(dim);
448 
449  const Matrix ZEROES(dim, dim, 0);
450 
451  // initialize random seed
452  srand((unsigned int)time(0));
453 
454  cout << endl << endl << "Testing diagonalize: n = " << dim;
455  cout << ", ntests = " << ntests << endl;
456  cout << setw(10) << "Test no." << setw(25) << "Max. rel. expm error";
457  cout << setw(25) << "Max. abs. P^-1*A*P-W" << endl << endl;
458 
459  for (Index i = 0; i < ntests; i++) {
460  // Generate a matrix that does not have complex answers...
461  random_fill_matrix_symmetric(A, 10, false);
462 
463  // Use the two expm methods to test that diagonalize works
464  matrix_exp(F1, A, 10);
465  matrix_exp2(F2, A);
466 
467  Numeric err1 = 0.0, err2 = 0.0;
468  err1 = get_maximum_error(F1, F2, true);
469 
470  // diagonalize directly to test that P^-1*A*P gives a diagonal matrix
471  diagonalize(P, Wr, Wi, A);
472 
473  // P^-1*A*P
474  inv(tmp1, P);
475  mult(tmp2, tmp1, A);
476  mult(tmp1, tmp2, P);
477 
478  // Minus W as diagonal matrix
479  for (Index j = 0; j < dim; j++) {
480  tmp1(j, j) -= Wr[j];
481  }
482 
483  err2 = get_maximum_error(ZEROES, tmp1, false);
484 
485  cout << setw(10) << i << setw(25) << err1 << setw(25) << err2 << endl;
486  }
487 }
488 
490  ComplexMatrix A(dim, dim), tmp1(dim, dim), tmp2(dim, dim), P(dim, dim);
491  ComplexVector W(dim);
492 
493  const ComplexMatrix ZEROES(dim, dim, 0);
494 
495  // initialize random seed
496  srand((unsigned int)time(0));
497 
498  cout << endl << endl << "Testing diagonalize: n = " << dim;
499  cout << ", ntests = " << ntests << endl;
500  cout << setw(10) << "Test no.";
501  cout << setw(25) << "Max. abs. P^-1*A*P-W" << endl << endl;
502 
503  for (Index i = 0; i < ntests; i++) {
504  // Generate a matrix that does not have complex answers...
505  random_fill_matrix_symmetric(A, 10, false);
506 
507  Numeric err = 0.0;
508 
509  // diagonalize directly to test that P^-1*A*P gives a diagonal matrix
510  diagonalize(P, W, A);
511 
512  // P^-1*A*P
513  inv(tmp1, P);
514  mult(tmp2, tmp1, A);
515  mult(tmp1, tmp2, P);
516 
517  // Minus W as diagonal matrix
518  for (Index j = 0; j < dim; j++) {
519  tmp1(j, j) -= W[j];
520  }
521 
522  err = get_maximum_error(ZEROES, tmp1, false);
523 
524  cout << setw(10) << i << setw(25) << err << endl;
525  }
526 }
527 
528 void test_matrix_exp_propmat(Index nruns, Index ndiffs) {
529  for (Index i = 0; i < nruns; i++) {
530  Numeric a = ((Numeric)rand() / RAND_MAX), b = ((Numeric)rand() / RAND_MAX),
531  c = ((Numeric)rand() / RAND_MAX), d = ((Numeric)rand() / RAND_MAX),
532  u = ((Numeric)rand() / RAND_MAX), v = ((Numeric)rand() / RAND_MAX),
533  w = ((Numeric)rand() / RAND_MAX);
534 
535  Matrix A(4, 4), F0(4, 4), F1(4, 4), F2(4, 4);
536 
537  A(0, 0) = A(1, 1) = A(2, 2) = A(3, 3) = a;
538  A(0, 1) = A(1, 0) = b;
539  A(0, 2) = A(2, 0) = c;
540  A(0, 3) = A(3, 0) = d;
541  A(1, 2) = u;
542  A(2, 1) = -A(1, 2);
543  A(1, 3) = v;
544  A(3, 1) = -A(1, 3);
545  A(2, 3) = w;
546  A(3, 2) = -A(2, 3);
547 
548  A *= -1.0;
549 
550  Tensor3 dF1(ndiffs, 4, 4), dF2(ndiffs, 4, 4), dF3(ndiffs, 4, 4),
551  dF4(ndiffs, 4, 4), dA1(ndiffs, 4, 4), dA2(ndiffs, 4, 4);
552 
553  Numeric da1 = ((Numeric)rand() / RAND_MAX),
554  db1 = ((Numeric)rand() / RAND_MAX),
555  dc1 = ((Numeric)rand() / RAND_MAX),
556  dd1 = ((Numeric)rand() / RAND_MAX),
557  du1 = ((Numeric)rand() / RAND_MAX),
558  dv1 = ((Numeric)rand() / RAND_MAX),
559  dw1 = ((Numeric)rand() / RAND_MAX);
560 
561  Numeric da2 = ((Numeric)rand() / RAND_MAX),
562  db2 = ((Numeric)rand() / RAND_MAX),
563  dc2 = ((Numeric)rand() / RAND_MAX),
564  dd2 = ((Numeric)rand() / RAND_MAX),
565  du2 = ((Numeric)rand() / RAND_MAX),
566  dv2 = ((Numeric)rand() / RAND_MAX),
567  dw2 = ((Numeric)rand() / RAND_MAX);
568 
569  dA1(joker, 0, 0) = dA1(joker, 1, 1) = dA1(joker, 2, 2) = dA1(joker, 3, 3) =
570  da1;
571  dA1(joker, 0, 1) = dA1(joker, 1, 0) = db1;
572  dA1(joker, 0, 2) = dA1(joker, 2, 0) = dc1;
573  dA1(joker, 0, 3) = dA1(joker, 3, 0) = dd1;
574  dA1(joker, 1, 2) = du1;
575  dA1(joker, 2, 1) = -du1;
576  dA1(joker, 1, 3) = dv1;
577  dA1(joker, 3, 1) = -dv1;
578  dA1(joker, 2, 3) = dw1;
579  dA1(joker, 3, 2) = -dw1;
580 
581  dA1 *= -1.0;
582 
583  dA2(joker, 0, 0) = dA2(joker, 1, 1) = dA2(joker, 2, 2) = dA2(joker, 3, 3) =
584  da2;
585  dA2(joker, 0, 1) = dA2(joker, 1, 0) = db2;
586  dA2(joker, 0, 2) = dA2(joker, 2, 0) = dc2;
587  dA2(joker, 0, 3) = dA2(joker, 3, 0) = dd2;
588  dA2(joker, 1, 2) = du2;
589  dA2(joker, 2, 1) = -du2;
590  dA2(joker, 1, 3) = dv2;
591  dA2(joker, 3, 1) = -dv2;
592  dA2(joker, 2, 3) = dw2;
593  dA2(joker, 3, 2) = -dw2;
594 
595  dA2 *= -1.0;
596 
597  special_matrix_exp_and_dmatrix_exp_dx_for_rt(F1, dF1, dF2, A, dA1, dA2, 10);
598 
600  F1, dF1, dF2, A, dA1, dA2);
602  F2, dF3, dF4, A, dA1, dA2);
603 
604  matrix_exp(F0, A);
607 
608  std::cout << MapToEigen(F1) << "\n\n";
609  std::cout << MapToEigen(F2) << "\n\n";
610  std::cout << MapToEigen(F1) - MapToEigen(F2) << "\n\n";
611  std::cout << MapToEigen(dF1(0, joker, joker)) << "\n\n";
612  std::cout << MapToEigen(dF3(0, joker, joker)) << "\n\n";
613  std::cout << MapToEigen(dF3(0, joker, joker)) -
614  MapToEigen(dF1(0, joker, joker))
615  << "\n\n";
616  std::cout << MapToEigen(dF2(0, joker, joker)) << "\n\n";
617  std::cout << MapToEigen(dF4(0, joker, joker)) << "\n\n";
618  std::cout << MapToEigen(dF4(0, joker, joker)) -
619  MapToEigen(dF2(0, joker, joker))
620  << "\n\n";
621  }
622 }
623 
624 int main(void) {
625  // test_lusolve4D();
626  // test_inv( 20, 1000 );
627  // test_solve_linear_system( 20, 1000, false );
628  // test_matrix_exp1D();
629  // test_real_diagonalize(20,100);
630  // test_complex_diagonalize(20,100);
631  test_matrix_exp_propmat(100000, 10);
632  return (0);
633 }
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
random_fill_matrix_symmetric
void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive)
Generate random, symmetric matrix.
Definition: test_utils.cc:156
test_matrix_exp_propmat
void test_matrix_exp_propmat(Index nruns, Index ndiffs)
Definition: test_linalg.cc:528
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
test_utils.h
Utility functions for testing.
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
id_mat
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:2240
joker
const Joker joker
inv
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
ComplexVector
The ComplexVector class.
Definition: complex.h:573
abs
#define abs(x)
Definition: legacy_continua.cc:20626
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
random_fill_vector
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
x0
#define x0
Definition: lineshapemodel.h:413
test_inv
void test_inv(Index ntests, Index dim, bool verbose=false)
Test matrix inversion.
Definition: test_linalg.cc:319
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
test_matrix_exp3D
void test_matrix_exp3D(void)
Test for the matrix exponential function (3D matrix)
Definition: test_linalg.cc:417
Array< Index >
MapToEigen
ComplexConstMatrixViewMap MapToEigen(const ConstComplexMatrixView &A)
Definition: complex.cc:1655
test_real_diagonalize
void test_real_diagonalize(Index ntests, Index dim)
Definition: test_linalg.cc:444
test_matrix_exp4D
void test_matrix_exp4D(void)
Test for the matrix exponential function (4D matrix)
Definition: test_linalg.cc:359
test_lusolve4D
void test_lusolve4D(void)
Definition: test_linalg.cc:106
diagonalize
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:245
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.
ARTS::Var::time
Time time(Workspace &ws) noexcept
Definition: autoarts.h:7043
test_matrix_exp1D
void test_matrix_exp1D(void)
Test for the matrix exponential function (3D matrix)
Definition: test_linalg.cc:397
test_lusolve1D
void test_lusolve1D(void)
Definition: test_linalg.cc:49
Range
The range class.
Definition: matpackI.h:160
q
#define q
Definition: legacy_continua.cc:21712
test_solve_linear_system
void test_solve_linear_system(Index ntests, Index dim, bool verbose)
Test ludcmp and lubacksub by solving a linear system of equations.
Definition: test_linalg.cc:252
get_maximum_error
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
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
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
test_complex_diagonalize
void test_complex_diagonalize(Index ntests, Index dim)
Definition: test_linalg.cc:489
F0
G0 G2 FVC Y DV F0
Definition: arts_api_classes.cc:218
ludcmp
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:56
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
main
int main(void)
Definition: test_linalg.cc:624
random_fill_matrix_pos_def
void random_fill_matrix_pos_def(MatrixView A, Numeric range, bool positive)
Generate random, positive definite matrix.
Definition: test_utils.cc:181
lubacksub
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:91
solve
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
Definition: covariance_matrix.cc:594