28 using namespace Eigen;
35 TripletEntry(
int r_,
int c_,
double x_) : r(r_), c(c_), x(x_) {}
39 bool operator()(
const TripletEntry& e1,
const TripletEntry& e2)
const 41 return e1.c < e2.c || (e1.c == e2.c && e1.r < e2.r);
45 template<
class T1,
class T2,
class Pred = std::less<T1> >
47 bool operator()(
const std::pair<T1,T2>& left,
const std::pair<T1,T2>& right) {
48 return Pred()(left.first, right.first);
53 template <
class MatrixType>
55 _rowBlockIndices(rbi,rbi+rb),
56 _colBlockIndices(cbi,cbi+cb),
57 _blockCols(cb), _hasStorage(hasStorage)
61 template <
class MatrixType>
63 _blockCols(0), _hasStorage(true)
67 template <
class MatrixType>
70 # pragma omp parallel for default (shared) if (_blockCols.size() > 100) 72 for (
int i=0; i < static_cast<int>(_blockCols.size()); ++i) {
75 if (_hasStorage && dealloc)
80 if (_hasStorage && dealloc)
81 _blockCols[i].
clear();
85 template <
class MatrixType>
91 template <
class MatrixType>
95 if (it==_blockCols[c].end()){
96 if (!_hasStorage && ! alloc )
99 int rb=rowsOfBlock(r);
100 int cb=colsOfBlock(c);
103 std::pair < typename SparseBlockMatrix<MatrixType>::IntBlockMap::iterator,
bool> result
104 =_blockCols[c].insert(std::make_pair(r,_block)); (void) result;
105 assert (result.second);
113 template <
class MatrixType>
116 if (it==_blockCols[c].end())
122 template <
class MatrixType>
125 for (
size_t i=0; i<_blockCols.size(); ++i){
128 ret->
_blockCols[i].insert(std::make_pair(it->first, b));
136 template <
class MatrixType>
137 template <
class MatrixTransposedType>
148 for (
size_t i=0; i<_rowBlockIndices.size(); ++i){
152 for (
size_t i=0; i<_colBlockIndices.size(); ++i){
158 for (
size_t i=0; i<_blockCols.size(); ++i){
168 template <
class MatrixType>
171 dest=
new SparseBlockMatrix(&_rowBlockIndices[0], &_colBlockIndices[0], _rowBlockIndices.size(), _colBlockIndices.size());
179 for (
size_t i=0; i<_rowBlockIndices.size(); ++i){
183 for (
size_t i=0; i<_colBlockIndices.size(); ++i){
188 for (
size_t i=0; i<_blockCols.size(); ++i){
198 template <
class MatrixType>
199 template <
class MatrixResultType,
class MatrixFactorType >
204 for (
size_t i=0; i<_colBlockIndices.size(); ++i){
213 for (
size_t i=0; i<M->
_blockCols.size(); ++i){
219 while(rbt!=_blockCols[it->first].end()){
234 template <
class MatrixType>
237 dest=
new double [_rowBlockIndices[_rowBlockIndices.size()-1] ];
238 memset(dest,0, _rowBlockIndices[_rowBlockIndices.size()-1]*
sizeof(double));
242 Eigen::Map<VectorXd> destVec(dest, rows());
243 const Eigen::Map<const VectorXd> srcVec(src, cols());
245 for (
size_t i=0; i<_blockCols.size(); ++i){
246 int srcOffset = i ? _colBlockIndices[i-1] : 0;
250 int destOffset = it->first ? _rowBlockIndices[it->first - 1] : 0;
257 template <
class MatrixType>
261 dest=
new double [_rowBlockIndices[_rowBlockIndices.size()-1] ];
262 memset(dest,0, _rowBlockIndices[_rowBlockIndices.size()-1]*
sizeof(double));
266 Eigen::Map<VectorXd> destVec(dest, rows());
267 const Eigen::Map<const VectorXd> srcVec(src, cols());
269 for (
size_t i=0; i<_blockCols.size(); ++i){
270 int srcOffset = colBaseOfBlock(i);
273 int destOffset = rowBaseOfBlock(it->first);
274 if (destOffset > srcOffset)
278 if (destOffset < srcOffset)
284 template <
class MatrixType>
289 dest=
new double [ destSize ];
290 memset(dest,0, destSize*
sizeof(
double));
294 Eigen::Map<VectorXd> destVec(dest, destSize);
295 Eigen::Map<const VectorXd> srcVec(src, rows());
298 # pragma omp parallel for default (shared) schedule(dynamic, 10) 300 for (
int i=0; i < static_cast<int>(_blockCols.size()); ++i){
301 int destOffset = colBaseOfBlock(i);
303 it!=_blockCols[i].end();
306 int srcOffset = rowBaseOfBlock(it->first);
314 template <
class MatrixType>
316 for (
size_t i=0; i<_blockCols.size(); ++i){
324 template <
class MatrixType>
329 rowIdx[0] = rowsOfBlock(rmin);
330 for (
int i=1; i<m; ++i){
331 rowIdx[i]=rowIdx[i-1]+rowsOfBlock(rmin+i);
335 colIdx[0] = colsOfBlock(cmin);
336 for (
int i=1; i<n; ++i){
337 colIdx[i]=colIdx[i-1]+colsOfBlock(cmin+i);
340 for (
int i=0; i<n; ++i){
343 if (it->first >= rmin && it->first < rmax){
345 s->
_blockCols[i].insert(std::make_pair(it->first-rmin, b));
353 template <
class MatrixType>
356 for (
size_t i=0; i<_blockCols.size(); ++i)
357 count+=_blockCols[i].size();
361 template <
class MatrixType>
363 if (MatrixType::SizeAtCompileTime != Eigen::Dynamic) {
364 size_t nnz = nonZeroBlocks() * MatrixType::SizeAtCompileTime;
368 for (
size_t i=0; i<_blockCols.size(); ++i){
378 template <
class MatrixType>
379 std::ostream& operator << (std::ostream& os, const SparseBlockMatrix<MatrixType>& m){
380 os <<
"RBI: " << m.rowBlockIndices().size();
381 for (
size_t i=0; i<m.rowBlockIndices().size(); ++i)
382 os <<
" " << m.rowBlockIndices()[i];
384 os <<
"CBI: " << m.colBlockIndices().size();
385 for (
size_t i=0; i<m.colBlockIndices().size(); ++i)
386 os <<
" " << m.colBlockIndices()[i];
389 for (
size_t i=0; i<m.blockCols().size(); ++i){
392 os <<
"BLOCK: " << it->first <<
" " << i << std::endl;
393 os << *b << std::endl;
399 template <
class MatrixType>
402 size_t n=_rowBlockIndices.size();
404 int blockSizes[_rowBlockIndices.size()];
405 blockSizes[0]=_rowBlockIndices[0];
406 for (
size_t i=1; i<n; ++i){
407 blockSizes[i]=_rowBlockIndices[i]-_rowBlockIndices[i-1];
410 int pBlockIndices[_rowBlockIndices.size()];
411 for (
size_t i=0; i<n; ++i){
412 pBlockIndices[pinv[i]]=blockSizes[i];
414 for (
size_t i=1; i<n; ++i){
415 pBlockIndices[i]+=pBlockIndices[i-1];
425 for (
size_t i=0; i<n; ++i){
434 for (
size_t i=0; i<n; ++i){
438 it!=_blockCols[i].end(); ++it){
439 int pj=pinv[it->first];
444 if (! upperTriangle || pj<=pi) {
445 b=dest->
block(pj,pi,
true);
450 b=dest->
block(pi,pj,
true);
464 template <
class MatrixType>
467 assert(Cx &&
"Target destination is NULL");
468 double* CxStart = Cx;
469 for (
size_t i=0; i<_blockCols.size(); ++i){
470 int cstart=i ? _colBlockIndices[i-1] : 0;
471 int csize=colsOfBlock(i);
472 for (
int c=0; c<csize; ++c) {
475 int rstart=it->first ? _rowBlockIndices[it->first-1] : 0;
477 int elemsToCopy = b->
rows();
478 if (upperTriangle && rstart == cstart)
480 memcpy(Cx, b->data() + c*b->
rows(), elemsToCopy *
sizeof(double));
489 template <
class MatrixType>
492 assert(Cp && Ci && Cx &&
"Target destination is NULL");
494 for (
size_t i=0; i<_blockCols.size(); ++i){
495 int cstart=i ? _colBlockIndices[i-1] : 0;
496 int csize=colsOfBlock(i);
497 for (
int c=0; c<csize; ++c) {
501 int rstart=it->first ? _rowBlockIndices[it->first-1] : 0;
503 int elemsToCopy = b->
rows();
504 if (upperTriangle && rstart == cstart)
506 for (
int r=0; r<elemsToCopy; ++r){
519 template <
class MatrixType>
522 int n = _colBlockIndices.size();
523 int nzMax = (int)nonZeroBlocks();
526 ms.
m = _rowBlockIndices.size();
531 for (
int i = 0; i < static_cast<int>(_blockCols.size()); ++i){
535 const int& r = it->first;
547 template <
class MatrixType>
550 std::string name = filename;
551 std::string::size_type lastDot = name.find_last_of(
'.');
552 if (lastDot != std::string::npos)
553 name = name.substr(0, lastDot);
555 std::vector<TripletEntry> entries;
556 for (
size_t i = 0; i<_blockCols.size(); ++i){
559 const int& r = it->first;
560 const MatrixType& m = *(it->second);
561 for (
int cc = 0; cc < m.cols(); ++cc)
562 for (
int rr = 0; rr < m.rows(); ++rr) {
563 int aux_r = rowBaseOfBlock(r) + rr;
564 int aux_c = colBaseOfBlock(c) + cc;
565 entries.push_back(TripletEntry(aux_r, aux_c, m(rr, cc)));
566 if (upperTriangle && r != c) {
567 entries.push_back(TripletEntry(aux_c, aux_r, m(rr, cc)));
573 int nz = entries.size();
574 std::sort(entries.begin(), entries.end(), TripletColSort());
576 std::ofstream fout(filename);
577 fout <<
"# name: " << name << std::endl;
578 fout <<
"# type: sparse matrix" << std::endl;
579 fout <<
"# nnz: " << nz << std::endl;
580 fout <<
"# rows: " << rows() << std::endl;
581 fout <<
"# columns: " << cols() << std::endl;
582 fout << std::setprecision(9) << std::fixed << std::endl;
584 for (std::vector<TripletEntry>::const_iterator it = entries.begin(); it != entries.end(); ++it) {
585 const TripletEntry& entry = *it;
586 fout << entry.r+1 <<
" " << entry.c+1 <<
" " << entry.x << std::endl;
591 template <
class MatrixType>
594 blockCCS.
blockCols().resize(blockCols().size());
596 for (
size_t i = 0; i < blockCols().size(); ++i) {
600 dest.reserve(row.size());
601 for (
typename IntBlockMap::const_iterator it = row.begin(); it != row.end(); ++it) {
609 template <
class MatrixType>
613 blockCCS.
blockCols().resize(_rowBlockIndices.size());
615 for (
size_t i = 0; i < blockCols().size(); ++i) {
617 for (
typename IntBlockMap::const_iterator it = row.begin(); it != row.end(); ++it) {
626 template <
class MatrixType>
631 typedef std::pair<int, MatrixType*> SparseColumnPair;
633 for (
size_t i = 0; i < hashMatrix.
blockCols().size(); ++i) {
635 HashSparseColumn& column = hashMatrix.
blockCols()[i];
636 if (column.size() == 0)
638 std::vector<SparseColumnPair> sparseRowSorted;
639 sparseRowSorted.reserve(column.size());
640 for (
typename HashSparseColumn::const_iterator it = column.begin(); it != column.end(); ++it)
641 sparseRowSorted.push_back(*it);
642 std::sort(sparseRowSorted.begin(), sparseRowSorted.end(), CmpPairFirst<int, MatrixType*>());
644 HashSparseColumn aux;
648 destColumnMap.insert(sparseRowSorted[0]);
649 for (
size_t j = 1; j < sparseRowSorted.size(); ++j) {
652 destColumnMap.insert(hint, sparseRowSorted[j]);
std::vector< int > _rowBlockIndices
vector of the indices of the blocks along the rows.
std::vector< int > _colBlockIndices
int cols() const
columns of the matrix
int * Ap
column pointers for A, of size n+1
void axpy(const MatrixType &A, const Eigen::Map< const Eigen::VectorXd > &x, int xoff, Eigen::Map< Eigen::VectorXd > &y, int yoff)
SparseMatrixBlock * block(int r, int c, bool alloc=false)
returns the block at location r,c. if alloc=true he block is created if it does not exist ...
representing the structure of a matrix in column compressed structure (only the upper triangular part...
void atxpy(const MatrixType &A, const Eigen::Map< const Eigen::VectorXd > &x, int xoff, Eigen::Map< Eigen::VectorXd > &y, int yoff)
void clear(bool dealloc=false)
this zeroes all the blocks. If dealloc=true the blocks are removed from memory
MatrixType SparseMatrixBlock
this is the type of the elementary block, it is an Eigen::Matrix.
std::vector< RowBlock > SparseColumn
const std::vector< SparseColumn > & blockCols() const
the block matrices per block-column
std::vector< IntBlockMap > _blockCols
int m
A is m-by-n. m must be >= 0.
const std::vector< SparseColumn > & blockCols() const
the block matrices per block-column
Sparse matrix which uses blocks.
std::map< int, SparseMatrixBlock *> IntBlockMap
void alloc(int n_, int nz)
std::tr1::unordered_map< int, MatrixType * > SparseColumn
const std::vector< IntBlockMap > & blockCols() const
the block matrices per block-column
int * Aii
row indices of A, of size nz = Ap [n]
bool transpose(SparseBlockMatrix< MatrixTransposedType > *&dest) const
transposes a block matrix, The transposed type should match the argument false on failure ...
Sparse matrix which uses blocks based on hash structures.
void swap(scoped_ptr< T > &, scoped_ptr< T > &)
Sparse matrix which uses blocks.
int rows() const
rows of the matrix