Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
MatrixFuncs.cpp
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////////////////
2 //
3 // File: MatrixFuncs.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: 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 
39 
40 namespace Nektar
41 {
42 
43  unsigned int BandedMatrixFuncs::GetRequiredStorageSize(unsigned int totalRows, unsigned int totalColumns,
44  unsigned int subDiags, unsigned int superDiags)
45  {
46  return CalculateNumberOfRows(totalRows, subDiags, superDiags)*totalColumns;
47  }
48 
49  unsigned int BandedMatrixFuncs::CalculateNumberOfDiags(unsigned int totalRows, unsigned int diags)
50  {
51  if( diags != std::numeric_limits<unsigned int>::max() )
52  {
53  return diags;
54  }
55  else if( totalRows > 0 )
56  {
57  return totalRows-1;
58  }
59  else
60  {
61  return 0;
62  }
63  }
64 
65  unsigned int BandedMatrixFuncs::CalculateNumberOfRows(unsigned int totalRows, unsigned int subDiags, unsigned int superDiags)
66  {
67  return CalculateNumberOfDiags(totalRows, subDiags) + CalculateNumberOfDiags(totalRows, superDiags) + 1;
68  }
69 
70  unsigned int BandedMatrixFuncs::CalculateIndex(unsigned int totalRows,
71  unsigned int totalColumns,
72  unsigned int row, unsigned int column,
73  unsigned int sub, unsigned int super)
74  {
75  if( (column <= row && (row - column) <= CalculateNumberOfDiags(totalRows, sub)) ||
76  (column > row && (column - row) <= CalculateNumberOfDiags(totalRows, super)) )
77  {
78  unsigned int elementRow = CalculateNumberOfDiags(totalRows, super)+row-column;
79  unsigned int elementColumn = column;
80 
81  return elementRow + elementColumn*CalculateNumberOfRows(totalRows, sub, super);
82  }
83  else
84  {
85  return std::numeric_limits<unsigned int>::max();
86  }
87  }
88 
89 
90  boost::tuples::tuple<unsigned int, unsigned int>
91  BandedMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
92  const unsigned int curRow, const unsigned int curColumn)
93  {
94  unsigned int nextRow = curRow;
95  unsigned int nextColumn = curColumn;
96 
97  //if( transpose == 'N' )
98  //{
99  // if( (column <= row && (row - column) <= data.GetNumberOfSubDiagonals(totalRows)) ||
100  // (column > row && (column - row) <= data.GetNumberOfSuperDiagonals(totalRows)) )
101  // {
102  // unsigned int arrayColumns = totalColumns;
103 
104  // unsigned int elementRow = data.GetNumberOfSuperDiagonals(totalRows)+row-column;
105  // unsigned int elementColumn = column;
106 
107  // return elementRow + elementColumn*CalculateNumberOfRows(totalRows, data);
108  // }
109  // else
110  // {
111  // return boost::optional<unsigned int>();
112  // }
113  //}
114  //else
115  //{
116  //}
117  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
118  }
119 
120 
121  unsigned int FullMatrixFuncs::GetRequiredStorageSize(unsigned int rows, unsigned int columns)
122  {
123  return rows*columns;
124  }
125 
126  unsigned int FullMatrixFuncs::CalculateIndex(unsigned int totalRows, unsigned int totalColumns, unsigned int curRow, unsigned int curColumn)
127  {
128  return curColumn*totalRows + curRow;
129  }
130 
131 
132  boost::tuples::tuple<unsigned int, unsigned int>
133  FullMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
134  const unsigned int curRow, const unsigned int curColumn)
135  {
136  unsigned int nextRow = curRow;
137  unsigned int nextColumn = curColumn;
138 
139  if( nextRow < totalRows )
140  {
141  ++nextRow;
142  }
143 
144  if( nextRow >= totalRows )
145  {
146  nextRow = 0;
147  ++nextColumn;
148  }
149 
150  if( nextColumn >= totalColumns )
151  {
152  nextRow = std::numeric_limits<unsigned int>::max();
153  nextColumn = std::numeric_limits<unsigned int>::max();
154  }
155 
156  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
157  }
158 
159 
160  unsigned int TriangularMatrixFuncs::GetRequiredStorageSize(unsigned int rows, unsigned int columns)
161  {
162  ASSERTL0(rows==columns, "Triangular matrices must be square.");
163  return rows*(rows+1)/2;
164  }
165 
166  unsigned int UpperTriangularMatrixFuncs::CalculateIndex(unsigned int curRow, unsigned int curColumn)
167  {
168  if( curRow <= curColumn )
169  {
170  return curRow + curColumn*(curColumn+1)/2;
171  }
172  else
173  {
174  return std::numeric_limits<unsigned int>::max();
175  }
176  }
177 
178  boost::tuples::tuple<unsigned int, unsigned int>
179  UpperTriangularMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
180  const unsigned int curRow, const unsigned int curColumn)
181  {
182  ASSERTL1(totalRows == totalColumns, "Triangular matrices must be square.");
183  ASSERTL1(curRow < totalRows, "Attemping to iterate through an element on row " +
184  boost::lexical_cast<std::string>(curRow) + " of a (" +
185  boost::lexical_cast<std::string>(totalRows) + ", " +
186  boost::lexical_cast<std::string>(totalColumns) + " upper triangular matrix.");
187  ASSERTL1(curColumn < totalColumns, "Attemping to iterate through an element on row " +
188  boost::lexical_cast<std::string>(curColumn) + " of a (" +
189  boost::lexical_cast<std::string>(totalRows) + ", " +
190  boost::lexical_cast<std::string>(totalColumns) + " upper triangular matrix.");
191  ASSERTL1(curRow <= curColumn, "Attemping to iterate through element (" +
192  boost::lexical_cast<std::string>(curRow) + ", " +
193  boost::lexical_cast<std::string>(curColumn) + ") of a (" +
194  boost::lexical_cast<std::string>(totalRows) + ", " +
195  boost::lexical_cast<std::string>(totalColumns) + " upper triangular matrix.");
196 
197  unsigned int nextRow = curRow;
198  unsigned int nextColumn = curColumn;
199 
200  if( nextRow <= nextColumn )
201  {
202  ++nextRow;
203  }
204 
205  if( nextRow > nextColumn )
206  {
207  ++nextColumn;
208  nextRow = 0;
209  }
210 
211  if( nextColumn >= totalColumns )
212  {
213  nextRow = std::numeric_limits<unsigned int>::max();
214  nextColumn = std::numeric_limits<unsigned int>::max();
215  }
216 
217  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
218  }
219 
220 
221  unsigned int LowerTriangularMatrixFuncs::CalculateIndex(unsigned int totalColumns, unsigned int curRow, unsigned int curColumn)
222  {
223  if( curRow >= curColumn )
224  {
225  return curRow + (2*totalColumns - curColumn - 1)*(curColumn)/2;
226  }
227  else
228  {
229  return std::numeric_limits<unsigned int>::max();
230  }
231  }
232 
233  boost::tuples::tuple<unsigned int, unsigned int>
234  LowerTriangularMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
235  const unsigned int curRow, const unsigned int curColumn,
236  char transpose)
237  {
238  ASSERTL1(totalRows == totalColumns, "Triangular matrices must be square.");
239  ASSERTL1(curRow < totalRows, "Attemping to iterate through an element on row " +
240  boost::lexical_cast<std::string>(curRow) + " of a (" +
241  boost::lexical_cast<std::string>(totalRows) + ", " +
242  boost::lexical_cast<std::string>(totalColumns) + " lower triangular matrix.");
243  ASSERTL1(curColumn < totalColumns, "Attemping to iterate through an element on row " +
244  boost::lexical_cast<std::string>(curColumn) + " of a (" +
245  boost::lexical_cast<std::string>(totalRows) + ", " +
246  boost::lexical_cast<std::string>(totalColumns) + " lower triangular matrix.");
247  ASSERTL1(curRow >= curColumn, "Attemping to iterate through element (" +
248  boost::lexical_cast<std::string>(curRow) + ", " +
249  boost::lexical_cast<std::string>(curColumn) + ") of a (" +
250  boost::lexical_cast<std::string>(totalRows) + ", " +
251  boost::lexical_cast<std::string>(totalColumns) + " lower triangular matrix.");
252 
253  if( transpose == 'T' )
254  {
256  totalColumns, totalRows, curColumn, curRow);
257  }
258 
259  unsigned int nextRow = curRow;
260  unsigned int nextColumn = curColumn;
261 
262  if( nextRow < totalRows )
263  {
264  ++nextRow;
265  }
266 
267  if( nextRow >= totalRows )
268  {
269  ++nextColumn;
270  nextRow = nextColumn;
271  }
272 
273  if( nextColumn >= totalColumns )
274  {
275  nextRow = std::numeric_limits<unsigned int>::max();
276  nextColumn = std::numeric_limits<unsigned int>::max();
277  }
278 
279  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
280  }
281 
282 
283  unsigned int SymmetricMatrixFuncs::CalculateIndex(unsigned int curRow, unsigned int curColumn)
284  {
285  if( curRow <= curColumn )
286  {
287  return curRow + curColumn*(curColumn+1)/2;
288  }
289  else
290  {
291  return curColumn + curRow*(curRow + 1)/2;
292  }
293  }
294 
295  boost::tuples::tuple<unsigned int, unsigned int>
296  SymmetricMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
297  const unsigned int curRow, const unsigned int curColumn)
298  {
299  ASSERTL1(totalRows == totalColumns, "Symmetric matrices must be square.");
300  ASSERTL1(curRow < totalRows, "Attemping to iterate through an element on row " +
301  boost::lexical_cast<std::string>(curRow) + " of a (" +
302  boost::lexical_cast<std::string>(totalRows) + ", " +
303  boost::lexical_cast<std::string>(totalColumns) + " symmetric matrix.");
304  ASSERTL1(curColumn < totalColumns, "Attemping to iterate through an element on row " +
305  boost::lexical_cast<std::string>(curColumn) + " of a (" +
306  boost::lexical_cast<std::string>(totalRows) + ", " +
307  boost::lexical_cast<std::string>(totalColumns) + " symmetric matrix.");
308 
309  unsigned int nextRow = curRow;
310  unsigned int nextColumn = curColumn;
311 
312  if( nextRow < totalRows )
313  {
314  ++nextRow;
315  }
316 
317  if( nextRow >= totalRows )
318  {
319  nextRow = 0;
320  ++nextColumn;
321  }
322 
323  if( nextColumn >= totalColumns )
324  {
325  nextRow = std::numeric_limits<unsigned int>::max();
326  nextColumn = std::numeric_limits<unsigned int>::max();
327  }
328 
329  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
330  }
331 
332  boost::tuples::tuple<unsigned int, unsigned int>
333  DiagonalMatrixFuncs::Advance(const unsigned int totalRows, const unsigned int totalColumns,
334  const unsigned int curRow, const unsigned int curColumn)
335  {
336  ASSERTL0(curRow == curColumn, "Iteration of a diagonal matrix is only valid along the diagonal.");
337 
338  unsigned int nextRow = curRow;
339  unsigned int nextColumn = curColumn;
340 
341  if( nextRow < totalRows )
342  {
343  ++nextRow;
344  ++nextColumn;
345  }
346 
347  if( nextRow >= totalRows )
348  {
349  nextRow = std::numeric_limits<unsigned int>::max();
350  nextColumn = std::numeric_limits<unsigned int>::max();
351  }
352 
353  return boost::tuples::tuple<unsigned int, unsigned int>(nextRow, nextColumn);
354  }
355 
356  unsigned int DiagonalMatrixFuncs::GetRequiredStorageSize(unsigned int rows, unsigned int columns)
357  {
358  ASSERTL0(rows==columns, "Diagonal matrices must be square.");
359  return rows;
360  }
361 
362  unsigned int DiagonalMatrixFuncs::CalculateIndex(unsigned int row, unsigned int col)
363  {
364  if( row == col )
365  {
366  return row;
367  }
368  else
369  {
370  return std::numeric_limits<unsigned int>::max();
371  }
372  }
373 
374  unsigned int TriangularBandedMatrixFuncs::GetRequiredStorageSize(unsigned int rows, unsigned int columns,
375  unsigned int nSubSuperDiags)
376  {
377  ASSERTL0(rows==columns, "Triangular matrices must be square.");
378  return (nSubSuperDiags+1)*columns;
379  }
380 
381  unsigned int SymmetricBandedMatrixFuncs::CalculateIndex(unsigned int curRow, unsigned int curColumn,
382  unsigned int nSuperDiags)
383  {
384  if( curRow <= curColumn )
385  {
386  if( (curColumn - curRow) <= nSuperDiags )
387  {
388  unsigned int elementRow = nSuperDiags - (curColumn - curRow);
389  unsigned int elementColumn = curColumn;
390 
391  return elementRow + elementColumn*(nSuperDiags+1);
392  }
393  else
394  {
395  return std::numeric_limits<unsigned int>::max();
396  }
397  }
398  else
399  {
400  return CalculateIndex(curColumn,curRow,nSuperDiags);
401  }
402 
403  }
404 }
405 
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static unsigned int CalculateIndex(unsigned int totalRows, unsigned int totalColumns, unsigned int row, unsigned int column, unsigned int sub, unsigned int super)
Definition: MatrixFuncs.cpp:70
static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn, unsigned int nSuperDiags)
static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn)
static unsigned int CalculateNumberOfDiags(unsigned int totalRows, unsigned int diags)
Definition: MatrixFuncs.cpp:49
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)
static unsigned int CalculateIndex(unsigned int totalColumns, unsigned int curRow, unsigned int curColumn)
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns, unsigned int nSubSuperDiags)
static unsigned int CalculateIndex(unsigned int curRow, unsigned int curColumn)
static unsigned int GetRequiredStorageSize(unsigned int rows, unsigned int columns)
static unsigned int CalculateIndex(unsigned int row, unsigned int col)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn, char transpose= 'N')
static unsigned int CalculateIndex(unsigned int totalRows, unsigned int totalColumns, unsigned int curRow, unsigned int curColumn)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
static unsigned int CalculateNumberOfRows(unsigned int totalRows, unsigned int subDiags, unsigned int superDiags)
Definition: MatrixFuncs.cpp:65
static boost::tuples::tuple< unsigned int, unsigned int > Advance(const unsigned int totalRows, const unsigned int totalColumns, const unsigned int curRow, const unsigned int curColumn)
Definition: MatrixFuncs.cpp:91
static unsigned int GetRequiredStorageSize(unsigned int totalRows, unsigned int totalColumns, unsigned int subDiags, unsigned int superDiags)
Calculates and returns the storage size required.
Definition: MatrixFuncs.cpp:43