1 #ifndef VIENNACL_COMPRESSED_MATRIX_HPP_
2 #define VIENNACL_COMPRESSED_MATRIX_HPP_
54 template <
typename CPU_MATRIX,
typename SCALARTYPE,
unsigned int ALIGNMENT>
55 void copy(
const CPU_MATRIX & cpu_matrix,
60 if ( cpu_matrix.size1() > 0 && cpu_matrix.size2() > 0 )
64 for (
typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
65 row_it != cpu_matrix.end1();
68 std::size_t entries_per_row = 0;
69 for (
typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
70 col_it != row_it.end();
75 num_entries += viennacl::tools::roundUpToNextMultiple<std::size_t>(entries_per_row, ALIGNMENT);
84 std::vector<cl_uint> row_buffer(cpu_matrix.size1() + 1);
85 std::vector<cl_uint> col_buffer(num_entries);
86 std::vector<SCALARTYPE> elements(num_entries);
88 std::size_t row_index = 0;
89 std::size_t data_index = 0;
91 for (
typename CPU_MATRIX::const_iterator1 row_it = cpu_matrix.begin1();
92 row_it != cpu_matrix.end1();
95 row_buffer[row_index] = data_index;
98 for (
typename CPU_MATRIX::const_iterator2 col_it = row_it.begin();
99 col_it != row_it.end();
102 col_buffer[data_index] =
static_cast<std::size_t
>(col_it.index2());
103 elements[data_index] = *col_it;
106 data_index = viennacl::tools::roundUpToNextMultiple<std::size_t>(data_index, ALIGNMENT);
108 row_buffer[row_index] = data_index;
110 gpu_matrix.
set(&row_buffer[0],
126 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
127 void copy(
const std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix,
133 #ifdef VIENNACL_HAVE_EIGEN
134 template <
typename SCALARTYPE,
int flags,
unsigned int ALIGNMENT>
135 void copy(
const Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix,
136 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix)
138 std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(eigen_matrix.rows());
140 for (
int k=0; k < eigen_matrix.outerSize(); ++k)
141 for (
typename Eigen::SparseMatrix<SCALARTYPE, flags>::InnerIterator it(eigen_matrix, k); it; ++it)
142 stl_matrix[it.row()][it.col()] = it.value();
144 copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix), gpu_matrix);
149 #ifdef VIENNACL_HAVE_MTL4
150 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
151 void copy(
const mtl::compressed2D<SCALARTYPE> & cpu_matrix,
152 compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix)
154 typedef mtl::compressed2D<SCALARTYPE> MatrixType;
156 std::vector< std::map<unsigned int, SCALARTYPE> > stl_matrix(cpu_matrix.num_rows());
158 using mtl::traits::range_generator;
159 using mtl::traits::range::min;
162 typedef typename min<range_generator<mtl::tag::row, MatrixType>,
163 range_generator<mtl::tag::col, MatrixType> >::type range_type;
167 typedef typename range_type::type c_type;
169 typedef typename mtl::traits::range_generator<mtl::tag::nz, c_type>::type ic_type;
172 typename mtl::traits::row<MatrixType>::type row(cpu_matrix);
173 typename mtl::traits::col<MatrixType>::type col(cpu_matrix);
174 typename mtl::traits::const_value<MatrixType>::type value(cpu_matrix);
177 for (c_type cursor(my_range.begin(cpu_matrix)), cend(my_range.end(cpu_matrix)); cursor != cend; ++cursor)
178 for (ic_type icursor(mtl::begin<mtl::tag::nz>(cursor)), icend(mtl::end<mtl::tag::nz>(cursor)); icursor != icend; ++icursor)
179 stl_matrix[row(*icursor)][col(*icursor)] = value(*icursor);
181 copy(tools::const_sparse_matrix_adapter<SCALARTYPE>(stl_matrix), gpu_matrix);
203 template <
typename CPU_MATRIX,
typename SCALARTYPE,
unsigned int ALIGNMENT>
205 CPU_MATRIX & cpu_matrix )
207 if ( gpu_matrix.
size1() > 0 && gpu_matrix.
size2() > 0 )
209 cpu_matrix.resize(gpu_matrix.
size1(), gpu_matrix.
size2(),
false);
212 std::vector<cl_uint> row_buffer(gpu_matrix.
size1() + 1);
213 std::vector<cl_uint> col_buffer(gpu_matrix.
nnz());
214 std::vector<SCALARTYPE> elements(gpu_matrix.
nnz());
228 std::size_t data_index = 0;
229 for (std::size_t row = 1; row <= gpu_matrix.
size1(); ++row)
231 while (data_index < row_buffer[row])
233 if (col_buffer[data_index] >= gpu_matrix.
size2())
235 std::cerr <<
"ViennaCL encountered invalid data at colbuffer[" << data_index <<
"]: " << col_buffer[data_index] << std::endl;
239 if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
240 cpu_matrix(row-1, col_buffer[data_index]) = elements[data_index];
253 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
255 std::vector< std::map<unsigned int, SCALARTYPE> > & cpu_matrix)
258 copy(gpu_matrix, temp);
262 #ifdef VIENNACL_HAVE_EIGEN
263 template <
typename SCALARTYPE,
int flags,
unsigned int ALIGNMENT>
264 void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix,
265 Eigen::SparseMatrix<SCALARTYPE, flags> & eigen_matrix)
267 if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
269 assert(static_cast<unsigned int>(eigen_matrix.rows()) >= gpu_matrix.size1()
270 &&
static_cast<unsigned int>(eigen_matrix.cols()) >= gpu_matrix.size2()
271 &&
"Provided Eigen compressed matrix is too small!");
274 std::vector<cl_uint> row_buffer(gpu_matrix.size1() + 1);
275 std::vector<cl_uint> col_buffer(gpu_matrix.nnz());
276 std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
280 CL_TRUE, 0,
sizeof(cl_uint)*(gpu_matrix.size1() + 1), &(row_buffer[0]), 0, NULL, NULL);
283 CL_TRUE, 0,
sizeof(cl_uint)*gpu_matrix.nnz(), &(col_buffer[0]), 0, NULL, NULL);
286 CL_TRUE, 0,
sizeof(SCALARTYPE)*gpu_matrix.nnz(), &(elements[0]), 0, NULL, NULL);
290 eigen_matrix.setZero();
291 eigen_matrix.startFill();
292 std::size_t data_index = 0;
293 for (std::size_t row = 1; row <= gpu_matrix.size1(); ++row)
295 while (data_index < row_buffer[row])
297 assert(col_buffer[data_index] < gpu_matrix.size2() &&
"ViennaCL encountered invalid data at col_buffer");
298 if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
299 eigen_matrix.fill(row-1, col_buffer[data_index]) = elements[data_index];
303 eigen_matrix.endFill();
310 #ifdef VIENNACL_HAVE_MTL4
311 template <
typename SCALARTYPE,
unsigned int ALIGNMENT>
312 void copy(compressed_matrix<SCALARTYPE, ALIGNMENT> & gpu_matrix,
313 mtl::compressed2D<SCALARTYPE> & mtl4_matrix)
315 if ( gpu_matrix.size1() > 0 && gpu_matrix.size2() > 0 )
317 assert(mtl4_matrix.num_rows() >= gpu_matrix.size1()
318 && mtl4_matrix.num_cols() >= gpu_matrix.size2()
319 &&
"Provided MTL4 compressed matrix is too small!");
322 std::vector<unsigned int> row_buffer(gpu_matrix.size1() + 1);
323 std::vector<unsigned int> col_buffer(gpu_matrix.nnz());
324 std::vector<SCALARTYPE> elements(gpu_matrix.nnz());
328 CL_TRUE, 0,
sizeof(cl_uint)*(gpu_matrix.size1() + 1), &(row_buffer[0]), 0, NULL, NULL);
331 CL_TRUE, 0,
sizeof(cl_uint)*gpu_matrix.nnz(), &(col_buffer[0]), 0, NULL, NULL);
334 CL_TRUE, 0,
sizeof(SCALARTYPE)*gpu_matrix.nnz(), &(elements[0]), 0, NULL, NULL);
341 mtl::matrix::inserter< mtl::compressed2D<SCALARTYPE> > ins(mtl4_matrix);
342 std::size_t data_index = 0;
343 for (std::size_t row = 1; row <= gpu_matrix.size1(); ++row)
345 while (data_index < row_buffer[row])
347 assert(col_buffer[data_index] < gpu_matrix.size2() &&
"ViennaCL encountered invalid data at col_buffer");
348 if (elements[data_index] != static_cast<SCALARTYPE>(0.0))
349 ins(row-1, col_buffer[data_index]) <<
typename mtl::Collection< mtl::compressed2D<SCALARTYPE> >::value_type(elements[data_index]);
367 template<
class SCALARTYPE,
unsigned int ALIGNMENT >
374 compressed_matrix() : _rows(0), _cols(0), _nonzeros(0) { viennacl::linalg::kernels::compressed_matrix<SCALARTYPE, ALIGNMENT>::init(); }
383 _rows(rows), _cols(cols), _nonzeros(nonzeros)
385 viennacl::linalg::kernels::compressed_matrix<SCALARTYPE, ALIGNMENT>::init();
397 std::size_t rows, std::size_t cols, std::size_t nonzeros) :
398 _rows(rows), _cols(cols), _nonzeros(nonzeros)
400 _row_buffer = mem_row_buffer;
402 _col_buffer = mem_col_buffer;
404 _elements = mem_elements;
418 void set(cl_uint * row_jumper,
419 cl_uint * col_buffer,
420 SCALARTYPE * elements,
423 std::size_t nonzeros)
426 assert(nonzeros > 0);
431 _nonzeros = nonzeros;
439 if (new_nonzeros > _nonzeros)
452 _nonzeros = new_nonzeros;
462 void resize(std::size_t new_size1, std::size_t new_size2,
bool preserve =
true)
464 assert(new_size1 > 0 && new_size2 > 0);
467 if (new_size1 != _rows || new_size2 != _cols)
469 std::vector<std::map<unsigned int, SCALARTYPE> > stl_sparse_matrix;
471 stl_sparse_matrix.resize(_rows);
473 if (preserve && _rows > 0)
476 stl_sparse_matrix.resize(new_size1);
479 if (new_size2 < _cols && _rows > 0)
481 for (
size_t i=0; i<stl_sparse_matrix.size(); ++i)
483 std::list<unsigned int> to_delete;
484 for (
typename std::map<unsigned int, SCALARTYPE>::iterator it = stl_sparse_matrix[i].begin();
485 it != stl_sparse_matrix[i].end();
488 if (it->first >= new_size2)
489 to_delete.push_back(it->first);
492 for (std::list<unsigned int>::iterator it = to_delete.begin(); it != to_delete.end(); ++it)
493 stl_sparse_matrix[i].erase(*it);
497 copy(stl_sparse_matrix, *
this);
505 const std::size_t &
size1()
const {
return _rows; }
507 const std::size_t &
size2()
const {
return _cols; }
509 const std::size_t &
nnz()
const {
return _nonzeros; }
528 std::size_t _nonzeros;