00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
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 "vnl/algo/vnl_svd.h"
00024
#include <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
00043 typedef LinearSystemWrapper::Float
Float;
00044
00045
00046 typedef LinearSystemWrapper SuperClass;
00047
00048
00049 typedef vnl_matrix<Float> MatrixRepresentation;
00050
00051
00052 typedef std::vector< MatrixRepresentation* >
MatrixHolder;
00053
00054
00055 LinearSystemWrapperDenseVNL() :
LinearSystemWrapper(), m_Matrices(0), m_Vectors(0), m_Solutions(0) {}
00056
virtual ~LinearSystemWrapperDenseVNL();
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 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
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 }}
00110
00111
#endif // #ifndef __itkFEMLinearSystemWrapperDenseVNL_h
00112
00113