ARTS  2.2.66
lin_alg.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 
30 #include "lin_alg.h"
31 
32 /*===========================================================================
33  === External declarations
34  ===========================================================================*/
35 
36 #include <stdexcept>
37 #include <cmath>
38 #include "arts.h"
39 #include "make_vector.h"
40 #include "array.h"
41 #include "logic.h"
42 
43 
44 
46 
55  ArrayOfIndex& indx,
56  ConstMatrixView A)
57 {
58  Index imax = 0;
59  const Numeric TINY=1.0e-20;
60  Numeric big, dum, sum, temp, d;
61 
62  LU = A;
63  const Index dim = A.nrows();
64  assert(is_size(A,dim,dim)); //check, if A is quadratic
65 
66  Vector vv(dim); //stores implicit scaling of each row
67  d = 1.0;
68  for (Index i=0; i<dim; i++)
69  {
70  indx[i]= i;
71  big = 0.0;
72  for (Index j=0; j<dim; j++)
73  {
74  if ((temp = abs(LU(i,j))) > big)
75  big = temp;
76  }
77  if (big == 0.)
78  throw runtime_error("ludcmp: Matrix is singular");
79  vv[i] = 1.0/big; // save scaling
80  }
81 
82  for (Index j=0; j<dim; j++)
83  {
84  for (Index i=0; i<j; i++)
85  {
86  sum = LU(i,j);
87  for (Index k=0; k<i; k++)
88  sum -= LU(i,k)*LU(k,j);
89  LU(i,j) = sum;
90  }
91  big = 0.0;
92  for( Index i=j; i<dim; i++)
93  {
94  sum = LU(i,j);
95  for (Index k=0; k<j; k++)
96  sum -= LU(i,k)*LU(k,j);
97  LU(i,j) = sum;
98  if( (dum = vv[i]*fabs(sum)) >= big)
99  {
100  big = dum;
101  imax = i;
102  }
103  }
104  if (j!=imax)
105  {
106  for(Index k=0; k<dim; k++)
107  {
108  dum = LU(imax,k);
109  LU(imax,k) = LU(j,k);
110  LU(j,k) = dum;
111  }
112  d = -d;
113  vv[imax] = vv[j];
114  indx[j] = imax;
115  indx[imax] =j;
116  }
117 
118  if(LU(j,j) == 0.0) LU(j,j) = TINY;
119 
120  if (j != dim)
121  {
122  dum=1.0/LU(j,j);
123  for (Index i=j+1; i<dim; i++)
124  LU(i,j) *=dum;
125  }
126  }
127 }
128 
129 
130 
132 
143  ConstMatrixView LU,
144  ConstVectorView b,
145  const ArrayOfIndex& indx)
146 {
147  Index dim = LU.nrows();
148 
149  /* Check if the dimensions of the input matrix and vectors agree and if LU is a quadratic matrix.*/
150  assert(is_size(LU, dim, dim));
151  assert(is_size(b, dim));
152  assert(is_size(indx, dim));
153 
154 
155  for(Index i=0; i<dim; i++)
156  {
157  x[indx[i]] = b[i];
158  }
159 
160  for (Index i=0; i<dim; i++)
161  {
162  Numeric sum = x[i];
163  for (Index j=0; j<=i-1; j++)
164  sum -= LU(i,j)*x[j];
165  x[i] = sum;
166  }
167 
168  for(Index i=dim-1; i>=0; i--)
169  {
170  Numeric sum = x[i];
171  for (Index j=i+1; j<dim; j++)
172  sum -= LU(i,j)*x[j];
173  x[i] = sum/LU(i,i);
174  }
175 }
176 
177 
178 
180 
195  MatrixView F,
196  ConstMatrixView A,
197  const Index& q )
198 {
199  const Index n = A.ncols();
200 
201  /* Check if A and F are a quadratic and of the same dimension. */
202  assert( is_size(A,n,n) );
203  assert( is_size(F,n,n) );
204 
205  Numeric A_norm_inf, c;
206  Numeric j;
207  Matrix D(n,n), N(n,n), X(n,n), cX(n,n,0.0), B(n,n,0.0);
208  Vector N_col_vec(n,0.), F_col_vec(n,0.);
209 
210  A_norm_inf = norm_inf(A);
211 
212  // This formular is derived in the book by Golub and Van Loan.
213  j = 1 + floor(1./log(2.)*log(A_norm_inf));
214 
215  if(j<0) j=0.;
216  Index j_index = (Index)(j);
217 
218  // Scale matrix
219  F = A;
220  F /= pow(2,j);
221 
222  /* The higher q the more accurate is the computation,
223  see user guide for accuracy */
224  // Index q = 8;
225  Numeric q_n = (Numeric)(q);
226  id_mat(D);
227  id_mat(N);
228  id_mat(X);
229  c = 1.;
230 
231  for(Index k=0; k<q; k++)
232  {
233  Numeric k_n = (Numeric)(k+1);
234  c *= (q_n-k_n+1)/((2*q_n-k_n+1)*k_n);
235  mult(B,F,X); // X = F * X
236  X = B;
237  cX = X;
238  cX *= c; // cX = X*c
239  N += cX; // N = N + X*c
240  cX *= pow(-1,k_n); // cX = (-1)^k*c*X
241  D += cX; // D = D + (-1)^k*c*X
242  }
243 
244  /*Solve the equation system DF=N for F using LU decomposition,
245  use the backsubstitution routine for columns of N*/
246 
247  /* Now use X for the LU decomposition matrix of D.*/
248  ArrayOfIndex indx(n);
249 
250  ludcmp(X, indx, D);
251 
252  for(Index i=0; i<n; i++)
253  {
254  N_col_vec = N(joker,i); // extract column vectors of N
255  lubacksub(F_col_vec, X, N_col_vec, indx);
256  F(joker,i) = F_col_vec; // construct F matrix from column vectors
257  }
258 
259  /* The square of F gives the result. */
260  for(Index k=0; k<j_index; k++)
261  {
262  mult(B,F,F); // F = F^2
263  F = B;
264  }
265 }
266 
267 
269 
278 {
279  Numeric norm_inf = 0;
280 
281  for(Index j=0; j<A.nrows(); j++)
282  {
283  Numeric row_sum = 0;
284  //Calculate the row sum for all rows
285  for(Index i=0; i<A.ncols(); i++)
286  row_sum += abs(A(i,j));
287  //Pick out the row with the highest row sum
288  if( norm_inf < row_sum)
289  norm_inf = row_sum;
290  }
291  return norm_inf;
292 }
293 
294 
296 
300 {
301 
302  const Index n = I.ncols();
303  assert(n == I.nrows());
304 
305  I = 0;
306  for(Index i=0; i<n; i++)
307  I(i,i) = 1.;
308 }
309 
310 
311 
321 {
322  const Index dim = A.nrows();
323  assert(dim == A.ncols());
324 
325  if(dim == 3)
326  return A(0,0)*A(1,1)*A(2,2) + A(0,1)*A(1,2)*A(2,0) +
327  A(0,2)*A(1,0)*A(2,1) - A(0,2)*A(1,1)*A(2,0) -
328  A(0,1)*A(1,0)*A(2,2) - A(0,0)*A(1,2)*A(2,1);
329  else if(dim == 2)
330  return A(0,0) * A(1,1) - A(0,1) * A(1,0);
331  else if(dim == 1)
332  return A(0,0);
333 
334  Numeric ret_val = 0.;
335 
336  for(Index j = 0; j < dim; j++)
337  {
338  Matrix temp(dim-1,dim-1);
339  for(Index I = 1; I < dim; I++)
340  for(Index J = 0; J < dim; J++)
341  {
342  if(J < j)
343  temp(I-1,J) = A(I,J);
344  else if(J > j)
345  temp(I-1,J-1) = A(I,J);
346  }
347 
348  Numeric tempNum = det(temp);
349 
350  ret_val += ((j%2 == 0)?-1.:1.) * tempNum * A(0,j);
351  }
352  return ret_val;
353 }
354 
355 
356 
371 void linreg(
372  Vector& p,
373  ConstVectorView x,
374  ConstVectorView y )
375 {
376  const Index n = x.nelem();
377 
378  assert( y.nelem() == n );
379 
380  p.resize(2);
381 
382  // Basic algorithm found at e.g.
383  // http://en.wikipedia.org/wiki/Simple_linear_regression
384  // The basic algorithm is as follows:
385  /*
386  Numeric s1=0, s2=0, s3=0, s4=0;
387  for( Index i=0; i<n; i++ )
388  {
389  s1 += x[i] * y[i];
390  s2 += x[i];
391  s3 += y[i];
392  s4 += x[i] * x[i];
393  }
394 
395  p[1] = ( s1 - (s2*s3)/n ) / ( s4 - s2*s2/n );
396  p[0] = s3/n - p[1]*s2/n;
397  */
398 
399  // A version abit more numerical stable:
400  // Mean value of x is removed before the fit: x' = (x-mean(x))
401  // This corresponds to that s2 in version above becomes 0
402  // y = a + b*x'
403  // p[1] = b
404  // p[0] = a - p[1]*mean(x)
405 
406  Numeric s1=0, xm=0, s3=0, s4=0;
407 
408  for( Index i=0; i<n; i++ )
409  { xm += x[i]/Numeric(n); }
410 
411  for( Index i=0; i<n; i++ )
412  {
413  const Numeric xv = x[i] - xm;
414  s1 += xv * y[i];
415  s3 += y[i];
416  s4 += xv * xv;
417  }
418 
419  p[1] = s1 / s4;
420  p[0] = s3/Numeric(n) - p[1]*xm;
421 }
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
#define N
Definition: rng.cc:195
The VectorView class.
Definition: matpackI.h:372
void matrix_exp(MatrixView F, ConstMatrixView A, const Index &q)
General exponential of a Matrix.
Definition: lin_alg.cc:194
void lubacksub(VectorView x, ConstMatrixView LU, ConstVectorView b, const ArrayOfIndex &indx)
LU backsubstitution.
Definition: lin_alg.cc:142
#define q
Definition: continua.cc:21469
The Vector class.
Definition: matpackI.h:556
The MatrixView class.
Definition: matpackI.h:679
Linear algebra functions.
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
This file contains the definition of Array.
void mult(VectorView y, const ConstMatrixView &M, const ConstVectorView &x)
Matrix Vector multiplication.
Definition: matpackI.cc:1648
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:838
The global header file for ARTS.
void ludcmp(MatrixView LU, ArrayOfIndex &indx, ConstMatrixView A)
LU decomposition.
Definition: lin_alg.cc:54
#define abs(x)
Definition: continua.cc:20458
const Joker joker
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
The Matrix class.
Definition: matpackI.h:788
Header file for logic.cc.
This can be used to make arrays out of anything.
Definition: array.h:40
The class MakeVector is a special kind of Vector that can be initialized explicitly from one or more ...
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
A constant view of a Vector.
Definition: matpackI.h:292
#define temp
Definition: continua.cc:20773
A constant view of a Matrix.
Definition: matpackI.h:596
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:91
void id_mat(MatrixView I)
Identity Matrix.
Definition: lin_alg.cc:299
void linreg(Vector &p, ConstVectorView x, ConstVectorView y)
Definition: lin_alg.cc:371
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:832
Numeric det(ConstMatrixView A)
Definition: lin_alg.cc:320
Numeric norm_inf(ConstMatrixView A)
Maximum absolute row sum norm.
Definition: lin_alg.cc:277