ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMLinearSystemWrapper.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
00016  *
00017  *=========================================================================*/
00018 
00019 #ifndef __itkFEMLinearSystemWrapper_h
00020 #define __itkFEMLinearSystemWrapper_h
00021 
00022 #include "itkMacro.h"
00023 #include "itkFEMSolution.h"
00024 #include "itkFEMException.h"
00025 
00026 #include <vector>
00027 #include <typeinfo>
00028 #include <string>
00029 
00030 namespace itk
00031 {
00032 namespace fem
00033 {
00053 class LinearSystemWrapper : public Solution
00054 {
00055 public:
00056   typedef LinearSystemWrapper Self;
00057   typedef Solution            Superclass;
00058   typedef Self *              Pointer;
00059   typedef const Self *        ConstPointer;
00060 
00061   typedef std::vector<unsigned int> ColumnArray;
00062 
00067   LinearSystemWrapper() :
00068     m_Order(0), m_NumberOfMatrices(1), m_NumberOfVectors(1), m_NumberOfSolutions(1)
00069   {
00070   }
00071   /* , m_PrimaryMatrixSetupFunction(0), m_PrimaryVectorSetupFunction(0),
00072     m_PrimarySolutionSetupFunction(0) {} */
00074 
00079   virtual ~LinearSystemWrapper()
00080   {
00081   }
00082 
00087   virtual void Clean(void);
00088 
00094   void SetSystemOrder(unsigned int N)
00095   {
00096     m_Order = N;
00097   }
00098 
00102   unsigned int GetSystemOrder() const
00103   {
00104     return m_Order;
00105   }
00106 
00111   void SetNumberOfMatrices(unsigned int nMatrices)
00112   {
00113     m_NumberOfMatrices = nMatrices;
00114   }
00115 
00116   /*
00117    * Set the maximum number of entries permitted in a matrix
00118    * \param matrixIndex index of matrix to set value for
00119    * \param maxNonZeros maximum number of entries allowed in matrix
00120    * \note in general this function does nothing, however it may
00121    *       redefined by the derived wrapper if necessary
00122    */
00123   // virtual void SetMaximumNonZeroValuesInMatrix(unsigned int maxNonZeroValues)
00124   // = 0;
00125 
00129   unsigned int GetNumberOfMatrices() const
00130   {
00131     return m_NumberOfMatrices;
00132   }
00133 
00138   void SetNumberOfVectors(unsigned int nVectors)
00139   {
00140     m_NumberOfVectors = nVectors;
00141   }
00142 
00146   unsigned int GetNumberOfVectors() const
00147   {
00148     return m_NumberOfVectors;
00149   }
00150 
00155   void SetNumberOfSolutions(unsigned int nSolutions)
00156   {
00157     m_NumberOfSolutions = nSolutions;
00158   }
00159 
00163   unsigned int GetNumberOfSolutions() const
00164   {
00165     return m_NumberOfSolutions;
00166   }
00167 
00175   virtual void InitializeMatrix(unsigned int matrixIndex = 0) = 0;
00176 
00181   virtual bool IsMatrixInitialized(unsigned int matrixIndex = 0) = 0;
00182 
00187   virtual void DestroyMatrix(unsigned int matrixIndex = 0) = 0;
00188 
00195   virtual void InitializeVector(unsigned int vectorIndex = 0) = 0;
00196 
00201   virtual bool IsVectorInitialized(unsigned int vectorIndex = 0) = 0;
00202 
00207   virtual void DestroyVector(unsigned int vectorIndex = 0) = 0;
00208 
00215   virtual void InitializeSolution(unsigned int solutionIndex = 0) = 0;
00216 
00221   virtual bool IsSolutionInitialized(unsigned int solutionIndex = 0) = 0;
00222 
00226   virtual void DestroySolution(unsigned int solutionIndex = 0) = 0;
00227 
00234   virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex = 0) const = 0;
00235 
00243   virtual void SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex = 0) = 0;
00244 
00252   virtual void AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex = 0) = 0;
00253 
00266   virtual void GetColumnsOfNonZeroMatrixElementsInRow(unsigned int row, ColumnArray & cols,
00267                                                       unsigned int matrixIndex = 0);
00268 
00274   virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex = 0) const = 0;
00275 
00282   virtual void SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex = 0) = 0;
00283 
00290   virtual void AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex = 0) = 0;
00291 
00299   virtual void SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex = 0) = 0;
00300 
00308   virtual void AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex = 0) = 0;
00309 
00317   virtual void Solve(void) = 0;
00318 
00324   virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2) = 0;
00325 
00333   virtual void CopyMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);
00334 
00340   virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2) = 0;
00341 
00347   virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2) = 0;
00348 
00354   virtual void ScaleMatrix(Float scale, unsigned int matrixIndex = 0);
00355 
00361   void ScaleVector(Float scale, unsigned int vectorIndex = 0);
00362 
00368   void ScaleSolution(Float scale, unsigned int solutionIndex = 0);
00369 
00376   virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex,
00377                                     unsigned int rightMatrixIndex) = 0;
00378 
00385   virtual void AddMatrixMatrix(unsigned int matrixIndex1, unsigned int matrixIndex2);
00386 
00393   virtual void AddVectorVector(unsigned int vectorIndex1, unsigned int vectorIndex2);
00394 
00401   virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00402 
00409   virtual void MultiplyMatrixSolution(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int solutionIndex);
00410 
00416   virtual void CopySolution2Vector(unsigned int solutionIndex, unsigned int vectorIndex) = 0;
00417 
00423   virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex) = 0;
00424 
00430   virtual void CopyVector(unsigned int vectorSource, unsigned int vectorDestination);
00431 
00438   virtual void OptimizeMatrixStorage(unsigned int matrixIndex, unsigned int tempMatrixIndex);
00439 
00445   virtual void ReverseCuthillMckeeOrdering(ColumnArray & newNumbering, unsigned int matrixIndex = 0);
00446 
00447 protected:
00448 
00450   unsigned int m_Order;
00451 
00455   unsigned int m_NumberOfMatrices;
00456 
00460   unsigned int m_NumberOfVectors;
00461 
00465   unsigned int m_NumberOfSolutions;
00466 
00467   /*
00468    * Function used to prepare primary matrix for numerical solving
00469    */
00470   // void (*m_PrimaryMatrixSetupFunction)(LinearSystemWrapper *lsw);
00471 
00472   /*
00473    * Function used to prepare primary vector for numerical solving
00474    */
00475   /* void (*m_PrimaryVectorSetupFunction)(LinearSystemWrapper *lsw);*/
00476 
00477   /*
00478    * Function used to prepare primary matrix for numerical solving
00479    */
00480   /* void (*m_PrimarySolutionSetupFunction)(LinearSystemWrapper *lsw); */
00481 private:
00482 
00486   void CuthillMckeeOrdering(ColumnArray & newNumbering, int startingRow, unsigned int matrixIndex = 0);
00487 
00488   void FollowConnectionsCuthillMckeeOrdering(unsigned int rowNumber, ColumnArray & rowDegree,
00489                                              ColumnArray & newNumbering, unsigned int nextRowNumber,
00490                                              unsigned int matrixIndex = 0);
00491 
00493   LinearSystemWrapper(const LinearSystemWrapper &);
00494 
00496   const LinearSystemWrapper & operator=(const LinearSystemWrapper &);
00497 
00498 };
00499 
00500 class ITK_ABI_EXPORT FEMExceptionLinearSystem : public FEMException
00501 {
00502 public:
00508   FEMExceptionLinearSystem(const char *file, unsigned int lineNumber, std::string location, std::string moreDescription);
00509 
00511   virtual ~FEMExceptionLinearSystem()
00512   throw ( )
00513   {
00514   }
00515 
00517   itkTypeMacro(FEMExceptionLinearSystem, FEMException);
00518 };
00519 
00520 class ITK_ABI_EXPORT FEMExceptionLinearSystemBounds : public FEMException
00521 {
00522 public:
00528   FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location,
00529                                  std::string moreDescription,
00530                                  unsigned int index1);
00531 
00536   FEMExceptionLinearSystemBounds(const char *file, unsigned int lineNumber, std::string location,
00537                                  std::string moreDescription, unsigned int index1,
00538                                  unsigned int index2);
00539 
00541   virtual ~FEMExceptionLinearSystemBounds()
00542   throw ( )
00543   {
00544   }
00545 
00547   itkTypeMacro(FEMExceptionLinearSystem, FEMException);
00548 };
00549 }
00550 }  // end namespace itk::fem
00551 
00552 #endif // #ifndef __itkFEMLinearSystemWrapper_h
00553