itkFEMLinearSystemWrapperDenseVNL.h Source File

00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkFEMLinearSystemWrapperDenseVNL.h,v $
00005   Language:  C++
00006   Date:      $Date: 2003/02/26 20:25:53 $
00007   Version:   $Revision: 1.5 $
00008 
00009   Copyright (c) 2002 Insight Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even 
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR 
00014      PURPOSE.  See the above copyright notices for more information.
00015 
00016 =========================================================================*/
00017 
00018 #ifndef __itkFEMLinearSystemWrapperDenseVNL_h
00019 #define __itkFEMLinearSystemWrapperDenseVNL_h
00020 #include "itkFEMLinearSystemWrapper.h"
00021 #include "vnl/vnl_matrix.h"
00022 #include "vnl/vnl_vector.h"
00023 #include "vxl/vnl/algo/vnl_svd.h"
00024 #include <vxl/vnl/algo/vnl_lsqr.h>
00025 #include <vector>
00026 
00027 
00028 namespace itk {
00029 namespace fem {
00030 
00031 
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   virtual ~LinearSystemWrapperDenseVNL();
00057 
00058   /* memory management routines */
00059   virtual void  InitializeMatrix(unsigned int matrixIndex);
00060   virtual bool  IsMatrixInitialized(unsigned int matrixIndex);
00061   virtual void  DestroyMatrix(unsigned int matrixIndex);
00062   virtual void  InitializeVector(unsigned int vectorIndex);
00063   virtual bool  IsVectorInitialized(unsigned int vectorIndex);
00064   virtual void  DestroyVector(unsigned int vectorIndex);
00065   virtual void  InitializeSolution(unsigned int solutionIndex);
00066   virtual bool  IsSolutionInitialized(unsigned int solutionIndex);
00067   virtual void  DestroySolution(unsigned int solutionIndex);
00068   virtual void  SetMaximumNonZeroValuesInMatrix(unsigned int, unsigned int) {}
00069 
00070 
00071   /* assembly & solving routines */
00072   virtual Float GetMatrixValue(unsigned int i, unsigned int j, unsigned int matrixIndex) const { return (*((*m_Matrices)[matrixIndex]))(i,j); }
00073   virtual void  SetMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) { (*((*m_Matrices)[matrixIndex]))(i,j) =  value; }
00074   virtual void  AddMatrixValue(unsigned int i, unsigned int j, Float value, unsigned int matrixIndex) { (*((*m_Matrices)[matrixIndex]))(i,j) += value; }
00075   virtual Float GetVectorValue(unsigned int i, unsigned int vectorIndex) const { return (* ( (*m_Vectors)[vectorIndex] ) )[i]; }
00076   virtual void  SetVectorValue(unsigned int i, Float value, unsigned int vectorIndex) { (*((*m_Vectors)[vectorIndex]))(i) =  value; }
00077   virtual void  AddVectorValue(unsigned int i, Float value, unsigned int vectorIndex) { (*((*m_Vectors)[vectorIndex]))(i) += value; }
00078   virtual Float GetSolutionValue(unsigned int i, unsigned int solutionIndex) const;
00079   virtual void  SetSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) { (*((*m_Solutions)[solutionIndex]))(i) =  value; }
00080   virtual void  AddSolutionValue(unsigned int i, Float value, unsigned int solutionIndex) { (*((*m_Solutions)[solutionIndex]))(i) += value; }
00081   virtual void  Solve(void);
00082 
00083   /* matrix & vector manipulation routines */
00084   virtual void  ScaleMatrix(Float scale, unsigned int matrixIndex);
00085   virtual void  ScaleVector(Float scale, unsigned int vectorIndex);
00086   virtual void  ScaleSolution(Float scale, unsigned int solutionIndex);
00087   virtual void  SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00088   virtual void  SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00089   virtual void  SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00090   virtual void  CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00091   virtual void  CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00092   virtual void  MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);
00093   virtual void  MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00094 
00095 private:
00096 
00098   //std::vector< vnl_sparse_matrix<Float>* > *m_Matrices;
00099   MatrixHolder *m_Matrices;
00100 
00102   std::vector< vnl_vector<Float>* > *m_Vectors;
00103 
00105   std::vector< vnl_vector<Float>* > *m_Solutions;
00106 
00107 };
00108 
00109 }} // end namespace itk::fem
00110 
00111 #endif // #ifndef __itkFEMLinearSystemWrapperDenseVNL_h
00112 
00113