|
Skylark (Sketching Library)
0.1
|
00001 #ifndef SKYLARK_HASH_TRANSFORM_ELEMENTAL_HPP 00002 #define SKYLARK_HASH_TRANSFORM_ELEMENTAL_HPP 00003 00004 #include "../utility/get_communicator.hpp" 00005 00006 namespace skylark { namespace sketch { 00007 00011 template <typename ValueType, 00012 template <typename> class IdxDistributionType, 00013 template <typename> class ValueDistribution> 00014 struct hash_transform_t < 00015 elem::Matrix<ValueType>, 00016 elem::Matrix<ValueType>, 00017 IdxDistributionType, 00018 ValueDistribution > : 00019 public hash_transform_data_t<IdxDistributionType, 00020 ValueDistribution> { 00021 00022 // Typedef matrix and distribution types so that we can use them regularly 00023 typedef ValueType value_type; 00024 typedef elem::Matrix<value_type> matrix_type; 00025 typedef elem::Matrix<value_type> output_matrix_type; 00026 typedef IdxDistributionType<size_t> idx_distribution_type; 00027 typedef ValueDistribution<value_type> value_distribution_type; 00028 typedef hash_transform_data_t<IdxDistributionType, 00029 ValueDistribution> data_type; 00033 hash_transform_t (int N, int S, base::context_t& context) : 00034 data_type (N, S, context) { 00035 00036 } 00037 00041 hash_transform_t (const hash_transform_t<matrix_type, 00042 output_matrix_type, 00043 IdxDistributionType, 00044 ValueDistribution>& other) : 00045 data_type(other) {} 00046 00050 hash_transform_t (const data_type& other_data) : 00051 data_type(other_data) {} 00052 00056 template <typename Dimension> 00057 void apply (const matrix_type& A, output_matrix_type& sketch_of_A, 00058 Dimension dimension) const { 00059 try { 00060 apply_impl(A, sketch_of_A, dimension); 00061 } catch (std::logic_error e) { 00062 SKYLARK_THROW_EXCEPTION ( 00063 base::elemental_exception() 00064 << base::error_msg(e.what()) ); 00065 } catch(boost::mpi::exception e) { 00066 SKYLARK_THROW_EXCEPTION ( 00067 base::mpi_exception() 00068 << base::error_msg(e.what()) ); 00069 } 00070 } 00071 00072 int get_N() const { return this->_N; } 00073 int get_S() const { return this->_S; } 00075 const sketch_transform_data_t* get_data() const { return this; } 00076 00077 private: 00078 00083 void apply_impl (const matrix_type& A, 00084 output_matrix_type& sketch_of_A, 00085 skylark::sketch::columnwise_tag) const { 00086 00087 elem::Zero(sketch_of_A); 00088 00089 // Construct Pi * A (directly on the fly) 00090 for (size_t row_idx = 0; row_idx < A.Height(); row_idx++) { 00091 00092 size_t new_row_idx = data_type::row_idx[row_idx]; 00093 value_type scale_factor = data_type::row_value[row_idx]; 00094 00095 for(size_t col_idx = 0; col_idx < A.Width(); col_idx++) { 00096 value_type value = scale_factor * A.Get(row_idx, col_idx); 00097 sketch_of_A.Update(new_row_idx, col_idx, value); 00098 } 00099 } 00100 } 00101 00106 void apply_impl (const matrix_type& A, 00107 output_matrix_type& sketch_of_A, 00108 skylark::sketch::rowwise_tag) const { 00109 00110 elem::Zero(sketch_of_A); 00111 00112 // Construct Pi * A (directly on the fly) 00113 for (size_t col_idx = 0; col_idx < A.Width(); col_idx++) { 00114 00115 size_t new_col_idx = data_type::row_idx[col_idx]; 00116 value_type scale_factor = data_type::row_value[col_idx]; 00117 00118 for(size_t row_idx = 0; row_idx < A.Height(); row_idx++) { 00119 value_type value = scale_factor * A.Get(row_idx, col_idx); 00120 sketch_of_A.Update(row_idx, new_col_idx, value); 00121 } 00122 } 00123 } 00124 }; 00125 00129 template <typename ValueType, 00130 template <typename> class IdxDistributionType, 00131 template <typename> class ValueDistribution> 00132 struct hash_transform_t < 00133 base::sparse_matrix_t<ValueType>, 00134 elem::Matrix<ValueType>, 00135 IdxDistributionType, 00136 ValueDistribution > : 00137 public hash_transform_data_t<IdxDistributionType, 00138 ValueDistribution> { 00139 00140 // Typedef matrix and distribution types so that we can use them regularly 00141 typedef ValueType value_type; 00142 typedef base::sparse_matrix_t<ValueType> matrix_type; 00143 typedef elem::Matrix<value_type> output_matrix_type; 00144 typedef IdxDistributionType<size_t> idx_distribution_type; 00145 typedef ValueDistribution<value_type> value_distribution_type; 00146 typedef hash_transform_data_t<IdxDistributionType, 00147 ValueDistribution> data_type; 00151 hash_transform_t (int N, int S, base::context_t& context) : 00152 data_type (N, S, context) { 00153 00154 } 00155 00159 hash_transform_t (const hash_transform_t<matrix_type, 00160 output_matrix_type, 00161 IdxDistributionType, 00162 ValueDistribution>& other) : 00163 data_type(other) {} 00164 00168 hash_transform_t (const data_type& other_data) : 00169 data_type(other_data) {} 00170 00174 template <typename Dimension> 00175 void apply (const matrix_type& A, output_matrix_type& sketch_of_A, 00176 Dimension dimension) const { 00177 try { 00178 apply_impl(A, sketch_of_A, dimension); 00179 } catch (std::logic_error e) { 00180 SKYLARK_THROW_EXCEPTION ( 00181 base::elemental_exception() 00182 << base::error_msg(e.what()) ); 00183 } catch(boost::mpi::exception e) { 00184 SKYLARK_THROW_EXCEPTION ( 00185 base::mpi_exception() 00186 << base::error_msg(e.what()) ); 00187 } 00188 } 00189 00190 int get_N() const { return this->_N; } 00191 int get_S() const { return this->_S; } 00193 const sketch_transform_data_t* get_data() const { return this; } 00194 00195 private: 00196 00201 void apply_impl (const matrix_type& A, 00202 output_matrix_type& sketch_of_A, 00203 skylark::sketch::columnwise_tag) const { 00204 00205 elem::Zero(sketch_of_A); 00206 00207 double *SA = sketch_of_A.Buffer(); 00208 int ld = sketch_of_A.LDim(); 00209 00210 const int* indptr = A.indptr(); 00211 const int* indices = A.indices(); 00212 const value_type* values = A.locked_values(); 00213 00214 # if SKYLARK_HAVE_OPENMP 00215 # pragma omp parallel for 00216 # endif 00217 for(int col = 0; col < A.width(); col++) { 00218 for (int j = indptr[col]; j < indptr[col + 1]; j++) { 00219 int row = indices[j]; 00220 value_type val = values[j]; 00221 SA[col * ld + data_type::row_idx[row]] += 00222 data_type::row_value[row] * val; 00223 } 00224 } 00225 } 00226 00231 void apply_impl (const matrix_type& A, 00232 output_matrix_type& sketch_of_A, 00233 skylark::sketch::rowwise_tag) const { 00234 00235 elem::Zero(sketch_of_A); 00236 00237 double *SA = sketch_of_A.Buffer(); 00238 int ld = sketch_of_A.LDim(); 00239 00240 const int* indptr = A.indptr(); 00241 const int* indices = A.indices(); 00242 const value_type* values = A.locked_values(); 00243 00244 for(int col = 0; col < A.width(); col++) { 00245 # if SKYLARK_HAVE_OPENMP 00246 # pragma omp parallel for 00247 # endif 00248 for (int j = indptr[col]; j < indptr[col + 1]; j++) { 00249 int row = indices[j]; 00250 value_type val = values[j]; 00251 SA[data_type::row_idx[col] * ld + row] += 00252 data_type::row_value[col] * val; 00253 } 00254 } 00255 00256 } 00257 }; 00258 00262 template <typename ValueType, 00263 elem::Distribution ColDist, 00264 elem::Distribution RowDist, 00265 template <typename> class IdxDistributionType, 00266 template <typename> class ValueDistribution> 00267 struct hash_transform_t < 00268 elem::DistMatrix<ValueType, ColDist, RowDist>, 00269 elem::Matrix<ValueType>, 00270 IdxDistributionType, 00271 ValueDistribution > : 00272 public hash_transform_data_t<IdxDistributionType, 00273 ValueDistribution> { 00274 00275 // Typedef matrix and distribution types so that we can use them regularly 00276 typedef ValueType value_type; 00277 typedef elem::DistMatrix<value_type, ColDist, RowDist> matrix_type; 00278 typedef elem::Matrix<value_type> output_matrix_type; 00279 typedef IdxDistributionType<size_t> idx_distribution_type; 00280 typedef ValueDistribution<value_type> value_distribution_type; 00281 typedef hash_transform_data_t<IdxDistributionType, 00282 ValueDistribution> data_type; 00286 hash_transform_t (int N, int S, base::context_t& context) : 00287 data_type (N, S, context) { 00288 00289 } 00290 00294 hash_transform_t (const hash_transform_t<matrix_type, 00295 output_matrix_type, 00296 IdxDistributionType, 00297 ValueDistribution>& other) : 00298 data_type(other) {} 00299 00303 hash_transform_t (const data_type& other_data) : 00304 data_type(other_data) {} 00305 00309 template <typename Dimension> 00310 void apply (const matrix_type& A, output_matrix_type& sketch_of_A, 00311 Dimension dimension) const { 00312 try { 00313 apply_impl(A, sketch_of_A, dimension); 00314 } catch (std::logic_error e) { 00315 SKYLARK_THROW_EXCEPTION ( 00316 base::elemental_exception() 00317 << base::error_msg(e.what()) ); 00318 } catch(boost::mpi::exception e) { 00319 SKYLARK_THROW_EXCEPTION ( 00320 base::mpi_exception() 00321 << base::error_msg(e.what()) ); 00322 } 00323 } 00324 00325 private: 00330 void apply_impl (const matrix_type& A, 00331 output_matrix_type& sketch_of_A, 00332 skylark::sketch::columnwise_tag) const { 00333 00334 // TODO this implementation is not communication efficient. 00335 // Sketching a nxd matrix to sxd will communicate O(sdP) 00336 // doubles, when you can sometime communicate less: 00337 // For [MC,MR] or [MR,MC] you need O(sd sqrt(P)). 00338 // For [*, VC/VR] you need only O(sd). 00339 00340 // Create space to hold local part of SA 00341 elem::Matrix<value_type> SA_part (sketch_of_A.Height(), 00342 sketch_of_A.Width(), 00343 sketch_of_A.LDim()); 00344 00345 elem::Zero(SA_part); 00346 00347 // Construct Pi * A (directly on the fly) 00348 for (size_t j = 0; j < A.LocalHeight(); j++) { 00349 00350 size_t row_idx = A.ColShift() + A.ColStride() * j; 00351 size_t new_row_idx = data_type::row_idx[row_idx]; 00352 value_type scale_factor = data_type::row_value[row_idx]; 00353 00354 for(size_t i = 0; i < A.LocalWidth(); i++) { 00355 size_t col_idx = A.RowShift() + A.RowStride() * i; 00356 value_type value = scale_factor * A.GetLocal(j, i); 00357 SA_part.Update(new_row_idx, col_idx, value); 00358 } 00359 } 00360 00361 // Pull everything to rank-0 00362 boost::mpi::reduce (utility::get_communicator(A), 00363 SA_part.LockedBuffer(), 00364 SA_part.MemorySize(), 00365 sketch_of_A.Buffer(), 00366 std::plus<value_type>(), 00367 0); 00368 } 00369 00374 void apply_impl (const matrix_type& A, 00375 output_matrix_type& sketch_of_A, 00376 skylark::sketch::rowwise_tag) const { 00377 00378 // TODO this implementation is not communication efficient. 00379 // Sketching a nxd matrix to sxd will communicate O(sdP) 00380 // doubles, when you can sometime communicate less: 00381 // For [MC,MR] or [MR,MC] you need O(sd sqrt(P)). 00382 // For [VC/VR, *] you need only O(sd). 00383 00384 // Create space to hold local part of SA 00385 elem::Matrix<value_type> SA_part (sketch_of_A.Height(), 00386 sketch_of_A.Width(), 00387 sketch_of_A.LDim()); 00388 00389 elem::Zero(SA_part); 00390 00391 // Construct A * Pi (directly on the fly) 00392 for (size_t j = 0; j < A.LocalWidth(); ++j) { 00393 00394 size_t col_idx = A.RowShift() + A.RowStride() * j; 00395 size_t new_col_idx = data_type::row_idx[col_idx]; 00396 value_type scale_factor = data_type::row_value[col_idx]; 00397 00398 for(size_t i = 0; i < A.LocalHeight(); ++i) { 00399 size_t row_idx = A.ColShift() + A.ColStride() * i; 00400 value_type value = scale_factor * A.GetLocal(i, j); 00401 SA_part.Update(row_idx, new_col_idx, value); 00402 } 00403 } 00404 00405 // Pull everything to rank-0 00406 boost::mpi::reduce (skylark::utility::get_communicator(A), 00407 SA_part.LockedBuffer(), 00408 SA_part.MemorySize(), 00409 sketch_of_A.Buffer(), 00410 std::plus<value_type>(), 00411 0); 00412 } 00413 }; 00414 00418 template <typename ValueType, 00419 elem::Distribution ColDist, 00420 elem::Distribution RowDist, 00421 template <typename> class IdxDistributionType, 00422 template <typename> class ValueDistribution> 00423 struct hash_transform_t < 00424 elem::DistMatrix<ValueType, ColDist, RowDist>, 00425 elem::DistMatrix<ValueType, elem::STAR, elem::STAR>, 00426 IdxDistributionType, 00427 ValueDistribution > : 00428 public hash_transform_data_t<IdxDistributionType, 00429 ValueDistribution>, 00430 virtual public sketch_transform_t<elem::DistMatrix<ValueType, ColDist, 00431 RowDist>, 00432 elem::DistMatrix<ValueType, 00433 elem::STAR, 00434 elem::STAR> > { 00435 00436 // Typedef matrix and distribution types so that we can use them regularly 00437 typedef ValueType value_type; 00438 typedef elem::DistMatrix<value_type, ColDist, RowDist> matrix_type; 00439 typedef elem::DistMatrix<value_type, elem::STAR, elem::STAR> output_matrix_type; 00440 typedef IdxDistributionType<size_t> idx_distribution_type; 00441 typedef ValueDistribution<value_type> value_distribution_type; 00442 typedef hash_transform_data_t<IdxDistributionType, 00443 ValueDistribution> data_type; 00447 hash_transform_t (int N, int S, base::context_t& context) : 00448 data_type (N, S, context) { 00449 00450 } 00451 00455 hash_transform_t (const hash_transform_t<matrix_type, 00456 output_matrix_type, 00457 IdxDistributionType, 00458 ValueDistribution>& other) : 00459 data_type(other) {} 00460 00464 hash_transform_t (const data_type& other_data) : 00465 data_type(other_data) {} 00466 00471 void apply (const matrix_type& A, output_matrix_type& sketch_of_A, 00472 columnwise_tag dimension) const { 00473 try { 00474 apply_impl(A, sketch_of_A, dimension); 00475 } catch (std::logic_error e) { 00476 SKYLARK_THROW_EXCEPTION ( 00477 base::elemental_exception() 00478 << base::error_msg(e.what()) ); 00479 } catch(boost::mpi::exception e) { 00480 SKYLARK_THROW_EXCEPTION ( 00481 base::mpi_exception() 00482 << base::error_msg(e.what()) ); 00483 } 00484 } 00485 00490 void apply (const matrix_type& A, output_matrix_type& sketch_of_A, 00491 rowwise_tag dimension) const { 00492 try { 00493 apply_impl(A, sketch_of_A, dimension); 00494 } catch (std::logic_error e) { 00495 SKYLARK_THROW_EXCEPTION ( 00496 base::elemental_exception() 00497 << base::error_msg(e.what()) ); 00498 } catch(boost::mpi::exception e) { 00499 SKYLARK_THROW_EXCEPTION ( 00500 base::mpi_exception() 00501 << base::error_msg(e.what()) ); 00502 } 00503 } 00504 00505 int get_N() const { return this->_N; } 00506 int get_S() const { return this->_S; } 00508 const sketch_transform_data_t* get_data() const { return this; } 00509 00510 private: 00515 void apply_impl (const matrix_type& A, 00516 output_matrix_type& sketch_of_A, 00517 skylark::sketch::columnwise_tag) const { 00518 00519 // TODO this implementation is not communication efficient. 00520 // Sketching a nxd matrix to sxd will communicate O(sdP) 00521 // doubles, when you can sometime communicate less: 00522 // For [MC,MR] or [MR,MC] you need O(sd sqrt(P)). 00523 // For [*, VC/VR] you need only O(sd). 00524 00525 // Create space to hold local part of SA 00526 elem::Matrix<value_type> SA_part (sketch_of_A.Height(), 00527 sketch_of_A.Width(), 00528 sketch_of_A.LDim()); 00529 00530 elem::Zero(SA_part); 00531 00532 // Construct Pi * A (directly on the fly) 00533 for (size_t j = 0; j < A.LocalHeight(); j++) { 00534 00535 size_t row_idx = A.ColShift() + A.ColStride() * j; 00536 size_t new_row_idx = data_type::row_idx[row_idx]; 00537 value_type scale_factor = data_type::row_value[row_idx]; 00538 00539 for(size_t i = 0; i < A.LocalWidth(); i++) { 00540 size_t col_idx = A.RowShift() + A.RowStride() * i; 00541 value_type value = scale_factor * A.GetLocal(j, i); 00542 SA_part.Update(new_row_idx, col_idx, value); 00543 } 00544 } 00545 00546 boost::mpi::all_reduce (utility::get_communicator(A), 00547 SA_part.LockedBuffer(), 00548 SA_part.MemorySize(), 00549 sketch_of_A.Buffer(), 00550 std::plus<value_type>()); 00551 } 00552 00557 void apply_impl (const matrix_type& A, 00558 output_matrix_type& sketch_of_A, 00559 skylark::sketch::rowwise_tag) const { 00560 00561 // TODO this implementation is not communication efficient. 00562 // Sketching a nxd matrix to sxd will communicate O(sdP) 00563 // doubles, when you can sometime communicate less: 00564 // For [MC,MR] or [MR,MC] you need O(sd sqrt(P)). 00565 // For [VC/VR, *] you need only O(sd). 00566 00567 // Create space to hold local part of SA 00568 elem::Matrix<value_type> SA_part (sketch_of_A.Height(), 00569 sketch_of_A.Width(), 00570 sketch_of_A.LDim()); 00571 00572 elem::Zero(SA_part); 00573 00574 // Construct A * Pi (directly on the fly) 00575 for (size_t j = 0; j < A.LocalWidth(); ++j) { 00576 00577 size_t col_idx = A.RowShift() + A.RowStride() * j; 00578 size_t new_col_idx = data_type::row_idx[col_idx]; 00579 value_type scale_factor = data_type::row_value[col_idx]; 00580 00581 for(size_t i = 0; i < A.LocalHeight(); ++i) { 00582 size_t row_idx = A.ColShift() + A.ColStride() * i; 00583 value_type value = scale_factor * A.GetLocal(i, j); 00584 SA_part.Update(row_idx, new_col_idx, value); 00585 } 00586 } 00587 00588 // Pull everything to rank-0 00589 boost::mpi::all_reduce (utility::get_communicator(A), 00590 SA_part.LockedBuffer(), 00591 SA_part.MemorySize(), 00592 sketch_of_A.Buffer(), 00593 std::plus<value_type>()); 00594 } 00595 }; 00596 00597 } } 00599 #endif // SKYLARK_HASH_TRANSFORM_ELEMENTAL_HPP
1.7.6.1