Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MatrixFuncs.h
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MatrixFuncs.h
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: Matrix functions that depend on storage policy. Putting
33 // methods in these separate classes makes it easier to use them from
34 // normal, scaled, and block matrices.
35 //
36 ///////////////////////////////////////////////////////////////////////////////
37 
38 #ifndef NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
39 #define NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
40 
41 #include <boost/lexical_cast.hpp>
42 #include <boost/tuple/tuple.hpp>
43 #include <limits>
44 #include <LibUtilities/LinearAlgebra/Lapack.hpp>
45 #include <LibUtilities/BasicUtils/SharedArray.hpp>
46 #include <LibUtilities/LibUtilitiesDeclspec.h>
47 
48 namespace Nektar
49 {
50  struct LIB_UTILITIES_EXPORT BandedMatrixFuncs
51  {
52  /// \brief Calculates and returns the storage size required.
53  ///
54  /// This method assumes that the matrix will be used with LU factorizationa and
55  /// allocates additional storage as appropriate.
56  static unsigned int GetRequiredStorageSize(unsigned int totalRows, unsigned int totalColumns,
57  unsigned int subDiags, unsigned int superDiags);
58 
59  static unsigned int CalculateNumberOfDiags(unsigned int totalRows, unsigned int diags);
60 
61  static unsigned int CalculateNumberOfRows(unsigned int totalRows, unsigned int subDiags, unsigned int superDiags);
62 
63  static unsigned int CalculateIndex(unsigned int totalRows,
64  unsigned int totalColumns,
65  unsigned int row, unsigned int column,
66  unsigned int sub, unsigned int super);
67 
68 
69  static boost::tuples::tuple<unsigned int, unsigned int>
70  Advance(const unsigned int totalRows, const unsigned int totalColumns,
71  const unsigned int curRow, const unsigned int curColumn);
72  };
73 
74  struct LIB_UTILITIES_EXPORT FullMatrixFuncs
75  {
76 
77  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns);
78  static unsigned int CalculateIndex(unsigned int totalRows, unsigned int totalColumns, unsigned int curRow, unsigned int curColumn);
79 
80 
81  static boost::tuples::tuple<unsigned int, unsigned int>
82  Advance(const unsigned int totalRows, const unsigned int totalColumns,
83  const unsigned int curRow, const unsigned int curColumn);
84 
85  template<typename DataType>
86  static void Invert(unsigned int rows, unsigned int columns,
87  Array<OneD, DataType>& data,
88  const char transpose)
89  {
90 #ifdef NEKTAR_USING_BLAS
91  ASSERTL0(rows==columns, "Only square matrices can be inverted.");
92  ASSERTL0(transpose=='N', "Only untransposed matrices may be inverted.");
93 
94  int m = rows;
95  int n = columns;
96  int info = 0;
97  Array<OneD, int> ipivot(n);
98  Array<OneD, DataType> work(n);
99 
100  Lapack::Dgetrf(m, n, data.get(), m, ipivot.get(), info);
101 
102  if( info < 0 )
103  {
104  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgetrf";
105  ASSERTL0(false, message.c_str());
106  }
107  else if( info > 0 )
108  {
109  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dgetrf";
110  ASSERTL0(false, message.c_str());
111  }
112 
113  Lapack::Dgetri(n, data.get(), n, ipivot.get(),
114  work.get(), n, info);
115 
116  if( info < 0 )
117  {
118  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dgetri";
119  ASSERTL0(false, message.c_str());
120  }
121  else if( info > 0 )
122  {
123  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dgetri";
124  ASSERTL0(false, message.c_str());
125  }
126 
127  #else
128  // error Full matrix inversion not supported without blas.
129  BOOST_STATIC_ASSERT(sizeof(DataType) == 0);
130  #endif
131  }
132 
133  static void EigenSolve(unsigned int n,
134  const Array<OneD, const double>& A,
135  Array<OneD, NekDouble> &EigValReal,
136  Array<OneD, NekDouble> &EigValImag,
137  Array<OneD, NekDouble> &EigVecs = NullNekDouble1DArray)
138  {
139  int lda = n,info = 0;
140  NekDouble dum;
141  int lwork = 3*lda;
142  Array<OneD,NekDouble> work(3*lda);
143  char uplo = 'N';
144 
145  if(EigVecs != NullNekDouble1DArray) // calculate Right Eigen Vectors
146  {
147  char lrev = 'V';
148  Lapack::Dgeev(uplo,lrev,lda, A.get(),lda,
149  EigValReal.get(),
150  EigValImag.get(),
151  &dum,1,
152  EigVecs.get(),lda,
153  &work[0],lwork,info);
154  }
155  else
156  {
157  char lrev = 'N';
158  Lapack::Dgeev(uplo,lrev,lda,
159  A.get(),lda,
160  EigValReal.get(),
161  EigValImag.get(),
162  &dum,1,&dum,1,
163  &work[0],lwork,info);
164  }
165  ASSERTL0(info == 0,"Info is not zero");
166 
167  }
168  };
169 
170  struct LIB_UTILITIES_EXPORT TriangularMatrixFuncs
171  {
172  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns);
173  };
174 
175  struct LIB_UTILITIES_EXPORT UpperTriangularMatrixFuncs : public TriangularMatrixFuncs
176  {
177  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn);
178 
179  static boost::tuples::tuple<unsigned int, unsigned int>
180  Advance(const unsigned int totalRows, const unsigned int totalColumns,
181  const unsigned int curRow, const unsigned int curColumn);
182  };
183 
184 
185  struct LIB_UTILITIES_EXPORT LowerTriangularMatrixFuncs : public TriangularMatrixFuncs
186  {
187  static unsigned int CalculateIndex(unsigned int totalColumns, unsigned int curRow, unsigned int curColumn);
188 
189  static boost::tuples::tuple<unsigned int, unsigned int>
190  Advance(const unsigned int totalRows, const unsigned int totalColumns,
191  const unsigned int curRow, const unsigned int curColumn,
192  char transpose = 'N');
193  };
194 
195  /// \internal
196  /// Symmetric matrices use upper triangular packed storage.
197  struct LIB_UTILITIES_EXPORT SymmetricMatrixFuncs : private TriangularMatrixFuncs
198  {
199  using TriangularMatrixFuncs::GetRequiredStorageSize;
200 
201  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn);
202 
203  static boost::tuples::tuple<unsigned int, unsigned int>
204  Advance(const unsigned int totalRows, const unsigned int totalColumns,
205  const unsigned int curRow, const unsigned int curColumn);
206  };
207 
208  struct LIB_UTILITIES_EXPORT DiagonalMatrixFuncs
209  {
210  static boost::tuples::tuple<unsigned int, unsigned int>
211  Advance(const unsigned int totalRows, const unsigned int totalColumns,
212  const unsigned int curRow, const unsigned int curColumn);
213 
214  template<typename DataType>
215  static void Invert(unsigned int rows, unsigned int columns,
216  Array<OneD, DataType>& data)
217  {
218  ASSERTL0(rows==columns, "Only square matrices can be inverted.");
219  for(unsigned int i = 0; i < rows; ++i)
220  {
221  data[i] = 1.0/data[i];
222  }
223  }
224 
225  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns);
226 
227  static unsigned int CalculateIndex(unsigned int row, unsigned int col);
228  };
229 
230 
231  struct LIB_UTILITIES_EXPORT TriangularBandedMatrixFuncs
232  {
233  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns,
234  unsigned int nSubSuperDiags);
235  };
236 
237  struct LIB_UTILITIES_EXPORT UpperTriangularBandedMatrixFuncs : public TriangularBandedMatrixFuncs
238  {
239  };
240 
241  struct LIB_UTILITIES_EXPORT LowerTriangularBandedMatrixFuncs : public TriangularBandedMatrixFuncs
242  {
243  };
244 
245  /// \internal
246  /// Symmetric banded matrices use upper triangular banded packed storage.
247  struct LIB_UTILITIES_EXPORT SymmetricBandedMatrixFuncs : private TriangularBandedMatrixFuncs
248  {
249  using TriangularBandedMatrixFuncs::GetRequiredStorageSize;
250 
251  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn,
252  unsigned int nSuperDiags);
253  };
254 
255 }
256 
257 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2