ReFRACtor
linear_algebra.h
Go to the documentation of this file.
1 #ifndef LINEAR_ALGEBRA_H
2 #define LINEAR_ALGEBRA_H
3 #include <blitz/array.h>
4 #include <limits>
5 
6 namespace FullPhysics {
9 
11 inline double sqr(double x) {return x * x;}
12 
13 //-----------------------------------------------------------------------
18 //-----------------------------------------------------------------------
19 
20 template<class T, int D> blitz::Array<T, D>
21 to_fortran(const blitz::Array<T, D> &In)
22 {
23  bool is_ok = In.isStorageContiguous();
24  for(int i = 0; i < D; ++i)
25  if(In.ordering(i) != i ||
26  !In.isRankStoredAscending(i))
27  is_ok = false;
28  if(is_ok)
29  return In;
30  else {
31  blitz::Array<T, D> res(In.shape(), blitz::ColumnMajorArray<D>());
32  res = In;
33  return res;
34  }
35 }
36 
37 //-----------------------------------------------------------------------
41 //-----------------------------------------------------------------------
42 
43 template<class T, int D> blitz::Array<T, D>
44 to_c_order(const blitz::Array<T, D> &In)
45 {
46  bool is_ok = In.isStorageContiguous();
47  for(int i = 0; i < D; ++i)
48  if(In.ordering(i) != D - i - 1||
49  !In.isRankStoredAscending(i))
50  is_ok = false;
51  if(is_ok)
52  return In;
53  else {
54  blitz::Array<T, D> res(In.shape());
55  res = In;
56  return res;
57  }
58 }
59 
60 //-----------------------------------------------------------------------
67 //-----------------------------------------------------------------------
68 
69 template<class T, int D> blitz::Array<T, D>
70 to_c_order_const(const blitz::Array<T, D> &In)
71 {
72  bool is_ok = In.isStorageContiguous();
73  for(int i = 0; i < D; ++i)
74  if(In.ordering(i) != D - i - 1||
75  !In.isRankStoredAscending(i))
76  is_ok = false;
77  if(is_ok)
78  return In;
79  else {
80  blitz::Array<T, D> res(In.shape());
81  res = In;
82  return res;
83  }
84 }
85 
86 //-----------------------------------------------------------------------
89 //-----------------------------------------------------------------------
90 
91 template<class T> void fortran_resize(blitz::Array<T, 1>& M, int sz1)
92 {
93  const int D = 1;
94  bool is_ok = M.isStorageContiguous();
95  for(int i = 0; i < D; ++i)
96  if(M.ordering(i) != i ||
97  !M.isRankStoredAscending(i))
98  is_ok = false;
99  if(is_ok)
100  M.resize(sz1);
101  else
102  M.reference(blitz::Array<T, D>(sz1, blitz::ColumnMajorArray<D>()));
103 }
104 
105 //-----------------------------------------------------------------------
108 //-----------------------------------------------------------------------
109 
110 template<class T> void fortran_resize(blitz::Array<T, 2>& M, int sz1, int sz2)
111 {
112  const int D = 2;
113  bool is_ok = M.isStorageContiguous();
114  for(int i = 0; i < D; ++i)
115  if(M.ordering(i) != i ||
116  !M.isRankStoredAscending(i))
117  is_ok = false;
118  if(is_ok)
119  M.resize(sz1, sz2);
120  else
121  M.reference(blitz::Array<T, D>(sz1, sz2, blitz::ColumnMajorArray<D>()));
122 }
123 
124 //-----------------------------------------------------------------------
127 //-----------------------------------------------------------------------
128 
129 template<class T> void fortran_resize(blitz::Array<T, 3>& M, int sz1, int sz2,
130  int sz3)
131 {
132  const int D = 3;
133  bool is_ok = M.isStorageContiguous();
134  for(int i = 0; i < D; ++i)
135  if(M.ordering(i) != i ||
136  !M.isRankStoredAscending(i))
137  is_ok = false;
138  if(is_ok)
139  M.resize(sz1, sz2, sz3);
140  else
141  M.reference(blitz::Array<T, D>(sz1, sz2, sz3,
142  blitz::ColumnMajorArray<D>()));
143 }
144 
145 //-----------------------------------------------------------------------
153 //-----------------------------------------------------------------------
154 
155 template<class T, int D> blitz::Array<T, D>
156 to_fortran_const(const blitz::Array<T, D> &In)
157 {
158  bool is_ok = In.isStorageContiguous();
159  for(int i = 0; i < D; ++i)
160  if(In.ordering(i) != i ||
161  !In.isRankStoredAscending(i))
162  is_ok = false;
163  if(is_ok)
164  return In;
165  else {
166  blitz::Array<T, D> res(In.shape(), blitz::ColumnMajorArray<D>());
167  res = In;
168  return res;
169  }
170 }
171 
172 blitz::Array<double, 1> solve_least_squares(const blitz::Array<double, 2>& A,
173  const blitz::Array<double, 1>& B,
174  double Rcond = 1e-12);
175 
176 blitz::Array<double, 1> solve_least_squares_qr(const blitz::Array<double, 2>& A,
177  const blitz::Array<double, 1>& B);
178 
179 void svd(const blitz::Array<double, 2>& A, blitz::Array<double, 1>& S,
180  blitz::Array<double, 2>& U, blitz::Array<double, 2>& VT);
181 
182 blitz::Array<double, 2> cholesky_decomposition(const blitz::Array<double, 2>& A);
183 
184 blitz::Array<double, 2> generalized_inverse(const blitz::Array<double, 2>& A,
185  double Rcond = std::numeric_limits<double>::epsilon());
186 blitz::Array<double, 2> generalized_inverse(const blitz::Array<double, 2>& A,
187  const blitz::Array<bool, 1>& Zero_unused,
188  double Rcond = std::numeric_limits<double>::epsilon());
189 
190 
191 blitz::Array<double, 2> multifit_covar(const blitz::Array<double, 2>& J,
192  double Eps_rel = std::numeric_limits<double>::epsilon());
193 bool multifit_test_gradient(const blitz::Array<double, 1>& g,
194  double Eps_abs = std::numeric_limits<double>::epsilon());
195 bool multifit_test_delta(const blitz::Array<double, 1>& dx, const blitz::Array<double, 1>& x,
196  double Eps_abs = std::numeric_limits<double>::epsilon(),
197  double Eps_rel = std::numeric_limits<double>::epsilon());
198 
200 }
201 
202 #endif
FullPhysics::svd
void svd(const blitz::Array< double, 2 > &A, blitz::Array< double, 1 > &S, blitz::Array< double, 2 > &U, blitz::Array< double, 2 > &VT)
Compute the SVD decomposition of a matrix A = U * SIGMA * V^T.
Definition: linear_algebra.cc:103
FullPhysics::solve_least_squares
blitz::Array< double, 1 > solve_least_squares(const blitz::Array< double, 2 > &A, const blitz::Array< double, 1 > &B, double Rcond=1e-12)
This solves the least squares system A*x = b, returning the minimum norm solution.
Definition: linear_algebra.cc:27
FullPhysics::to_fortran_const
blitz::Array< T, D > to_fortran_const(const blitz::Array< T, D > &In)
Ensure that a given blitz::Array is contiguous, not reversed, and in fortran ColumnMajorArray format...
Definition: linear_algebra.h:156
FullPhysics::fortran_resize
void fortran_resize(blitz::Array< T, 1 > &M, int sz1)
If matrix is in Fortran order, resize it.
Definition: linear_algebra.h:91
FullPhysics::to_c_order_const
blitz::Array< T, D > to_c_order_const(const blitz::Array< T, D > &In)
Ensure that a given blitz::Array is contiguous, not reversed, and in C RowMajorArray format...
Definition: linear_algebra.h:70
FullPhysics::cholesky_decomposition
blitz::Array< double, 2 > cholesky_decomposition(const blitz::Array< double, 2 > &A)
This calculates the Cholesky Decompostion of A so that A = L L^T.
Definition: linear_algebra.cc:212
FullPhysics::units::J
const Unit J("J", N *m)
FullPhysics::units::A
const Unit A("A", 1.0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0)
FullPhysics::solve_least_squares_qr
blitz::Array< double, 1 > solve_least_squares_qr(const blitz::Array< double, 2 > &A, const blitz::Array< double, 1 > &B)
This finds the least squares solution to the overdetermined system A x = b where the matrix A has mor...
Definition: linear_algebra.cc:66
FullPhysics::multifit_covar
blitz::Array< double, 2 > multifit_covar(const blitz::Array< double, 2 > &J, double Eps_rel=std::numeric_limits< double >::epsilon())
Definition: linear_algebra.cc:226
FullPhysics::multifit_test_gradient
bool multifit_test_gradient(const blitz::Array< double, 1 > &g, double Eps_abs=std::numeric_limits< double >::epsilon())
Definition: linear_algebra.cc:234
FullPhysics::generalized_inverse
blitz::Array< double, 2 > generalized_inverse(const blitz::Array< double, 2 > &A, double Rcond=std::numeric_limits< double >::epsilon())
This returns the generalized inverse of the given matrix.
Definition: linear_algebra.cc:136
FullPhysics::to_fortran
blitz::Array< T, D > to_fortran(const blitz::Array< T, D > &In)
Ensure that a given blitz::Array is contiguous, not reversed, and in fortran ColumnMajorArray format...
Definition: linear_algebra.h:21
FullPhysics::multifit_test_delta
bool multifit_test_delta(const blitz::Array< double, 1 > &dx, const blitz::Array< double, 1 > &x, double Eps_abs=std::numeric_limits< double >::epsilon(), double Eps_rel=std::numeric_limits< double >::epsilon())
Definition: linear_algebra.cc:239
FullPhysics::sqr
double sqr(double x)
Definition: linear_algebra.h:11
FullPhysics::units::g
const Unit g("g", 1e-3 *kg)
FullPhysics
Contains classes to abstract away details in various Spurr Radiative Transfer software.
Definition: doxygen_python.h:1
FullPhysics::to_c_order
blitz::Array< T, D > to_c_order(const blitz::Array< T, D > &In)
Ensure that a given blitz::Array is contiguous, not reversed, and in C RowMajorArray format...
Definition: linear_algebra.h:44
Copyright © 2017, California Institute of Technology.
ALL RIGHTS RESERVED.
U.S. Government Sponsorship acknowledged.
Generated Fri Aug 24 2018 15:44:11