ITK
4.1.0
Insight Segmentation and Registration Toolkit
|
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