ARTS 2.5.10 (git: 2f1c442c)
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 <stdexcept>
39#include <vector>
40
41#include "array.h"
42#include "arts.h"
43#include "arts_omp.h"
44#include "lapack.h"
45#include "logic.h"
46#include "matpackIII.h"
47
49
58 // Assert that A is quadratic.
59 const Index n = A.nrows();
60 ARTS_ASSERT(is_size(LU, n, n));
61
62 int n_int, info;
63 auto ipiv = std::vector<int>(n);
64
65 LU = transpose(A);
66
67 // Standard case: The arts matrix is not transposed, the leading
68 // dimension is the row stride of the matrix.
69 n_int = (int)n;
70
71 // Compute LU decomposition using LAPACK dgetrf_.
72 lapack::dgetrf_(&n_int, &n_int, LU.mdata, &n_int, ipiv.data(), &info);
73
74 // Copy pivot array to pivot vector.
75 for (Index i = 0; i < n; i++) {
76 indx[i] = ipiv[i];
77 }
78}
79
81
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 std::vector<int> ipiv(n);
112 std::vector<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
121 &trans, &n_int, &one, LU.mdata, &n_int, ipiv.data(), rhs.data(), &n_int, &info);
122
123 for (Index i = 0; i < n; i++) {
124 x[i] = rhs[i];
125 }
126}
127
129
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
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 auto ipiv = std::vector<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.data(), &info);
181
182 // Invert matrix.
183 int lwork = n_int;
184 auto work = std::vector<double>(lwork);
185
186 lapack::dgetri_(&n_int, LU.mdata, &n_int, ipiv.data(), work.data(), &lwork, &info);
187
188 // Check for success.
189 ARTS_USER_ERROR_IF(info not_eq 0,
190 "Error inverting matrix: Matrix not of full rank.");
191 Ainv = LU;
192}
193
195 // A must be a square matrix.
196 ARTS_ASSERT(A.ncols() == A.nrows());
197
198 Index n = A.ncols();
199
200 // Workdata
201 Ainv = A;
202 int n_int = int(n), lwork = n_int, info;
203 std::vector<int> ipiv(n);
204 std::vector<Complex> work(lwork);
205
206 // Compute LU decomposition using LAPACK dgetrf_.
208 &n_int, &n_int, Ainv.get_c_array(), &n_int, ipiv.data(), &info);
209 lapack::zgetri_(&n_int,
210 Ainv.get_c_array(),
211 &n_int,
212 ipiv.data(),
213 work.data(),
214 &lwork,
215 &info);
216
217 // Check for success.
218 ARTS_USER_ERROR_IF(info not_eq 0,
219 "Error inverting matrix: Matrix not of full rank.");
220}
221
223
241 VectorView WR,
242 VectorView WI,
243 ConstMatrixView A) {
244 Index n = A.ncols();
245
246 // A must be a square matrix.
247 ARTS_ASSERT(n == A.nrows());
248 ARTS_ASSERT(n == WR.nelem());
249 ARTS_ASSERT(n == WI.nelem());
250 ARTS_ASSERT(n == P.nrows());
251 ARTS_ASSERT(n == P.ncols());
252
253 Matrix A_tmp = A;
254 Matrix P2 = P;
255 Vector WR2 = WR;
256 Vector WI2 = WI;
257
258 // Integers
259 int LDA, LDA_L, LDA_R, n_int, info;
260 n_int = (int)n;
261 LDA = LDA_L = LDA_R = (int)A.mcr.get_extent();
262
263 // We want to calculate RP not LP
264 char l_eig = 'N', r_eig = 'V';
265
266 // Work matrix
267 int lwork = 2 * n_int + n_int * n_int;
268 auto work = std::vector<double>(lwork);
269 auto rwork = std::vector<double>(2 * n_int);
270
271 // Memory references
272 double* adata = A_tmp.mdata;
273 double* rpdata = P2.mdata;
274 double* wrdata = WR2.mdata;
275 double* widata = WI2.mdata;
276
277 // Main calculations. Note that errors in the output is ignored
278 lapack::dgeev_(&l_eig,
279 &r_eig,
280 &n_int,
281 adata,
282 &LDA,
283 wrdata,
284 widata,
285 nullptr,
286 &LDA_L,
287 rpdata,
288 &LDA_R,
289 work.data(),
290 &lwork,
291 rwork.data(),
292 &info);
293
294 // Re-order. This can be done better?
295 for (Index i = 0; i < n; i++)
296 for (Index j = 0; j < n; j++) P(j, i) = P2(i, j);
297 WI = WI2;
298 WR = WR2;
299}
300
302
315 const ConstComplexMatrixView A) {
316 Index n = A.ncols();
317
318 // A must be a square matrix.
319 ARTS_ASSERT(n == A.nrows());
320 ARTS_ASSERT(n == W.nelem());
321 ARTS_ASSERT(n == P.nrows());
322 ARTS_ASSERT(n == P.ncols());
323
324 ComplexMatrix A_tmp = A;
325
326 // Integers
327 int LDA = int(A.get_column_extent()), LDA_L = int(A.get_column_extent()),
328 LDA_R = int(A.get_column_extent()), n_int = int(n), info;
329
330 // We want to calculate RP not LP
331 char l_eig = 'N', r_eig = 'V';
332
333 // Work matrix
334 int lwork = 2 * n_int + n_int * n_int;
335 ComplexVector work(lwork);
336 ComplexVector lpdata(0);
337 Vector rwork(2 * n_int);
338
339 // Main calculations. Note that errors in the output is ignored
340 lapack::zgeev_(&l_eig,
341 &r_eig,
342 &n_int,
343 A_tmp.get_c_array(),
344 &LDA,
345 W.get_c_array(),
346 lpdata.get_c_array(),
347 &LDA_L,
348 P.get_c_array(),
349 &LDA_R,
350 work.get_c_array(),
351 &lwork,
352 rwork.get_c_array(),
353 &info);
354
355 for (Index i = 0; i < n; i++)
356 for (Index j = 0; j <= i; j++) std::swap(P(j, i), P(i, j));
357}
358
360
375 const Index n = A.ncols();
376
377 /* Check if A and F are a quadratic and of the same dimension. */
378 ARTS_ASSERT(is_size(A, n, n));
379 ARTS_ASSERT(is_size(F, n, n));
380
381 Numeric A_norm_inf, c;
382 Numeric j;
383 Matrix D(n, n), N(n, n), X(n, n), cX(n, n, 0.0), B(n, n, 0.0);
384 Vector N_col_vec(n, 0.), F_col_vec(n, 0.);
385
386 A_norm_inf = norm_inf(A);
387
388 // This formular is derived in the book by Golub and Van Loan.
389 j = 1 + floor(1. / log(2.) * log(A_norm_inf));
390
391 if (j < 0) j = 0.;
392 auto j_index = (Index)(j);
393
394 // Scale matrix
395 F = A;
396 F /= pow(2, j);
397
398 /* The higher q the more accurate is the computation,
399 see user guide for accuracy */
400 // Index q = 8;
401 auto q_n = (Numeric)(q);
402 id_mat(D);
403 id_mat(N);
404 id_mat(X);
405 c = 1.;
406
407 for (Index k = 0; k < q; k++) {
408 auto k_n = (Numeric)(k + 1);
409 c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
410 mult(B, F, X); // X = F * X
411 X = B;
412 cX = X;
413 cX *= c; // cX = X*c
414 N += cX; // N = N + X*c
415 cX *= pow(-1, k_n); // cX = (-1)^k*c*X
416 D += cX; // D = D + (-1)^k*c*X
417 }
418
419 /*Solve the equation system DF=N for F using LU decomposition,
420 use the backsubstitution routine for columns of N*/
421
422 /* Now use X for the LU decomposition matrix of D.*/
423 ArrayOfIndex indx(n);
424
425 ludcmp(X, indx, D);
426
427 for (Index i = 0; i < n; i++) {
428 N_col_vec = N(joker, i); // extract column vectors of N
429 lubacksub(F_col_vec, X, N_col_vec, indx);
430 F(joker, i) = F_col_vec; // construct F matrix from column vectors
431 }
432
433 /* The square of F gives the result. */
434 for (Index k = 0; k < j_index; k++) {
435 mult(B, F, F); // F = F^2
436 F = B;
437 }
438}
439
441
450 Numeric norm_inf = 0;
451
452 for (Index j = 0; j < A.nrows(); j++) {
453 Numeric row_sum = 0;
454 //Calculate the row sum for all rows
455 for (Index i = 0; i < A.ncols(); i++) row_sum += abs(A(i, j));
456 //Pick out the row with the highest row sum
457 if (norm_inf < row_sum) norm_inf = row_sum;
458 }
459 return norm_inf;
460}
461
463
467 const Index n = I.ncols();
468 ARTS_ASSERT(n == I.nrows());
469
470 I = 0;
471 for (Index i = 0; i < n; i++) I(i, i) = 1.;
472}
473
483 const Index dim = A.nrows();
484 ARTS_ASSERT(dim == A.ncols());
485
486 if (dim == 3)
487 return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +
488 A(0, 2) * A(1, 0) * A(2, 1) - A(0, 2) * A(1, 1) * A(2, 0) -
489 A(0, 1) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 2) * A(2, 1);
490 if (dim == 2) return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
491 if (dim == 1) return A(0, 0);
492
493 Numeric ret_val = 0.;
494
495 for (Index j = 0; j < dim; j++) {
496 Matrix temp(dim - 1, dim - 1);
497 for (Index I = 1; I < dim; I++)
498 for (Index J = 0; J < dim; J++) {
499 if (J < j)
500 temp(I - 1, J) = A(I, J);
501 else if (J > j)
502 temp(I - 1, J - 1) = A(I, J);
503 }
504
505 Numeric tempNum = det(temp);
506
507 ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
508 }
509 return ret_val;
510}
511
527 const Index n = x.nelem();
528
529 ARTS_ASSERT(y.nelem() == n);
530
531 p.resize(2);
532
533 // Basic algorithm found at e.g.
534 // http://en.wikipedia.org/wiki/Simple_linear_regression
535 // The basic algorithm is as follows:
536 /*
537 Numeric s1=0, s2=0, s3=0, s4=0;
538 for( Index i=0; i<n; i++ )
539 {
540 s1 += x[i] * y[i];
541 s2 += x[i];
542 s3 += y[i];
543 s4 += x[i] * x[i];
544 }
545
546 p[1] = ( s1 - (s2*s3)/n ) / ( s4 - s2*s2/n );
547 p[0] = s3/n - p[1]*s2/n;
548 */
549
550 // A version abit more numerical stable:
551 // Mean value of x is removed before the fit: x' = (x-mean(x))
552 // This corresponds to that s2 in version above becomes 0
553 // y = a + b*x'
554 // p[1] = b
555 // p[0] = a - p[1]*mean(x)
556
557 Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
558
559 for (Index i = 0; i < n; i++) {
560 xm += x[i] / Numeric(n);
561 }
562
563 for (Index i = 0; i < n; i++) {
564 const Numeric xv = x[i] - xm;
565 s1 += xv * y[i];
566 s3 += y[i];
567 s4 += xv * xv;
568 }
569
570 p[1] = s1 / s4;
571 p[0] = s3 / Numeric(n) - p[1] * xm;
572}
573
577 bool residual) noexcept {
578 // Size of the problem
579 const Index n = x.nelem();
580 Matrix AT, ATA(n, n);
581 Vector ATy(n);
582
583 // Solver
584 AT = transpose(A);
585 mult(ATA, AT, A);
586 mult(ATy, AT, y);
587 solve(x, ATA, ATy);
588
589 // Residual
590 if (residual) {
591 Vector r(n);
592 mult(r, ATA, x);
593 r -= ATy;
594 return r * r;
595 }
596
597 return 0;
598}
This file contains the definition of Array.
The global header file for ARTS.
Header file for helper functions for OpenMP.
The ComplexMatrixView class.
The ComplexMatrix class.
The ComplexVectorView class.
The ComplexVector class.
A constant view of a ComplexMatrix.
Index ncols() const noexcept
Index nrows() const noexcept
Index get_column_extent() const
Get the extent of the underlying data.
Complex * get_c_array() const noexcept
Complex * get_c_array() const noexcept
Index nelem() const noexcept
A constant view of a Matrix.
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1176
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1174
A constant view of a Vector.
Definition: matpackI.h:521
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:660
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Definition: matpackI.h:633
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The Matrix class.
Definition: matpackI.h:1285
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:343
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:345
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define DEBUG_ONLY(...)
Definition: debug.h:71
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
Interface for the LAPACK library.
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:526
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:466
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Definition: lin_alg.cc:138
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
Definition: lin_alg.cc:449
void diagonalize(MatrixView P, VectorView WR, VectorView WI, ConstMatrixView A)
Matrix Diagonalization.
Definition: lin_alg.cc:240
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:91
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:374
void ludcmp(Matrix &LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:57
Numeric det(ConstMatrixView A)
Definition: lin_alg.cc:482
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:167
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:574
Linear algebra functions.
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:82
Header file for logic.cc.
void abs(Sparse &A, const Sparse &B)
Absolute value of sparse matrix elements.
Definition: matpackII.cc:394
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
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
const Joker joker
void zgetri_(int *n, std::complex< double > *A, int *lda, int *ipiv, std::complex< double > *work, int *lwork, int *info)
void zgeev_(char *jobvl, char *jobvr, int *n, std::complex< double > *A, int *lda, std::complex< double > *W, std::complex< double > *VL, int *ldvl, std::complex< double > *VR, int *ldvr, std::complex< double > *work, int *lwork, double *rwork, int *info)
void dgetri_(int *n, double *A, int *lda, int *ipiv, double *work, int *lwork, int *info)
Matrix inversion.
void dgeev_(char *jobvl, char *jobvr, int *n, double *A, int *lda, double *WR, double *WI, double *VL, int *ldvl, double *VR, int *ldvr, double *work, int *lwork, double *rwork, int *info)
void zgetrf_(int *m, int *n, std::complex< double > *A, int *lda, int *ipiv, int *info)
void dgetrs_(char *trans, int *n, int *nrhs, double *A, int *lda, int *ipiv, double *b, int *ldb, int *info)
Solve linear system of equations.
void dgetrf_(int *m, int *n, double *A, int *lda, int *ipiv, int *info)
LU decomposition.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
#define N
Definition: rng.cc:164
#define c
#define b