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 __VNLSparseLUSolverTraits_h 19 #define __VNLSparseLUSolverTraits_h 20 21 #include "vnl/vnl_vector.h" 22 #include "vnl/vnl_sparse_matrix.h" 23 #include "vnl/algo/vnl_sparse_lu.h" 24 25 /** \class VNLSparseLUSolverTraits 26 * \brief Generic interface for sparse LU solver. 27 * 28 * This generic interface (common to several sparse solvers), allow to 29 * interchange solver solutions when dealing with sparse linear systems. See 30 * itk::ParameterizationQuadEdgeMeshFilter for reference. 31 * 32 * It internally uses the VNL library to represent and deal with vectors 33 * (vnl_vector) and sparse matrices (vnl_sparse_matrix). The solver by itself 34 * is made of a sparse LU decomposition followed by solving upper triangular 35 * system. see vnl_sparse_lu for more details on the method used. 36 * 37 * \ingroup ITKCommon 38 * 39 * \sa VNLIterativeSparseSolverTraits:w 40 */ 41 template< typename T = double > 42 class VNLSparseLUSolverTraits 43 { 44 public: 45 using ValueType = T; 46 using MatrixType = vnl_sparse_matrix< ValueType >; 47 using VectorType = vnl_vector< ValueType >; 48 using SolverType = vnl_sparse_lu; 49 50 /** \return false (it is not a direct solver, it is an iterative solver) */ IsDirectSolver()51 static bool IsDirectSolver() 52 { 53 return true; 54 } 55 56 /** \brief initialize a square sparse matrix of size iN x iN */ InitializeSparseMatrix(const unsigned int & iN)57 static MatrixType InitializeSparseMatrix(const unsigned int & iN) 58 { 59 return MatrixType(iN, iN); 60 } 61 62 /** \brief initialize a sparse matrix of size iRow x iCol */ InitializeSparseMatrix(const unsigned int & iRow,const unsigned int & iCol)63 static MatrixType InitializeSparseMatrix(const unsigned int & iRow, const unsigned int& iCol) 64 { 65 return MatrixType(iRow, iCol); 66 } 67 68 /** \brief initialize a vector of size iN */ InitializeVector(const unsigned int & iN)69 static VectorType InitializeVector(const unsigned int & iN) 70 { 71 return VectorType(iN); 72 } 73 74 /** \brief iA[iR][iC] = iV */ FillMatrix(MatrixType & iA,const unsigned int & iR,const unsigned int & iC,const ValueType & iV)75 static void FillMatrix(MatrixType & iA, const unsigned int & iR, const unsigned int & iC, const ValueType & iV) 76 { 77 iA(iR, iC) = iV; 78 } 79 80 /** \brief iA[iR][iC] += iV */ AddToMatrix(MatrixType & iA,const unsigned int & iR,const unsigned int & iC,const ValueType & iV)81 static void AddToMatrix(MatrixType & iA, const unsigned int & iR, const unsigned int & iC, const ValueType & iV) 82 { 83 iA(iR, iC) += iV; 84 } 85 86 /** \brief Solve the linear system \f$ iA \cdot oX = iB \f$ */ Solve(const MatrixType & iA,const VectorType & iB,VectorType & oX)87 static bool Solve(const MatrixType & iA, const VectorType & iB, VectorType & oX) 88 { 89 SolverType solver( iA ); 90 Solve( solver, iB, oX ); 91 92 return true; 93 } 94 95 /** \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)96 static bool Solve(const MatrixType & iA, 97 const VectorType & iBx, const VectorType & iBy, const VectorType & iBz, 98 VectorType & oX, VectorType & oY, VectorType & oZ ) 99 { 100 SolverType solver( iA ); 101 Solve( solver, iBx, iBy, iBz, oX, oY, oZ ); 102 103 return true; 104 } 105 106 /** \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)107 static bool Solve(const MatrixType & iA, 108 const VectorType & iBx, const VectorType & iBy, 109 VectorType & oX, VectorType & oY) 110 { 111 SolverType solver( iA ); 112 Solve( solver, iBx, iBy, oX, oY ); 113 114 return true; 115 } 116 117 /** \brief Solve the linear system \f$ iA \cdot oX = iB \f$ factoring the internal matrix if needed */ Solve(SolverType & solver,const VectorType & iB,VectorType & oX)118 static void Solve( SolverType & solver, const VectorType & iB, VectorType & oX ) 119 { 120 oX = solver.solve( iB ); 121 } 122 123 /** \brief Solve the linear systems: \f$ iA \cdot oX = iBx \f$, \f$ iA \cdot oY = iBy \f$, \f$ iA \cdot oZ = iBz \f$ factoring the internal matrix if needed */ Solve(SolverType & solver,const VectorType & iBx,const VectorType & iBy,const VectorType & iBz,VectorType & oX,VectorType & oY,VectorType & oZ)124 static void Solve( SolverType & solver, 125 const VectorType & iBx, const VectorType & iBy, const VectorType & iBz, 126 VectorType & oX, VectorType & oY, VectorType & oZ ) 127 { 128 oX = solver.solve( iBx ); 129 oY = solver.solve( iBy ); 130 oZ = solver.solve( iBz ); 131 } 132 133 /** \brief Solve the linear systems: \f$ iA \cdot oX = iBx \f$, \f$ iA \cdot oY = iBy \f$ factoring the internal matrix if needed */ Solve(SolverType & solver,const VectorType & iBx,const VectorType & iBy,VectorType & oX,VectorType & oY)134 static void Solve( SolverType & solver, 135 const VectorType & iBx, const VectorType & iBy, 136 VectorType & oX, VectorType & oY ) 137 { 138 oX = solver.solve( iBx ); 139 oY = solver.solve( iBy ); 140 } 141 }; 142 143 #endif 144