Symaptic: os/ossrv/ossrv_pub/boost_apis/boost/numeric/ublas/lu.hpp@bde4ae8d615e

     1 //

     2 //  Copyright (c) 2000-2002

     3 //  Joerg Walter, Mathias Koch

     4 //

     5 //  Permission to use, copy, modify, distribute and sell this software

     6 //  and its documentation for any purpose is hereby granted without fee,

     7 //  provided that the above copyright notice appear in all copies and

     8 //  that both that copyright notice and this permission notice appear

     9 //  in supporting documentation.  The authors make no representations

    10 //  about the suitability of this software for any purpose.

    11 //  It is provided "as is" without express or implied warranty.

    12 //

    13 //  The authors gratefully acknowledge the support of

    14 //  GeNeSys mbH & Co. KG in producing this work.

    15 //

    17 #ifndef _BOOST_UBLAS_LU_

    18 #define _BOOST_UBLAS_LU_

    20 #include <boost/numeric/ublas/operation.hpp>

    21 #include <boost/numeric/ublas/vector_proxy.hpp>

    22 #include <boost/numeric/ublas/matrix_proxy.hpp>

    23 #include <boost/numeric/ublas/vector.hpp>

    24 #include <boost/numeric/ublas/triangular.hpp>

    26 // LU factorizations in the spirit of LAPACK and Golub & van Loan

    28 namespace boost { namespace numeric { namespace ublas {

    30     template<class T = std::size_t, class A = unbounded_array<T> >

    31     class permutation_matrix:

    32         public vector<T, A> {

    33     public:

    34         typedef vector<T, A> vector_type;

    35         typedef typename vector_type::size_type size_type;

    37         // Construction and destruction

    38         BOOST_UBLAS_INLINE

    39         permutation_matrix (size_type size):

    40             vector<T, A> (size) {

    41             for (size_type i = 0; i < size; ++ i)

    42                 (*this) (i) = i;

    43         }

    44         BOOST_UBLAS_INLINE

    45         ~permutation_matrix () {}

    47         // Assignment

    48         BOOST_UBLAS_INLINE

    49         permutation_matrix &operator = (const permutation_matrix &m) {

    50             vector_type::operator = (m);

    51             return *this;

    52         }

    53     };

    55     template<class PM, class MV>

    56     BOOST_UBLAS_INLINE

    57     void swap_rows (const PM &pm, MV &mv, vector_tag) {

    58         typedef typename PM::size_type size_type;

    59         typedef typename MV::value_type value_type;

    61         size_type size = pm.size ();

    62         for (size_type i = 0; i < size; ++ i) {

    63             if (i != pm (i))

    64                 std::swap (mv (i), mv (pm (i)));

    65         }

    66     }

    67     template<class PM, class MV>

    68     BOOST_UBLAS_INLINE

    69     void swap_rows (const PM &pm, MV &mv, matrix_tag) {

    70         typedef typename PM::size_type size_type;

    71         typedef typename MV::value_type value_type;

    73         size_type size = pm.size ();

    74         for (size_type i = 0; i < size; ++ i) {

    75             if (i != pm (i))

    76                 row (mv, i).swap (row (mv, pm (i)));

    77         }

    78     }

    79     // Dispatcher

    80     template<class PM, class MV>

    81     BOOST_UBLAS_INLINE

    82     void swap_rows (const PM &pm, MV &mv) {

    83         swap_rows (pm, mv, typename MV::type_category ());

    84     }

    86     // LU factorization without pivoting

    87     template<class M>

    88     typename M::size_type lu_factorize (M &m) {

    89         typedef M matrix_type;

    90         typedef typename M::size_type size_type;

    91         typedef typename M::value_type value_type;

    93 #if BOOST_UBLAS_TYPE_CHECK

    94         matrix_type cm (m);

    95 #endif

    96         int singular = 0;

    97         size_type size1 = m.size1 ();

    98         size_type size2 = m.size2 ();

    99         size_type size = (std::min) (size1, size2);

   100         for (size_type i = 0; i < size; ++ i) {

   101             matrix_column<M> mci (column (m, i));

   102             matrix_row<M> mri (row (m, i));

   103             if (m (i, i) != value_type/*zero*/()) {

   104                 project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);

   105             } else if (singular == 0) {

   106                 singular = i + 1;

   107             }

   108             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (

   109                 outer_prod (project (mci, range (i + 1, size1)),

   110                             project (mri, range (i + 1, size2))));

   111         }

   112 #if BOOST_UBLAS_TYPE_CHECK

   113         BOOST_UBLAS_CHECK (singular != 0 ||

   114                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),

   115                                                                 triangular_adaptor<matrix_type, upper> (m)),

   116                                                           cm), internal_logic ());

   117 #endif

   118         return singular;

   119     }

   121     // LU factorization with partial pivoting

   122     template<class M, class PM>

   123     typename M::size_type lu_factorize (M &m, PM &pm) {

   124         typedef M matrix_type;

   125         typedef typename M::size_type size_type;

   126         typedef typename M::value_type value_type;

   128 #if BOOST_UBLAS_TYPE_CHECK

   129         matrix_type cm (m);

   130 #endif

   131         int singular = 0;

   132         size_type size1 = m.size1 ();

   133         size_type size2 = m.size2 ();

   134         size_type size = (std::min) (size1, size2);

   135         for (size_type i = 0; i < size; ++ i) {

   136             matrix_column<M> mci (column (m, i));

   137             matrix_row<M> mri (row (m, i));

   138             size_type i_norm_inf = i + index_norm_inf (project (mci, range (i, size1)));

   139             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());

   140             if (m (i_norm_inf, i) != value_type/*zero*/()) {

   141                 if (i_norm_inf != i) {

   142                     pm (i) = i_norm_inf;

   143                     row (m, i_norm_inf).swap (mri);

   144                 } else {

   145                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());

   146                 }

   147                 project (mci, range (i + 1, size1)) *= value_type (1) / m (i, i);

   148             } else if (singular == 0) {

   149                 singular = i + 1;

   150             }

   151             project (m, range (i + 1, size1), range (i + 1, size2)).minus_assign (

   152                 outer_prod (project (mci, range (i + 1, size1)),

   153                             project (mri, range (i + 1, size2))));

   154         }

   155 #if BOOST_UBLAS_TYPE_CHECK

   156         swap_rows (pm, cm);

   157         BOOST_UBLAS_CHECK (singular != 0 ||

   158                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),

   159                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());

   160 #endif

   161         return singular;

   162     }

   164     template<class M, class PM>

   165     typename M::size_type axpy_lu_factorize (M &m, PM &pm) {

   166         typedef M matrix_type;

   167         typedef typename M::size_type size_type;

   168         typedef typename M::value_type value_type;

   169         typedef vector<value_type> vector_type;

   171 #if BOOST_UBLAS_TYPE_CHECK

   172         matrix_type cm (m);

   173 #endif

   174         int singular = 0;

   175         size_type size1 = m.size1 ();

   176         size_type size2 = m.size2 ();

   177         size_type size = (std::min) (size1, size2);

   178 #ifndef BOOST_UBLAS_LU_WITH_INPLACE_SOLVE

   179         matrix_type mr (m);

   180         mr.assign (zero_matrix<value_type> (size1, size2));

   181         vector_type v (size1);

   182         for (size_type i = 0; i < size; ++ i) {

   183             matrix_range<matrix_type> lrr (project (mr, range (0, i), range (0, i)));

   184             vector_range<matrix_column<matrix_type> > urr (project (column (mr, i), range (0, i)));

   185             urr.assign (solve (lrr, project (column (m, i), range (0, i)), unit_lower_tag ()));

   186             project (v, range (i, size1)).assign (

   187                 project (column (m, i), range (i, size1)) -

   188                 axpy_prod<vector_type> (project (mr, range (i, size1), range (0, i)), urr));

   189             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));

   190             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());

   191             if (v (i_norm_inf) != value_type/*zero*/()) {

   192                 if (i_norm_inf != i) {

   193                     pm (i) = i_norm_inf;

   194                     std::swap (v (i_norm_inf), v (i));

   195                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));

   196                 } else {

   197                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());

   198                 }

   199                 project (column (mr, i), range (i + 1, size1)).assign (

   200                     project (v, range (i + 1, size1)) / v (i));

   201                 if (i_norm_inf != i) {

   202                     project (row (mr, i_norm_inf), range (0, i)).swap (project (row (mr, i), range (0, i)));

   203                 }

   204             } else if (singular == 0) {

   205                 singular = i + 1;

   206             }

   207             mr (i, i) = v (i);

   208         }

   209         m.assign (mr);

   210 #else

   211         matrix_type lr (m);

   212         matrix_type ur (m);

   213         lr.assign (identity_matrix<value_type> (size1, size2));

   214         ur.assign (zero_matrix<value_type> (size1, size2));

   215         vector_type v (size1);

   216         for (size_type i = 0; i < size; ++ i) {

   217             matrix_range<matrix_type> lrr (project (lr, range (0, i), range (0, i)));

   218             vector_range<matrix_column<matrix_type> > urr (project (column (ur, i), range (0, i)));

   219             urr.assign (project (column (m, i), range (0, i)));

   220             inplace_solve (lrr, urr, unit_lower_tag ());

   221             project (v, range (i, size1)).assign (

   222                 project (column (m, i), range (i, size1)) -

   223                 axpy_prod<vector_type> (project (lr, range (i, size1), range (0, i)), urr));

   224             size_type i_norm_inf = i + index_norm_inf (project (v, range (i, size1)));

   225             BOOST_UBLAS_CHECK (i_norm_inf < size1, external_logic ());

   226             if (v (i_norm_inf) != value_type/*zero*/()) {

   227                 if (i_norm_inf != i) {

   228                     pm (i) = i_norm_inf;

   229                     std::swap (v (i_norm_inf), v (i));

   230                     project (row (m, i_norm_inf), range (i + 1, size2)).swap (project (row (m, i), range (i + 1, size2)));

   231                 } else {

   232                     BOOST_UBLAS_CHECK (pm (i) == i_norm_inf, external_logic ());

   233                 }

   234                 project (column (lr, i), range (i + 1, size1)).assign (

   235                     project (v, range (i + 1, size1)) / v (i));

   236                 if (i_norm_inf != i) {

   237                     project (row (lr, i_norm_inf), range (0, i)).swap (project (row (lr, i), range (0, i)));

   238                 }

   239             } else if (singular == 0) {

   240                 singular = i + 1;

   241             }

   242             ur (i, i) = v (i);

   243         }

   244         m.assign (triangular_adaptor<matrix_type, strict_lower> (lr) +

   245                   triangular_adaptor<matrix_type, upper> (ur));

   246 #endif

   247 #if BOOST_UBLAS_TYPE_CHECK

   248         swap_rows (pm, cm);

   249         BOOST_UBLAS_CHECK (singular != 0 ||

   250                            detail::expression_type_check (prod (triangular_adaptor<matrix_type, unit_lower> (m),

   251                                                                 triangular_adaptor<matrix_type, upper> (m)), cm), internal_logic ());

   252 #endif

   253         return singular;

   254     }

   256     // LU substitution

   257     template<class M, class E>

   258     void lu_substitute (const M &m, vector_expression<E> &e) {

   259         typedef const M const_matrix_type;

   260         typedef vector<typename E::value_type> vector_type;

   262 #if BOOST_UBLAS_TYPE_CHECK

   263         vector_type cv1 (e);

   264 #endif

   265         inplace_solve (m, e, unit_lower_tag ());

   266 #if BOOST_UBLAS_TYPE_CHECK

   267         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cv1), internal_logic ());

   268         vector_type cv2 (e);

   269 #endif

   270         inplace_solve (m, e, upper_tag ());

   271 #if BOOST_UBLAS_TYPE_CHECK

   272         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cv2), internal_logic ());

   273 #endif

   274     }

   275     template<class M, class E>

   276     void lu_substitute (const M &m, matrix_expression<E> &e) {

   277         typedef const M const_matrix_type;

   278         typedef matrix<typename E::value_type> matrix_type;

   280 #if BOOST_UBLAS_TYPE_CHECK

   281         matrix_type cm1 (e);

   282 #endif

   283         inplace_solve (m, e, unit_lower_tag ());

   284 #if BOOST_UBLAS_TYPE_CHECK

   285         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, unit_lower> (m), e), cm1), internal_logic ());

   286         matrix_type cm2 (e);

   287 #endif

   288         inplace_solve (m, e, upper_tag ());

   289 #if BOOST_UBLAS_TYPE_CHECK

   290         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (triangular_adaptor<const_matrix_type, upper> (m), e), cm2), internal_logic ());

   291 #endif

   292     }

   293     template<class M, class PMT, class PMA, class MV>

   294     void lu_substitute (const M &m, const permutation_matrix<PMT, PMA> &pm, MV &mv) {

   295         swap_rows (pm, mv);

   296         lu_substitute (m, mv);

   297     }

   298     template<class E, class M>

   299     void lu_substitute (vector_expression<E> &e, const M &m) {

   300         typedef const M const_matrix_type;

   301         typedef vector<typename E::value_type> vector_type;

   303 #if BOOST_UBLAS_TYPE_CHECK

   304         vector_type cv1 (e);

   305 #endif

   306         inplace_solve (e, m, upper_tag ());

   307 #if BOOST_UBLAS_TYPE_CHECK

   308         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cv1), internal_logic ());

   309         vector_type cv2 (e);

   310 #endif

   311         inplace_solve (e, m, unit_lower_tag ());

   312 #if BOOST_UBLAS_TYPE_CHECK

   313         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cv2), internal_logic ());

   314 #endif

   315     }

   316     template<class E, class M>

   317     void lu_substitute (matrix_expression<E> &e, const M &m) {

   318         typedef const M const_matrix_type;

   319         typedef matrix<typename E::value_type> matrix_type;

   321 #if BOOST_UBLAS_TYPE_CHECK

   322         matrix_type cm1 (e);

   323 #endif

   324         inplace_solve (e, m, upper_tag ());

   325 #if BOOST_UBLAS_TYPE_CHECK

   326         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, upper> (m)), cm1), internal_logic ());

   327         matrix_type cm2 (e);

   328 #endif

   329         inplace_solve (e, m, unit_lower_tag ());

   330 #if BOOST_UBLAS_TYPE_CHECK

   331         BOOST_UBLAS_CHECK (detail::expression_type_check (prod (e, triangular_adaptor<const_matrix_type, unit_lower> (m)), cm2), internal_logic ());

   332 #endif

   333     }

   334     template<class MV, class M, class PMT, class PMA>

   335     void lu_substitute (MV &mv, const M &m, const permutation_matrix<PMT, PMA> &pm) {

   336         swap_rows (pm, mv);

   337         lu_substitute (mv, m);

   338     }

   340 }}}

   342 #endif

author	sl@SLION-WIN7.fritz.box
	Fri, 15 Jun 2012 03:10:57 +0200
changeset 0	bde4ae8d615e
permissions	-rw-r--r--