Rev Author Line No. Line
178 kaklik 1 #ifndef LINALG_H
2 #define LINALG_H
3  
4 #include <boost/numeric/ublas/banded.hpp>
5 #include <boost/numeric/ublas/matrix.hpp>
6 #include <boost/numeric/ublas/symmetric.hpp>
7 #include <boost/numeric/ublas/vector.hpp>
8 #include <boost/numeric/ublas/operation.hpp>
9 #include <boost/numeric/ublas/lu.hpp>
10 #include <boost/numeric/ublas/vector_proxy.hpp>
11 #include <cassert>
12 #include "mimasexception.h"
13  
14 namespace mimas {
15  
16 /** @defgroup linearAlgebra Linear Algebra
17 mimas offers support for linear algebra using
18 <A HREF="http://www.boost.org/libs/numeric/ublas/doc/index.htm">boost's
19 uBLAS</A>-library and the
20 <A HREF="http://www.netlib.org/lapack/lug/lapack_lug.html">LAPACK</A>-library.
21 @{ */
22 /** @defgroup ublas Mimas linear algebra operations.
23 This group contains Mimas operations, for vectors and matrices. All
24 functionality provided by \c mimas::vector formerly, should be provided
25 here.
26 @author Stuart Meikle (stu@stumeikle.org)
27 @date Sat Mar 1 2003 BST
28 @{ */
29 /** Compute unit vector.
30 Compute unit vector \f$\frac{\vec{x}}{|\vec{x}|}\f$.
31 @param x input vector (non-zero).
32 @return Vector scaled to length 1. */
33 template< typename T >
34 boost::numeric::ublas::vector< T > unit
35 ( const boost::numeric::ublas::vector< T > &x )
36 {
37 assert( boost::numeric::ublas::norm_2( x ) > 0 );
38 return x / boost::numeric::ublas::norm_2( x );
39 }
40  
41 /** Scalar cross product.
42 The scalar cross product for 2-dimensional vectors is
43 \f$a\times b:=a_1\,b_2-a_2\,b_1\f$. It is rotational invariant.
44  
45 For scalar product see \c boost::inner_prod.
46 @param a first vector.
47 @param b second vector. */
48 template< typename T >
49 T scalar_cross_product( const boost::numeric::ublas::vector< T > &a,
50 const boost::numeric::ublas::vector< T > &b )
51 {
52 assert( a.size() == 2 );
53 assert( b.size() == 2 );
54 return a(0) * b(1) - a(1) * b(0);
55 }
56  
57 /** Cross product for 3-dimensional vectors.
58 Compute
59 \f$a\times b:=(a_2\,b_3-a_3\,b_2,a_3\,b_1-a_1\,b_3,
60 a_1\,b_2-a_2\,b_1)^\top\f$
61 @param a first vector.
62 @param b second vector.
63 @return cross product of \c a and \c b. */
64 template< typename T >
65 boost::numeric::ublas::vector< T > crossProduct
66 ( boost::numeric::ublas::vector< T > &a,
67 boost::numeric::ublas::vector< T > &b );
68  
69 /** Rodrigues' rotation formula
70  
71 Computes the rotation of a point around an axis of given vector,
72 and passing through the origin.
73  
74 @param u the axis of the rotation (must be a unit 3-D vector)
75 @param v the vector to rotate
76 @param angle the angle of the rotation
77 @author Julien Faucher (faucherj@gmail.com)
78 */
79 template< typename T >
80 boost::numeric::ublas::vector< T > rodrigues
81 ( boost::numeric::ublas::vector< T > const &u,
82 boost::numeric::ublas::vector< T > const &v,
83 double angle);
84  
85 /** Determinant of a matrix using the LU factorization.
86  
87 The decomposition used is the one provided by uBLAS.
88  
89 @param M the matrix to compute the determinant of.
90 @author Julien Faucher (faucherj@gmail.com)
91 */
92 template< typename T >
93 double determinant (boost::numeric::ublas::matrix< T > const &M);
94  
95 ///@}
96  
97 /** @defgroup lapack Lapack functions
98 Wrappers or bindings for Lapack functions are provided. The matrix- and
99 vector-classes of boost are used as datatypes.
100  
101 Note, that lapack uses matrices in <B>column-major</B> orientation. You
102 are forced to convert to
103 \code
104 boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major >
105 \endcode
106 The default orientation of boost matrices is row-major.
107  
108 @author Bala Amavasai (bala@amavasai.org)
109 @author Jan Wedekind (jan@wedesoft.de)
110 @todo Extend the solver-classes to complex values.
111 @date Fri May 31 12:07:52 UTC 2002
112 @{ */
113 /** Function-object for gglse.
114 @see gglse */
115 template< typename T >
116 struct gglse_t
117 {
118 ///
119 typedef boost::numeric::ublas::vector< T > Vector;
120 ///
121 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
122 /// Invoke lapack method.
123 Vector operator()( const Matrix &_A, const Matrix &_B,
124 const Vector &_c, const Vector &_d ) const;
125 };
126  
127 /** Solve the linear equality-constrained least squares (LSE) problem.
128 \f$A\f$, \f$\vec{c}\f$, \f$B\f$ and \f$\vec{d}\f$ are given and
129 \f$\vec{x}\f$ is
130 sought:
131  
132 minimize \f$\vec{x}:=\displaystyle\mathop{argmin}_{\vec{x}}|\vec{c}-A\,\vec{x}|\f$
133 subject to \f$B\,\vec{x}=\vec{d}\f$.
134  
135 See manpage <A HREF="man:/dgglse">dgglse(1)</A> for more documentation.
136 @see gglse_t */
137 template< typename T >
138 boost::numeric::ublas::vector< T > gglse
139 ( const boost::numeric::ublas::matrix
140 < T, boost::numeric::ublas::column_major > &A,
141 const boost::numeric::ublas::matrix
142 < T, boost::numeric::ublas::column_major > &B,
143 const boost::numeric::ublas::vector< T > &c,
144 const boost::numeric::ublas::vector< T > &d )
145 { return gglse_t< T >()( A, B, c, d ); }
146  
147 /** Function-object for inv.
148 @see inv */
149 template< typename T >
150 struct inv_t
151 {
152 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
153 /// Invoke lapack method.
154 Matrix operator()( const Matrix &A ) const;
155 };
156  
157 /** Compute matrix inverse by using LU factorization.
158 The matrix \f$A\in\mathbf{K}^{n\times n}\f$ is given and the inverse
159 \f$A^{-1}\f$ has to be computed.
160  
161 See manpages <A HREF="man:/dgetrf">dgetrf(1)</A> and
162 <A HREF="man:/dgetri">dgetri(1)</A> for more documentation.
163 @see inv_t */
164 template< typename T >
165 boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > inv
166 ( const boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > &A )
167 { return inv_t< T >()( A ); }
168  
169 /** Function-object for gels.
170 @see gels */
171 template< typename T >
172 struct gels_t
173 {
174 ///
175 typedef boost::numeric::ublas::vector< T > Vector;
176 ///
177 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
178 /// Invoke lapack method.
179 Vector operator()( const Matrix &_A, const Vector &_b ) const;
180 };
181  
182 /** solve overdetermined linear systems using QR factorization.
183 Solve the minimum least square problem. \f$A\f$ and \f$\vec{b}\f$ are given
184 and \f$\vec{x}\f$ is the desired solution for
185 \f[\vec{x}:=\displaystyle\mathop{argmin}_{\vec{x}}|\vec{b}-A\,\vec{x}|\f]
186  
187 This is much faster than computing
188 \f$\vec{x}:=(A^\top\,A)^{-1}\,A^\top\,\vec{b}\f$ with the matrix-classes.
189  
190 See manpage <A HREF="man:/dgels">dgels(1)</A> for more documentation. The
191 parametrisation for underdetermined systems wasn't working as expected and
192 was not included in the interface therefore.
193 @see gels_t */
194 template< typename T >
195 boost::numeric::ublas::vector< T > gels
196 ( const boost::numeric::ublas::matrix
197 < T, boost::numeric::ublas::column_major > &A,
198 const boost::numeric::ublas::vector< T > &b )
199 { return gels_t< T >()( A, b ); }
200  
201 /** Function-object for gesvd.
202 @see gesvd */
203 template< typename T >
204 struct gesvd_t
205 {
206 ///
207 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
208 ///
209 typedef boost::numeric::ublas::banded_matrix< T > BandedMatrix;
210 /// Invoke lapack-method.
211 BandedMatrix operator()( const Matrix &_A, Matrix *_U, Matrix *_Vt )
212 throw(mimasexception);
213 };
214  
215 /** Compute the singular value decomposition (SVD) of a matrix A.
216 \f$A\f$ is given. \f$U\f$, \f$\Sigma\f$ and \f$V^\top\f$ have to be
217 determined, such that:
218 \f$A=U\,\Sigma\,V^\top\f$ and
219 \f$\Sigma=\mathop{diag}(\sigma_1,\sigma_2,\ldots,\sigma_n)\f$ diagonal
220 matrix with the singular values \f$\sigma_1,\sigma_2,\ldots,\sigma_n\f$
221 (\f$\sigma_i\ge 0\f$) in descending order. The singular values are sorted
222 in descending order.
223  
224 Note that the singular value decomposition doesn't necessarily converge.
225 In the later case an exception is thrown.
226  
227 See \c mimas::gesdd, which is a faster algorithm.
228  
229 See manpage <A HREF="man:/dgesvd">dgesvd(1)</A> for more documentation.
230  
231 @param A Matrix to be decomposed.
232 @param U Returns matrix with left hand singular vectors \f$U\f$ (may be
233 \c NULL, if it is of no interest).
234 @param Vt Returns transposed matrix with right hand singular vectors
235 \f$V^\top\f$ (may be \c NULL, if it is of no interest).
236 @return Diagonal matrix \f$\Sigma\f$.
237  
238 @todo Does not converge on AMD64-processor.
239 @see gesdd
240 @see gesvd_t */
241 template< typename T >
242 boost::numeric::ublas::banded_matrix< T > gesvd
243 ( const boost::numeric::ublas::matrix
244 < T, boost::numeric::ublas::column_major > &A,
245 boost::numeric::ublas::matrix
246 < T, boost::numeric::ublas::column_major > *U,
247 boost::numeric::ublas::matrix
248 < T, boost::numeric::ublas::column_major > *Vt )
249 { return gesvd_t< T >()( A, U, Vt ); }
250  
251 /** Function-wrapper for gesdd.
252 @see gesdd */
253 template< typename T >
254 struct gesdd_t
255 {
256 ///
257 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
258 ///
259 typedef boost::numeric::ublas::banded_matrix< T > BandedMatrix;
260 /// Invoke lapack-method.
261 BandedMatrix operator()( const Matrix &_A, Matrix *_U, Matrix *_Vt )
262 throw (mimasexception);
263 };
264  
265 /** Compute the singular value decomposition (SVD) of a matrix A (fast).
266 Same as \c mimas::gesvd but faster algorithm (divide-and-conquer
267 approach).
268  
269 Note that the singular value decomposition doesn't necessarily converge.
270 In the later case an exception is thrown.
271  
272 See manpage <A HREF="man:/dgesdd">dgesdd(1)</A> for more documentation.
273  
274 Note, that, in contrast to the interface of \c mimas::gesvd, \c U
275 and \c Vt must both be defined or both be \c NULL!
276 @param A Matrix to be decomposed.
277 @param U Returns matrix with left hand singular vectors \f$U\f$ (may be
278 \c NULL, if it is of no interest).
279 @param Vt Returns transposed matrix with right hand singular vectors
280 `\f$V^\top\f$ (may be \c NULL, if it is of no interest).
281 @return Diagonal matrix \f$\Sigma\f$.
282 @see gesvd
283 @see gesdd_t */
284 template< typename T >
285 boost::numeric::ublas::banded_matrix< T > gesdd
286 ( const boost::numeric::ublas::matrix
287 < T, boost::numeric::ublas::column_major > &A,
288 boost::numeric::ublas::matrix
289 < T, boost::numeric::ublas::column_major > *U,
290 boost::numeric::ublas::matrix
291 < T, boost::numeric::ublas::column_major > *Vt )
292 { return gesdd_t< T >()( A, U, Vt ); }
293  
294 /** Function-wrapper for syev.
295 @see syev */
296 template< typename T >
297 struct syev_t
298 {
299 ///
300 typedef boost::numeric::ublas::matrix< T, boost::numeric::ublas::column_major > Matrix;
301 ///
302 typedef boost::numeric::ublas::symmetric_matrix< T, boost::numeric::ublas::upper > SymmetricMatrix;
303 ///
304 typedef boost::numeric::ublas::diagonal_matrix< T > EigenValues;
305 /// Invoke lapack-method.
306 EigenValues operator()( const SymmetricMatrix &_A, Matrix *_E ) const
307 throw (mimasexception);
308 };
309  
310 /** Compute all eigenvalues and, optionally, eigenvectors of a real symmetric matrix A.
311 The eigenvalues and eigenvectors of a symmetric matrix
312 \f$A\in\mathbf{K}^{n\times n}\f$ can be computed.
313  
314 The eigenvectors are orthogonal by nature and normalised by definition.
315 If the matrix \f$E\f$ contains the eigenvectors of \f$A\f$, then the
316 eigentransform can be written as (\f$E\f$ is orthonormal \f$\Rightarrow\f$
317 \f$E^{-1}=E^\top\f$):
318 \f[A=E\,\Lambda\,E^\top\f]
319  
320 Where \f$\Lambda=\mathop{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n)\f$ is
321 a diagonal matrix with the eigenvalues
322 \f$\lambda_1,\lambda_2,\ldots,\lambda_n\f$ in ascending order.
323  
324 See manpage <A HREF="man:/dsyev">dsyev(1)</A> for more documentation.
325  
326 Compute eigenvalues and, optionally, the eigenvectors of A.
327 @param A Symmetric matrix (with upper storage) to compute
328 eigentransform for.
329 @param E Pointer to matrix, where the eigenvectors have to be written
330 to (may be \c NULL, if it is of no interest).
331 @return Diagonal matrix with eigenvalues in ascending order.
332 @see syev_t */
333 template< typename T >
334 boost::numeric::ublas::diagonal_matrix< T > syev
335 ( const boost::numeric::ublas::symmetric_matrix
336 < T, boost::numeric::ublas::upper > &A,
337 boost::numeric::ublas::matrix
338 < T, boost::numeric::ublas::column_major > *E )
339 { return syev_t< T >()( A, E ); }
340  
341  
342 /** Function-object for pptrf.
343 @see pptrf */
344 template< typename T >
345 struct pptrf_t
346 {
347 ///
348 typedef boost::numeric::ublas::triangular_matrix
349 < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > TriangularMatrix;
350 ///
351 typedef boost::numeric::ublas::symmetric_matrix
352 < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > SymmetricMatrix;
353  
354 /// Invoke lapack-method.
355 TriangularMatrix operator()( const SymmetricMatrix &_A )
356 throw(mimasexception);
357 };
358  
359 /** Compute the Cholesky decomposition of a symmetric matrix A.
360  
361 Given a symmetric positive-definite matrix\f$A\in\mathbf{K}^{n\times n}\f$,
362 the Cholesky decomposition returns the upper triangular matrix
363 \f$U\in\mathbf{K}^{n\times n}\f$ such that:
364 \f[A=U^\top\,U\f]
365  
366 @param A Symmetric matrix (with upper and column-based storage) to compute
367 Cholesky decomposition for.
368 @return Upper triangular matrix (with column-based storage).
369 @author Julien Faucher (faucherj@gmail.com)
370 @see pptrf_t */
371 template< typename T >
372 boost::numeric::ublas::triangular_matrix
373 < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > pptrf
374 ( const boost::numeric::ublas::symmetric_matrix
375 < T, boost::numeric::ublas::upper, boost::numeric::ublas::column_major > &A)
376 { return pptrf_t< T >()( A ); }
377  
378 ///@}
379  
380 ///@}
381  
382 };
383  
384 #include "linalg.tcc"
385  
386 #endif
387