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