Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 
38 #include <LibUtilities/LinearAlgebra/MatrixFuncs.h>
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  {
255  return UpperTriangularMatrixFuncs::Advance(
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 
Generated on Sun Mar 15 2015 00:11:32 for Nektar++ by   doxygen 1.8.1.2