Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties 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>
47 
48 namespace Nektar
49 {
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 
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,
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,
135  Array<OneD, NekDouble> &EigValReal,
136  Array<OneD, NekDouble> &EigValImag,
138  {
139  int lda = n,info = 0;
140  NekDouble dum;
141  char uplo = 'N';
142 
143  if(EigVecs != NullNekDouble1DArray) // calculate Right Eigen Vectors
144  {
145  int lwork = 4*lda;
146  Array<OneD,NekDouble> work(4*lda);
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  int lwork = 3*lda;
158  Array<OneD,NekDouble> work(3*lda);
159  char lrev = 'N';
160  Lapack::Dgeev(uplo,lrev,lda,
161  A.get(),lda,
162  EigValReal.get(),
163  EigValImag.get(),
164  &dum,1,&dum,1,
165  &work[0],lwork,info);
166  }
167  ASSERTL0(info == 0,"Info is not zero");
168 
169  }
170  };
171 
173  {
174  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns);
175  };
176 
178  {
179  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn);
180 
181  static boost::tuples::tuple<unsigned int, unsigned int>
182  Advance(const unsigned int totalRows, const unsigned int totalColumns,
183  const unsigned int curRow, const unsigned int curColumn);
184  };
185 
186 
188  {
189  static unsigned int CalculateIndex(unsigned int totalColumns, unsigned int curRow, unsigned int curColumn);
190 
191  static boost::tuples::tuple<unsigned int, unsigned int>
192  Advance(const unsigned int totalRows, const unsigned int totalColumns,
193  const unsigned int curRow, const unsigned int curColumn,
194  char transpose = 'N');
195  };
196 
197  /// \internal
198  /// Symmetric matrices use upper triangular packed storage.
200  {
202 
203  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn);
204 
205  template<typename DataType>
206  static void Invert(unsigned int rows, unsigned int columns,
207  Array<OneD, DataType>& data)
208  {
209 #ifdef NEKTAR_USING_BLAS
210  ASSERTL0(rows==columns, "Only square matrices can be inverted.");
211 
212  int n = columns;
213  int info = 0;
214  Array<OneD, int> ipivot(n);
215  Array<OneD, DataType> work(n);
216 
217  Lapack::Dsptrf('U', n, data.get(), ipivot.get(), info);
218 
219  if( info < 0 )
220  {
221  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dsptrf";
222  ASSERTL0(false, message.c_str());
223  }
224  else if( info > 0 )
225  {
226  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dsptrf";
227  ASSERTL0(false, message.c_str());
228  }
229 
230  Lapack::Dsptri('U', n, data.get(), ipivot.get(),
231  work.get(), info);
232 
233  if( info < 0 )
234  {
235  std::string message = "ERROR: The " + boost::lexical_cast<std::string>(-info) + "th parameter had an illegal parameter for dsptri";
236  ASSERTL0(false, message.c_str());
237  }
238  else if( info > 0 )
239  {
240  std::string message = "ERROR: Element u_" + boost::lexical_cast<std::string>(info) + boost::lexical_cast<std::string>(info) + " is 0 from dsptri";
241  ASSERTL0(false, message.c_str());
242  }
243 
244 #else
245  // error Full matrix inversion not supported without blas.
246  BOOST_STATIC_ASSERT(sizeof(DataType) == 0);
247 #endif
248  }
249 
250  static boost::tuples::tuple<unsigned int, unsigned int>
251  Advance(const unsigned int totalRows, const unsigned int totalColumns,
252  const unsigned int curRow, const unsigned int curColumn);
253  };
254 
256  {
257  static boost::tuples::tuple<unsigned int, unsigned int>
258  Advance(const unsigned int totalRows, const unsigned int totalColumns,
259  const unsigned int curRow, const unsigned int curColumn);
260 
261  template<typename DataType>
262  static void Invert(unsigned int rows, unsigned int columns,
263  Array<OneD, DataType>& data)
264  {
265  ASSERTL0(rows==columns, "Only square matrices can be inverted.");
266  for(unsigned int i = 0; i < rows; ++i)
267  {
268  data[i] = 1.0/data[i];
269  }
270  }
271 
272  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns);
273 
274  static unsigned int CalculateIndex(unsigned int row, unsigned int col);
275  };
276 
277 
279  {
280  static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns,
281  unsigned int nSubSuperDiags);
282  };
283 
285  {
286  };
287 
289  {
290  };
291 
292  /// \internal
293  /// Symmetric banded matrices use upper triangular banded packed storage.
295  {
297 
298  static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn,
299  unsigned int nSuperDiags);
300  };
301 
302 }
303 
304 #endif //NEKTAR_LIB_UTILITIES_LINEAR_ALGEBRA_MATRIX_FUNCS_H
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static Array< OneD, NekDouble > NullNekDouble1DArray
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data, const char transpose)
Definition: MatrixFuncs.h:86
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)
#define LIB_UTILITIES_EXPORT
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns, unsigned int nSubSuperDiags)
double NekDouble
static void EigenSolve(unsigned int n, const Array< OneD, const double > &A, Array< OneD, NekDouble > &EigValReal, Array< OneD, NekDouble > &EigValImag, Array< OneD, NekDouble > &EigVecs=NullNekDouble1DArray)
Definition: MatrixFuncs.h:133
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:206
static void Invert(unsigned int rows, unsigned int columns, Array< OneD, DataType > &data)
Definition: MatrixFuncs.h:262