IncompleteLUT.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_INCOMPLETE_LUT_H
11 #define EIGEN_INCOMPLETE_LUT_H
12 
13 namespace Eigen {
14 
43 template <typename _Scalar>
44 class IncompleteLUT : internal::noncopyable
45 {
46  typedef _Scalar Scalar;
47  typedef typename NumTraits<Scalar>::Real RealScalar;
51  typedef typename FactorType::Index Index;
52 
53  public:
55 
57  : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
58  m_analysisIsOk(false), m_factorizationIsOk(false), m_isInitialized(false)
59  {}
60 
61  template<typename MatrixType>
62  IncompleteLUT(const MatrixType& mat, RealScalar droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
63  : m_droptol(droptol),m_fillfactor(fillfactor),
64  m_analysisIsOk(false),m_factorizationIsOk(false),m_isInitialized(false)
65  {
66  eigen_assert(fillfactor != 0);
67  compute(mat);
68  }
69 
70  Index rows() const { return m_lu.rows(); }
71 
72  Index cols() const { return m_lu.cols(); }
73 
80  {
81  eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
82  return m_info;
83  }
84 
85  template<typename MatrixType>
86  void analyzePattern(const MatrixType& amat);
87 
88  template<typename MatrixType>
89  void factorize(const MatrixType& amat);
90 
96  template<typename MatrixType>
98  {
99  analyzePattern(amat);
100  factorize(amat);
101  eigen_assert(m_factorizationIsOk == true);
102  m_isInitialized = true;
103  return *this;
104  }
105 
106  void setDroptol(RealScalar droptol);
107  void setFillfactor(int fillfactor);
108 
109  template<typename Rhs, typename Dest>
110  void _solve(const Rhs& b, Dest& x) const
111  {
112  x = m_Pinv * b;
113  x = m_lu.template triangularView<UnitLower>().solve(x);
114  x = m_lu.template triangularView<Upper>().solve(x);
115  x = m_P * x;
116  }
117 
118  template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs>
119  solve(const MatrixBase<Rhs>& b) const
120  {
121  eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
122  eigen_assert(cols()==b.rows()
123  && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
124  return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived());
125  }
126 
127 protected:
128 
129  template <typename VectorV, typename VectorI>
130  int QuickSplit(VectorV &row, VectorI &ind, int ncut);
131 
132 
134  struct keep_diag {
135  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
136  {
137  return row!=col;
138  }
139  };
140 
141 protected:
142 
143  FactorType m_lu;
144  RealScalar m_droptol;
145  int m_fillfactor;
146  bool m_analysisIsOk;
147  bool m_factorizationIsOk;
148  bool m_isInitialized;
149  ComputationInfo m_info;
150  PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation
151  PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // Inverse permutation
152 };
153 
158 template<typename Scalar>
159 void IncompleteLUT<Scalar>::setDroptol(RealScalar droptol)
160 {
161  this->m_droptol = droptol;
162 }
163 
168 template<typename Scalar>
170 {
171  this->m_fillfactor = fillfactor;
172 }
173 
174 
184 template <typename Scalar>
185 template <typename VectorV, typename VectorI>
186 int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut)
187 {
188  using std::swap;
189  int mid;
190  int n = row.size(); /* length of the vector */
191  int first, last ;
192 
193  ncut--; /* to fit the zero-based indices */
194  first = 0;
195  last = n-1;
196  if (ncut < first || ncut > last ) return 0;
197 
198  do {
199  mid = first;
200  RealScalar abskey = std::abs(row(mid));
201  for (int j = first + 1; j <= last; j++) {
202  if ( std::abs(row(j)) > abskey) {
203  ++mid;
204  swap(row(mid), row(j));
205  swap(ind(mid), ind(j));
206  }
207  }
208  /* Interchange for the pivot element */
209  swap(row(mid), row(first));
210  swap(ind(mid), ind(first));
211 
212  if (mid > ncut) last = mid - 1;
213  else if (mid < ncut ) first = mid + 1;
214  } while (mid != ncut );
215 
216  return 0; /* mid is equal to ncut */
217 }
218 
219 template <typename Scalar>
220 template<typename _MatrixType>
221 void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
222 {
223  // Compute the Fill-reducing permutation
225  SparseMatrix<Scalar,ColMajor, Index> mat2 = amat.transpose();
226  // Symmetrize the pattern
227  // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
228  // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered...
229  SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1;
230  AtA.prune(keep_diag());
231  internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering...
232 
233  m_Pinv = m_P.inverse(); // ... and the inverse permutation
234 
235  m_analysisIsOk = true;
236 }
237 
238 template <typename Scalar>
239 template<typename _MatrixType>
240 void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
241 {
242  using std::sqrt;
243  using std::swap;
244  using std::abs;
245 
246  eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
247  int n = amat.cols(); // Size of the matrix
248  m_lu.resize(n,n);
249  // Declare Working vectors and variables
250  Vector u(n) ; // real values of the row -- maximum size is n --
251  VectorXi ju(n); // column position of the values in u -- maximum size is n
252  VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
253 
254  // Apply the fill-reducing permutation
255  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
256  SparseMatrix<Scalar,RowMajor, Index> mat;
257  mat = amat.twistedBy(m_Pinv);
258 
259  // Initialization
260  jr.fill(-1);
261  ju.fill(0);
262  u.fill(0);
263 
264  // number of largest elements to keep in each row:
265  int fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
266  if (fill_in > n) fill_in = n;
267 
268  // number of largest nonzero elements to keep in the L and the U part of the current row:
269  int nnzL = fill_in/2;
270  int nnzU = nnzL;
271  m_lu.reserve(n * (nnzL + nnzU + 1));
272 
273  // global loop over the rows of the sparse matrix
274  for (int ii = 0; ii < n; ii++)
275  {
276  // 1 - copy the lower and the upper part of the row i of mat in the working vector u
277 
278  int sizeu = 1; // number of nonzero elements in the upper part of the current row
279  int sizel = 0; // number of nonzero elements in the lower part of the current row
280  ju(ii) = ii;
281  u(ii) = 0;
282  jr(ii) = ii;
283  RealScalar rownorm = 0;
284 
285  typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
286  for (; j_it; ++j_it)
287  {
288  int k = j_it.index();
289  if (k < ii)
290  {
291  // copy the lower part
292  ju(sizel) = k;
293  u(sizel) = j_it.value();
294  jr(k) = sizel;
295  ++sizel;
296  }
297  else if (k == ii)
298  {
299  u(ii) = j_it.value();
300  }
301  else
302  {
303  // copy the upper part
304  int jpos = ii + sizeu;
305  ju(jpos) = k;
306  u(jpos) = j_it.value();
307  jr(k) = jpos;
308  ++sizeu;
309  }
310  rownorm += internal::abs2(j_it.value());
311  }
312 
313  // 2 - detect possible zero row
314  if(rownorm==0)
315  {
316  m_info = NumericalIssue;
317  return;
318  }
319  // Take the 2-norm of the current row as a relative tolerance
320  rownorm = sqrt(rownorm);
321 
322  // 3 - eliminate the previous nonzero rows
323  int jj = 0;
324  int len = 0;
325  while (jj < sizel)
326  {
327  // In order to eliminate in the correct order,
328  // we must select first the smallest column index among ju(jj:sizel)
329  int k;
330  int minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
331  k += jj;
332  if (minrow != ju(jj))
333  {
334  // swap the two locations
335  int j = ju(jj);
336  swap(ju(jj), ju(k));
337  jr(minrow) = jj; jr(j) = k;
338  swap(u(jj), u(k));
339  }
340  // Reset this location
341  jr(minrow) = -1;
342 
343  // Start elimination
344  typename FactorType::InnerIterator ki_it(m_lu, minrow);
345  while (ki_it && ki_it.index() < minrow) ++ki_it;
346  eigen_internal_assert(ki_it && ki_it.col()==minrow);
347  Scalar fact = u(jj) / ki_it.value();
348 
349  // drop too small elements
350  if(abs(fact) <= m_droptol)
351  {
352  jj++;
353  continue;
354  }
355 
356  // linear combination of the current row ii and the row minrow
357  ++ki_it;
358  for (; ki_it; ++ki_it)
359  {
360  Scalar prod = fact * ki_it.value();
361  int j = ki_it.index();
362  int jpos = jr(j);
363  if (jpos == -1) // fill-in element
364  {
365  int newpos;
366  if (j >= ii) // dealing with the upper part
367  {
368  newpos = ii + sizeu;
369  sizeu++;
370  eigen_internal_assert(sizeu<=n);
371  }
372  else // dealing with the lower part
373  {
374  newpos = sizel;
375  sizel++;
376  eigen_internal_assert(sizel<=ii);
377  }
378  ju(newpos) = j;
379  u(newpos) = -prod;
380  jr(j) = newpos;
381  }
382  else
383  u(jpos) -= prod;
384  }
385  // store the pivot element
386  u(len) = fact;
387  ju(len) = minrow;
388  ++len;
389 
390  jj++;
391  } // end of the elimination on the row ii
392 
393  // reset the upper part of the pointer jr to zero
394  for(int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
395 
396  // 4 - partially sort and insert the elements in the m_lu matrix
397 
398  // sort the L-part of the row
399  sizel = len;
400  len = (std::min)(sizel, nnzL);
401  typename Vector::SegmentReturnType ul(u.segment(0, sizel));
402  typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
403  QuickSplit(ul, jul, len);
404 
405  // store the largest m_fill elements of the L part
406  m_lu.startVec(ii);
407  for(int k = 0; k < len; k++)
408  m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
409 
410  // store the diagonal element
411  // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization)
412  if (u(ii) == Scalar(0))
413  u(ii) = sqrt(m_droptol) * rownorm;
414  m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
415 
416  // sort the U-part of the row
417  // apply the dropping rule first
418  len = 0;
419  for(int k = 1; k < sizeu; k++)
420  {
421  if(abs(u(ii+k)) > m_droptol * rownorm )
422  {
423  ++len;
424  u(ii + len) = u(ii + k);
425  ju(ii + len) = ju(ii + k);
426  }
427  }
428  sizeu = len + 1; // +1 to take into account the diagonal element
429  len = (std::min)(sizeu, nnzU);
430  typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
431  typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
432  QuickSplit(uu, juu, len);
433 
434  // store the largest elements of the U part
435  for(int k = ii + 1; k < ii + len; k++)
436  m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
437  }
438 
439  m_lu.finalize();
440  m_lu.makeCompressed();
441 
442  m_factorizationIsOk = true;
443  m_info = Success;
444 }
445 
446 namespace internal {
447 
448 template<typename _MatrixType, typename Rhs>
449 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
450  : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
451 {
452  typedef IncompleteLUT<_MatrixType> Dec;
453  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
454 
455  template<typename Dest> void evalTo(Dest& dst) const
456  {
457  dec()._solve(rhs(),dst);
458  }
459 };
460 
461 } // end namespace internal
462 
463 } // end namespace Eigen
464 
465 #endif // EIGEN_INCOMPLETE_LUT_H
466