ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkFEMLinearSystemWrapperDenseVNL.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 __itkFEMLinearSystemWrapperDenseVNL_h
00020 #define __itkFEMLinearSystemWrapperDenseVNL_h
00021 #include "itkFEMLinearSystemWrapper.h"
00022 #include "vnl/vnl_matrix.h"
00023 #include "vnl/vnl_vector.h"
00024 #include "vnl/algo/vnl_svd.h"
00025 #include <vector>
00026 
00027 namespace itk
00028 {
00029 namespace fem
00030 {
00038 class LinearSystemWrapperDenseVNL : public LinearSystemWrapper
00039 {
00040 public:
00041 
00042   /* values stored in matrices & vectors */
00043   typedef LinearSystemWrapper::Float Float;
00044 
00045   /* superclass */
00046   typedef LinearSystemWrapper SuperClass;
00047 
00048   /* matrix typedef */
00049   typedef vnl_matrix<Float> MatrixRepresentation;
00050 
00051   /* matrix holder typedef */
00052   typedef std::vector<MatrixRepresentation *> MatrixHolder;
00053 
00054   /* constructor & destructor */
00055   LinearSystemWrapperDenseVNL() : LinearSystemWrapper(), m_Matrices(0), m_Vectors(0), m_Solutions(0)
00056   {
00057   }
00058   virtual ~LinearSystemWrapperDenseVNL();
00059 
00060   /* memory management routines */
00061   virtual void  InitializeMatrix(unsigned int matrixIndex);
00062 
00063   virtual bool  IsMatrixInitialized(unsigned int matrixIndex);
00064 
00065   virtual void  DestroyMatrix(unsigned int matrixIndex);
00066 
00067   virtual void  InitializeVector(unsigned int vectorIndex);
00068 
00069   virtual bool  IsVectorInitialized(unsigned int vectorIndex);
00070 
00071   virtual void  DestroyVector(unsigned int vectorIndex);
00072 
00073   virtual void  InitializeSolution(unsigned int solutionIndex);
00074 
00075   virtual bool  IsSolutionInitialized(unsigned int solutionIndex);
00076 
00077   virtual void  DestroySolution(unsigned int solutionIndex);
00078 
00079   virtual void  SetMaximumNonZeroValuesInMatrix(unsigned int, unsigned int)
00080   {
00081   }
00082 
00083   /* assembly & solving routines */
00084   virtual Float GetMatrixValue(unsigned int i, unsigned int j,
00085                                unsigned int matrixIndex) const
00086   {
00087     return ( *( ( *m_Matrices )[matrixIndex] ) )(i, j);
00088   }
00089   virtual void  SetMatrixValue(unsigned int i, unsigned int j, Float value,
00090                                unsigned int matrixIndex)
00091   {
00092     ( *( ( *m_Matrices )[matrixIndex] ) )(i, j) =  value;
00093   }
00094   virtual void  AddMatrixValue(unsigned int i, unsigned int j, Float value,
00095                                unsigned int matrixIndex)
00096   {
00097     ( *( ( *m_Matrices )[matrixIndex] ) )(i, j) += value;
00098   }
00099   virtual Float GetVectorValue(unsigned int i,
00100                                unsigned int vectorIndex) const
00101   {
00102     return ( *( ( *m_Vectors )[vectorIndex] ) )[i];
00103   }
00104   virtual void  SetVectorValue(unsigned int i, Float value,
00105                                unsigned int vectorIndex)
00106   {
00107     ( *( ( *m_Vectors )[vectorIndex] ) )(i) =  value;
00108   }
00109   virtual void  AddVectorValue(unsigned int i, Float value,
00110                                unsigned int vectorIndex)
00111   {
00112     ( *( ( *m_Vectors )[vectorIndex] ) )(i) += value;
00113   }
00114   virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
00115 
00116   virtual void  SetSolutionValue(unsigned int i, Float value,
00117                                  unsigned int solutionIndex)
00118   {
00119     ( *( ( *m_Solutions )[solutionIndex] ) )(i) =  value;
00120   }
00121   virtual void  AddSolutionValue(unsigned int i, Float value,
00122                                  unsigned int solutionIndex)
00123   {
00124     ( *( ( *m_Solutions )[solutionIndex] ) )(i) += value;
00125   }
00126   virtual void  Solve(void);
00127 
00128   /* matrix & vector manipulation routines */
00129   virtual void  ScaleMatrix(Float scale, unsigned int matrixIndex);
00130 
00131   virtual void  ScaleVector(Float scale, unsigned int vectorIndex);
00132 
00133   virtual void  ScaleSolution(Float scale, unsigned int solutionIndex);
00134 
00135   virtual void  SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00136 
00137   virtual void  SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00138 
00139   virtual void  SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00140 
00141   virtual void  CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00142 
00143   virtual void  CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00144 
00145   virtual void  MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex,
00146                                      unsigned int rightMatrixIndex);
00147 
00148   virtual void  MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00149 
00150 private:
00151 
00153   // std::vector< vnl_sparse_matrix<Float>* > *m_Matrices;
00154   MatrixHolder *m_Matrices;
00155 
00157   std::vector<vnl_vector<Float> *> *m_Vectors;
00158 
00160   std::vector<vnl_vector<Float> *> *m_Solutions;
00161 };
00162 }
00163 }  // end namespace itk::fem
00164 
00165 #endif
00166