00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkFEMLinearSystemWrapperVNL_h
00019 #define __itkFEMLinearSystemWrapperVNL_h
00020 #include "itkFEMLinearSystemWrapper.h"
00021 #include "vnl/vnl_sparse_matrix.h"
00022 #include "vnl/vnl_vector.h"
00023 #include <vnl/vnl_sparse_matrix_linear_system.h>
00024 #include <vnl/algo/vnl_lsqr.h>
00025 #include <vector>
00026
00027
00028 namespace itk {
00029 namespace fem {
00030
00031
00038 class LinearSystemWrapperVNL : public LinearSystemWrapper
00039 {
00040 public:
00041
00042
00043 typedef LinearSystemWrapper::Float Float;
00044
00045
00046 typedef LinearSystemWrapper SuperClass;
00047
00048
00049 typedef vnl_sparse_matrix<Float> MatrixRepresentation;
00050
00051
00052 typedef std::vector< MatrixRepresentation* > MatrixHolder;
00053
00054
00055 LinearSystemWrapperVNL() : LinearSystemWrapper(), m_Matrices(0), m_Vectors(0), m_Solutions(0) {}
00056 virtual ~LinearSystemWrapperVNL();
00057
00058
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
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
00084 virtual void ScaleMatrix(Float scale, unsigned int matrixIndex);
00085 virtual void SwapMatrices(unsigned int matrixIndex1, unsigned int matrixIndex2);
00086 virtual void SwapVectors(unsigned int vectorIndex1, unsigned int vectorIndex2);
00087 virtual void SwapSolutions(unsigned int solutionIndex1, unsigned int solutionIndex2);
00088 virtual void CopySolution2Vector(unsigned solutionIndex, unsigned int vectorIndex);
00089 virtual void CopyVector2Solution(unsigned int vectorIndex, unsigned int solutionIndex);
00090 virtual void MultiplyMatrixMatrix(unsigned int resultMatrixIndex, unsigned int leftMatrixIndex, unsigned int rightMatrixIndex);
00091 virtual void MultiplyMatrixVector(unsigned int resultVectorIndex, unsigned int matrixIndex, unsigned int vectorIndex);
00092
00093 private:
00094
00096
00097 MatrixHolder *m_Matrices;
00098
00100 std::vector< vnl_vector<Float>* > *m_Vectors;
00101
00103 std::vector< vnl_vector<Float>* > *m_Solutions;
00104
00105 };
00106
00107 }}
00108
00109 #endif
00110