ARTS 2.5.9 (git: 825fa5f2)
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#include "minimize.h"
40
41// #include "gui/plot.h"
42
43using std::abs;
44using std::cout;
45using std::endl;
46using std::setw;
47
49
53void test_lusolve1D(void) {
54 Matrix a(1, 1);
55 ArrayOfIndex indx(1);
56 Matrix orig(1, 1);
57 Matrix b(1, 1);
58
59 /* Assign test-matrix element. */
60 a(0, 0) = 3;
61
62 /* ------------------------------------------------------------------------
63 Test the function ludcmp.
64 ----------------------------------------------------------------------- */
65 cout << "\n LU decomposition test \n";
66 cout << "initial matrix: \n";
67
68 cout << " " << a(0, 0) << endl;
69
70 /* input: Test-matrix a,
71 output: Decomposed matrix b (includes upper and lower triangle, cp.
72 Numerical Recipies)
73 and index which gives information about pivoting. */
74 ludcmp(b, indx, a);
75
76 cout << "\n after decomposition: ";
77 cout << b(0, 0) << endl;
78
79 /* Seperate b into the two triangular matrices. */
80 Matrix l(1, 1, 0.0);
81 Matrix u(1, 1, 0.0);
82 Matrix lu(1, 1, 0.0);
83
84 l(0, 0) = 1.0;
85 u(0, 0) = b(0, 0);
86
87 /*-------------------------------------------------------------------
88 end of ludcmp test
89 ------------------------------------------------------------------*/
90 /*--------------------------------------------------------------------
91 test backsubstitution routine lubacksub
92 -------------------------------------------------------------------*/
93
94 Vector c(1);
95 c[0] = 6;
96 cout << indx[0] << " " << c[0] << endl;
97
98 Vector x(1);
99 lubacksub(x, b, c, indx);
100
101 cout << "\n solution vector x: ";
102 cout << x[0] << endl;
103}
104
106
110void test_lusolve4D(void) {
111 Matrix a(4, 4);
112 ArrayOfIndex indx(4);
113 Matrix orig(4, 4);
114 Matrix b(4, 4);
115
116 /* Assign test-matrix elements. */
117
118 a(0, 0) = 1;
119 a(0, 1) = 3;
120 a(0, 2) = 5;
121 a(0, 3) = 6;
122 a(1, 0) = 2;
123 a(1, 1) = 3;
124 a(1, 2) = 4;
125 a(1, 3) = 4;
126 a(2, 0) = 1;
127 a(2, 1) = 2;
128 a(2, 2) = 5;
129 a(2, 3) = 1;
130 a(3, 0) = 7;
131 a(3, 1) = 2;
132 a(3, 2) = 4;
133 a(3, 3) = 3;
134
135 /* ------------------------------------------------------------------------
136 Test the function ludcmp.
137 ----------------------------------------------------------------------- */
138
139 cout << "\n LU decomposition test \n";
140 cout << "initial matrix: \n";
141 for (Index i = 0; i < 4; i++) {
142 cout << "\n";
143 for (Index j = 0; j < 4; j++) cout << " " << a(i, j);
144 }
145 cout << "\n";
146
147 /* input: Test-matrix a,
148 output: Decomposed matrix b (includes upper and lower triangle, cp.
149 Numerical Recipies)
150 and index which gives information about pivoting. */
151 ludcmp(b, indx, a);
152
153 cout << "\n after decomposition";
154 for (Index i = 0; i < 4; i++) {
155 cout << "\n";
156 for (Index j = 0; j < 4; j++) cout << " " << b(i, j);
157 }
158 cout << "\n";
159
160 /* Seperate b into the two triangular matrices. */
161 Matrix l(4, 4, 0.0);
162 Matrix u(4, 4, 0.0);
163 Matrix lu(4, 4, 0.0);
164
165 for (Index i = 0; i < 4; i++) l(i, i) = 1.0;
166 l(1, 0) = b(1, 0);
167 l(2, Range(0, 2)) = b(2, Range(0, 2));
168 l(3, Range(0, 3)) = b(3, Range(0, 3));
169
170 cout << "\n Matrix L";
171 for (Index i = 0; i < 4; i++) {
172 cout << "\n";
173 for (Index j = 0; j < 4; j++) cout << " " << l(i, j);
174 }
175 cout << "\n";
176
177 u(0, Range(0, 4)) = b(0, Range(0, 4));
178 u(1, Range(1, 3)) = b(1, Range(1, 3));
179 u(2, Range(2, 2)) = b(2, Range(2, 2));
180 u(3, Range(3, 1)) = b(3, Range(3, 1));
181
182 cout << "\n Matrix U";
183 for (Index i = 0; i < 4; i++) {
184 cout << "\n";
185 for (Index j = 0; j < 4; j++) cout << " " << u(i, j);
186 }
187 cout << "\n";
188
189 /* Test, if LU = a. */
190 mult(lu, l, u);
191
192 cout << "\n product L*U";
193 for (Index i = 0; i < 4; i++) {
194 cout << "\n";
195 for (Index j = 0; j < 4; j++) cout << " " << lu(i, j);
196 }
197 cout << "\n";
198
199 /*-------------------------------------------------------------------
200 end of ludcmp test
201 ------------------------------------------------------------------*/
202
203 /*--------------------------------------------------------------------
204 test backsubstitution routine lubacksub
205 -------------------------------------------------------------------*/
206
207 Vector c(4);
208 c[0] = 2;
209 c[1] = 5;
210 c[2] = 6;
211 c[3] = 7;
212
213 cout << "\n vector indx";
214 for (Index i = 0; i < 4; i++) {
215 cout << "\n";
216 cout << indx[i] << " " << c[i];
217 }
218 cout << endl;
219
220 Vector x(4);
221 lubacksub(x, b, c, indx);
222
223 cout << "\n solution vector x" << endl;
224 for (Index i = 0; i < 4; i++) {
225 cout << "\n";
226 cout << x[i];
227 }
228 cout << endl;
229
230 cout << "\n test solution LU*x";
231 Vector y(4);
232 mult(y, a, x);
233 for (Index i = 0; i < 4; i++) {
234 cout << "\n";
235 cout << y[i];
236 }
237 cout << "\n";
238}
239
241
256void test_solve_linear_system(Index ntests, Index dim, bool verbose) {
257 Matrix A(dim, dim);
258 Matrix LU(dim, dim);
259 Vector x0(dim);
260 Vector x(dim);
261 Vector b(dim);
262 ArrayOfIndex indx(dim);
263
264 // initialize random seed
265 srand((unsigned int)time(0));
266
267 cout << endl << endl << "Testing linear system solution: n = " << dim;
268 cout << ", ntests = " << ntests << endl;
269 cout << endl << setw(10) << "Test no." << setw(20) << "lubacksub(...)";
270 cout << setw(20) << "solve(...)" << endl << endl;
271
272 for (Index i = 0; i < ntests; i++) {
273 // Generate linear system, make sure the determinant
274 // is non-zero.
275 random_fill_matrix_pos_def(A, 10, false);
276 random_fill_vector(x0, 10, false);
277 mult(b, A, x0);
278
279 // Test ludcmp/lubacksub.
280 ludcmp(LU, indx, A);
281 lubacksub(x, LU, b, indx);
282
283 Numeric err = 0.0;
284 err = get_maximum_error(x, x0, true);
285
286 cout << setw(10) << i << setw(20) << err;
287
288 // Test solve.
289 solve(x, A, b);
290
291 err = get_maximum_error(x, x0, true);
292
293 cout << setw(20) << err << endl;
294
295 if (verbose) {
296 cout << endl;
297 cout << "A:" << endl << A << endl << endl;
298 cout << "x0:" << endl << x0 << endl << endl;
299 cout << "x:" << endl << x << endl << endl;
300 cout << "Permutation Vector:" << endl << indx << endl;
301 cout << endl;
302 }
303 }
304}
305
307
323void test_inv(Index ntests, Index dim, bool verbose = false) {
324 Matrix A(dim, dim);
325 Matrix Ainv(dim, dim);
326 Matrix I0(dim, dim);
327 id_mat(I0);
328 Matrix I(dim, dim);
329
330 // initialize random seed
331 srand((unsigned int)time(0));
332
333 cout << endl << endl << "Testing matrix inversion: n = " << dim;
334 cout << ", ntests = " << ntests << endl << endl;
335 cout << setw(10) << "Test no." << setw(20) << "Max. rel. error" << endl
336 << endl;
337
338 for (Index i = 0; i < ntests; i++) {
339 // Generate random matrix, make sure the determinant
340 // is non-zero.
341 random_fill_matrix_pos_def(A, 10, false);
342
343 inv(Ainv, A);
344 mult(I, Ainv, A);
345
346 Numeric err = get_maximum_error(I, I0, false);
347 // Print results.
348 cout << setw(10) << i << setw(20) << err << endl;
349
350 if (verbose) {
351 cout << endl;
352 cout << "A:" << endl << A << endl << endl;
353 cout << "Ainv:" << endl << Ainv << endl << endl;
354 cout << "A*Ainv:" << endl << I << endl << endl;
355 }
356 }
357}
358
360
364 Matrix A(4, 4);
365 Matrix F(4, 4);
366 A(0, 0) = 1;
367 A(0, 1) = 3;
368 A(0, 2) = 5;
369 A(0, 3) = 6;
370 A(1, 0) = 2;
371 A(1, 1) = 3;
372 A(1, 2) = 4;
373 A(1, 3) = 4;
374 A(2, 0) = 1;
375 A(2, 1) = 2;
376 A(2, 2) = 5;
377 A(2, 3) = 1;
378 A(3, 0) = 7;
379 A(3, 1) = 2;
380 A(3, 2) = 4;
381 A(3, 3) = 3;
382
383 /* set parameter for accuracy */
384 Index q = 8;
385
386 /*execute matrix exponential function*/
387 matrix_exp(F, A, q);
388
389 cout << "\n Exponential of Matrix K";
390 for (Index i = 0; i < 4; i++) {
391 cout << "\n";
392 for (Index j = 0; j < 4; j++) cout << " " << F(i, j);
393 }
394 cout << "\n";
395}
396
398
402 Matrix A(1, 1);
403 Matrix F(1, 1);
404 A(0, 0) = 5;
405
406 /* set parameter for accuracy */
407 Index q = 8;
408
409 /*execute matrix exponential function*/
410 matrix_exp(F, A, q);
411
412 cout << "\n Exponential of Matrix A:\n";
413 cout << F(0, 0);
414 cout << "\n";
415}
416
418
422 Matrix A(3, 3);
423 Matrix F(3, 3);
424 A(0, 0) = 1;
425 A(0, 1) = 3;
426 A(0, 2) = 5;
427 A(1, 0) = 2;
428 A(1, 1) = 3;
429 A(1, 2) = 4;
430 A(2, 0) = 1;
431 A(2, 1) = 2;
432 A(2, 2) = 5;
433
434 /* set parameter for accuracy */
435 Index q = 8;
436
437 /*execute matrix exponential function*/
438 matrix_exp(F, A, q);
439
440 cout << "\n Exponential of Matrix A";
441 for (Index i = 0; i < 3; i++) {
442 cout << "\n";
443 for (Index j = 0; j < 3; j++) cout << " " << F(i, j);
444 }
445 cout << "\n";
446}
447
449 Matrix A(dim, dim), F1(dim, dim), F2(dim, dim), tmp1(dim, dim),
450 tmp2(dim, dim), P(dim, dim);
451 Vector Wr(dim), Wi(dim);
452
453 const Matrix ZEROES(dim, dim, 0);
454
455 // initialize random seed
456 srand((unsigned int)time(0));
457
458 cout << endl << endl << "Testing diagonalize: n = " << dim;
459 cout << ", ntests = " << ntests << endl;
460 cout << setw(10) << "Test no." << setw(25) << "Max. rel. expm error";
461 cout << setw(25) << "Max. abs. P^-1*A*P-W" << endl << endl;
462
463 for (Index i = 0; i < ntests; i++) {
464 // Generate a matrix that does not have complex answers...
465 random_fill_matrix_symmetric(A, 10, false);
466
467 // Use the two expm methods to test that diagonalize works
468 matrix_exp(F1, A, 10);
469 matrix_exp(F2, A);
470
471 Numeric err1 = 0.0, err2 = 0.0;
472 err1 = get_maximum_error(F1, F2, true);
473
474 // diagonalize directly to test that P^-1*A*P gives a diagonal matrix
475 diagonalize(P, Wr, Wi, A);
476
477 // P^-1*A*P
478 inv(tmp1, P);
479 mult(tmp2, tmp1, A);
480 mult(tmp1, tmp2, P);
481
482 // Minus W as diagonal matrix
483 for (Index j = 0; j < dim; j++) {
484 tmp1(j, j) -= Wr[j];
485 }
486
487 err2 = get_maximum_error(ZEROES, tmp1, false);
488
489 cout << setw(10) << i << setw(25) << err1 << setw(25) << err2 << endl;
490 }
491}
492
494 ComplexMatrix A(dim, dim), tmp1(dim, dim), tmp2(dim, dim), P(dim, dim);
495 ComplexVector W(dim);
496
497 const ComplexMatrix ZEROES(dim, dim, 0);
498
499 // initialize random seed
500 srand((unsigned int)time(0));
501
502 cout << endl << endl << "Testing diagonalize: n = " << dim;
503 cout << ", ntests = " << ntests << endl;
504 cout << setw(10) << "Test no.";
505 cout << setw(25) << "Max. abs. P^-1*A*P-W" << endl << endl;
506
507 for (Index i = 0; i < ntests; i++) {
508 // Generate a matrix that does not have complex answers...
509 random_fill_matrix_symmetric(A, 10, false);
510
511 Numeric err = 0.0;
512
513 // diagonalize directly to test that P^-1*A*P gives a diagonal matrix
514 diagonalize(P, W, A);
515
516 // P^-1*A*P
517 inv(tmp1, P);
518 mult(tmp2, tmp1, A);
519 mult(tmp1, tmp2, P);
520
521 // Minus W as diagonal matrix
522 for (Index j = 0; j < dim; j++) {
523 tmp1(j, j) -= W[j];
524 }
525
526 err = get_maximum_error(ZEROES, tmp1, false);
527
528 cout << setw(10) << i << setw(25) << err << endl;
529 }
530}
531
532int main() {
533 // test_lusolve4D();
534 // test_inv( 20, 1000 );
535 // test_solve_linear_system( 20, 1000, false );
536 // test_matrix_exp1D();
537 // test_real_diagonalize(20,100);
539 return (0);
540}
This file contains the definition of Array.
The ComplexMatrix class.
The ComplexVector class.
The Matrix class.
Definition: matpackI.h:1285
The range class.
Definition: matpackI.h:160
The Vector class.
Definition: matpackI.h:910
void mult(MatrixView C, ConstMatrixView A, const Block &B)
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:250
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:92
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:396
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:57
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:168
Linear algebra functions.
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
void test_inv(Index ntests, Index dim, bool verbose=false)
Test matrix inversion.
Definition: test_linalg.cc:323
void test_matrix_exp3D(void)
Test for the matrix exponential function (3D matrix)
Definition: test_linalg.cc:421
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:256
void test_real_diagonalize(Index ntests, Index dim)
Definition: test_linalg.cc:448
void test_lusolve1D(void)
Definition: test_linalg.cc:53
void test_complex_diagonalize(Index ntests, Index dim)
Definition: test_linalg.cc:493
void test_matrix_exp1D(void)
Test for the matrix exponential function (3D matrix)
Definition: test_linalg.cc:401
void test_matrix_exp4D(void)
Test for the matrix exponential function (4D matrix)
Definition: test_linalg.cc:363
int main()
Definition: test_linalg.cc:532
void test_lusolve4D(void)
Definition: test_linalg.cc:110
void random_fill_matrix_symmetric(MatrixView A, Numeric range, bool positive)
Generate random, symmetric matrix.
Definition: test_utils.cc:156
void random_fill_matrix_pos_def(MatrixView A, Numeric range, bool positive)
Generate random, positive definite matrix.
Definition: test_utils.cc:181
Numeric get_maximum_error(ConstVectorView v1, ConstVectorView v2, bool relative)
Maximum element-wise error of two vectors.
Definition: test_utils.cc:297
void random_fill_vector(VectorView v, Numeric range, bool positive)
Fill vector with random values.
Definition: test_utils.cc:230
Utility functions for testing.
#define u
#define a
#define c
#define b