00001
00002 #include <sys/time.h>
00003 #include <unistd.h>
00004
00005 #include <iostream>
00006
00007 #include "mtl/mtl.h"
00008 #include "mtl/matrix.h"
00009 #include "mtl/mtl_algo.h"
00010
00011 using std::cout;
00012 using std::endl;
00013
00014 using namespace mtl;
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 template <class Matrix>
00032 void do_lu(Matrix& LU)
00033 {
00034 typedef typename Matrix::size_type size_type;
00035
00036 typename Matrix::iterator r = LU.begin();
00037 typename Matrix::iterator r_end = LU.end();
00038
00039 ++r;
00040 while (r != r_end) {
00041 size_type i = r.index();
00042 typename Matrix::iterator rk = LU.begin();
00043 typename Matrix::iterator rk_end = LU.end();
00044 typename Matrix::OneD::iterator c_start = (*r).begin();
00045 typename Matrix::OneD::iterator c_end = (*r).end();
00046 typename Matrix::OneD::iterator c;
00047 typename Matrix::OneD::iterator kk;
00048 while (c_start.index() < i) {
00049 size_type k = c_start.index();
00050 while (rk.index() < k) ++rk;
00051 c = c_start;
00052 kk = (*rk).begin();
00053 while (kk.index() < k) ++kk;
00054 (*c) /= (*kk);
00055 ++c; ++kk;
00056 while (kk != (*rk).end()) {
00057
00058 LU(r.index(),kk.index()) -= (*c_start) * (*kk);
00059 ++kk;
00060 }
00061 ++c_start;
00062 }
00063 ++r;
00064 }
00065 }
00066
00067 template <class Matrix>
00068 void spy(std::ostream& os, const Matrix& A)
00069 {
00070 const double eps = 1e-10;
00071 typename Matrix::size_type i,j,m,n;
00072 m = A.nrows();
00073 n = A.ncols();
00074 for (i=0; i<m; ++i) {
00075 for (j=0; j<n; ++j) {
00076 if (A(i,j) < -eps) {
00077 os << "-";
00078 } else if (A(i,j) > eps) {
00079 os << "+";
00080 } else {
00081 os << " ";
00082 }
00083 }
00084 os << "\n";
00085 }
00086 os << "\n";
00087 }
00088
00089 int main(int argc, char *argv[])
00090 {
00091
00092 typedef mtl::matrix< double,
00093 rectangle<>,
00094 mtl::array< mtl::compressed<> >,
00095 row_major >::type MY_STIMA;
00096
00097 typedef mtl::dense1D<double> MY_VEC;
00098
00099 const int COUNT= 2;
00100
00101 const int size = 20;
00102 const int sub1 = 1;
00103 const int sub2 = 10;
00104 const int sup1 = 1;
00105 const int sup2 = 10;
00106
00107 timeval t[7];
00108
00109
00110 MY_STIMA M(size,size);
00111 MY_VEC x(size);
00112 MY_VEC y(size);
00113 set(y,0.0);
00114
00115 for (int i=sup2; i<size; ++i) M(i-sup2,i ) =-1.0;
00116 for (int i=sup1; i<size; ++i) M(i-sup1,i ) =-1.0;
00117 for (int i=0; i<size; ++i) M(i ,i ) = 4.0;
00118 for (int i=sub1; i<size; ++i) M(i ,i-sub1) =-1.0;
00119 for (int i=sub2; i<size; ++i) M(i ,i-sub2) =-1.0;
00120
00121
00122 if (size < 72) {
00123 cout << "matrix spy of input matrix\n";
00124 spy(cout, M);
00125 }
00126
00127 for (int i=0; i<x.size(); ++i) x[i] = 1.0;
00128
00129 gettimeofday( &t[0], NULL );
00130 for (int i=0; i<COUNT; ++i) {
00131 mult(M, x, y);
00132 }
00133 gettimeofday( &t[1], NULL );
00134
00135
00136 do_lu(M);
00137 gettimeofday( &t[2], NULL );
00138
00139 if (size < 72) {
00140 cout << "matrix spy after LU decomposition\n";
00141 spy(cout, M);
00142 }
00143
00144
00145 cout << "start solver\n";
00146
00147
00148 gettimeofday( &t[3], NULL );
00149 tri_solve(triangle_view<MY_STIMA, unit_lower>::type (M), y);
00150 print_vector(y);
00151 gettimeofday( &t[4], NULL );
00152 tri_solve(triangle_view<MY_STIMA, upper>::type (M), y);
00153 print_vector(y);
00154 gettimeofday( &t[5], NULL );
00155
00156 add(scaled(y, -1.0), x);
00157 cout << "L2 error : " << (two_norm(x)) << "\n";
00158
00159 cout << "### running times " << size << endl;
00160 double t1,t2;
00161 t1=(double(t[0].tv_sec)+t[0].tv_usec/1000000.0);
00162 t2=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
00163 cout << "### axpy_prod : " << (t2-t1)/COUNT << " s " << endl;
00164 t1=(double(t[1].tv_sec)+t[1].tv_usec/1000000.0);
00165 t2=(double(t[2].tv_sec)+t[2].tv_usec/1000000.0);
00166 cout << "### compute LU : " << (t2-t1) << " s " << endl;
00167 t1=(double(t[3].tv_sec)+t[3].tv_usec/1000000.0);
00168 t2=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
00169 cout << "### solving unit lower system : " << (t2-t1) << " s " << endl;
00170 t1=(double(t[4].tv_sec)+t[4].tv_usec/1000000.0);
00171 t2=(double(t[5].tv_sec)+t[5].tv_usec/1000000.0);
00172 cout << "### solving upper system : " << (t2-t1) << " s " << endl;
00173
00174 return EXIT_SUCCESS;
00175 };