1 // 2 // BAGEL - Brilliantly Advanced General Electronic Structure Library 3 // Filename: blocksparsematrix.h 4 // Copyright (C) 2014 Toru Shiozaki 5 // 6 // Author: Shane Parker <shane.parker@u.northwestern.edu> 7 // Maintainer: Shiozaki Group 8 // 9 // This file is part of the BAGEL package. 10 // 11 // This program is free software: you can redistribute it and/or modify 12 // it under the terms of the GNU General Public License as published by 13 // the Free Software Foundation, either version 3 of the License, or 14 // (at your option) any later version. 15 // 16 // This program is distributed in the hope that it will be useful, 17 // but WITHOUT ANY WARRANTY; without even the implied warranty of 18 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 19 // GNU General Public License for more details. 20 // 21 // You should have received a copy of the GNU General Public License 22 // along with this program. If not, see <http://www.gnu.org/licenses/>. 23 // 24 25 #ifndef __math_blocksparsematrix_h 26 #define __math_blocksparsematrix_h 27 28 #include <src/util/math/matrix.h> 29 30 namespace bagel { 31 32 /// Class for blockwise sparse matrices 33 class BlockSparseMatrix { 34 protected: 35 /// Store blocks as separate Matrix objects with offsets 36 std::map<std::pair<size_t, size_t>, std::shared_ptr<Matrix>> data_; 37 38 int ndim_; 39 int mdim_; 40 41 public: 42 /// A single block (so actually a dense matrix) 43 BlockSparseMatrix(std::shared_ptr<Matrix> m); 44 /// Pre-determined blocking structure 45 BlockSparseMatrix(const int n, const int m, std::map<std::pair<size_t, size_t>, std::shared_ptr<Matrix>> d); 46 // for completeness, a constructor that accepts a matrix and a threshold should probably exist at some point 47 data()48 std::map<std::pair<size_t, size_t>, std::shared_ptr<Matrix>>& data() { return data_; } data()49 const std::map<std::pair<size_t, size_t>, std::shared_ptr<Matrix>>& data() const { return data_; } 50 51 /// get block based on starting location. this is the best way to alter the data get_block(size_t nstart,size_t mstart)52 std::shared_ptr<Matrix> get_block(size_t nstart, size_t mstart) { return (data_.find({nstart, mstart})!=data_.end() ? data_[{nstart,mstart}] : nullptr); } 53 /// const get block based on starting location. get_block(size_t nstart,size_t mstart)54 std::shared_ptr<const Matrix> get_block(size_t nstart, size_t mstart) const { return (data_.find({nstart, mstart})!=data_.end() ? data_.at({nstart,mstart}) : nullptr); } 55 ndim()56 int ndim() const { return ndim_; } mdim()57 int mdim() const { return mdim_; } size()58 int size() const { return accumulate(data_.begin(), data_.end(), 0, [] (int x, std::pair<std::pair<size_t, size_t>, std::shared_ptr<Matrix>> p) { return x + p.second->size(); }); } 59 60 /// element-wise read-only: is slow and not recommended. element(const int n,const int m)61 double element(const int n, const int m) const { 62 auto iter = std::find_if(data_.begin(), data_.end(), [&n, &m] (std::pair<std::pair<size_t, size_t>, std::shared_ptr<Matrix>> p) 63 { return (n >= p.first.first && n < p.first.first + p.second->ndim()) && (m >= p.first.second && m < p.first.second + p.second->mdim()); }); 64 if (iter!=data_.end()) 65 return iter->second->element(n - iter->first.first, m - iter->first.second); 66 else 67 return 0.0; 68 } 69 70 /// returns a VectorB containing the diagonal elements 71 std::shared_ptr<VectorB> diagonal() const; 72 }; 73 74 void mat_block_multiply(bool Atrans, bool Btrans, double alpha, const Matrix& A, const BlockSparseMatrix& B, double beta, Matrix& C); 75 76 } 77 78 #endif 79