Skip to content
Newer
Older
100644 194 lines (178 sloc) 6.81 KB
167081e [spin.cpp] fix signs of magnetic field terms; use the flux!
mlenz authored Feb 2, 2009
1 // From http://www.guwi17.de/ublas/examples/sparse_io.cpp
2 #include <boost/numeric/ublas/vector.hpp>
3 #include <boost/numeric/ublas/matrix_sparse.hpp>
4 #include <boost/numeric/ublas/io.hpp>
5 #include <boost/numeric/ublas/operation.hpp>
6
7 #include <iostream>
8
9 namespace boost { namespace numeric { namespace ublas {
10
11 namespace io {
12
13 template<class ME>
14 class sparse_ioproxy {
15 public:
16 typedef ME matrix_expression_type;
17 typedef sparse_ioproxy<ME> self_type;
18 sparse_ioproxy(matrix_expression_type & me) : me(me) {}
19
20 matrix_expression_type & operator() () const {
21 return me;
22 }
23 matrix_expression_type & operator() () {
24 return me;
25 }
26
27 private:
28 sparse_ioproxy () {}
29 matrix_expression_type &me;
30 };
31
32 template < class ME >
33 sparse_ioproxy<ME> sparse(ME& me) {
34 return sparse_ioproxy<ME>(me);
35 }
36
37 template<class E, class T, class ME>
38 std::basic_ostream<E, T> &
39 output (std::basic_ostream<E, T> &os,
40 const sparse_ioproxy<ME> m,
41 row_major_tag) {
42
43 typedef typename ME::size_type size_type;
44 typedef ME ex_type;
45 size_type size1 = m () .size1 ();
46 size_type size2 = m () .size2 ();
47 std::basic_ostringstream<E, T, std::allocator<E> > s;
48 s.flags (os.flags ());
49 s.imbue (os.getloc ());
50 s.precision (os.precision ());
51 s << '[' << size1 << ',' << size2 << "]";
52 s << "(";
53 typename ex_type::const_iterator1 i = m().begin1();
54 typename ex_type::const_iterator1 n = m().end1();
55 bool one_back = (i != n);
56 for (; i != n; ++i) {
57 typename ex_type::const_iterator2 ri=i.begin();
58 typename ex_type::const_iterator2 rn=i.end();
59 bool nonempty = (ri != rn);
60 if (nonempty) {
61 s << "(" << ri.index1() << "," << ri.index2()
62 << ":" << *ri << ")";
63 ++ri;
64 }
65 for (; ri != rn; ++ri) {
66 s << ",(" << ri.index1() << "," << ri.index2()
67 << ":" << *ri << ")";
68 }
69 if (nonempty) { s << ';'; }
70 }
71 if (one_back) s.seekp(-1, std::ios_base::cur);
72 s << ')';
73 return os << s.str ().c_str ();
74 }
75
76 template<class E, class T, class ME>
77 std::basic_ostream<E, T> &
78 output (std::basic_ostream<E, T> &os,
79 const sparse_ioproxy<ME> m,
80 column_major_tag) {
81
82 typedef typename ME::size_type size_type;
83 typedef ME ex_type;
84 size_type size1 = m () .size1 ();
85 size_type size2 = m () .size2 ();
86 std::basic_ostringstream<E, T, std::allocator<E> > s;
87 s.flags (os.flags ());
88 s.imbue (os.getloc ());
89 s.precision (os.precision ());
90 s << '[' << size1 << ',' << size2 << "]";
91 s << "(";
92 typename ex_type::const_iterator2 i = m().begin2();
93 typename ex_type::const_iterator2 n = m().end2();
94 bool one_back = (i != n);
95 for (; i != n; ++i) {
96 typename ex_type::const_iterator1 ri=i.begin();
97 typename ex_type::const_iterator1 rn=i.end();
98 bool nonempty = (ri != rn);
99 if (nonempty) {
100 s << "(" << ri.index1() << "," << ri.index2()
101 << ":" << *ri << ")";
102 ++ri;
103 }
104 for (; ri != rn; ++ri) {
105 s << ",(" << ri.index1() << "," << ri.index2()
106 << ":" << *ri << ")";
107 }
108 if (nonempty) { s << ';'; }
109 }
110 if (one_back) s.seekp(-1, std::ios_base::cur);
111 s << ')';
112 return os << s.str ().c_str ();
113 }
114
115 template<class E, class T, class ME>
116 std::basic_ostream<E, T> &
117 operator << (std::basic_ostream<E, T> &os,
118 const sparse_ioproxy<ME> m) {
119 return output(os, m, typename ME::orientation_category());
120 }
121
122 // note: the orientation of the data must match the orientation
123 // of the matrix (otherwise push_back() will fail)
124 template<class E, class T, class ME>
125 std::basic_istream<E, T> &
126 operator >> (std::basic_istream<E, T> &is,
127 sparse_ioproxy<ME> m) {
128
129 typedef typename ME::size_type size_type;
130 typedef typename ME::value_type value_type;
131 typedef ME ex_type;
132 E ch;
133 size_type size1, size2;
134 if (is >> ch && ch != '[') {
135 is.putback (ch);
136 is.setstate (std::ios_base::failbit);
137 } else if (is >> size1 >> ch && ch != ',') {
138 is.putback (ch);
139 is.setstate (std::ios_base::failbit);
140 } else if (is >> size2 >> ch && ch != ']') {
141 is.putback (ch);
142 is.setstate (std::ios_base::failbit);
143 } else if (! is.fail ()) {
144 ME s(size1, size2);
145 if (is >> ch && ch != '(') {
146 is.putback (ch);
147 is.setstate (std::ios_base::failbit);
148 } else {
149 while ( ! is.fail() ) {
150 is >> ch;
151 if (ch == ')') {
152 break;
153 } else if (ch != '(') {
154 is.putback (ch);
155 is.setstate (std::ios_base::failbit);
156 } else {
157 // read triplet
158 size_type i,j;
159 value_type a;
160 is >> i >> ch;
161 if (is.fail() || ch != ',') {
162 is.putback (ch);
163 is.setstate (std::ios_base::failbit);
164 }
165 is >> j >> ch;
166 if (is.fail() || ch != ':') {
167 is.putback (ch);
168 is.setstate (std::ios_base::failbit);
169 }
170 is >> a >> ch;
171 if (is.fail() || ch != ')') {
172 is.putback (ch);
173 is.setstate (std::ios_base::failbit);
174 }
175 s.push_back(i,j,a);
176 is >> ch;
177 if ( ch == ')' ) {
178 break;
179 } else if ( ch != ',' && ch != ';' ) {
180 is.putback (ch);
181 is.setstate (std::ios_base::failbit);
182 }
183 }
184 }
185 }
186 if (! is.fail ())
187 m().assign_temporary (s);
188 }
189 return is;
190 }
191 }
192
193 }}}
Something went wrong with that request. Please try again.