Main Page | Namespace List | Data Structures | File List | Namespace Members | Data Fields | Globals

ex_residual.cpp

Go to the documentation of this file.
00001 // -*- c++ -*- //
00002 /** \file main.cpp  \brief main routines */
00003 /***************************************************************************
00004  -   begin                : 2003-09-05
00005  -   copyright            : (C) 2003 by Gunter Winkler
00006  -   email                : guwi17@gmx.de
00007  ***************************************************************************/
00008 
00009 /***************************************************************************
00010  *                                                                         *
00011  *   This program is free software; you can redistribute it and/or modify  *
00012  *   it under the terms of the GNU General Public License as published by  *
00013  *   the Free Software Foundation; either version 2 of the License, or     *
00014  *   (at your option) any later version.                                   *
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     /** \brief compute residual r = b - Ax
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 

Generated on Wed Oct 1 14:41:00 2003 for Sample Code by doxygen 1.3.2