1 /*=========================================================================
2 *
3 * Copyright Insight Software Consortium
4 *
5 * Licensed under the Apache License, Version 2.0 (the "License");
6 * you may not use this file except in compliance with the License.
7 * You may obtain a copy of the License at
8 *
9 * http://www.apache.org/licenses/LICENSE-2.0.txt
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 *
17 *=========================================================================*/
18
19 #include "itkFEMLinearSystemWrapperDenseVNL.h"
20
21 namespace itk
22 {
23 namespace fem
24 {
InitializeMatrix(unsigned int matrixIndex)25 void LinearSystemWrapperDenseVNL::InitializeMatrix(unsigned int matrixIndex)
26 {
27 // allocate if necessary
28 if( m_Matrices == nullptr )
29 {
30 m_Matrices = new MatrixHolder(m_NumberOfMatrices);
31 }
32
33 // out with old, in with new
34 delete ( *m_Matrices )[matrixIndex];
35
36 ( *m_Matrices )[matrixIndex] = new MatrixRepresentation( this->GetSystemOrder(), this->GetSystemOrder() );
37 ( *m_Matrices )[matrixIndex]->fill(0.0);
38 }
39
IsMatrixInitialized(unsigned int matrixIndex)40 bool LinearSystemWrapperDenseVNL::IsMatrixInitialized(unsigned int matrixIndex)
41 {
42 if( !m_Matrices )
43 {
44 return false;
45 }
46 if( !( ( *m_Matrices )[matrixIndex] ) )
47 {
48 return false;
49 }
50
51 return true;
52 }
53
DestroyMatrix(unsigned int matrixIndex)54 void LinearSystemWrapperDenseVNL::DestroyMatrix(unsigned int matrixIndex)
55 {
56 if( m_Matrices )
57 {
58 delete ( *m_Matrices )[matrixIndex];
59 ( *m_Matrices )[matrixIndex] = nullptr;
60 }
61 }
62
InitializeVector(unsigned int vectorIndex)63 void LinearSystemWrapperDenseVNL::InitializeVector(unsigned int vectorIndex)
64 {
65 // allocate if necessary
66 if( m_Vectors == nullptr )
67 {
68 m_Vectors = new std::vector<vnl_vector<Float> *>(m_NumberOfVectors);
69 }
70
71 // out with old, in with new
72 delete ( *m_Vectors )[vectorIndex];
73
74 ( *m_Vectors )[vectorIndex] = new vnl_vector<Float>( this->GetSystemOrder() );
75 ( *m_Vectors )[vectorIndex]->fill(0.0);
76 }
77
IsVectorInitialized(unsigned int vectorIndex)78 bool LinearSystemWrapperDenseVNL::IsVectorInitialized(unsigned int vectorIndex)
79 {
80 if( !m_Vectors )
81 {
82 return false;
83 }
84 if( !( *m_Vectors )[vectorIndex] )
85 {
86 return false;
87 }
88
89 return true;
90 }
91
DestroyVector(unsigned int vectorIndex)92 void LinearSystemWrapperDenseVNL::DestroyVector(unsigned int vectorIndex)
93 {
94 if( m_Vectors )
95 {
96 delete ( *m_Vectors )[vectorIndex];
97 ( *m_Vectors )[vectorIndex] = nullptr;
98 }
99 }
100
InitializeSolution(unsigned int solutionIndex)101 void LinearSystemWrapperDenseVNL::InitializeSolution(unsigned int solutionIndex)
102 {
103 // allocate if necessary
104 if( m_Solutions == nullptr )
105 {
106 m_Solutions = new std::vector<vnl_vector<Float> *>(m_NumberOfSolutions);
107 }
108
109 // out with old, in with new
110 delete ( *m_Solutions )[solutionIndex];
111
112 ( *m_Solutions )[solutionIndex] = new vnl_vector<Float>( this->GetSystemOrder() );
113 ( *m_Solutions )[solutionIndex]->fill(0.0);
114 }
115
IsSolutionInitialized(unsigned int solutionIndex)116 bool LinearSystemWrapperDenseVNL::IsSolutionInitialized(unsigned int solutionIndex)
117 {
118 if( !m_Solutions )
119 {
120 return false;
121 }
122 if( !( *m_Solutions )[solutionIndex] )
123 {
124 return false;
125 }
126
127 return true;
128 }
129
DestroySolution(unsigned int solutionIndex)130 void LinearSystemWrapperDenseVNL::DestroySolution(unsigned int solutionIndex)
131 {
132 if( m_Solutions )
133 {
134 delete ( *m_Solutions )[solutionIndex];
135 ( *m_Solutions )[solutionIndex] = nullptr;
136 }
137 }
138
GetSolutionValue(unsigned int i,unsigned int solutionIndex) const139 LinearSystemWrapperDenseVNL::Float LinearSystemWrapperDenseVNL::GetSolutionValue(unsigned int i,
140 unsigned int solutionIndex) const
141 {
142 if( m_Solutions == nullptr )
143 {
144 return 0.0;
145 }
146 if( ( ( *m_Solutions )[solutionIndex] )->size() <= i )
147 {
148 return 0.0;
149 }
150 else
151 {
152 return ( *( ( *m_Solutions )[solutionIndex] ) )(i);
153 }
154 }
155
Solve()156 void LinearSystemWrapperDenseVNL::Solve()
157 {
158 if( ( m_Matrices->size() == 0 ) || ( m_Vectors->size() == 0 ) || ( m_Solutions->size() == 0 ) )
159 {
160 throw FEMException(
161 __FILE__,
162 __LINE__,
163 "FEM error!");
164 }
165
166 /* use functions to make sure that zero based matrix, vector, & index store
167 final system to solve */
168 /*
169 if (m_PrimaryMatrixSetupFunction != nullptr) (*m_PrimaryMatrixSetupFunction)(static_cast<Superclass*>(this));
170 if (m_PrimaryVectorSetupFunction != nullptr) (*m_PrimaryVectorSetupFunction)(static_cast<Superclass*>(this));
171 if (m_PrimarySolutionSetupFunction != nullptr) (*m_PrimarySolutionSetupFunction)(static_cast<Superclass*>(this));
172 */
173
174 /**
175 * Solve the system of linear equation and store the result in
176 * m_Solutions(0).
177 * Here we use the SVD method.
178 */
179
180 vnl_svd<Float> svd( ( *( ( *m_Matrices )[0] ) ) );
181 ( *( ( *m_Solutions )[0] ) ) = svd.solve( ( *( ( *m_Vectors )[0] ) ) );
182 }
183
SwapMatrices(unsigned int MatrixIndex1,unsigned int MatrixIndex2)184 void LinearSystemWrapperDenseVNL::SwapMatrices(unsigned int MatrixIndex1, unsigned int MatrixIndex2)
185 {
186 vnl_matrix<Float> *tmp;
187 tmp = ( *m_Matrices )[MatrixIndex1];
188 ( *m_Matrices )[MatrixIndex1] = ( *m_Matrices )[MatrixIndex2];
189 ( *m_Matrices )[MatrixIndex2] = tmp;
190 }
191
SwapVectors(unsigned int VectorIndex1,unsigned int VectorIndex2)192 void LinearSystemWrapperDenseVNL::SwapVectors(unsigned int VectorIndex1, unsigned int VectorIndex2)
193 {
194 vnl_vector<Float> tmp;
195 tmp = *( *m_Vectors )[VectorIndex1];
196 *( *m_Vectors )[VectorIndex1] = *( *m_Vectors )[VectorIndex2];
197 *( *m_Vectors )[VectorIndex2] = tmp;
198 }
199
SwapSolutions(unsigned int SolutionIndex1,unsigned int SolutionIndex2)200 void LinearSystemWrapperDenseVNL::SwapSolutions(unsigned int SolutionIndex1, unsigned int SolutionIndex2)
201 {
202 vnl_vector<Float> *tmp;
203 tmp = ( *m_Solutions )[SolutionIndex1];
204 ( *m_Solutions )[SolutionIndex1] = ( *m_Solutions )[SolutionIndex2];
205 ( *m_Solutions )[SolutionIndex2] = tmp;
206 }
207
CopySolution2Vector(unsigned int SolutionIndex,unsigned int VectorIndex)208 void LinearSystemWrapperDenseVNL::CopySolution2Vector(unsigned int SolutionIndex, unsigned int VectorIndex)
209 {
210 delete ( *m_Vectors )[VectorIndex];
211 ( *m_Vectors )[VectorIndex] = new vnl_vector<Float>( *( ( *m_Solutions )[SolutionIndex] ) );
212 }
213
CopyVector2Solution(unsigned int VectorIndex,unsigned int SolutionIndex)214 void LinearSystemWrapperDenseVNL::CopyVector2Solution(unsigned int VectorIndex, unsigned int SolutionIndex)
215 {
216 delete ( *m_Solutions )[SolutionIndex];
217 ( *m_Solutions )[SolutionIndex] = new vnl_vector<Float>( *( ( *m_Vectors )[VectorIndex] ) );
218 }
219
MultiplyMatrixMatrix(unsigned int ResultMatrixIndex,unsigned int LeftMatrixIndex,unsigned int RightMatrixIndex)220 void LinearSystemWrapperDenseVNL::MultiplyMatrixMatrix(unsigned int ResultMatrixIndex,
221 unsigned int LeftMatrixIndex,
222 unsigned int RightMatrixIndex)
223 {
224 // delete (*m_Matrices)[ResultMatrixIndex];
225 // (*m_Matrices)[ResultMatrixIndex] = new vnl_matrix<Float>(
226 // this->GetSystemOrder(), this->GetSystemOrder() );
227
228 ( *( ( *m_Matrices )[ResultMatrixIndex] ) ) = ( *( ( *m_Matrices )[LeftMatrixIndex] ) )
229 * ( *( ( *m_Matrices )[RightMatrixIndex] ) );
230 }
231
MultiplyMatrixVector(unsigned int resultVectorIndex,unsigned int matrixIndex,unsigned int vectorIndex)232 void LinearSystemWrapperDenseVNL::MultiplyMatrixVector(unsigned int resultVectorIndex,
233 unsigned int matrixIndex,
234 unsigned int vectorIndex)
235 {
236 // delete (*m_Matrices)[ResultMatrixIndex];
237 // (*m_Matrices)[ResultMatrixIndex] = new vnl_matrix<Float>(
238 // this->GetSystemOrder(), this->GetSystemOrder() );
239 ( *( *m_Vectors )[resultVectorIndex] ) = ( *( *m_Vectors )[vectorIndex] );
240 ( *m_Vectors )[resultVectorIndex]->pre_multiply( *( ( *m_Matrices )[matrixIndex] ) );
241 }
242
ScaleMatrix(Float scale,unsigned int matrixIndex)243 void LinearSystemWrapperDenseVNL::ScaleMatrix(Float scale, unsigned int matrixIndex)
244 {
245 ( *( ( *m_Matrices )[matrixIndex] ) ) = ( *( ( *m_Matrices )[matrixIndex] ) ) * scale;
246 }
247
ScaleVector(Float scale,unsigned int vectorIndex)248 void LinearSystemWrapperDenseVNL::ScaleVector(Float scale, unsigned int vectorIndex)
249 {
250 ( *( *m_Vectors )[vectorIndex] ) = ( *( *m_Vectors )[vectorIndex] ) * scale;
251 }
252
ScaleSolution(Float scale,unsigned int solutionIndex)253 void LinearSystemWrapperDenseVNL::ScaleSolution(Float scale, unsigned int solutionIndex)
254 {
255 ( *( *m_Solutions )[solutionIndex] ) = ( *( *m_Solutions )[solutionIndex] ) * scale;
256 }
257
~LinearSystemWrapperDenseVNL()258 LinearSystemWrapperDenseVNL::~LinearSystemWrapperDenseVNL()
259 {
260 unsigned int i;
261 for( i = 0; i < m_NumberOfMatrices; i++ )
262 {
263 this->DestroyMatrix(i);
264 }
265 for( i = 0; i < m_NumberOfVectors; i++ )
266 {
267 this->DestroyVector(i);
268 }
269 for( i = 0; i < m_NumberOfSolutions; i++ )
270 {
271 this->DestroySolution(i);
272 }
273
274 delete m_Matrices;
275 delete m_Vectors;
276 delete m_Solutions;
277 }
278
279 } // end namespace fem
280 } // end namespace itk
281