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
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
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;
00057 (*c) /= (*kk);
00058 ++c; ++kk;
00059 while (kk != rk.end()) {
00060
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 {
00084 os << " ";
00085 }
00086 }
00087 os << "\n";
00088 }
00089 os << "\n";
00090 }
00091
00092
00093
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
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128 }
00129 }
00130
00131
00132
00133
00134
00135
00136
00137
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
00165
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
00177
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
00203
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
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
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
00264
00265 cout << "start solver\n";
00266
00267
00268 gettimeofday( &t[3], NULL );
00269 inplace_solve(triangular_adaptor<const MY_STIMA, unit_lower>(M), y, unit_lower_tag());
00270
00271 gettimeofday( &t[4], NULL );
00272 inplace_solve(M, y, upper_tag());
00273
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 };