|
ViennaCL - The Vienna Computing Library
1.4.2
|
00001 #ifndef VIENNACL_COMPRESSED_MATRIX_HPP_ 00002 #define VIENNACL_COMPRESSED_MATRIX_HPP_ 00003 00004 /* ========================================================================= 00005 Copyright (c) 2010-2013, Institute for Microelectronics, 00006 Institute for Analysis and Scientific Computing, 00007 TU Wien. 00008 Portions of this software are copyright by UChicago Argonne, LLC. 00009 00010 ----------------- 00011 ViennaCL - The Vienna Computing Library 00012 ----------------- 00013 00014 Project Head: Karl Rupp rupp@iue.tuwien.ac.at 00015 00016 (A list of authors and contributors can be found in the PDF manual) 00017 00018 License: MIT (X11), see file LICENSE in the base directory 00019 ============================================================================= */ 00020 00025 #include <vector> 00026 #include <list> 00027 #include <map> 00028 #include "viennacl/forwards.h" 00029 #include "viennacl/vector.hpp" 00030 00031 #include "viennacl/linalg/sparse_matrix_operations.hpp" 00032 00033 #include "viennacl/tools/tools.hpp" 00034 #include "viennacl/tools/entry_proxy.hpp" 00035 00036 namespace viennacl 00037 { 00038 namespace detail 00039 { 00040 template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT> 00041 void copy_impl(const CPU_MATRIX & cpu_matrix, 00042 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix, 00043 std::size_t nonzeros) 00044 { 00045 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1); 00046 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), nonzeros); 00047 std::vector<SCALARTYPE> elements(nonzeros); 00048 00049 std::size_t row_index = 0; 00050 std::size_t data_index = 0; 00051 00052 for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1(); 00053 row_it != cpu_matrix.end1(); 00054 ++row_it) 00055 { 00056 row_buffer.set(row_index, data_index); 00057 ++row_index; 00058 00059 for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin(); 00060 col_it != row_it.end(); 00061 ++col_it) 00062 { 00063 col_buffer.set(data_index, col_it.index2()); 00064 elements[data_index] = *col_it; 00065 ++data_index; 00066 } 00067 data_index = viennacl::tools::roundUpToNextMultiple<std::size_t>(data_index, ALIGNMENT); //take care of alignment 00068 } 00069 row_buffer.set(row_index, data_index); 00070 00071 gpu_matrix.set(row_buffer.get(), 00072 col_buffer.get(), 00073 &elements[0], 00074 cpu_matrix.size1(), 00075 cpu_matrix.size2(), 00076 nonzeros); 00077 } 00078 } 00079 00080 //provide copy-operation: 00095 template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT> 00096 void copy(const CPU_MATRIX & cpu_matrix, 00097 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix ) 00098 { 00099 assert( (gpu_matrix.size1() == 0 || cpu_matrix.size1() == gpu_matrix.size1()) && bool("Size mismatch") ); 00100 assert( (gpu_matrix.size2() == 0 || cpu_matrix.size2() == gpu_matrix.size2()) && bool("Size mismatch") ); 00101 00102 //std::cout << "copy for (" << cpu_matrix.size1() << ", " << cpu_matrix.size2() << ", " << cpu_matrix.nnz() << ")" << std::endl; 00103 00104 if ( cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0 ) 00105 { 00106 //determine nonzeros: 00107 long num_entries = 0; 00108 for (typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1(); 00109 row_it != cpu_matrix.end1(); 00110 ++row_it) 00111 { 00112 std::size_t entries_per_row = 0; 00113 for (typename CPU_MATRIX::const_iterator2 col_it = row_it.begin(); 00114 col_it != row_it.end(); 00115 ++col_it) 00116 { 00117 ++entries_per_row; 00118 } 00119 num_entries += viennacl::tools::roundUpToNextMultiple<std::size_t>(entries_per_row, ALIGNMENT); 00120 } 00121 00122 if (num_entries == 0) //we copy an empty matrix 00123 num_entries = 1; 00124 00125 //set up matrix entries: 00126 detail::copy_impl(cpu_matrix, gpu_matrix, num_entries); 00127 } 00128 } 00129 00130 00131 //adapted for std::vector< std::map < > > argument: 00137 template <typename SizeType, typename SCALARTYPE, unsigned int ALIGNMENT> 00138 void copy(const std::vector< std::map<SizeType, SCALARTYPE> > & cpu_matrix, 00139 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix ) 00140 { 00141 std::size_t nonzeros = 0; 00142 std::size_t max_col = 0; 00143 for (std::size_t i=0; i<cpu_matrix.size(); ++i) 00144 { 00145 nonzeros += cpu_matrix[i].size(); 00146 if (cpu_matrix[i].size() > 0) 00147 max_col = std::max<std::size_t>(max_col, (cpu_matrix[i].rbegin())->first); 00148 } 00149 00150 viennacl::detail::copy_impl(tools::const_sparse_matrix_adapter<SCALARTYPE, SizeType>(cpu_matrix, cpu_matrix.size(), max_col + 1), 00151 gpu_matrix, 00152 nonzeros); 00153 } 00154 00155 #ifdef VIENNACL_WITH_UBLAS 00156 template <typename ScalarType, typename F, std::size_t IB, typename IA, typename TA, unsigned int ALIGNMENT> 00157 void copy(const boost::numeric::ublas::compressed_matrix<ScalarType, F, IB, IA, TA> & ublas_matrix, 00158 viennacl::compressed_matrix<ScalarType, ALIGNMENT> & gpu_matrix) 00159 { 00160 //we just need to copy the CSR arrays: 00161 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), ublas_matrix.size1() + 1); 00162 for (std::size_t i=0; i<=ublas_matrix.size1(); ++i) 00163 row_buffer.set(i, ublas_matrix.index1_data()[i]); 00164 00165 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), ublas_matrix.nnz()); 00166 for (std::size_t i=0; i<ublas_matrix.nnz(); ++i) 00167 col_buffer.set(i, ublas_matrix.index2_data()[i]); 00168 00169 gpu_matrix.set(row_buffer.get(), 00170 col_buffer.get(), 00171 &(ublas_matrix.value_data()[0]), 00172 ublas_matrix.size1(), 00173 ublas_matrix.size2(), 00174 ublas_matrix.nnz()); 00175 00176 } 00177 #endif 00178 00179 #ifdef VIENNACL_WITH_EIGEN 00180 template <typename SCALARTYPE, int flags, unsigned int ALIGNMENT> 00181 void copy(const Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix, 00182 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix) 00183 { 00184 std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(eigen_matrix.rows()); 00185 00186 for (int k=0; k < eigen_matrix.outerSize(); ++k) 00187 for (typename Eigen::SparseMatrix<SCALARTYPE, flags>::InnerIterator it(eigen_matrix, k); it; ++it) 00188 stl_matrix[it.row()][it.col()] = it.value(); 00189 00190 copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix, eigen_matrix.rows(), eigen_matrix.cols()), gpu_matrix); 00191 } 00192 #endif 00193 00194 00195 #ifdef VIENNACL_WITH_MTL4 00196 template <typename SCALARTYPE, unsigned int ALIGNMENT> 00197 void copy(const mtl::compressed2D<SCALARTYPE> & cpu_matrix, 00198 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix) 00199 { 00200 typedef mtl::compressed2D<SCALARTYPE> MatrixType; 00201 00202 std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(cpu_matrix.num_rows()); 00203 00204 using mtl::traits::range_generator; 00205 using mtl::traits::range::min; 00206 00207 // Choose between row and column traversal 00208 typedef typename min<range_generator<mtl::tag::row, MatrixType>, 00209 range_generator<mtl::tag::col, MatrixType> >::type range_type; 00210 range_type my_range; 00211 00212 // Type of outer cursor 00213 typedef typename range_type::type c_type; 00214 // Type of inner cursor 00215 typedef typename mtl::traits::range_generator<mtl::tag::nz, c_type>::type ic_type; 00216 00217 // Define the property maps 00218 typename mtl::traits::row<MatrixType>::type row(cpu_matrix); 00219 typename mtl::traits::col<MatrixType>::type col(cpu_matrix); 00220 typename mtl::traits::const_value<MatrixType>::type value(cpu_matrix); 00221 00222 // Now iterate over the matrix 00223 for (c_type cursor(my_range.begin(cpu_matrix)), cend(my_range.end(cpu_matrix)); cursor != cend; ++cursor) 00224 for (ic_type icursor(mtl::begin<mtl::tag::nz>(cursor)), icend(mtl::end<mtl::tag::nz>(cursor)); icursor != icend; ++icursor) 00225 stl_matrix[row(*icursor)][col(*icursor)] = value(*icursor); 00226 00227 copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix, cpu_matrix.num_rows(), cpu_matrix.num_cols()), gpu_matrix); 00228 } 00229 #endif 00230 00231 00232 00233 00234 00235 00236 00237 // 00238 // gpu to cpu: 00239 // 00249 template <typename CPU_MATRIX, typename SCALARTYPE, unsigned int ALIGNMENT> 00250 void copy(const compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix, 00251 CPU_MATRIX & cpu_matrix ) 00252 { 00253 assert( (cpu_matrix.size1() == 0 || cpu_matrix.size1() == gpu_matrix.size1()) && bool("Size mismatch") ); 00254 assert( (cpu_matrix.size2() == 0 || cpu_matrix.size2() == gpu_matrix.size2()) && bool("Size mismatch") ); 00255 00256 if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 ) 00257 { 00258 if (cpu_matrix.size1() == 0 || cpu_matrix.size2() == 0) 00259 cpu_matrix.resize(gpu_matrix.size1(), gpu_matrix.size2(), false); 00260 00261 //get raw data from memory: 00262 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), cpu_matrix.size1() + 1); 00263 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz()); 00264 std::vector<SCALARTYPE> elements(gpu_matrix.nnz()); 00265 00266 //std::cout << "GPU->CPU, nonzeros: " << gpu_matrix.nnz() << std::endl; 00267 00268 viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get()); 00269 viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get()); 00270 viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0])); 00271 00272 //fill the cpu_matrix: 00273 std::size_t data_index = 0; 00274 for (std::size_t row = 1; row <= gpu_matrix.size1(); ++row) 00275 { 00276 while (data_index < row_buffer[row]) 00277 { 00278 if (col_buffer[data_index] >= gpu_matrix.size2()) 00279 { 00280 std::cerr << "ViennaCL encountered invalid data at colbuffer[" << data_index << "]: " << col_buffer[data_index] << std::endl; 00281 return; 00282 } 00283 00284 if (elements[data_index] != static_cast<SCALARTYPE>(0.0)) 00285 cpu_matrix(row-1, col_buffer[data_index]) = elements[data_index]; 00286 ++data_index; 00287 } 00288 } 00289 } 00290 } 00291 00292 00298 template <typename SCALARTYPE, unsigned int ALIGNMENT> 00299 void copy(const compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix, 00300 std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix) 00301 { 00302 tools::sparse_matrix_adapter<SCALARTYPE> temp(cpu_matrix, cpu_matrix.size(), cpu_matrix.size()); 00303 copy(gpu_matrix, temp); 00304 } 00305 00306 #ifdef VIENNACL_WITH_UBLAS 00307 template <typename ScalarType, unsigned int ALIGNMENT, typename F, std::size_t IB, typename IA, typename TA> 00308 void copy(viennacl::compressed_matrix<ScalarType, ALIGNMENT> const & gpu_matrix, 00309 boost::numeric::ublas::compressed_matrix<ScalarType> & ublas_matrix) 00310 { 00311 assert( (ublas_matrix.size1() == 0 || ublas_matrix.size1() == gpu_matrix.size1()) && bool("Size mismatch") ); 00312 assert( (ublas_matrix.size2() == 0 || ublas_matrix.size2() == gpu_matrix.size2()) && bool("Size mismatch") ); 00313 00314 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1); 00315 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz()); 00316 00317 viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get()); 00318 viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get()); 00319 00320 if (ublas_matrix.size1() == 0 || ublas_matrix.size2() == 0) 00321 ublas_matrix.resize(gpu_matrix.size1(), gpu_matrix.size2(), false); 00322 00323 ublas_matrix.clear(); 00324 ublas_matrix.reserve(gpu_matrix.nnz()); 00325 00326 ublas_matrix.set_filled(gpu_matrix.size1() + 1, gpu_matrix.nnz()); 00327 00328 for (std::size_t i=0; i<ublas_matrix.size1() + 1; ++i) 00329 ublas_matrix.index1_data()[i] = row_buffer[i]; 00330 00331 for (std::size_t i=0; i<ublas_matrix.nnz(); ++i) 00332 ublas_matrix.index2_data()[i] = col_buffer[i]; 00333 00334 viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(ScalarType) * gpu_matrix.nnz(), &(ublas_matrix.value_data()[0])); 00335 00336 } 00337 #endif 00338 00339 #ifdef VIENNACL_WITH_EIGEN 00340 template <typename SCALARTYPE, int flags, unsigned int ALIGNMENT> 00341 void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix, 00342 Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix) 00343 { 00344 if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 ) 00345 { 00346 assert(static_cast<unsigned int>(eigen_matrix.rows()) >= gpu_matrix.size1() 00347 && static_cast<unsigned int>(eigen_matrix.cols()) >= gpu_matrix.size2() 00348 && bool("Provided Eigen compressed matrix is too small!")); 00349 00350 //get raw data from memory: 00351 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1); 00352 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz()); 00353 std::vector<SCALARTYPE> elements(gpu_matrix.nnz()); 00354 00355 viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get()); 00356 viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get()); 00357 viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0])); 00358 00359 eigen_matrix.setZero(); 00360 std::size_t data_index = 0; 00361 for (std::size_t row = 1; row <= gpu_matrix.size1(); ++row) 00362 { 00363 while (data_index < row_buffer[row]) 00364 { 00365 assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer")); 00366 if (elements[data_index] != static_cast<SCALARTYPE>(0.0)) 00367 eigen_matrix.insert(row-1, col_buffer[data_index]) = elements[data_index]; 00368 ++data_index; 00369 } 00370 } 00371 } 00372 } 00373 #endif 00374 00375 00376 00377 #ifdef VIENNACL_WITH_MTL4 00378 template <typename SCALARTYPE, unsigned int ALIGNMENT> 00379 void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix, 00380 mtl::compressed2D<SCALARTYPE> & mtl4_matrix) 00381 { 00382 if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 ) 00383 { 00384 assert(mtl4_matrix.num_rows() >= gpu_matrix.size1() 00385 && mtl4_matrix.num_cols() >= gpu_matrix.size2() 00386 && bool("Provided MTL4 compressed matrix is too small!")); 00387 00388 //get raw data from memory: 00389 viennacl::backend::typesafe_host_array<unsigned int> row_buffer(gpu_matrix.handle1(), gpu_matrix.size1() + 1); 00390 viennacl::backend::typesafe_host_array<unsigned int> col_buffer(gpu_matrix.handle2(), gpu_matrix.nnz()); 00391 std::vector<SCALARTYPE> elements(gpu_matrix.nnz()); 00392 00393 viennacl::backend::memory_read(gpu_matrix.handle1(), 0, row_buffer.raw_size(), row_buffer.get()); 00394 viennacl::backend::memory_read(gpu_matrix.handle2(), 0, col_buffer.raw_size(), col_buffer.get()); 00395 viennacl::backend::memory_read(gpu_matrix.handle(), 0, sizeof(SCALARTYPE)* gpu_matrix.nnz(), &(elements[0])); 00396 00397 //set_to_zero(mtl4_matrix); 00398 //mtl4_matrix.change_dim(gpu_matrix.size1(), gpu_matrix.size2()); 00399 00400 mtl::matrix::inserter< mtl::compressed2D<SCALARTYPE> > ins(mtl4_matrix); 00401 std::size_t data_index = 0; 00402 for (std::size_t row = 1; row <= gpu_matrix.size1(); ++row) 00403 { 00404 while (data_index < row_buffer[row]) 00405 { 00406 assert(col_buffer[data_index] < gpu_matrix.size2() && bool("ViennaCL encountered invalid data at col_buffer")); 00407 if (elements[data_index] != static_cast<SCALARTYPE>(0.0)) 00408 ins(row-1, col_buffer[data_index]) << typename mtl::Collection< mtl::compressed2D<SCALARTYPE> >::value_type(elements[data_index]); 00409 ++data_index; 00410 } 00411 } 00412 } 00413 } 00414 #endif 00415 00416 00417 00418 00419 00421 00426 template<class SCALARTYPE, unsigned int ALIGNMENT /* see VCLForwards.h */> 00427 class compressed_matrix 00428 { 00429 public: 00430 typedef viennacl::backend::mem_handle handle_type; 00431 typedef scalar<typename viennacl::tools::CHECK_SCALAR_TEMPLATE_ARGUMENT<SCALARTYPE>::ResultType> value_type; 00432 typedef vcl_size_t size_type; 00433 00435 compressed_matrix() : rows_(0), cols_(0), nonzeros_(0) {} 00436 00443 explicit compressed_matrix(std::size_t rows, std::size_t cols, std::size_t nonzeros = 0) : rows_(rows), cols_(cols), nonzeros_(nonzeros) 00444 { 00445 if (rows > 0) 00446 { 00447 viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * (rows + 1)); 00448 } 00449 if (nonzeros > 0) 00450 { 00451 viennacl::backend::memory_create(col_buffer_, viennacl::backend::typesafe_host_array<unsigned int>().element_size() * nonzeros); 00452 viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * nonzeros); 00453 } 00454 } 00455 00456 #ifdef VIENNACL_WITH_OPENCL 00457 explicit compressed_matrix(cl_mem mem_row_buffer, cl_mem mem_col_buffer, cl_mem mem_elements, 00458 std::size_t rows, std::size_t cols, std::size_t nonzeros) : 00459 rows_(rows), cols_(cols), nonzeros_(nonzeros) 00460 { 00461 row_buffer_.switch_active_handle_id(viennacl::OPENCL_MEMORY); 00462 row_buffer_.opencl_handle() = mem_row_buffer; 00463 row_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed. 00464 row_buffer_.raw_size(sizeof(cl_uint) * (rows + 1)); 00465 00466 col_buffer_.switch_active_handle_id(viennacl::OPENCL_MEMORY); 00467 col_buffer_.opencl_handle() = mem_col_buffer; 00468 col_buffer_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed. 00469 col_buffer_.raw_size(sizeof(cl_uint) * nonzeros); 00470 00471 elements_.switch_active_handle_id(viennacl::OPENCL_MEMORY); 00472 elements_.opencl_handle() = mem_elements; 00473 elements_.opencl_handle().inc(); //prevents that the user-provided memory is deleted once the matrix object is destroyed. 00474 elements_.raw_size(sizeof(SCALARTYPE) * nonzeros); 00475 } 00476 #endif 00477 00478 00480 compressed_matrix & operator=(compressed_matrix const & other) 00481 { 00482 assert( (rows_ == 0 || rows_ == other.size1()) && bool("Size mismatch") ); 00483 assert( (cols_ == 0 || cols_ == other.size2()) && bool("Size mismatch") ); 00484 00485 rows_ = other.size1(); 00486 cols_ = other.size2(); 00487 nonzeros_ = other.nnz(); 00488 00489 viennacl::backend::typesafe_memory_copy<unsigned int>(other.row_buffer_, row_buffer_); 00490 viennacl::backend::typesafe_memory_copy<unsigned int>(other.col_buffer_, col_buffer_); 00491 viennacl::backend::typesafe_memory_copy<SCALARTYPE>(other.elements_, elements_); 00492 00493 return *this; 00494 } 00495 00496 00506 void set(const void * row_jumper, 00507 const void * col_buffer, 00508 const SCALARTYPE * elements, 00509 std::size_t rows, 00510 std::size_t cols, 00511 std::size_t nonzeros) 00512 { 00513 assert( (rows > 0) && bool("Error in compressed_matrix::set(): Number of rows must be larger than zero!")); 00514 assert( (cols > 0) && bool("Error in compressed_matrix::set(): Number of columns must be larger than zero!")); 00515 assert( (nonzeros > 0) && bool("Error in compressed_matrix::set(): Number of nonzeros must be larger than zero!")); 00516 //std::cout << "Setting memory: " << cols + 1 << ", " << nonzeros << std::endl; 00517 00518 //row_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY); 00519 viennacl::backend::memory_create(row_buffer_, viennacl::backend::typesafe_host_array<unsigned int>(row_buffer_).element_size() * (rows + 1), row_jumper); 00520 00521 //col_buffer_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY); 00522 viennacl::backend::memory_create(col_buffer_, viennacl::backend::typesafe_host_array<unsigned int>(col_buffer_).element_size() * nonzeros, col_buffer); 00523 00524 //elements_.switch_active_handle_id(viennacl::backend::OPENCL_MEMORY); 00525 viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * nonzeros, elements); 00526 00527 nonzeros_ = nonzeros; 00528 rows_ = rows; 00529 cols_ = cols; 00530 } 00531 00533 void reserve(std::size_t new_nonzeros) 00534 { 00535 if (new_nonzeros > nonzeros_) 00536 { 00537 handle_type col_buffer_old; 00538 handle_type elements_old; 00539 viennacl::backend::memory_shallow_copy(col_buffer_, col_buffer_old); 00540 viennacl::backend::memory_shallow_copy(elements_, elements_old); 00541 00542 viennacl::backend::typesafe_host_array<unsigned int> size_deducer(col_buffer_); 00543 viennacl::backend::memory_create(col_buffer_, size_deducer.element_size() * new_nonzeros); 00544 viennacl::backend::memory_create(elements_, sizeof(SCALARTYPE) * new_nonzeros); 00545 00546 viennacl::backend::memory_copy(col_buffer_old, col_buffer_, 0, 0, size_deducer.element_size() * nonzeros_); 00547 viennacl::backend::memory_copy(elements_old, elements_, 0, 0, sizeof(SCALARTYPE)* nonzeros_); 00548 00549 nonzeros_ = new_nonzeros; 00550 } 00551 } 00552 00559 void resize(std::size_t new_size1, std::size_t new_size2, bool preserve = true) 00560 { 00561 assert(new_size1 > 0 && new_size2 > 0 && bool("Cannot resize to zero size!")); 00562 00563 if (new_size1 != rows_ || new_size2 != cols_) 00564 { 00565 std::vector<std::map<unsigned int, SCALARTYPE> > stl_sparse_matrix; 00566 if (rows_ > 0) 00567 { 00568 if (preserve) 00569 { 00570 stl_sparse_matrix.resize(rows_); 00571 viennacl::copy(*this, stl_sparse_matrix); 00572 } else 00573 stl_sparse_matrix[0][0] = 0; 00574 } else { 00575 stl_sparse_matrix.resize(new_size1); 00576 stl_sparse_matrix[0][0] = 0; //enforces nonzero array sizes if matrix was initially empty 00577 } 00578 00579 stl_sparse_matrix.resize(new_size1); 00580 00581 //discard entries with column index larger than new_size2 00582 if (new_size2 < cols_ && rows_ > 0) 00583 { 00584 for (std::size_t i=0; i<stl_sparse_matrix.size(); ++i) 00585 { 00586 std::list<unsigned int> to_delete; 00587 for (typename std::map<unsigned int, SCALARTYPE>::iterator it = stl_sparse_matrix[i].begin(); 00588 it != stl_sparse_matrix[i].end(); 00589 ++it) 00590 { 00591 if (it->first >= new_size2) 00592 to_delete.push_back(it->first); 00593 } 00594 00595 for (std::list<unsigned int>::iterator it = to_delete.begin(); it != to_delete.end(); ++it) 00596 stl_sparse_matrix[i].erase(*it); 00597 } 00598 } 00599 00600 viennacl::copy(stl_sparse_matrix, *this); 00601 00602 rows_ = new_size1; 00603 cols_ = new_size2; 00604 } 00605 } 00606 00608 entry_proxy<SCALARTYPE> operator()(std::size_t i, std::size_t j) 00609 { 00610 assert( (i < rows_) && (j < cols_) && bool("compressed_matrix access out of bounds!")); 00611 00612 std::size_t index = element_index(i, j); 00613 00614 // check for element in sparsity pattern 00615 if (index < nonzeros_) 00616 return entry_proxy<SCALARTYPE>(index, elements_); 00617 00618 // Element not found. Copying required. Very slow, but direct entry manipulation is painful anyway... 00619 std::vector< std::map<unsigned int, SCALARTYPE> > cpu_backup(rows_); 00620 tools::sparse_matrix_adapter<SCALARTYPE> adapted_cpu_backup(cpu_backup, rows_, cols_); 00621 viennacl::copy(*this, adapted_cpu_backup); 00622 cpu_backup[i][j] = 0.0; 00623 viennacl::copy(adapted_cpu_backup, *this); 00624 00625 index = element_index(i, j); 00626 00627 assert(index < nonzeros_); 00628 00629 return entry_proxy<SCALARTYPE>(index, elements_); 00630 } 00631 00633 const std::size_t & size1() const { return rows_; } 00635 const std::size_t & size2() const { return cols_; } 00637 const std::size_t & nnz() const { return nonzeros_; } 00638 00640 const handle_type & handle1() const { return row_buffer_; } 00642 const handle_type & handle2() const { return col_buffer_; } 00644 const handle_type & handle() const { return elements_; } 00645 00647 handle_type & handle1() { return row_buffer_; } 00649 handle_type & handle2() { return col_buffer_; } 00651 handle_type & handle() { return elements_; } 00652 00653 void switch_memory_domain(viennacl::memory_types new_domain) 00654 { 00655 viennacl::backend::switch_memory_domain<unsigned int>(row_buffer_, new_domain); 00656 viennacl::backend::switch_memory_domain<unsigned int>(col_buffer_, new_domain); 00657 viennacl::backend::switch_memory_domain<SCALARTYPE>(elements_, new_domain); 00658 } 00659 00660 viennacl::memory_types memory_domain() const 00661 { 00662 return row_buffer_.get_active_handle_id(); 00663 } 00664 00665 private: 00666 00667 std::size_t element_index(std::size_t i, std::size_t j) 00668 { 00669 //read row indices 00670 viennacl::backend::typesafe_host_array<unsigned int> row_indices(row_buffer_, 2); 00671 viennacl::backend::memory_read(row_buffer_, row_indices.element_size()*i, row_indices.element_size()*2, row_indices.get()); 00672 00673 //get column indices for row i: 00674 viennacl::backend::typesafe_host_array<unsigned int> col_indices(col_buffer_, row_indices[1] - row_indices[0]); 00675 viennacl::backend::memory_read(col_buffer_, col_indices.element_size()*row_indices[0], row_indices.element_size()*col_indices.size(), col_indices.get()); 00676 00677 //get entries for row i: 00678 viennacl::backend::typesafe_host_array<SCALARTYPE> row_entries(elements_, row_indices[1] - row_indices[0]); 00679 viennacl::backend::memory_read(elements_, sizeof(SCALARTYPE)*row_indices[0], sizeof(SCALARTYPE)*row_entries.size(), row_entries.get()); 00680 00681 for (std::size_t k=0; k<col_indices.size(); ++k) 00682 { 00683 if (col_indices[k] == j) 00684 return row_indices[0] + k; 00685 } 00686 00687 // if not found, return index past the end of the matrix (cf. matrix.end() in the spirit of the STL) 00688 return nonzeros_; 00689 } 00690 00691 // /** @brief Copy constructor is by now not available. */ 00692 //compressed_matrix(compressed_matrix const &); 00693 00694 00695 std::size_t rows_; 00696 std::size_t cols_; 00697 std::size_t nonzeros_; 00698 handle_type row_buffer_; 00699 handle_type col_buffer_; 00700 handle_type elements_; 00701 }; 00702 00703 00704 00705 00706 } 00707 00708 #endif
1.7.6.1