ceph/src/boost/libs/numeric/ublas/test/test_lu.cpp

   1 // Copyright 2008 Gunter Winkler <guwi17@gmx.de>
   2 // Distributed under the Boost Software License, Version 1.0. (See
   3 // accompanying file LICENSE_1_0.txt or copy at
   4 // http://www.boost.org/LICENSE_1_0.txt)
   5
   6 // switch automatic singular check off
   7 #define BOOST_UBLAS_TYPE_CHECK 0
   8
   9 #include <boost/numeric/ublas/io.hpp>
  10 #include <boost/numeric/ublas/lu.hpp>
  11 #include <boost/cstdlib.hpp>
  12
  13 #include "common/testhelper.hpp"
  14
  15 #include <iostream>
  16 #include <sstream>
  17
  18 using namespace boost::numeric::ublas;
  19 using std::string;
  20
  21 static const string matrix_IN = "[3,3]((1,2,2),(2,3,3),(3,4,6))\0";
  22 static const string matrix_LU = "[3,3]((3,4,6),(3.33333343e-01,6.66666627e-01,0),(6.66666687e-01,4.99999911e-01,-1))\0";
  23 static const string matrix_INV= "[3,3]((-3,2,-7.94728621e-08),(1.50000012,0,-5.00000060e-01),(4.99999911e-01,-1,5.00000060e-01))\0";
  24 static const string matrix_PM = "[3](2,2,2)";
  25
  26 int main () {
  27
  28   typedef float TYPE;
  29
  30   typedef matrix<TYPE> MATRIX;
  31
  32   MATRIX A;
  33   MATRIX LU;
  34   MATRIX INV;
  35
  36   {
  37     std::istringstream is(matrix_IN);
  38     is >> A;
  39   }
  40   {
  41     std::istringstream is(matrix_LU);
  42     is >> LU;
  43   }
  44   {
  45     std::istringstream is(matrix_INV);
  46     is >> INV;
  47   }
  48   permutation_matrix<>::vector_type temp;
  49   {
  50     std::istringstream is(matrix_PM);
  51     is >> temp;
  52   }
  53   permutation_matrix<> PM(temp);
  54
  55   permutation_matrix<> pm(3);
  56
  57   int result = lu_factorize<MATRIX, permutation_matrix<> >(A, pm);
  58
  59   assertTrue("factorization completed: ", 0 == result);
  60   assertTrue("LU factors are correct: ", compare(A, LU));
  61   assertTrue("permutation is correct: ", compare(pm, PM));
  62
  63   MATRIX B = identity_matrix<TYPE>(A.size2());
  64
  65   lu_substitute(A, pm, B);
  66
  67   assertTrue("inverse is correct: ", compare(B, INV));
  68
  69   return (getResults().second > 0) ? boost::exit_failure : boost::exit_success;
  70 }