Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lu.hpp
Go to the documentation of this file.
1 /* -*- mode: c++; coding: utf-8; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; show-trailing-whitespace: t -*- vim:fenc=utf-8:ft=tcl:et:sw=4:ts=4:sts=4
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2005-02-08
7 
8  Copyright (C) 2005,2006 EPFL
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 3.0 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __DenseLU_H
30 #define __DenseLU_H 1
31 
32 #include <boost/mpl/int.hpp>
33 #include <feel/feelcore/feel.hpp>
34 #include <feel/feelcore/traits.hpp>
35 #include <feel/feelalg/glas.hpp>
36 
37 namespace Feel
38 {
39 namespace mpl = boost::mpl;
40 
41 namespace details
42 {
43 template<typename Matrix>
44 inline typename Matrix::value_type
45 det( Matrix const& M, mpl::int_<1> )
46 {
47  return M( 0, 0 );
48 }
49 template<typename Matrix>
50 inline typename Matrix::value_type
51 det( Matrix const& M, mpl::int_<2> )
52 {
53  return M( 0, 0 )*M( 1, 1 )-M( 0, 1 )*M( 1, 0 );
54 }
55 template<typename Matrix>
56 inline typename Matrix::value_type
57 det( Matrix const& M, mpl::int_<3> )
58 {
59  return ( M( 0, 0 )*M( 1, 1 )*M( 2, 2 )-
60  M( 0, 0 )*M( 1, 2 )*M( 2, 1 )-
61  M( 1, 0 )*M( 0, 1 )*M( 2, 2 )+
62  M( 1, 0 )*M( 0, 2 )*M( 2, 1 )+
63  M( 2, 0 )*M( 0, 1 )*M( 1, 2 )-
64  M( 2, 0 )*M( 0, 2 )*M( 1, 1 ) );
65 }
66 
67 template<typename Matrix>
68 inline typename Matrix::value_type
69 det( Matrix const& M, mpl::int_<4> )
70 {
71  return ( M( 0, 0 ) * M( 1, 1 ) * M( 2, 2 ) * M( 3, 3 ) -
72  M( 0, 0 ) * M( 1, 1 ) * M( 2, 3 ) * M( 3, 2 ) -
73  M( 0, 0 ) * M( 2, 1 ) * M( 1, 2 ) * M( 3, 3 ) +
74  M( 0, 0 ) * M( 2, 1 ) * M( 1, 3 ) * M( 3, 2 ) +
75  M( 0, 0 ) * M( 3, 1 ) * M( 1, 2 ) * M( 2, 3 ) -
76  M( 0, 0 ) * M( 3, 1 ) * M( 1, 3 ) * M( 2, 2 ) -
77  M( 1, 0 ) * M( 0, 1 ) * M( 2, 2 ) * M( 3, 3 ) +
78  M( 1, 0 ) * M( 0, 1 ) * M( 2, 3 ) * M( 3, 2 ) +
79  M( 1, 0 ) * M( 2, 1 ) * M( 0, 2 ) * M( 3, 3 ) -
80  M( 1, 0 ) * M( 2, 1 ) * M( 0, 3 ) * M( 3, 2 ) -
81  M( 1, 0 ) * M( 3, 1 ) * M( 0, 2 ) * M( 2, 3 ) +
82  M( 1, 0 ) * M( 3, 1 ) * M( 0, 3 ) * M( 2, 2 ) +
83  M( 2, 0 ) * M( 0, 1 ) * M( 1, 2 ) * M( 3, 3 ) -
84  M( 2, 0 ) * M( 0, 1 ) * M( 1, 3 ) * M( 3, 2 ) -
85  M( 2, 0 ) * M( 1, 1 ) * M( 0, 2 ) * M( 3, 3 ) +
86  M( 2, 0 ) * M( 1, 1 ) * M( 0, 3 ) * M( 3, 2 ) +
87  M( 2, 0 ) * M( 3, 1 ) * M( 0, 2 ) * M( 1, 3 ) -
88  M( 2, 0 ) * M( 3, 1 ) * M( 0, 3 ) * M( 1, 2 ) -
89  M( 3, 0 ) * M( 0, 1 ) * M( 1, 2 ) * M( 2, 3 ) +
90  M( 3, 0 ) * M( 0, 1 ) * M( 1, 3 ) * M( 2, 2 ) +
91  M( 3, 0 ) * M( 1, 1 ) * M( 0, 2 ) * M( 2, 3 ) -
92  M( 3, 0 ) * M( 1, 1 ) * M( 0, 3 ) * M( 2, 2 ) -
93  M( 3, 0 ) * M( 2, 1 ) * M( 0, 2 ) * M( 1, 3 ) +
94  M( 3, 0 ) * M( 2, 1 ) * M( 0, 3 ) * M( 1, 2 ) );
95 
96 }
97 
98 
99 template<typename Matrix>
100 inline void
101 inverse( Matrix const& M, Matrix& Minv, mpl::int_<1> )
102 {
103  Minv( 0, 0 ) = 1.0/M( 0, 0 );
104 }
105 template<typename Matrix>
106 inline void
107 inverse( Matrix const& M, Matrix& Minv, mpl::int_<2> )
108 {
109  typename Matrix::value_type J = M( 0, 0 )*M( 1, 1 )-M( 0, 1 )*M( 1, 0 );
110  Minv( 0, 0 ) = M( 1, 1 )/J;
111  Minv( 1, 0 ) = -M( 1, 0 )/J;
112  Minv( 0, 1 ) = -M( 0, 1 )/J;
113  Minv( 1, 1 ) = M( 0, 0 )/J;
114 }
115 template<typename Matrix>
116 inline void
117 inverse( Matrix const& M, Matrix& Minv, mpl::int_<3> )
118 {
119  typedef typename Matrix::value_type value_type;
120 
121  value_type t4 = M( 0, 0 )*M( 1, 1 );
122  value_type t6 = M( 0, 0 )*M( 1, 2 );
123  value_type t8 = M( 0, 1 )*M( 1, 0 );
124  value_type t10 = M( 0, 2 )*M( 1, 0 );
125 
126  value_type t12 = M( 0, 1 )*M( 2, 0 );
127  value_type t14 = M( 0, 2 )*M( 2, 0 );
128  value_type t17 = 1.0/( M( 2, 2 )*t4-t6*M( 2, 1 )-t8*M( 2, 2 )+t10*M( 2, 1 )+t12*M( 1, 2 )-t14*M( 1, 1 ) );
129 
130  Minv( 0, 0 ) = ( M( 1, 1 )*M( 2, 2 )-M( 1, 2 )*M( 2, 1 ) )*t17;
131  Minv( 0, 1 ) = -( M( 0, 1 )*M( 2, 2 )-M( 0, 2 )*M( 2, 1 ) )*t17;
132  Minv( 0, 2 ) = -( -M( 0, 1 )*M( 1, 2 )+M( 0, 2 )*M( 1, 1 ) )*t17;
133  Minv( 1, 0 ) = -( M( 1, 0 )*M( 2, 2 )-M( 1, 2 )*M( 2, 0 ) )*t17;
134  Minv( 1, 1 ) = ( M( 0, 0 )*M( 2, 2 )-t14 )*t17;
135  Minv( 1, 2 ) = -( t6-t10 )*t17;
136  Minv( 2, 0 ) = -( -M( 1, 0 )*M( 2, 1 )+M( 1, 1 )*M( 2, 0 ) )*t17;
137  Minv( 2, 1 ) = -( M( 0, 0 )*M( 2, 1 )-t12 )*t17;
138  Minv( 2, 2 ) = ( t4-t8 )*t17;
139 }
140 
141 template<typename Matrix>
142 inline void
143 inverse( Matrix const& __restrict__ M, Matrix& __restrict__ Minv, typename Matrix::value_type const& J, mpl::int_<1> )
144 {
145  Minv( 0, 0 ) = typename Matrix::value_type( 1.0 )/J;
146 }
147 template<typename Matrix>
148 inline void
149 inverse( Matrix const& __restrict__ M, Matrix& __restrict__ Minv, typename Matrix::value_type const& J, mpl::int_<2> )
150 {
151  Minv( 0, 0 ) = M( 1, 1 )/J;
152  Minv( 1, 0 ) = -M( 1, 0 )/J;
153  Minv( 0, 1 ) = -M( 0, 1 )/J;
154  Minv( 1, 1 ) = M( 0, 0 )/J;
155 }
156 template<typename Matrix>
157 inline void
158 inverse( Matrix const& __restrict__ M, Matrix& __restrict__ Minv, typename Matrix::value_type const& J, mpl::int_<3> )
159 {
160  typedef typename Matrix::value_type value_type;
161 
162  value_type t4 = M( 0, 0 )*M( 1, 1 );
163  value_type t6 = M( 0, 0 )*M( 1, 2 );
164  value_type t8 = M( 0, 1 )*M( 1, 0 );
165  value_type t10 = M( 0, 2 )*M( 1, 0 );
166 
167  value_type t12 = M( 0, 1 )*M( 2, 0 );
168  value_type t14 = M( 0, 2 )*M( 2, 0 );
169 
170  Minv( 0, 0 ) = ( M( 1, 1 )*M( 2, 2 )-M( 1, 2 )*M( 2, 1 ) )/J;
171  Minv( 0, 1 ) = -( M( 0, 1 )*M( 2, 2 )-M( 0, 2 )*M( 2, 1 ) )/J;
172  Minv( 0, 2 ) = -( -M( 0, 1 )*M( 1, 2 )+M( 0, 2 )*M( 1, 1 ) )/J;
173  Minv( 1, 0 ) = -( M( 1, 0 )*M( 2, 2 )-M( 1, 2 )*M( 2, 0 ) )/J;
174  Minv( 1, 1 ) = ( M( 0, 0 )*M( 2, 2 )-t14 )/J;
175  Minv( 1, 2 ) = -( t6-t10 )/J;
176  Minv( 2, 0 ) = -( -M( 1, 0 )*M( 2, 1 )+M( 1, 1 )*M( 2, 0 ) )/J;
177  Minv( 2, 1 ) = -( M( 0, 0 )*M( 2, 1 )-t12 )/J;
178  Minv( 2, 2 ) = ( t4-t8 )/J;
179 }
180 
181 
182 } // details
183 
184 template<int Dim, typename Matrix>
185 inline typename Matrix::value_type
186 det( Matrix const& M )
187 {
188  return details::det( M, mpl::int_<Dim>() );
189 }
190 
191 template<int Dim, typename Matrix>
192 inline void
193 inverse( Matrix const& M, Matrix& Minv )
194 {
195  details::inverse( M, Minv, mpl::int_<Dim>() );
196 }
197 
198 template<int Dim, typename Matrix>
199 inline void
200 inverse( Matrix const& __restrict__ M, Matrix& __restrict__ Minv, typename Matrix::value_type const& J )
201 {
202  details::inverse( M, Minv, J, mpl::int_<Dim>() );
203 }
204 
205 
218 template <typename MatrixType>
219 class LU
220 {
221 public :
222  //typedef Real value_type;
223  //typedef boost::numeric::ublas::matrix<value_type> matrix_type;
224  typedef typename MatrixType::value_type value_type;
225  typedef MatrixType matrix_type;
226  typedef boost::numeric::ublas::vector<value_type> vector_type;
227  typedef boost::numeric::ublas::vector<uint> vector_uint_type;
228 
234  LU ( const matrix_type &A )
235  :
236  __LU( A.size1(), A.size2() ),
237  m( A.size1() ),
238  n( A.size2() ),
239  pivsign( 1 ),
240  piv( A.size1() )
241 
242  {
243  __LU.assign( A );
244 
245  DVLOG(2) << "LU m = "<< m << "\n";
246  DVLOG(2) << "LU n = "<< n << "\n";
247 
248 
249  // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
250  for ( uint i = 0; i < m; i++ )
251  {
252  piv( i ) = i;
253  }
254 
255  pivsign = 1;
256 
257  //value_type *LUrowi = 0;;
258  vector_type LUcolj( m );
259 
260  // Outer loop.
261 
262  for ( uint j = 0; j < n; j++ )
263  {
264 
265  // Make a copy of the j-th column to localize references.
266  //matrix_column<matrix<value_type> > mc (m, j);
267  for ( uint i = 0; i < m; i++ )
268  {
269  LUcolj( i ) = __LU( i,j );
270  }
271 
272  // Apply previous transformations.
273 
274  for ( uint i = 0; i < m; i++ )
275  {
276  boost::numeric::ublas::matrix_row<matrix_type> LUrowi ( __LU, i );
277  //LUrowi = __LU(i);
278 
279  // Most of the time is spent in the following dot product.
280 
281  uint kmax = std::min( i,j );
282  value_type s = value_type( 0 );
283 
284  for ( uint k = 0; k < kmax; k++ )
285  {
286  s += LUrowi( k )*LUcolj( k );
287  }
288 
289  LUrowi( j ) = LUcolj( i ) -= s;
290  }
291 
292  // Find pivot and exchange if necessary.
293 
294  uint p = j;
295 
296  for ( uint i = j+1; i < m; i++ )
297  {
298  if ( math::abs( LUcolj( i ) ) > math::abs( LUcolj( p ) ) )
299  {
300  p = i;
301  }
302  }
303 
304  if ( p != j )
305  {
306  for ( uint k = 0; k < n; k++ )
307  {
308  value_type t = __LU( p,k );
309  __LU( p,k ) = __LU( j,k );
310  __LU( j,k ) = t;
311  }
312 
313  uint k = piv( p );
314  piv( p ) = piv( j );
315  piv( j ) = k;
316  pivsign = -pivsign;
317  }
318 
319  // Compute multipliers.
320 
321  if ( ( j < m ) && ( __LU( j,j ) != value_type( 0.0 ) ) )
322  {
323  for ( uint i = j+1; i < m; i++ )
324  {
325  __LU( i,j ) /= __LU( j,j );
326  }
327  }
328  }
329  }
330 
331 
338  {
339  for ( uint j = 0; j < n; j++ )
340  {
341  if ( __LU( j,j ) == value_type( 0.0 ) )
342  return 0;
343  }
344 
345  return 1;
346  }
347 
353  matrix_type getL ()
354  {
355  matrix_type L_( m,n );
356 
357  for ( uint i = 0; i < m; i++ )
358  {
359  for ( uint j = 0; j < n; j++ )
360  {
361  if ( i > j )
362  {
363  L_( i,j ) = __LU( i,j );
364  }
365 
366  else if ( i == j )
367  {
368  L_( i,j ) = 1.0;
369  }
370 
371  else
372  {
373  L_( i,j ) = 0.0;
374  }
375  }
376  }
377 
378  return L_;
379  }
380 
385  matrix_type getU ()
386  {
387  matrix_type U_( n,n );
388 
389  for ( uint i = 0; i < n; i++ )
390  {
391  for ( uint j = 0; j < n; j++ )
392  {
393  if ( i <= j )
394  {
395  U_( i,j ) = __LU( i,j );
396  }
397 
398  else
399  {
400  U_( i,j ) = 0.0;
401  }
402  }
403  }
404 
405  return U_;
406  }
407 
411  vector_uint_type getPivot ()
412  {
413  return piv;
414  }
415 
416 
420  value_type det ()
421  {
422  if ( m != n )
423  {
424  return value_type( 0 );
425  }
426 
427  value_type d = value_type( pivsign );
428 
429  //DVLOG(2) << "LU::det() d= " << pivsign << "\n";
430  for ( uint j = 0; j < n; j++ )
431  {
432  d *= __LU( j,j );
433  }
434 
435  return d;
436  }
437 
438  void inverse( matrix_type& __inv )
439  {
440  __inv = solve( ublas::identity_matrix<double>( m, m ) );
441 #if 0
442  vector_type __t( m );
443  __t = ZeroVector( __t.size() );
444 
445  for ( size_type i = 0; i < piv.size(); ++i )
446  {
447  __t[i] = value_type( 1 );
448  ublas::row( __inv, i ) = solve( __t );
449  __t[i] = value_type( 0 );
450  }
451 
452 #endif
453  }
459  matrix_type solve ( const matrix_type &B )
460  {
461 
462  /* Dimensions: A is mxn, X is nxk, B is mxk */
463 
464  if ( B.size1() != m )
465  {
466  return matrix_type( 0,0 );
467  }
468 
469  if ( !isNonsingular() )
470  {
471  return matrix_type( 0,0 );
472  }
473 
474  // Copy right hand side with pivoting
475  uint nx = B.size2();
476 
477 
478  matrix_type X ( permute_copy( B, piv, 0, nx-1 ) );
479 
480  // Solve L*Y = B(piv,:)
481  for ( uint k = 0; k < n; k++ )
482  {
483  for ( uint i = k+1; i < n; i++ )
484  {
485  for ( uint j = 0; j < nx; j++ )
486  {
487  X( i,j ) -= X( k,j )*__LU( i,k );
488  }
489  }
490  }
491 
492  // Solve U*X = Y;
493  for ( int k = ( int )n-1; k >= 0; k-- )
494  {
495  for ( uint j = 0; j < nx; j++ )
496  {
497  X( k,j ) /= __LU( k,k );
498  }
499 
500  for ( int i = 0; i < k; i++ )
501  {
502  for ( uint j = 0; j < nx; j++ )
503  {
504  X( i,j ) -= X( k,j )*__LU( i,k );
505  }
506  }
507  }
508 
509  return X;
510  }
511 
512 
522  vector_type solve ( const vector_type &b )
523  {
524 
525  /* Dimensions: A is mxn, X is nxk, B is mxk */
526 
527  if ( b.size() != m )
528  {
529  return vector_type();
530  }
531 
532  if ( !isNonsingular() )
533  {
534  return vector_type();
535  }
536 
537 
538  vector_type x = permute_copy( b, piv );
539 
540  // Solve L*Y = B(piv)
541  for ( uint k = 0; k < n; k++ )
542  {
543  for ( uint i = k+1; i < n; i++ )
544  {
545  x( i ) -= x( k )*__LU( i,k );
546  }
547  }
548 
549  // Solve U*X = Y;
550  for ( int k = ( int )n-1; k >= 0; k-- )
551  {
552  x( k ) /= __LU( k,k );
553 
554  for ( int i = 0; i < k; i++ )
555  {
556  x( i ) -= x( k )*__LU( i,k );
557  }
558  }
559 
560 
561  return x;
562  }
563 
564 private:
565  /* Array for internal storage of decomposition. */
566  matrix_type __LU;
567  size_type m, n;
568  int pivsign;
569  vector_uint_type piv;
570 
571 
572  matrix_type
573  permute_copy( const matrix_type &A, const vector_uint_type &piv, uint j0, uint j1 )
574  {
575  uint piv_length = piv.size();
576 
577  matrix_type X( piv_length, j1-j0+1 );
578 
579 
580  for ( uint i = 0; i < piv_length; i++ )
581  for ( uint j = j0; j <= j1; j++ )
582  X( i,j-j0 ) = A( piv( i ),j );
583 
584  return X;
585  }
586 
587  vector_type
588  permute_copy( const vector_type &A, const vector_uint_type &piv )
589  {
590  uint piv_length = piv.size();
591 
592  if ( piv_length != A.size() )
593  return vector_type();
594 
595  vector_type x( piv_length );
596 
597 
598  for ( uint i = 0; i < piv_length; i++ )
599  x( i ) = A( piv( i ) );
600 
601  return x;
602  }
603 
604 }; /* class LU */
605 
606 }
607 #endif /* __DenseLU_H */
vector_uint_type getPivot()
Definition: lu.hpp:411
Definition: solverlinear.hpp:33
LU(const matrix_type &A)
Definition: lu.hpp:234
Definition: lu.hpp:219
vector_type solve(const vector_type &b)
Definition: lu.hpp:522
value_type det()
Definition: lu.hpp:420
matrix_type getU()
Definition: lu.hpp:385
matrix_type getL()
Definition: lu.hpp:353
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
uint isNonsingular()
Definition: lu.hpp:337
matrix_type solve(const matrix_type &B)
Definition: lu.hpp:459

Generated on Sun Dec 22 2013 13:11:06 for Feel++ by doxygen 1.8.5