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

trisolve.cpp

Go to the documentation of this file.
00001 
00002 #include <sys/time.h>
00003 #include <unistd.h>
00004 
00005 #include <iostream>
00006 
00007 #include <boost/numeric/ublas/matrix.hpp>
00008 #include <boost/numeric/ublas/matrix_sparse.hpp>
00009 #include <boost/numeric/ublas/triangular.hpp>
00010 #include <boost/numeric/ublas/vector.hpp>
00011 #include <boost/numeric/ublas/io.hpp>
00012 #include <boost/numeric/ublas/operation.hpp>
00013 
00014 using std::cout;
00015 using std::endl;
00016 
00017 using namespace boost::numeric::ublas;
00018 
00019 /*
00020  * LU for matrix A
00021  *
00022  * for i = 2:n
00023  *   for k = 1:i-1
00024  *     A(i,k) = A(i,k) / A(k,k);
00025  *     for j = k+1:n
00026  *       A(i,j) = A(i,j) - A(i,k)*A(k,j)
00027  *     next j
00028  *   next k
00029  * next i
00030  *
00031  * we assume: all diagonal elements are nonzero
00032  */
00033     
00034 template <class Matrix>
00035 void do_lu(Matrix& LU)
00036 {
00037   typedef  typename Matrix::size_type   size_type;
00038 
00039   typename Matrix::iterator1 r = LU.begin1();
00040   typename Matrix::iterator1 r_end = LU.end1();
00041 
00042   ++r;
00043   while (r != r_end) {
00044         size_type  i = r.index1();
00045         typename Matrix::iterator1 rk = LU.begin1();
00046         typename Matrix::iterator1 rk_end = LU.end1();
00047         typename Matrix::iterator2 c_start = r.begin();
00048         typename Matrix::iterator2 c_end = r.end();
00049         typename Matrix::iterator2 c;
00050         typename Matrix::iterator2 kk;
00051         while (c_start.index2() < i) {
00052           size_type  k = c_start.index2();
00053           while (rk.index1() < k) ++rk;
00054           c = c_start;
00055           kk = rk.begin();
00056           while (kk.index2() < k) ++kk; // better use find2(...) ?
00057           (*c) /= (*kk);
00058           ++c; ++kk;
00059           while (kk != rk.end()) {
00060                 // insert elements if necessary
00061                 LU(c.index1(),kk.index2()) -= (*c_start) * (*kk);
00062                 ++kk;
00063           }
00064           ++c_start;
00065         }
00066         ++r;
00067   }
00068 }
00069 
00070 template <class Matrix>
00071 void spy(std::ostream& os, const Matrix& A)
00072 {
00073   const double eps = 1e-10;
00074   typename Matrix::size_type i,j,m,n;
00075   m = A.size1();
00076   n = A.size2();
00077   for (i=0; i<m; ++i) {
00078         for (j=0; j<n; ++j) {
00079           if (A(i,j) < -eps) {
00080                 os << "-";
00081           } else if (A(i,j) > eps) {
00082                 os << "+";
00083           } else /* if A(i,j) == 0 */ {
00084                 os << " ";
00085           }
00086         }
00087         os << "\n";
00088   }
00089   os << "\n";
00090 }
00091 
00092 // Sparse (proxy) case, 
00093 // row-major version, iterators must have correct .index2()
00094 template<class E1, class E2>
00095 BOOST_UBLAS_INLINE
00096 void inplace_solve (const matrix_expression<E1> &e1, E2 &e2,
00097                                         unit_lower_tag, vector_tag, unknown_storage_tag) {
00098   typedef BOOST_UBLAS_TYPENAME E2::size_type size_type;
00099   typedef BOOST_UBLAS_TYPENAME E2::difference_type difference_type;
00100   typedef BOOST_UBLAS_TYPENAME E2::value_type value_type;
00101 
00102   BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
00103   BOOST_UBLAS_CHECK (e1 ().size2 () == e2.size (), bad_size ());
00104   size_type size = e2.size ();
00105 
00106   typename E1::const_iterator1 row;
00107   typename E1::const_iterator1 row_end = e1().end1();
00108 
00109   for (row = e1().begin1(); row != row_end; ++ row) {
00110         size_type n = row.index1();
00111 
00112         typename E1::const_iterator2 col = row.begin();
00113         typename E1::const_iterator2 col_end = row.end();
00114         while ((col != col_end) && (col.index2() < n) ) {
00115           e2 (n) -= *col * e2 (col.index2());
00116           ++ col;
00117         }
00118         /* only if e1()(n,n) != 1 
00119         // assert( col.index2() == n )
00120 #ifndef BOOST_UBLAS_SINGULAR_CHECK
00121         BOOST_UBLAS_CHECK ( (*col) != value_type (), singular ());
00122 #else
00123         if ( (*col) == value_type ())
00124           singular ().raise ();
00125 #endif
00126         e2 (n) /= (*col); 
00127         */
00128   }
00129 }
00130 
00131 
00132 /* reverse iterators look very good, but have a problem:
00133    index calculation takes time, because their native index
00134    is shifted by one. ( revit.begin == it.end == 1 behind last! )
00135 */
00136 // Sparse (proxy) case, 
00137 // row-major version, iterators must have correct .index2()
00138 template<class E1, class E2>
00139 BOOST_UBLAS_INLINE
00140 void inplace_solve1 (const matrix_expression<E1> &e1, E2 &e2,
00141                                          upper_tag, vector_tag, unknown_storage_tag ) {
00142   typedef BOOST_UBLAS_TYPENAME E2::size_type size_type;
00143   typedef BOOST_UBLAS_TYPENAME E2::difference_type difference_type;
00144   typedef BOOST_UBLAS_TYPENAME E2::value_type value_type;
00145 
00146   BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
00147   BOOST_UBLAS_CHECK (e1 ().size2 () == e2.size (), bad_size ());
00148   size_type size = e2.size ();
00149 
00150   typename E1::const_iterator1 row = e1().end1();
00151   typename E1::const_iterator1 row_end = e1().begin1();
00152 
00153   while ( row != row_end ) {
00154         --row;
00155         size_type n = row.index1();
00156 
00157         typename E1::const_iterator2 col = row.end();
00158         typename E1::const_iterator2 col_end = row.begin();
00159         while ((col != col_end) && (col.index2() > n+1) ) {
00160           -- col;
00161           e2 (n) -= *col * e2 (col.index2());
00162         }
00163         --col;
00164         /* only if e1()(n,n) != 1 */
00165         // assert( col.index2() == n )
00166 #ifndef BOOST_UBLAS_SINGULAR_CHECK
00167         BOOST_UBLAS_CHECK ( (*col) != value_type (), singular ());
00168 #else
00169         if ( (*col) == value_type ())
00170           singular ().raise ();
00171 #endif
00172         e2 (n) /= (*col); 
00173   }
00174 }
00175 
00176 // Sparse (proxy) case, 
00177 // row-major version, iterators must have correct .index2()
00178 template<class E1, class E2>
00179 BOOST_UBLAS_INLINE
00180 void inplace_solve2 (const matrix_expression<E1> &e1, E2 &e2,
00181                                         upper_tag, vector_tag, unknown_storage_tag) {
00182   typedef BOOST_UBLAS_TYPENAME E2::size_type size_type;
00183   typedef BOOST_UBLAS_TYPENAME E2::difference_type difference_type;
00184   typedef BOOST_UBLAS_TYPENAME E2::value_type value_type;
00185 
00186   BOOST_UBLAS_CHECK (e1 ().size1 () == e1 ().size2 (), bad_size ());
00187   BOOST_UBLAS_CHECK (e1 ().size2 () == e2.size (), bad_size ());
00188   size_type size = e2.size ();
00189 
00190   typename E1::const_reverse_iterator1 row;
00191   typename E1::const_reverse_iterator1 row_end = e1().rend1();
00192 
00193   for (row = e1().rbegin1(); row != row_end; ++ row) {
00194         size_type n = row.index1();
00195 
00196         typename E1::const_reverse_iterator2 col = row.rbegin();
00197         typename E1::const_reverse_iterator2 col_end = row.rend();
00198         while ((col != col_end) && (col.index2() > n) ) {
00199           e2 (n) -= *col * e2 (col.index2());
00200           ++ col;
00201         }
00202         /* only if e1()(n,n) != 1  */
00203         // assert( col.index2() == n )
00204 #ifndef BOOST_UBLAS_SINGULAR_CHECK
00205         BOOST_UBLAS_CHECK ( (*col) != value_type (), singular ());
00206 #else
00207         if ( (*col) == value_type ())
00208           singular ().raise ();
00209 #endif
00210         e2 (n) /= (*col); 
00211   }
00212 }
00213 
00214 
00215 int main(int argc, char *argv[])
00216 {
00217 
00218   typedef sparse_matrix<double, column_major>  MY_STIMA;
00219   typedef vector<double>         MY_VEC;
00220 
00221   const int COUNT= 2;
00222 
00223   const int size = 2000;
00224   const int sub1 = 1;
00225   const int sub2 = 10;
00226   const int sup1 = 1;
00227   const int sup2 = 10;
00228 
00229   timeval t[7];
00230 
00231   MY_STIMA M(size,size,(sup2+1+sub2)*size);
00232   MY_VEC   x(size);
00233   MY_VEC   y(size);
00234 
00235   for (int i=sup2; i<M.size1(); ++i)  M(i-sup2,i     ) =-1.0;
00236   for (int i=sup1; i<M.size1(); ++i)  M(i-sup1,i     ) =-1.0;
00237   for (int i=0;    i<M.size1(); ++i)  M(i     ,i     ) = 4.0;
00238   for (int i=sub1; i<M.size1(); ++i)  M(i     ,i-sub1) =-1.0;
00239   for (int i=sub2; i<M.size1(); ++i)  M(i     ,i-sub2) =-1.0;
00240 
00241   //cout << "input M = " << M << "\n";
00242   if (size < 72) {
00243         cout << "matrix spy of input matrix\n";
00244         spy(cout, M);
00245   }
00246 
00247   for (int i=0; i<x.size(); ++i)  x(i) = 1.0;
00248 
00249   gettimeofday( &t[0], NULL );
00250   for (int i=0; i<COUNT; ++i) {
00251         axpy_prod(M, x, y, true);
00252   }
00253   gettimeofday( &t[1], NULL );
00254 
00255   /* uncomment LU if you want to check the results, but my LU is slow ... */
00256   do_lu(M);
00257   gettimeofday( &t[2], NULL );
00258 
00259   if (size < 72) {
00260         cout << "matrix spy after LU decomposition\n";
00261         spy(cout, M);
00262   }
00263   // cout << "output M = " << M << endl;
00264   
00265   cout << "start solver\n";
00266 
00267   // why can't I use inplace_solve(M, y, unit_low_tag()) ?
00268   gettimeofday( &t[3], NULL );
00269   inplace_solve(triangular_adaptor<const MY_STIMA, unit_lower>(M), y, unit_lower_tag());
00270   // inplace_solve(M, y, unit_lower_tag(), vector_tag(), unknown_storage_tag());
00271   gettimeofday( &t[4], NULL );
00272   inplace_solve(M, y, upper_tag());
00273   //inplace_solve1(M, y, upper_tag(), vector_tag(), unknown_storage_tag());
00274   gettimeofday( &t[5], NULL );
00275   
00276   cout << "L2 error : " << norm_2(x-y) << "\n";
00277 
00278   cout << "### running times " << size << endl;
00279   double t1,t2;
00280   t1=(double(t[0].tv_sec)+t[0].tv_usec/1000000.0);
00281   t2=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
00282   cout << "### axpy_prod   : " << (t2-t1)/COUNT << " s " << endl;
00283   t1=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
00284   t2=(double(t[2].tv_sec)+t[2].tv_usec/1000000.0);
00285   cout << "### compute LU  : " << (t2-t1) << " s " << endl;
00286   t1=(double(t[3].tv_sec)+t[3].tv_usec/1000000.0);
00287   t2=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
00288   cout << "### solving unit lower system : " << (t2-t1) << " s " << endl;
00289   t1=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
00290   t2=(double(t[5].tv_sec)+t[5].tv_usec/1000000.0);
00291   cout << "### solving upper system      : " << (t2-t1) << " s " << endl;
00292 
00293   return EXIT_SUCCESS;
00294 };

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