Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
StorageSmvBsr.hpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: StorageSmvBsr.hpp
4 //
5 // For more information, please see: http://www.nektar.info
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: 0-based sparse BSR storage class with own unrolled and
33 // LibSMV multiply kernels.
34 //
35 ///////////////////////////////////////////////////////////////////////////////
36 
37 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STORAGE_SMV_BSR_HPP
38 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STORAGE_SMV_BSR_HPP
39 
40 #include <map>
41 #include <vector>
42 #include <utility>
43 #include <fstream>
44 
47 
48 #include <boost/call_traits.hpp>
49 
50 
51 namespace Nektar
52 {
53  /*
54  * Zero-based BSR (Block Sparse Row) storage class with its sparse
55  * multiply kernels built upon dense multiply kernels provided by
56  * LibSMV library. It also includes its own unrolled multiply kernels
57  * up to 4x4 matrices, provided 'just in case' LibSMV library
58  * is not available. When matrix is larger than either one
59  * supported by LibSMV library or 4x4 (whichever is greater), the
60  * multiply kernel calls dense dgemv from BLAS.
61  *
62  * The BSR sparse format assumes sparse matrix is a CSR collection of
63  * dense square blocks of same size. In contrast with Nist BSR class
64  * this one uses zero-based storage. The constructor takes input matrix in
65  * block coordinate storage (BCO) sparse format.
66  *
67  */
68 
69  template<typename T>
70  class StorageSmvBsr
71  {
72 
73  public:
74  typedef T DataType;
78 
79  typedef void (*MultiplyKernel)(const double*, const double*, double*);
80 
81  /// \internal
82  /// \brief Forward iterator through nonzero (double) elements of the matrix
83  /// that mimics forward iteration of BCOMatType.
84  /// It's a dirty hack, not a real iterator in the C++ sense.
86  {
87  struct IterType
88  {
89  CoordType first; //< (row, column)
90  DataType second; //< value
91  IndexType nnzindex; //< index of this nnz entry
92  IndexType storageindex;//< offset of this nnz entry in the storage
93  };
94 
95  public:
98  IndexType end,
99  IndexType blkDim,
100  const DataVectorType& val,
101  const IndexVectorType& indx,
102  const IndexVectorType& pntr);
103  const_iterator(const const_iterator& src);
104  ~const_iterator();
105 
108  const IterType& operator*();
109  const IterType* operator->();
110  const bool operator==(const const_iterator& rhs);
111  const bool operator!=(const const_iterator& rhs);
112 
113  // one way conversion: iterator -> const_iterator
114  // operator const_iterator<T const, Tag>() const;
115 
116  private:
117  void forward();
119 
125  const DataVectorType& m_val;
126  const IndexVectorType& m_indx;
127  const IndexVectorType& m_pntr;
128  };
129 
130 
131  public:
132  // Constructs zero-based BSR sparse matrix based on input BCO storage
134  const IndexType blkCols,
135  const IndexType blkDim,
136  const BCOMatType& bcoMat,
137  const MatrixStorage matType = eFULL);
138 
139  // Copy constructor
141 
143 
144  LIB_UTILITIES_EXPORT const IndexType GetRows() const;
149  LIB_UTILITIES_EXPORT const DataType GetFillInRatio() const;
150  LIB_UTILITIES_EXPORT const size_t GetMemoryUsage(IndexType nnz, IndexType nRows) const;
151 
154 
155  LIB_UTILITIES_EXPORT const typename boost::call_traits<DataType>::const_reference
156  GetValue(IndexType row, IndexType column) const;
157 
158  LIB_UTILITIES_EXPORT void Multiply(const DataType* in,
159  DataType* out);
160  LIB_UTILITIES_EXPORT void Multiply(const DataVectorType &in,
161  DataVectorType &out);
162  LIB_UTILITIES_EXPORT void MultiplyLight(const DataVectorType &in,
163  DataVectorType &out);
164 
165 
166  protected:
167 
168  // converts input vector of BCO blocks
169  // to the internal zero-based BSR representation
170  void processBcoInput(
171  const IndexType blkRows,
172  const IndexType blkColumns,
173  const IndexType blkDim,
174  const BCOMatType& bcoMat);
175 
176 
177  void Multiply_1x1(const int mb, const int kb, const double* val,
178  const int* bindx, const int* bpntrb, const int* bpntre,
179  const double* b, double* c);
180 
181  void Multiply_2x2(const int mb, const int kb, const double* val,
182  const int* bindx, const int* bpntrb, const int* bpntre,
183  const double* b, double* c);
184 
185  void Multiply_3x3(const int mb, const int kb, const double* val,
186  const int* bindx, const int* bpntrb, const int* bpntre,
187  const double* b, double* c);
188 
189  void Multiply_4x4(const int mb, const int kb, const double* val,
190  const int* bindx, const int* bpntrb, const int* bpntre,
191  const double* b, double* c);
192 
193  void Multiply_generic(const int mb, const int kb, const double* val,
194  const int* bindx, const int* bpntrb, const int* bpntre,
195  const double* b, double* c);
196 
197 #ifdef NEKTAR_USING_SMV
198  void Multiply_libsmv(const int mb, const int kb, const double* val,
199  const int* bindx, const int* bpntrb, const int* bpntre,
200  const double* b, double* c);
201 #endif
202 
203 
204  // interface to lowest level LibSMV multiply kernels
206 
208 
209  IndexType m_blkRows; // number of block rows
210  IndexType m_blkCols; // number of block columns
211  IndexType m_blkDim; // rank of each block
212 
213  IndexType m_bnnz; //< number of nonzero blocks
214  IndexType m_nnz; //< number of factual nonzero entries in the sparse matrix
215 
216  DataVectorType m_val; // values of non-zero entries
217  IndexVectorType m_indx; // column indices of non-zero entries
218  IndexVectorType m_pntr; // m_pntr(i) contains index in m_val of first non-zero element in row i
219 
220  private:
221 
222  };
223 
224 
225 
226 } // namespace
227 
228 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_STORAGE_SMV_BSR_HPP
Array< OneD, DataType > DataVectorType
const size_t GetMemoryUsage(IndexType nnz, IndexType nRows) const
Array< OneD, const DataType > ConstDataVectorType
const IndexType GetColumns() const
IndexVectorType m_indx
StorageSmvBsr(const IndexType blkRows, const IndexType blkCols, const IndexType blkDim, const BCOMatType &bcoMat, const MatrixStorage matType=eFULL)
const_iterator begin() const
void Multiply_1x1(const int mb, const int kb, const double *val, const int *bindx, const int *bpntrb, const int *bpntre, const double *b, double *c)
Zero-based CSR multiply. Essentially this is slightly modified copy-paste from NIST Sparse Blas 0...
Array< OneD, IndexType > IndexVectorType
const IndexType GetNumNonZeroEntries() const
void Multiply(const DataType *in, DataType *out)
const IndexType GetRows() const
void(* MultiplyKernel)(const double *, const double *, double *)
void Multiply_4x4(const int mb, const int kb, const double *val, const int *bindx, const int *bpntrb, const int *bpntre, const double *b, double *c)
Zero-based BSR multiply unrolled for 4x4 blocks.
std::map< CoordType, BCOEntryType > BCOMatType
#define LIB_UTILITIES_EXPORT
const boost::call_traits< DataType >::const_reference GetValue(IndexType row, IndexType column) const
void Multiply_2x2(const int mb, const int kb, const double *val, const int *bindx, const int *bpntrb, const int *bpntre, const double *b, double *c)
Zero-based BSR multiply unrolled for 2x2 blocks. Essentially this is slightly optimised copy-paste fr...
const IndexType GetNumStoredDoubles() const
const IndexType GetBlkSize() const
unsigned int IndexType
void processBcoInput(const IndexType blkRows, const IndexType blkColumns, const IndexType blkDim, const BCOMatType &bcoMat)
const bool operator!=(const const_iterator &rhs)
IndexVectorType m_pntr
void Multiply_generic(const int mb, const int kb, const double *val, const int *bindx, const int *bpntrb, const int *bpntre, const double *b, double *c)
Generic zero-based BSR multiply for higher matrix ranks.
CoordType storageIndexToFullCoord(IndexType storageIndex)
const_iterator(MatrixStorage matType, IndexType begin, IndexType end, IndexType blkDim, const DataVectorType &val, const IndexVectorType &indx, const IndexVectorType &pntr)
const bool operator==(const const_iterator &rhs)
const DataType GetFillInRatio() const
const_iterator end() const
MultiplyKernel m_mvKernel
void Multiply_3x3(const int mb, const int kb, const double *val, const int *bindx, const int *bpntrb, const int *bpntre, const double *b, double *c)
Zero-based BSR multiply unrolled for 3x3 blocks.
std::pair< IndexType, IndexType > CoordType
1D Array of constant elements with garbage collection and bounds checking.
Definition: SharedArray.hpp:59
void MultiplyLight(const DataVectorType &in, DataVectorType &out)