ITK  4.4.0
Insight Segmentation and Registration Toolkit
VNLIterativeSparseSolverTraits.h
Go to the documentation of this file.
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 #ifndef __VNLIterativeSparseSolverTraits_h
19 #define __VNLIterativeSparseSolverTraits_h
20 
21 #include "vnl/vnl_vector.h"
22 #include "vnl/vnl_sparse_matrix.h"
23 #include "vnl/vnl_sparse_matrix_linear_system.h"
24 #include "vnl/algo/vnl_lsqr.h"
25 
42 template< typename T = double >
44 {
45 public:
46  typedef T ValueType;
47  typedef vnl_sparse_matrix< ValueType > MatrixType;
48  typedef vnl_vector< ValueType > VectorType;
49  typedef vnl_lsqr SolverType;
50 
52  static bool IsDirectSolver()
53  {
54  return false;
55  }
56 
58  static MatrixType InitializeSparseMatrix(const unsigned int & iN)
59  {
60  return MatrixType(iN, iN);
61  }
62 
64  static MatrixType InitializeSparseMatrix(const unsigned int & iRow, const unsigned int& iCol)
65  {
66  return MatrixType(iRow, iCol);
67  }
68 
70  static VectorType InitializeVector(const unsigned int & iN)
71  {
72  return VectorType(iN);
73  }
74 
76  static void FillMatrix(MatrixType & iA, const unsigned int & iR, const unsigned int & iC, const ValueType & iV)
77  {
78  iA(iR, iC) = iV;
79  }
80 
82  static void AddToMatrix(MatrixType & iA, const unsigned int & iR, const unsigned int & iC, const ValueType & iV)
83  {
84  iA(iR, iC) += iV;
85  }
86 
88  static bool Solve(const MatrixType & iA, const VectorType & iB, VectorType & oX)
89  {
90  typedef vnl_sparse_matrix_linear_system< ValueType > SparseLinearSystemType;
91  SparseLinearSystemType system(iA, iB);
92 
93  SolverType solver(system);
94  return solver.minimize(oX);
95  }
96 
98  static bool Solve(const MatrixType & iA,
99  const VectorType & iBx, const VectorType & iBy, const VectorType & iBz,
100  VectorType & oX, VectorType & oY, VectorType & oZ )
101  {
102  bool result1 = Solve(iA, iBx, 100000, oX);
103  bool result2 = Solve(iA, iBy, 100000, oY);
104  bool result3 = Solve(iA, iBz, 100000, oZ);
105 
106  return ( result1 && result2 && result3 );
107  }
108 
110  static bool Solve(const MatrixType & iA,
111  const VectorType & iBx, const VectorType & iBy,
112  VectorType & oX, VectorType & oY)
113  {
114  bool result1 = Solve(iA, iBx, oX);
115  bool result2 = Solve(iA, iBy, oY);
116 
117  return ( result1 && result2 );
118  }
119 
121  static bool Solve(const MatrixType & iA,
122  const VectorType & iB,
123  const long & iNbIter,
124  VectorType & oX)
125  {
126  typedef vnl_sparse_matrix_linear_system< ValueType > SparseLinearSystemType;
127  SparseLinearSystemType system(iA, iB);
128 
129  SolverType solver(system);
130  solver.set_max_iterations(iNbIter);
131  return solver.minimize(oX);
132  }
133 };
134 
135 #endif
136