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 26 /** \class VNLIterativeSparseSolverTraits 27 * \brief Generic interface for iterative sparse linear solver. 28 * 29 * This generic interface (common to several sparse solvers), allow to 30 * interchange solver solutions when dealing with sparse linear system. See 31 * itk::ParameterizationQuadEdgeMeshFilter for reference. 32 * 33 * It internally uses the VNL library to represent and deal with vectors 34 * (vnl_vector) and sparse matrices (vnl_sparse_matrix). The solver by itself 35 * is an iterative one see vnl_sparse_matrix_linear_system for more details on 36 * the method used. 37 * 38 * \ingroup ITKCommon 39 * 40 * \sa VNLSparseLUSolverTraits 41 */ 42 template< typename T = double > 43 class VNLIterativeSparseSolverTraits 44 { 45 public: 46 using ValueType = T; 47 using MatrixType = vnl_sparse_matrix< ValueType >; 48 using VectorType = vnl_vector< ValueType >; 49 using SolverType = vnl_lsqr; 50 51 /** \return false (it is not a direct solver, it is an iterative solver) */ IsDirectSolver()52 static bool IsDirectSolver() 53 { 54 return false; 55 } 56 57 /** \brief initialize a square sparse matrix of size iN x iN */ InitializeSparseMatrix(const unsigned int & iN)58 static MatrixType InitializeSparseMatrix(const unsigned int & iN) 59 { 60 return MatrixType(iN, iN); 61 } 62 63 /** \brief initialize a sparse matrix of size iRow x iCol */ InitializeSparseMatrix(const unsigned int & iRow,const unsigned int & iCol)64 static MatrixType InitializeSparseMatrix(const unsigned int & iRow, const unsigned int& iCol) 65 { 66 return MatrixType(iRow, iCol); 67 } 68 69 /** \brief initialize a vector of size iN */ InitializeVector(const unsigned int & iN)70 static VectorType InitializeVector(const unsigned int & iN) 71 { 72 return VectorType(iN); 73 } 74 75 /** \brief iA[iR][iC] = iV */ FillMatrix(MatrixType & iA,const unsigned int & iR,const unsigned int & iC,const ValueType & iV)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 81 /** \brief iA[iR][iC] += iV */ AddToMatrix(MatrixType & iA,const unsigned int & iR,const unsigned int & iC,const ValueType & iV)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 87 /** \brief Solve the linear system \f$ iA \cdot oX = iB \f$ */ Solve(const MatrixType & iA,const VectorType & iB,VectorType & oX)88 static bool Solve(const MatrixType & iA, const VectorType & iB, VectorType & oX) 89 { 90 using SparseLinearSystemType = vnl_sparse_matrix_linear_system< ValueType >; 91 SparseLinearSystemType system(iA, iB); 92 93 SolverType solver(system); 94 return solver.minimize(oX); 95 } 96 97 /** \brief Solve the linear systems: \f$ iA \cdot oX = iBx \f$, \f$ iA \cdot oY = iBy \f$, \f$ iA \cdot oZ = iBz \f$ */ Solve(const MatrixType & iA,const VectorType & iBx,const VectorType & iBy,const VectorType & iBz,VectorType & oX,VectorType & oY,VectorType & oZ)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 109 /** \brief Solve the linear systems: \f$ iA \cdot oX = iBx \f$, \f$ iA \cdot oY = iBy \f$ */ Solve(const MatrixType & iA,const VectorType & iBx,const VectorType & iBy,VectorType & oX,VectorType & oY)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 120 /** \brief Solve the linear systems: \f$ iA \cdot oX = iBx \f$ in N iterations */ Solve(const MatrixType & iA,const VectorType & iB,const long & iNbIter,VectorType & oX)121 static bool Solve(const MatrixType & iA, 122 const VectorType & iB, 123 const long & iNbIter, 124 VectorType & oX) 125 { 126 using SparseLinearSystemType = vnl_sparse_matrix_linear_system< ValueType >; 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