ITK  4.3.0
Insight Segmentation and Registration Toolkit
itkFEMLinearSystemWrapper.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright Insight Software Consortium
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 
19 #ifndef __itkFEMLinearSystemWrapper_h
20 #define __itkFEMLinearSystemWrapper_h
21 
22 #include "itkMacro.h"
23 #include "itkFEMSolution.h"
24 #include "itkFEMException.h"
25 
26 #include <vector>
27 #include <typeinfo>
28 #include <string>
29 
30 namespace itk
31 {
32 namespace fem
33 {
54 {
55 public:
58  typedef Self * Pointer;
59  typedef const Self * ConstPointer;
60 
61  typedef std::vector<unsigned int> ColumnArray;
62 
69  {
70  }
71  /* , m_PrimaryMatrixSetupFunction(0), m_PrimaryVectorSetupFunction(0),
72  m_PrimarySolutionSetupFunction(0) {} */
74 
80  {
81  }
82 
87  virtual void Clean(void);
88 
94  void SetSystemOrder(unsigned int N)
95  {
96  m_Order = N;
97  }
98 
102  unsigned int GetSystemOrder() const
103  {
104  return m_Order;
105  }
106 
111  void SetNumberOfMatrices(unsigned int nMatrices)
112  {
113  m_NumberOfMatrices = nMatrices;
114  }
115 
116  /*
117  * Set the maximum number of entries permitted in a matrix
118  * \param matrixIndex index of matrix to set value for
119  * \param maxNonZeros maximum number of entries allowed in matrix
120  * \note in general this function does nothing, however it may
121  * redefined by the derived wrapper if necessary
122  */
123  // virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues)
124  // = 0;
125 
129  unsigned int GetNumberOfMatrices() const
130  {
131  return m_NumberOfMatrices;
132  }
133 
138  void SetNumberOfVectors(unsigned int nVectors)
139  {
140  m_NumberOfVectors = nVectors;
141  }
142 
146  unsigned int GetNumberOfVectors() const
147  {
148  return m_NumberOfVectors;
149  }
150 
155  void SetNumberOfSolutions(unsigned int nSolutions)
156  {
157  m_NumberOfSolutions = nSolutions;
158  }
159 
163  unsigned int GetNumberOfSolutions() const
164  {
165  return m_NumberOfSolutions;
166  }
167 
175  virtual void InitializeMatrix(unsigned int matrixIndex = 0) = 0;
176 
181  virtual bool IsMatrixInitialized(unsigned int matrixIndex = 0) = 0;
182 
187  virtual void DestroyMatrix(unsigned int matrixIndex = 0) = 0;
188 
195  virtual void InitializeVector(unsigned int vectorIndex = 0) = 0;
196 
201  virtual bool IsVectorInitialized(unsigned int vectorIndex = 0) = 0;
202 
207  virtual void DestroyVector(unsigned int vectorIndex = 0) = 0;
208 
215  virtual void InitializeSolution(unsigned int solutionIndex = 0) = 0;
216 
221  virtual bool IsSolutionInitialized(unsigned int solutionIndex = 0) = 0;
222 
226  virtual void DestroySolution(unsigned int solutionIndex = 0) = 0;
227 
234  virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex = 0) const = 0;
235 
243  virtual void SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex = 0) = 0;
244 
252  virtual void AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex = 0) = 0;
253 
266  virtual void GetColumnsOfNonZeroMatrixElementsInRow(unsigned int row, ColumnArray & cols,
267  unsigned int matrixIndex = 0);
268 
274  virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex = 0) const = 0;
275 
282  virtual void SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex = 0) = 0;
283 
290  virtual void AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex = 0) = 0;
291 
299  virtual void SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex = 0) = 0;
300 
308  virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex = 0) = 0;
309 
317  virtual void Solve(void) = 0;
318 
324  virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2) = 0;
325 
333  virtual void CopyMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);
334 
340  virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2) = 0;
341 
347  virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2) = 0;
348 
354  virtual void ScaleMatrix(Float scale, unsigned int matrixIndex = 0);
355 
361  void ScaleVector(Float scale, unsigned int vectorIndex = 0);
362 
368  void ScaleSolution(Float scale, unsigned int solutionIndex = 0);
369 
376  virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex,
377  unsigned int rightMatrixIndex) = 0;
378 
385  virtual void AddMatrixMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);
386 
393  virtual void AddVectorVector(unsigned int vectorIndex1, unsigned int vectorIndex2);
394 
401  virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
402 
409  virtual void MultiplyMatrixSolution(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int solutionIndex);
410 
416  virtual void CopySolution2Vector(unsigned int solutionIndex, unsigned int vectorIndex) = 0;
417 
423  virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex) = 0;
424 
430  virtual void CopyVector(unsigned int vectorSource, unsigned int vectorDestination);
431 
438  virtual void OptimizeMatrixStorage(unsigned int matrixIndex, unsigned int tempMatrixIndex);
439 
445  virtual void ReverseCuthillMckeeOrdering(ColumnArray & newNumbering, unsigned int matrixIndex = 0);
446 
447 protected:
448 
450  unsigned int m_Order;
451 
455  unsigned int m_NumberOfMatrices;
456 
460  unsigned int m_NumberOfVectors;
461 
465  unsigned int m_NumberOfSolutions;
466 
467  /*
468  * Function used to prepare primary matrix for numerical solving
469  */
470  // void (*m_PrimaryMatrixSetupFunction)(LinearSystemWrapper *lsw);
471 
472  /*
473  * Function used to prepare primary vector for numerical solving
474  */
475  /* void (*m_PrimaryVectorSetupFunction)(LinearSystemWrapper *lsw);*/
476 
477  /*
478  * Function used to prepare primary matrix for numerical solving
479  */
480  /* void (*m_PrimarySolutionSetupFunction)(LinearSystemWrapper *lsw); */
481 
482 private:
483 
487  void CuthillMckeeOrdering(ColumnArray & newNumbering, int startingRow, unsigned int matrixIndex = 0);
488 
489  void FollowConnectionsCuthillMckeeOrdering(unsigned int rowNumber, ColumnArray & rowDegree,
490  ColumnArray & newNumbering, unsigned int nextRowNumber,
491  unsigned int matrixIndex = 0);
492 
495 
498 
499 };
500 
501 class ITK_ABI_EXPORT FEMExceptionLinearSystem : public FEMException
502 {
503 public:
509  FEMExceptionLinearSystem(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription);
510 
513  throw ( )
514  {
515  }
516 
518  itkTypeMacro(FEMExceptionLinearSystem, FEMException);
519 };
520 
521 class ITK_ABI_EXPORT FEMExceptionLinearSystemBounds : public FEMException
522 {
523 public:
529  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location,
530  std::string moreDescription,
531  unsigned int index1);
532 
537  FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location,
538  std::string moreDescription, unsigned int index1,
539  unsigned int index2);
540 
543  throw ( )
544  {
545  }
546 
548  itkTypeMacro(FEMExceptionLinearSystem, FEMException);
549 };
550 }
551 } // end namespace itk::fem
552 
553 #endif // #ifndef __itkFEMLinearSystemWrapper_h
554