00001 #include <iostream>
00002
00003
00004
00005
00006 #include <boost/numeric/ublas/matrix.hpp>
00007 #include <boost/numeric/ublas/matrix_sparse.hpp>
00008 #include <boost/numeric/ublas/vector.hpp>
00009 #include <boost/numeric/ublas/io.hpp>
00010
00011 using std::cout;
00012 using std::endl;
00013
00014 using namespace boost::numeric::ublas;
00015
00016 int main(int argc, char *argv[])
00017 {
00018
00019
00020 typedef compressed_matrix<double> MY_STIMA;
00021
00022 MY_STIMA M(4,4,16);
00023
00024 for (int i=0; i<M.size1(); ++i) M(i,i)=1.0;
00025 M(0,0)=3;
00026 M(0,1)=-1;
00027 M(0,2)=-2;
00028 M(1,0)=-1;
00029 M(1,1)=4;
00030 M(1,2)=-2;
00031 M(2,0)=-2;
00032 M(2,1)=-2;
00033 M(2,2)=9;
00034
00035 cout << "input M = " << M << "\n";
00036
00037 MY_STIMA::iterator1 r_it1 = M.begin1();
00038 MY_STIMA::iterator2 c_it1 = M.begin2();
00039 while ( (r_it1 != M.end1()) && (c_it1 != M.end2()) ) {
00040 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00041 MY_STIMA::iterator1 r_it2 = c_it1.begin();
00042 MY_STIMA::iterator2 c_it2 = r_it1.begin();
00043 #else
00044 MY_STIMA::iterator1 r_it2 = begin(c_it1, iterator2_tag());
00045 MY_STIMA::iterator2 c_it2 = begin(r_it1, iterator1_tag());
00046 #endif
00047 MY_STIMA::iterator2 c;
00048
00049 c = c_it2;
00050 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00051 while ( (r_it2 != c_it1.end()) && (c != r_it1.end()) ) {
00052 #else
00053 while ( (r_it2 != end(c_it1, iterator2_tag())) && (c != end(r_it1, iterator1_tag())) ) {
00054 #endif
00055 while (c.index2() < r_it2.index1()) {
00056 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00057 assert (c != r_it1.end());
00058 #else
00059 assert (c != end(r_it1, iterator1_tag()));
00060 #endif
00061 ++c;
00062 }
00063 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00064 assert (c != r_it1.end());
00065 #else
00066 assert (c != end(r_it1, iterator1_tag()));
00067 #endif
00068 while (c.index2() > r_it2.index1()) {
00069 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00070 assert (r_it2 != c_it1.end());
00071 #else
00072 assert (r_it2 != end(c_it1, iterator2_tag()));
00073 #endif
00074 ++r_it2;
00075 }
00076 #ifndef BOOST_UBLAS_NO_NESTED_CLASS_RELATION
00077 assert (r_it2 != c_it1.end());
00078 #else
00079 assert (r_it2 != end(c_it1, iterator2_tag()));
00080 #endif
00081 if (c.index2() == r_it2.index1()) {
00082 cout << r_it2.index1() << ", " << c.index2() << ": " << (*c) << " += " << (*r_it2);
00083 double temp = (*r_it2);
00084 (*c) += temp;
00085 cout << " --> " << (*c) << " ? \n";
00086 ++ r_it2;
00087 ++ c;
00088 }
00089 }
00090 ++ r_it1;
00091 ++ c_it1;
00092 }
00093
00094 cout << "output M = " << M << endl;
00095
00096 return EXIT_SUCCESS;
00097 }
00098