00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #include <iostream>
00019
00020 #include <boost/timer.hpp>
00021
00022 #include <boost/numeric/ublas/matrix.hpp>
00023 #include <boost/numeric/ublas/matrix_sparse.hpp>
00024 #include <boost/numeric/ublas/vector.hpp>
00025 #include <boost/numeric/ublas/io.hpp>
00026 #include <boost/numeric/ublas/operation.hpp>
00027
00028
00029 #ifndef NDEBUG
00030 const size_t default_size = 100;
00031 #else
00032 const size_t default_size = 1000;
00033 #endif
00034
00035 const size_t subsize = 3;
00036
00037 using namespace boost::numeric::ublas ;
00038
00039 using std::cout;
00040 using std::endl;
00041
00042 namespace boost { namespace numeric { namespace ublas {
00043
00044
00045
00046
00047 template<class V, class T1, class IA1, class TA1, class E2>
00048 void
00049 residual ( const compressed_matrix<T1, row_major, 0, IA1, TA1> & A,
00050 const vector_expression<E2> &x,
00051 const vector_expression<E2> &b,
00052 V &r, row_major_tag) {
00053
00054 BOOST_UBLAS_CHECK( x().size() >= A.size2(), external_logic() );
00055 BOOST_UBLAS_CHECK( b().size() >= A.size1(), external_logic() );
00056 BOOST_UBLAS_CHECK( r.size() >= A.size1(), external_logic() );
00057
00058 typedef typename V::size_type size_type;
00059 typedef typename V::value_type value_type;
00060
00061 for (size_type i = 0; i < A.size1 (); ++ i) {
00062 size_type begin = A.index1_data () [i];
00063 size_type end = A.index1_data () [i + 1];
00064 value_type t (b() (i));
00065 for (size_type j = begin; j < end; ++ j)
00066 t -= A.value_data () [j] * x () (A.index2_data () [j]);
00067 r (i) = t;
00068 }
00069 return;
00070 }
00071
00072 }}}
00073
00074 template <class MAT>
00075 void matrix_fill(MAT& A)
00076 {
00077 typedef typename MAT::size_type size_type;
00078
00079 for (size_type i=0; i<A.size1(); ++i) {
00080 if (i>=1) { A.push_back(i,i-1,-1.0); }
00081 A.push_back(i,i,4);
00082 if (i+1<A.size2()) { A.push_back(i,i+1,-1.0); }
00083 }
00084 }
00085
00086 int main(int argc, char *argv[])
00087 {
00088
00089 typedef size_t size_type;
00090
00091 size_type size = default_size;
00092
00093 if (argc > 1)
00094 size = ::atoi (argv [1]);
00095
00096 cout << "size " << size << endl;
00097
00098 compressed_matrix<double, row_major> A(size,size);
00099 matrix_fill(A);
00100
00101 vector<double> b( scalar_vector<double>(size, 2.0) );
00102 b(0)=3;
00103 b(size-1)=3;
00104
00105 vector<double> x( scalar_vector<double>(size, 0.99) );
00106
00107 {
00108 boost::timer t;
00109 t.restart();
00110 vector<double> r(size);
00111 r = prod(A, x) - b;
00112 cout << "r=Ax-b : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
00113 }
00114
00115 {
00116 boost::timer t;
00117 t.restart();
00118 vector<double> r(size);
00119 noalias(r) = b - prod(A, x);
00120 cout << "noalias(r)=b-Ax : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
00121 }
00122
00123 {
00124 boost::timer t;
00125 t.restart();
00126 vector<double> r(b - prod(A, x));
00127 cout << "vector<double> r(b - prod(A, x)) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
00128 }
00129
00130 {
00131 boost::timer t;
00132 t.restart();
00133 vector<double> r(b);
00134 axpy_prod(A,-x,r,false);
00135 cout << "axpy_prod(A,-x,r,false) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
00136 }
00137
00138 {
00139 boost::timer t;
00140 t.restart();
00141 vector<double> r(size);
00142 residual(A,x,b,r,row_major_tag());
00143 cout << "residual(A,x,b,r,row_major_tag()) : " << t.elapsed() << " (" << norm_2(r) << ")" << endl;
00144 }
00145
00146 return EXIT_SUCCESS;
00147
00148 }
00149