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