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

   1 #include <iostream>
   2 #include <stdlib.h>
   3 #include <cmath>
   4
   5 #include <boost/numeric/ublas/vector.hpp>
   6 #include <boost/numeric/ublas/matrix.hpp>
   7 #include <boost/numeric/ublas/matrix_sparse.hpp>
   8 #include <boost/numeric/ublas/triangular.hpp>
   9 #include <boost/numeric/ublas/io.hpp>
  10
  11 #include <boost/timer/timer.hpp>
  12
  13 namespace ublas  = boost::numeric::ublas;
  14
  15 template<class mat, class vec>
  16 double diff(const mat& A, const vec& x, const vec& b) {
  17   vec temp(prod(A, x) - b);
  18   double result = 0;
  19   for (typename vec::size_type i=0; i<temp.size(); ++i) {
  20     result += temp(i)*temp(i);
  21   }
  22   return std::sqrt(result);
  23 }
  24
  25 template<class mat, class vec>
  26 double diff(const vec& x, const mat& A, const vec& b) {
  27   return diff(trans(A), x, b);
  28 }
  29
  30 namespace ublas  = boost::numeric::ublas;
  31
  32
  33 int main() {
  34   const int n=7000;
  35 #if 1
  36   ublas::compressed_matrix<double, ublas::row_major>     mat_row_upp(n, n);
  37   ublas::compressed_matrix<double, ublas::column_major>  mat_col_upp(n, n);
  38   ublas::compressed_matrix<double, ublas::row_major>     mat_row_low(n, n);
  39   ublas::compressed_matrix<double, ublas::column_major>  mat_col_low(n, n);
  40 #else
  41   ublas::matrix<double, ublas::row_major>     mat_row_upp(n, n, 0);
  42   ublas::matrix<double, ublas::column_major>  mat_col_upp(n, n, 0);
  43   ublas::matrix<double, ublas::row_major>     mat_row_low(n, n, 0);
  44   ublas::matrix<double, ublas::column_major>  mat_col_low(n, n, 0);
  45 #endif
  46   ublas::vector<double>  b(n, 1);
  47
  48   std::cerr << "Constructing..." << std::endl;
  49   for (int i=0; i<n; ++i) {
  50     b(i) = std::rand() % 10;
  51     double main = -10 + std::rand() % 20 ;
  52     if (main == 0) main+=1;
  53     double side = -10 + std::rand() % 20 ;
  54     if (i-1>=0) {
  55       mat_row_low(i, i-1) = side;
  56     }
  57     mat_row_low(i, i) = main;
  58
  59     mat_col_low(i, i) = main;
  60     if (i+1<n) {
  61       mat_col_low(i+1, i) = side;
  62     }
  63
  64     mat_row_upp(i, i) = main;
  65     if (i+1<n) {
  66       mat_row_upp(i, i+1) = side;
  67     }
  68
  69     if (i-1>=0) {
  70       mat_col_upp(i-1, i) = side;
  71     }
  72     mat_col_upp(i, i) = main;
  73   }
  74
  75   std::cerr << "Starting..." << std::endl;
  76   {
  77     boost::timer::auto_cpu_timer t(std::cerr, "col_low x: %t sec CPU, %w sec real\n");
  78     ublas::vector<double>  x(b);
  79     ublas::inplace_solve(mat_col_low,  x, ublas::lower_tag());
  80     std::cerr << "delta: " << diff(mat_col_low, x, b) << "\n";
  81   }
  82   {
  83     boost::timer::auto_cpu_timer t(std::cerr, "row_low x: %t sec CPU, %w sec real\n");
  84     ublas::vector<double>  x(b);
  85     ublas::inplace_solve(mat_row_low, x, ublas::lower_tag());
  86     std::cerr << "delta: " << diff(mat_row_low, x, b) << "\n";
  87   }
  88
  89   {
  90     boost::timer::auto_cpu_timer t(std::cerr, "col_upp x: %t sec CPU, %w sec real\n");
  91     ublas::vector<double>  x(b);
  92     ublas::inplace_solve(mat_col_upp,  x, ublas::upper_tag());
  93     std::cerr << "delta: " << diff(mat_col_upp, x, b) << "\n";
  94   }
  95   {
  96     boost::timer::auto_cpu_timer t(std::cerr, "row_upp x: %t sec CPU, %w sec real\n");
  97     ublas::vector<double>  x(b);
  98     ublas::inplace_solve(mat_row_upp, x, ublas::upper_tag());
  99     std::cerr << "delta: " << diff(mat_row_upp, x, b) << "\n";
 100   }
 101
 102   {
 103     boost::timer::auto_cpu_timer t(std::cerr, "x col_low: %t sec CPU, %w sec real\n");
 104     ublas::vector<double>  x(b);
 105     ublas::inplace_solve(x, mat_col_low, ublas::lower_tag());
 106     std::cerr << "delta: " << diff(x, mat_col_low, b) << "\n";
 107   }
 108   {
 109     boost::timer::auto_cpu_timer t(std::cerr, "x row_low: %t sec CPU, %w sec real\n");
 110     ublas::vector<double>  x(b);
 111     ublas::inplace_solve(x, mat_row_low, ublas::lower_tag());
 112     std::cerr << "delta: " << diff(x, mat_row_low, b) << "\n";
 113   }
 114
 115   {
 116     boost::timer::auto_cpu_timer t(std::cerr, "x col_upp: %t sec CPU, %w sec real\n");
 117     ublas::vector<double>  x(b);
 118     ublas::inplace_solve(x, mat_col_upp, ublas::upper_tag());
 119     std::cerr << "delta: " << diff(x, mat_col_upp, b) << "\n";
 120   }
 121   {
 122     boost::timer::auto_cpu_timer t(std::cerr, "x row_upp: %t sec CPU, %w sec real\n");
 123     ublas::vector<double>  x(b);
 124     ublas::inplace_solve(x, mat_row_upp, ublas::upper_tag());
 125     std::cerr << "delta: " << diff(x, mat_row_upp, b) << "\n";
 126   }
 127
 128
 129 }