Mercurial > fem-fenics-eugenio

comparison src/PETSc_factory.cc @ 268:61830a4f9ab9

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Improve formatting
author Eugenio Gianniti <eugenio.gianniti@mail.polimi.it>
date Thu, 14 Aug 2014 12:26:55 +0200
parents 8fe68d94ab76
children
comparison
equal deleted inserted replaced
267:53039ac90368 268:61830a4f9ab9
19 #ifdef LATEST_DOLFIN 19 #ifdef LATEST_DOLFIN
20 #include <boost/mpi.hpp> 20 #include <boost/mpi.hpp>
21 #include <boost/serialization/vector.hpp> 21 #include <boost/serialization/vector.hpp>
22 #endif 22 #endif
23 23
24 femfenics_base_factory const& 24 femfenics_base_factory const &
25 PETSc_factory::instance (void) 25 PETSc_factory::instance (void)
26 { 26 {
27 static PETSc_factory const theinstance; 27 static PETSc_factory const theinstance;
28 return theinstance; 28 return theinstance;
29 } 29 }
30 30
31 octave_value 31 octave_value
32 PETSc_factory::matrix (dolfin::Matrix const& A) const 32 PETSc_factory::matrix (dolfin::Matrix const & A) const
33 { 33 {
34 octave_value retval; 34 octave_value retval;
35 35
36 unsigned const size = 36 unsigned const size =
37 #ifdef LATEST_DOLFIN 37 #ifdef LATEST_DOLFIN
38 dolfin::MPI::size (MPI_COMM_WORLD); 38 dolfin::MPI::size (MPI_COMM_WORLD);
39 #else 39 #else
40 dolfin::MPI::num_processes (); 40 dolfin::MPI::num_processes ();
41 #endif 41 #endif
42 42
43 if (size > 1) 43 if (size > 1)
44 retval = do_matrix_parallel (A); 44 { retval = do_matrix_parallel (A); }
45 else 45 else
46 retval = do_matrix_serial (A); 46 { retval = do_matrix_serial (A); }
47 47
48 return retval; 48 return retval;
49 } 49 }
50 50
51 octave_value 51 octave_value
52 PETSc_factory::vector (dolfin::Vector const& b) const 52 PETSc_factory::vector (dolfin::Vector const & b) const
53 { 53 {
54 octave_value retval; 54 octave_value retval;
55 55
56 unsigned const size = 56 unsigned const size =
57 #ifdef LATEST_DOLFIN 57 #ifdef LATEST_DOLFIN
58 dolfin::MPI::size (MPI_COMM_WORLD); 58 dolfin::MPI::size (MPI_COMM_WORLD);
59 #else 59 #else
60 dolfin::MPI::num_processes (); 60 dolfin::MPI::num_processes ();
61 #endif 61 #endif
62 62
63 if (size > 1) 63 if (size > 1)
64 retval = do_vector_parallel (b); 64 { retval = do_vector_parallel (b); }
65 else 65 else
66 retval = do_vector_serial (b); 66 { retval = do_vector_serial (b); }
67 67
68 return retval; 68 return retval;
69 } 69 }
70 70
71 void 71 void
72 PETSc_factory::add_to_arrays (Array <octave_idx_type> & ridx, 72 PETSc_factory::add_to_arrays (Array <octave_idx_type> & ridx,
73 Array <octave_idx_type> & cidx, 73 Array <octave_idx_type> & cidx,
74 Array <double> & values, 74 Array <double> & values,
75 std::vector <std::size_t> const& row, 75 std::vector <std::size_t> const & row,
76 std::vector <std::size_t> const& col, 76 std::vector <std::size_t> const & col,
77 std::vector <double> const& val) 77 std::vector <double> const & val)
78 { 78 {
79 octave_idx_type pos = cidx.numel (); 79 octave_idx_type pos = cidx.numel ();
80 80
81 dim_vector dims = cidx.dims (); 81 dim_vector dims = cidx.dims ();
82 if (dims(1) == 1) 82 if (dims(1) == 1)
83 dims(0) += col.size (); 83 { dims(0) += col.size (); }
84 else 84 else
85 dims(1) += col.size (); 85 { dims(1) += col.size (); }
86 ridx.resize (dims, 0); 86 ridx.resize (dims, 0);
87 cidx.resize (dims, 0); 87 cidx.resize (dims, 0);
88 values.resize (dims, 0); 88 values.resize (dims, 0);
89 89
90 std::vector <std::size_t>::const_iterator rit = row.begin (), 90 std::vector <std::size_t>::const_iterator rit = row.begin (),
91 cit = col.begin (); 91 cit = col.begin ();
92 std::vector <double>::const_iterator vit = val.begin (); 92 std::vector <double>::const_iterator vit = val.begin ();
93 93
94 while (cit != col.end ()) 94 while (cit != col.end ())
95 { 95 {
96 ridx(pos) = *rit; 96 ridx(pos) = *rit;
97 cidx(pos) = *cit; 97 cidx(pos) = *cit;
98 values(pos) = *vit; 98 values(pos) = *vit;
99 99
100 ++rit; 100 ++rit;
101 ++cit; 101 ++cit;
102 ++vit; 102 ++vit;
103 ++pos; 103 ++pos;
104 } 104 }
105 105
106 return; 106 return;
107 } 107 }
108 108
109 octave_value 109 octave_value
110 PETSc_factory::do_matrix_serial (dolfin::Matrix const& A) const 110 PETSc_factory::do_matrix_serial (dolfin::Matrix const & A) const
111 { 111 {
112 octave_value retval; 112 octave_value retval;
113 113
114 std::size_t nr = A.size (0), nc = A.size (1); 114 std::size_t nr = A.size (0), nc = A.size (1);
115 115
116 dim_vector dims (0, 1); 116 dim_vector dims (0, 1);
117 Array <octave_idx_type> 117 Array <octave_idx_type> ridx (dims, 0), cidx (dims, 0);
118 ridx (dims, 0), 118 Array<double> data (dims, 0);
119 cidx (dims, 0); 119
120 Array<double> data (dims, 0); 120 for (std::size_t i = 0; i < nr; ++i)
121 121 {
122 for (std::size_t i = 0; i < nr; ++i) 122 std::vector <double> data_tmp;
123 { 123 std::vector <std::size_t> cidx_tmp;
124 std::vector <double> data_tmp; 124 A.getrow (i, cidx_tmp, data_tmp);
125 std::vector <std::size_t> cidx_tmp; 125 std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i);
126 A.getrow (i, cidx_tmp, data_tmp); 126 add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp);
127 std::vector <std::size_t> ridx_tmp (cidx_tmp.size (), i); 127 }
128 add_to_arrays (ridx, cidx, data, ridx_tmp, cidx_tmp, data_tmp); 128
129 } 129 SparseMatrix sm (data, ridx, cidx, nr, nc);
130 130 retval = sm;
131 SparseMatrix sm (data, ridx, cidx, nr, nc); 131
132 retval = sm; 132 return retval;
133 133 }
134 return retval; 134
135 } 135 octave_value
136 136 PETSc_factory::do_matrix_parallel (dolfin::Matrix const & A) const
137 octave_value 137 {
138 PETSc_factory::do_matrix_parallel (dolfin::Matrix const& A) const 138 octave_value retval;
139 { 139
140 octave_value retval; 140 std::vector <std::size_t> rows, cols;
141 141 std::vector <double> vals;
142 std::vector <std::size_t> rows, cols; 142 std::size_t const nrows = A.size (0), ncols = A.size (1);
143 std::vector <double> vals; 143
144 std::size_t const nrows = A.size (0), ncols = A.size (1); 144 std::pair <std::size_t, std::size_t> const range = A.local_range (0);
145 145 for (std::size_t i = range.first; i < range.second; ++i)
146 std::pair <std::size_t, std::size_t> const range = A.local_range (0); 146 {
147 for (std::size_t i = range.first; i < range.second; ++i) 147 std::vector <std::size_t> cols_tmp;
148 { 148 std::vector <double> vals_tmp;
149 std::vector <std::size_t> cols_tmp; 149 A.getrow (i, cols_tmp, vals_tmp);
150 std::vector <double> vals_tmp; 150 cols.insert (cols.end (), cols_tmp.begin (), cols_tmp.end ());
151 A.getrow (i, cols_tmp, vals_tmp); 151 vals.insert (vals.end (), vals_tmp.begin (), vals_tmp.end ());
152 cols.insert (cols.end (), cols_tmp.begin (), cols_tmp.end ()); 152 for (std::size_t j = 0; j < cols_tmp.size (); ++j)
153 vals.insert (vals.end (), vals_tmp.begin (), vals_tmp.end ()); 153 { rows.push_back (i); }
154 for (std::size_t j = 0; j < cols_tmp.size (); ++j) 154 }
155 rows.push_back (i); 155
156 } 156 boost::mpi::communicator world;
157 157
158 boost::mpi::communicator world; 158 unsigned const size = world.size ();
159 159 unsigned const rank = world.rank ();
160 unsigned const size = world.size (); 160
161 unsigned const rank = world.rank (); 161 for (unsigned p = 1; p < size; ++p)
162 162 {
163 for (unsigned p = 1; p < size; ++p) 163 if (rank == p)
164 { 164 {
165 if (rank == p) 165 unsigned const tag = 3 * (p - 1);
166 { 166 world.send (0, tag, rows);
167 unsigned const tag = 3*(p-1); 167 world.send (0, tag + 1, cols);
168 world.send (0, tag, rows); 168 world.send (0, tag + 2, vals);
169 world.send (0, tag+1, cols); 169 }
170 world.send (0, tag+2, vals); 170 }
171 } 171
172 } 172 if (rank == 0)
173 173 {
174 if (rank == 0) 174 dim_vector dims (0, 1);
175 { 175 Array <octave_idx_type> ridx (dims, 0), cidx (dims, 0);
176 dim_vector dims (0, 1); 176 Array <double> data (dims, 0);
177 Array <octave_idx_type> ridx (dims, 0), cidx (dims, 0); 177 add_to_arrays (ridx, cidx, data, rows, cols, vals);
178 Array <double> data (dims, 0); 178
179 add_to_arrays (ridx, cidx, data, rows, cols, vals); 179 for (unsigned p = 1; p < size; ++p)
180 180 {
181 for (unsigned p = 1; p < size; ++p) 181 std::vector <double> new_vals;
182 { 182 std::vector <std::size_t> new_rows, new_cols;
183 std::vector <double> new_vals; 183 unsigned const tag = 3 * (p - 1);
184 std::vector <std::size_t> new_rows, new_cols; 184
185 unsigned const tag = 3*(p-1); 185 world.recv (p, tag, new_rows);
186 186 world.recv (p, tag + 1, new_cols);
187 world.recv (p, tag, new_rows); 187 world.recv (p, tag + 2, new_vals);
188 world.recv (p, tag+1, new_cols); 188
189 world.recv (p, tag+2, new_vals); 189 add_to_arrays (ridx, cidx, data, new_rows, new_cols, new_vals);
190 190 }
191 add_to_arrays (ridx, cidx, data, new_rows, new_cols, new_vals); 191
192 } 192 SparseMatrix sm (data, ridx, cidx, nrows, ncols);
193 193 retval = sm;
194 SparseMatrix sm (data, ridx, cidx, nrows, ncols); 194 }
195 retval = sm; 195 else
196 } 196 // Return an identity matrix just to avoid warnings or errors
197 else 197 {
198 // Return an identity matrix just to avoid warnings or errors 198 dim_vector dims (nrows, 1);
199 { 199 Array <octave_idx_type> ridx (dims), cidx (dims);
200 dim_vector dims (nrows, 1); 200 Array <double> data (dims, 1);
201 Array <octave_idx_type> ridx (dims), cidx (dims); 201
202 Array <double> data (dims, 1); 202 for (std::size_t i = 0; i < nrows; ++i)
203 203 {
204 for (std::size_t i = 0; i < nrows; ++i) 204 ridx.xelem (i) = i;
205 { 205 cidx.xelem (i) = i;
206 ridx.xelem (i) = i; 206 }
207 cidx.xelem (i) = i; 207
208 } 208 SparseMatrix sm (data, ridx, cidx, nrows, ncols);
209 209 retval = sm;
210 SparseMatrix sm (data, ridx, cidx, nrows, ncols); 210 }
211 retval = sm; 211
212 } 212 return retval;
213 213 }
214 return retval; 214
215 } 215 octave_value
216 216 PETSc_factory::do_vector_serial (dolfin::Vector const & b) const
217 octave_value 217 {
218 PETSc_factory::do_vector_serial (dolfin::Vector const& b) const 218 octave_value retval;
219 { 219
220 octave_value retval; 220 dim_vector dims;
221 221 dims.resize (2);
222 dim_vector dims; 222 dims(0) = b.size ();
223 dims.resize (2); 223 dims(1) = 1;
224 dims(0) = b.size (); 224 Array <double> myb (dims);
225 dims(1) = 1; 225
226 Array <double> myb (dims); 226 for (std::size_t i = 0; i < b.size (); ++i)
227 227 { myb.xelem (i) = b[i]; }
228 for (std::size_t i = 0; i < b.size (); ++i) 228
229 myb.xelem (i) = b[i]; 229 retval = myb;
230 230
231 retval = myb; 231 return retval;
232 232 }
233 return retval; 233
234 } 234 octave_value
235 235 PETSc_factory::do_vector_parallel (dolfin::Vector const & b) const
236 octave_value 236 {
237 PETSc_factory::do_vector_parallel (dolfin::Vector const& b) const 237 octave_value retval;
238 { 238
239 octave_value retval; 239 std::size_t const length = b.size ();
240 240 dim_vector dims (length, 1);
241 std::size_t const length = b.size (); 241 Array <double> vec (dims, 0);
242 dim_vector dims (length, 1); 242
243 Array <double> vec (dims, 0); 243 std::vector <double> entries;
244 244 b.gather_on_zero (entries);
245 std::vector <double> entries; 245
246 b.gather_on_zero (entries); 246 #ifdef LATEST_DOLFIN
247 247 if (dolfin::MPI::rank (MPI_COMM_WORLD) == 0)
248 #ifdef LATEST_DOLFIN
249 if (dolfin::MPI::rank (MPI_COMM_WORLD) == 0)
250 #else 248 #else
251 if (dolfin::MPI::process_number () == 0) 249 if (dolfin::MPI::process_number () == 0)
252 #endif 250 #endif
253 for (std::size_t i = 0; i < length; ++i) 251 for (std::size_t i = 0; i < length; ++i)
254 vec.xelem (i) = entries[i]; 252 { vec.xelem (i) = entries[i]; }
255 253
256 // On nodes other than zero, return an all zero vector of the right length 254 // On nodes other than zero, return an all zero vector of the right length
257 retval = vec; 255 retval = vec;
258 256
259 return retval; 257 return retval;
260 } 258 }