ARTS 2.5.4 (git: 31ce4f0e)
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
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 int* ipiv = new 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, &info);
73
74 // Copy pivot array to pivot vector.
75 for (Index i = 0; i < n; i++) {
76 indx[i] = ipiv[i];
77 }
78 delete[] ipiv;
79}
80
82
95 const ArrayOfIndex& indx) {
96 Index n = LU.nrows();
97
98 /* Check if the dimensions of the input matrix and vectors agree and if LU
99 is a quadratic matrix.*/
100 DEBUG_ONLY(Index column_stride = LU.mcr.get_stride());
101 DEBUG_ONLY(Index vec_stride = b.mrange.get_stride());
102
103 ARTS_ASSERT(is_size(LU, n, n));
104 ARTS_ASSERT(is_size(b, n));
105 ARTS_ASSERT(is_size(indx, n));
106 ARTS_ASSERT(column_stride == 1);
107 ARTS_ASSERT(vec_stride == 1);
108
109 char trans = 'N';
110 int n_int = (int)n;
111 int one = (int)1;
112 std::vector<int> ipiv(n);
113 std::vector<double> rhs(n);
114 int info;
115
116 for (Index i = 0; i < n; i++) {
117 ipiv[i] = (int)indx[i];
118 rhs[i] = b[i];
119 }
120
122 &trans, &n_int, &one, LU.mdata, &n_int, ipiv.data(), rhs.data(), &n_int, &info);
123
124 for (Index i = 0; i < n; i++) {
125 x[i] = rhs[i];
126 }
127}
128
130
140 Index n = A.ncols();
141
142 // Check dimensions of the system.
143 ARTS_ASSERT(n == A.nrows());
144 ARTS_ASSERT(n == x.nelem());
145 ARTS_ASSERT(n == b.nelem());
146
147 // Allocate matrix and index vector for the LU decomposition.
148 Matrix LU = Matrix(n, n);
149 ArrayOfIndex indx(n);
150
151 // Perform LU decomposition.
152 ludcmp(LU, indx, A);
153
154 // Solve the system using backsubstitution.
155 lubacksub(x, LU, b, indx);
156}
157
159
169 // A must be a square matrix.
170 ARTS_ASSERT(A.ncols() == A.nrows());
171
172 Index n = A.ncols();
173 Matrix LU(A);
174
175 int n_int, info;
176 int* ipiv = new int[n];
177
178 n_int = (int)n;
179
180 // Compute LU decomposition using LAPACK dgetrf_.
181 lapack::dgetrf_(&n_int, &n_int, LU.mdata, &n_int, ipiv, &info);
182
183 // Invert matrix.
184 int lwork = n_int;
185 double* work;
186
187 try {
188 work = new double[lwork];
189 } catch (std::bad_alloc& ba) {
191 "Error inverting matrix: Could not allocate workspace memory.");
192 }
193
194 lapack::dgetri_(&n_int, LU.mdata, &n_int, ipiv, work, &lwork, &info);
195 delete[] work;
196 delete[] ipiv;
197
198 // Check for success.
199 ARTS_USER_ERROR_IF(info not_eq 0,
200 "Error inverting matrix: Matrix not of full rank.");
201 Ainv = LU;
202}
203
205 // A must be a square matrix.
206 ARTS_ASSERT(A.ncols() == A.nrows());
207
208 Index n = A.ncols();
209
210 // Workdata
211 Ainv = A;
212 int n_int = int(n), lwork = n_int, info;
213 std::vector<int> ipiv(n);
214 ComplexVector work(lwork);
215
216 // Compute LU decomposition using LAPACK dgetrf_.
218 &n_int, &n_int, Ainv.get_c_array(), &n_int, ipiv.data(), &info);
219 lapack::zgetri_(&n_int,
220 Ainv.get_c_array(),
221 &n_int,
222 ipiv.data(),
223 work.get_c_array(),
224 &lwork,
225 &info);
226
227 // Check for success.
228 ARTS_USER_ERROR_IF(info not_eq 0,
229 "Error inverting matrix: Matrix not of full rank.");
230}
231
233
251 VectorView WR,
252 VectorView WI,
253 ConstMatrixView A) {
254 Index n = A.ncols();
255
256 // A must be a square matrix.
257 ARTS_ASSERT(n == A.nrows());
258 ARTS_ASSERT(n == WR.nelem());
259 ARTS_ASSERT(n == WI.nelem());
260 ARTS_ASSERT(n == P.nrows());
261 ARTS_ASSERT(n == P.ncols());
262
263 Matrix A_tmp = A;
264 Matrix P2 = P;
265 Vector WR2 = WR;
266 Vector WI2 = WI;
267
268 // Integers
269 int LDA, LDA_L, LDA_R, n_int, info;
270 n_int = (int)n;
271 LDA = LDA_L = LDA_R = (int)A.mcr.get_extent();
272
273 // We want to calculate RP not LP
274 char l_eig = 'N', r_eig = 'V';
275
276 // Work matrix
277 int lwork = 2 * n_int + n_int * n_int;
278 double *work, *rwork;
279 try {
280 rwork = new double[2 * n_int];
281 work = new double[lwork];
282 } catch (std::bad_alloc& ba) {
284 "Error diagonalizing: Could not allocate workspace memory.");
285 }
286
287 // Memory references
288 double* adata = A_tmp.mdata;
289 double* rpdata = P2.mdata;
290 auto* lpdata = new double[0]; //To not confuse the compiler
291 double* wrdata = WR2.mdata;
292 double* widata = WI2.mdata;
293
294 // Main calculations. Note that errors in the output is ignored
295 lapack::dgeev_(&l_eig,
296 &r_eig,
297 &n_int,
298 adata,
299 &LDA,
300 wrdata,
301 widata,
302 lpdata,
303 &LDA_L,
304 rpdata,
305 &LDA_R,
306 work,
307 &lwork,
308 rwork,
309 &info);
310
311 // Free memory. Can these be sent in to speed things up?
312 delete[] work;
313 delete[] rwork;
314 delete[] lpdata;
315
316 // Re-order. This can be done better?
317 for (Index i = 0; i < n; i++)
318 for (Index j = 0; j < n; j++) P(j, i) = P2(i, j);
319 WI = WI2;
320 WR = WR2;
321}
322
324
337 const ConstComplexMatrixView A) {
338 Index n = A.ncols();
339
340 // A must be a square matrix.
341 ARTS_ASSERT(n == A.nrows());
342 ARTS_ASSERT(n == W.nelem());
343 ARTS_ASSERT(n == P.nrows());
344 ARTS_ASSERT(n == P.ncols());
345
346 ComplexMatrix A_tmp = A;
347
348 // Integers
349 int LDA = int(A.get_column_extent()), LDA_L = int(A.get_column_extent()),
350 LDA_R = int(A.get_column_extent()), n_int = int(n), info;
351
352 // We want to calculate RP not LP
353 char l_eig = 'N', r_eig = 'V';
354
355 // Work matrix
356 int lwork = 2 * n_int + n_int * n_int;
357 ComplexVector work(lwork);
358 ComplexVector lpdata(0);
359 Vector rwork(2 * n_int);
360
361 // Main calculations. Note that errors in the output is ignored
362 lapack::zgeev_(&l_eig,
363 &r_eig,
364 &n_int,
365 A_tmp.get_c_array(),
366 &LDA,
367 W.get_c_array(),
368 lpdata.get_c_array(),
369 &LDA_L,
370 P.get_c_array(),
371 &LDA_R,
372 work.get_c_array(),
373 &lwork,
374 rwork.get_c_array(),
375 &info);
376
377 for (Index i = 0; i < n; i++)
378 for (Index j = 0; j <= i; j++) std::swap(P(j, i), P(i, j));
379}
380
382
397 const Index n = A.ncols();
398
399 /* Check if A and F are a quadratic and of the same dimension. */
400 ARTS_ASSERT(is_size(A, n, n));
401 ARTS_ASSERT(is_size(F, n, n));
402
403 Numeric A_norm_inf, c;
404 Numeric j;
405 Matrix D(n, n), N(n, n), X(n, n), cX(n, n, 0.0), B(n, n, 0.0);
406 Vector N_col_vec(n, 0.), F_col_vec(n, 0.);
407
408 A_norm_inf = norm_inf(A);
409
410 // This formular is derived in the book by Golub and Van Loan.
411 j = 1 + floor(1. / log(2.) * log(A_norm_inf));
412
413 if (j < 0) j = 0.;
414 auto j_index = (Index)(j);
415
416 // Scale matrix
417 F = A;
418 F /= pow(2, j);
419
420 /* The higher q the more accurate is the computation,
421 see user guide for accuracy */
422 // Index q = 8;
423 auto q_n = (Numeric)(q);
424 id_mat(D);
425 id_mat(N);
426 id_mat(X);
427 c = 1.;
428
429 for (Index k = 0; k < q; k++) {
430 auto k_n = (Numeric)(k + 1);
431 c *= (q_n - k_n + 1) / ((2 * q_n - k_n + 1) * k_n);
432 mult(B, F, X); // X = F * X
433 X = B;
434 cX = X;
435 cX *= c; // cX = X*c
436 N += cX; // N = N + X*c
437 cX *= pow(-1, k_n); // cX = (-1)^k*c*X
438 D += cX; // D = D + (-1)^k*c*X
439 }
440
441 /*Solve the equation system DF=N for F using LU decomposition,
442 use the backsubstitution routine for columns of N*/
443
444 /* Now use X for the LU decomposition matrix of D.*/
445 ArrayOfIndex indx(n);
446
447 ludcmp(X, indx, D);
448
449 for (Index i = 0; i < n; i++) {
450 N_col_vec = N(joker, i); // extract column vectors of N
451 lubacksub(F_col_vec, X, N_col_vec, indx);
452 F(joker, i) = F_col_vec; // construct F matrix from column vectors
453 }
454
455 /* The square of F gives the result. */
456 for (Index k = 0; k < j_index; k++) {
457 mult(B, F, F); // F = F^2
458 F = B;
459 }
460}
461
463
472 Numeric norm_inf = 0;
473
474 for (Index j = 0; j < A.nrows(); j++) {
475 Numeric row_sum = 0;
476 //Calculate the row sum for all rows
477 for (Index i = 0; i < A.ncols(); i++) row_sum += abs(A(i, j));
478 //Pick out the row with the highest row sum
479 if (norm_inf < row_sum) norm_inf = row_sum;
480 }
481 return norm_inf;
482}
483
485
489 const Index n = I.ncols();
490 ARTS_ASSERT(n == I.nrows());
491
492 I = 0;
493 for (Index i = 0; i < n; i++) I(i, i) = 1.;
494}
495
505 const Index dim = A.nrows();
506 ARTS_ASSERT(dim == A.ncols());
507
508 if (dim == 3)
509 return A(0, 0) * A(1, 1) * A(2, 2) + A(0, 1) * A(1, 2) * A(2, 0) +
510 A(0, 2) * A(1, 0) * A(2, 1) - A(0, 2) * A(1, 1) * A(2, 0) -
511 A(0, 1) * A(1, 0) * A(2, 2) - A(0, 0) * A(1, 2) * A(2, 1);
512 if (dim == 2) return A(0, 0) * A(1, 1) - A(0, 1) * A(1, 0);
513 if (dim == 1) return A(0, 0);
514
515 Numeric ret_val = 0.;
516
517 for (Index j = 0; j < dim; j++) {
518 Matrix temp(dim - 1, dim - 1);
519 for (Index I = 1; I < dim; I++)
520 for (Index J = 0; J < dim; J++) {
521 if (J < j)
522 temp(I - 1, J) = A(I, J);
523 else if (J > j)
524 temp(I - 1, J - 1) = A(I, J);
525 }
526
527 Numeric tempNum = det(temp);
528
529 ret_val += ((j % 2 == 0) ? -1. : 1.) * tempNum * A(0, j);
530 }
531 return ret_val;
532}
533
549 const Index n = x.nelem();
550
551 ARTS_ASSERT(y.nelem() == n);
552
553 p.resize(2);
554
555 // Basic algorithm found at e.g.
556 // http://en.wikipedia.org/wiki/Simple_linear_regression
557 // The basic algorithm is as follows:
558 /*
559 Numeric s1=0, s2=0, s3=0, s4=0;
560 for( Index i=0; i<n; i++ )
561 {
562 s1 += x[i] * y[i];
563 s2 += x[i];
564 s3 += y[i];
565 s4 += x[i] * x[i];
566 }
567
568 p[1] = ( s1 - (s2*s3)/n ) / ( s4 - s2*s2/n );
569 p[0] = s3/n - p[1]*s2/n;
570 */
571
572 // A version abit more numerical stable:
573 // Mean value of x is removed before the fit: x' = (x-mean(x))
574 // This corresponds to that s2 in version above becomes 0
575 // y = a + b*x'
576 // p[1] = b
577 // p[0] = a - p[1]*mean(x)
578
579 Numeric s1 = 0, xm = 0, s3 = 0, s4 = 0;
580
581 for (Index i = 0; i < n; i++) {
582 xm += x[i] / Numeric(n);
583 }
584
585 for (Index i = 0; i < n; i++) {
586 const Numeric xv = x[i] - xm;
587 s1 += xv * y[i];
588 s3 += y[i];
589 s4 += xv * xv;
590 }
591
592 p[1] = s1 / s4;
593 p[0] = s3 / Numeric(n) - p[1] * xm;
594}
595
599 bool residual) noexcept {
600 // Size of the problem
601 const Index n = x.nelem();
602 Matrix AT, ATA(n, n);
603 Vector ATy(n);
604
605 // Solver
606 AT = transpose(A);
607 mult(ATA, AT, A);
608 mult(ATy, AT, y);
609 solve(x, ATA, ATy);
610
611 // Residual
612 if (residual) {
613 Vector r(n);
614 mult(r, ATA, x);
615 r -= ATy;
616 return r * r;
617 }
618
619 return 0;
620}
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:1043
Index nrows() const noexcept
Definition: matpackI.h:1055
Index ncols() const noexcept
Definition: matpackI.h:1056
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:1152
Range mcr
The column range of mdata that is actually used.
Definition: matpackI.h:1150
A constant view of a Vector.
Definition: matpackI.h:512
Numeric * mdata
Pointer to the plain C array that holds the data.
Definition: matpackI.h:649
Numeric * get_c_array() const noexcept
Conversion to plain C-array, const-version.
Definition: matpackI.h:622
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:536
The MatrixView class.
Definition: matpackI.h:1164
The Matrix class.
Definition: matpackI.h:1261
constexpr Index get_extent() const noexcept
Returns the extent of the range.
Definition: matpackI.h:342
constexpr Index get_stride() const noexcept
Returns the stride of the range.
Definition: matpackI.h:344
The VectorView class.
Definition: matpackI.h:663
The Vector class.
Definition: matpackI.h:899
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:52
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
Interface for the LAPACK library.
#define abs(x)
#define q
#define temp
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:548
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:488
void solve(VectorView x, ConstMatrixView A, ConstVectorView b)
Solve a linear system.
Definition: lin_alg.cc:139
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
Definition: lin_alg.cc:471
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
Numeric det(ConstMatrixView A)
Definition: lin_alg.cc:504
void inv(MatrixView Ainv, ConstMatrixView A)
Matrix Inverse.
Definition: lin_alg.cc:168
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:596
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.
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.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
const Joker joker
constexpr Numeric k
Boltzmann constant convenience name [J/K].
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.
invlib::Matrix< ArtsMatrix > Matrix
invlib wrapper type for ARTS matrices.
Definition: oem.h:33
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
#define N
Definition: rng.cc:164
#define c
#define b