Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GlobalMatrix.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: GlobalMatrix.cpp
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: GlobalMatrix definition
33 //
34 ///////////////////////////////////////////////////////////////////////////////
35 
40 
41 #include <iomanip>
42 #include <fstream>
43 
44 
45 namespace Nektar
46 {
47  namespace MultiRegions
48  {
50  RegisterDefaultSolverInfo("GlobalMatrixStorageType","SmvBSR");
51  std::string GlobalMatrix::lookupIds[3] = {
53  "GlobalMatrixStorageType", "SmvBSR", MultiRegions::eSmvBSR)
54  };
55 
56 
57  /**
58  * @class GlobalMatrix
59  * This matrix is essentially a wrapper around a DNekSparseMat.
60  */
61 
62  /**
63  * Allocates a new DNekSparseMat object from the given specification.
64  * @param rows Number of rows in matrix.
65  * @param columns Number of columns in matrix.
66  * @param cooMat ?
67  */
70  unsigned int rows,
71  unsigned int columns,
72  const COOMatType &cooMat,
73  const MatrixStorage& matStorage):
74  m_smvbsrmatrix(),
75  m_rows(rows),
76  m_mulCallsCounter(0)
77  {
78  MatrixStorageType storageType = pSession->
79  GetSolverInfoAsEnum<MatrixStorageType>("GlobalMatrixStorageType");
80 
81  unsigned int brows, bcols;
82 
83  // Size of dense matrix sub-blocks
84  int block_size = 1;
85 
86  BCOMatType bcoMat;
87 
88  // assuming current sparse format allows
89  // block-sparse data representation
90 
91  if(pSession->DefinesParameter("SparseBlockSize"))
92  {
93  pSession->LoadParameter("SparseBlockSize", block_size);
94  ASSERTL1(block_size > 0,"SparseBlockSize parameter must to be positive");
95  }
96 
97  brows = rows / block_size + (rows % block_size > 0);
98  bcols = columns / block_size + (columns % block_size > 0);
99 
100  if (rows % block_size > 0) m_copyOp = true;
101 
102  if (m_copyOp)
103  {
104  m_tmpin = Array<OneD, NekDouble> (brows*block_size, 0.0);
105  m_tmpout = Array<OneD, NekDouble> (brows*block_size, 0.0);
106  }
107 
108  convertCooToBco(brows, bcols, block_size, cooMat, bcoMat);
109 
110  size_t matBytes;
111  switch(storageType)
112  {
113  case eSmvBSR:
114  {
115 
116  // Create zero-based Smv-multiply BSR sparse storage holder
120  brows, bcols, block_size, bcoMat, matStorage );
121 
122  // Create sparse matrix
124  AllocateSharedPtr( sparseStorage );
125 
126  matBytes = m_smvbsrmatrix->GetMemoryFootprint();
127 
128  }
129  break;
130 
131  default:
132  NEKERROR(ErrorUtil::efatal,"Unsupported sparse storage type chosen");
133  }
134 
135  cout << "Global matrix storage type: "
136  << MatrixStorageTypeMap[storageType] << endl;
137  std::cout << "Global matrix memory, bytes = " << matBytes;
138  if (matBytes/(1024*1024) > 0)
139  {
140  std::cout << " ("<< matBytes/(1024*1024) <<" MB)" << std::endl;
141  }
142  else
143  {
144  std::cout << " ("<< matBytes/1024 <<" KB)" << std::endl;
145  }
146  std::cout << "Sparse storage block size = " << block_size << std::endl;
147  }
148 
149  /**
150  * Performs a matrix-vector multiply using the sparse format-specific
151  * multiply routine.
152  * @param in Input vector.
153  * @param out Output vector.
154  */
155  void GlobalMatrix::Multiply(const Array<OneD,const NekDouble> &in,
156  Array<OneD, NekDouble> &out)
157  {
158  if (!m_copyOp)
159  {
160  if (m_smvbsrmatrix) m_smvbsrmatrix->Multiply(in,out);
161  }
162  else
163  {
164  // if block size makes the last row/column bigger, one needs
165  // using temporary storage for rhs and result vectors.
166  Vmath::Vcopy(m_rows, &in[0], 1, &m_tmpin[0], 1);
167 
169 
170  Vmath::Vcopy(m_rows, &m_tmpout[0], 1, &out[0], 1);
171  }
172 
174  }
175 
176  const unsigned long GlobalMatrix::GetMulCallsCounter() const
177  {
178  if (m_smvbsrmatrix) return m_smvbsrmatrix->GetMulCallsCounter();
179  return -1;
180  }
181 
182  const unsigned int GlobalMatrix::GetNumNonZeroEntries() const
183  {
184  if (m_smvbsrmatrix) return m_smvbsrmatrix->GetNumNonZeroEntries();
185  return -1;
186  }
187 
188 
189  } //end of namespace
190 } //end of namespace
191