IncompleteLUT.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_INCOMPLETE_LUT_H
26 #define EIGEN_INCOMPLETE_LUT_H
27 
28 namespace Eigen {
29 
53 template <typename _Scalar>
54 class IncompleteLUT : internal::noncopyable
55 {
56  typedef _Scalar Scalar;
57  typedef typename NumTraits<Scalar>::Real RealScalar;
61  typedef typename FactorType::Index Index;
62 
63  public:
65 
67  : m_droptol(NumTraits<Scalar>::dummy_precision()), m_fillfactor(10),
69  {}
70 
71  template<typename MatrixType>
72  IncompleteLUT(const MatrixType& mat, RealScalar droptol=NumTraits<Scalar>::dummy_precision(), int fillfactor = 10)
73  : m_droptol(droptol),m_fillfactor(fillfactor),
75  {
76  eigen_assert(fillfactor != 0);
77  compute(mat);
78  }
79 
80  Index rows() const { return m_lu.rows(); }
81 
82  Index cols() const { return m_lu.cols(); }
83 
90  {
91  eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
92  return m_info;
93  }
94 
95  template<typename MatrixType>
96  void analyzePattern(const MatrixType& amat);
97 
98  template<typename MatrixType>
99  void factorize(const MatrixType& amat);
100 
106  template<typename MatrixType>
108  {
109  analyzePattern(amat);
110  factorize(amat);
112  m_isInitialized = true;
113  return *this;
114  }
115 
116  void setDroptol(RealScalar droptol);
117  void setFillfactor(int fillfactor);
118 
119  template<typename Rhs, typename Dest>
120  void _solve(const Rhs& b, Dest& x) const
121  {
122  x = m_Pinv * b;
123  x = m_lu.template triangularView<UnitLower>().solve(x);
124  x = m_lu.template triangularView<Upper>().solve(x);
125  x = m_P * x;
126  }
127 
128  template<typename Rhs> inline const internal::solve_retval<IncompleteLUT, Rhs>
129  solve(const MatrixBase<Rhs>& b) const
130  {
131  eigen_assert(m_isInitialized && "IncompleteLUT is not initialized.");
132  eigen_assert(cols()==b.rows()
133  && "IncompleteLUT::solve(): invalid number of rows of the right hand side matrix b");
134  return internal::solve_retval<IncompleteLUT, Rhs>(*this, b.derived());
135  }
136 
137 protected:
138 
139  template <typename VectorV, typename VectorI>
140  int QuickSplit(VectorV &row, VectorI &ind, int ncut);
141 
142 
144  struct keep_diag {
145  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
146  {
147  return row!=col;
148  }
149  };
150 
151 protected:
152 
154  RealScalar m_droptol;
160  PermutationMatrix<Dynamic,Dynamic,Index> m_P; // Fill-reducing permutation
162 };
163 
168 template<typename Scalar>
169 void IncompleteLUT<Scalar>::setDroptol(RealScalar droptol)
170 {
171  this->m_droptol = droptol;
172 }
173 
178 template<typename Scalar>
180 {
181  this->m_fillfactor = fillfactor;
182 }
183 
184 
194 template <typename Scalar>
195 template <typename VectorV, typename VectorI>
196 int IncompleteLUT<Scalar>::QuickSplit(VectorV &row, VectorI &ind, int ncut)
197 {
198  using std::swap;
199  int mid;
200  int n = row.size(); /* length of the vector */
201  int first, last ;
202 
203  ncut--; /* to fit the zero-based indices */
204  first = 0;
205  last = n-1;
206  if (ncut < first || ncut > last ) return 0;
207 
208  do {
209  mid = first;
210  RealScalar abskey = std::abs(row(mid));
211  for (int j = first + 1; j <= last; j++) {
212  if ( std::abs(row(j)) > abskey) {
213  ++mid;
214  swap(row(mid), row(j));
215  swap(ind(mid), ind(j));
216  }
217  }
218  /* Interchange for the pivot element */
219  swap(row(mid), row(first));
220  swap(ind(mid), ind(first));
221 
222  if (mid > ncut) last = mid - 1;
223  else if (mid < ncut ) first = mid + 1;
224  } while (mid != ncut );
225 
226  return 0; /* mid is equal to ncut */
227 }
228 
229 template <typename Scalar>
230 template<typename _MatrixType>
231 void IncompleteLUT<Scalar>::analyzePattern(const _MatrixType& amat)
232 {
233  // Compute the Fill-reducing permutation
236  // Symmetrize the pattern
237  // FIXME for a matrix with nearly symmetric pattern, mat2+mat1 is the appropriate choice.
238  // on the other hand for a really non-symmetric pattern, mat2*mat1 should be prefered...
239  SparseMatrix<Scalar,ColMajor, Index> AtA = mat2 + mat1;
240  AtA.prune(keep_diag());
241  internal::minimum_degree_ordering<Scalar, Index>(AtA, m_P); // Then compute the AMD ordering...
242 
243  m_Pinv = m_P.inverse(); // ... and the inverse permutation
244 
245  m_analysisIsOk = true;
246 }
247 
248 template <typename Scalar>
249 template<typename _MatrixType>
250 void IncompleteLUT<Scalar>::factorize(const _MatrixType& amat)
251 {
252  using std::sqrt;
253  using std::swap;
254  using std::abs;
255 
256  eigen_assert((amat.rows() == amat.cols()) && "The factorization should be done on a square matrix");
257  int n = amat.cols(); // Size of the matrix
258  m_lu.resize(n,n);
259  // Declare Working vectors and variables
260  Vector u(n) ; // real values of the row -- maximum size is n --
261  VectorXi ju(n); // column position of the values in u -- maximum size is n
262  VectorXi jr(n); // Indicate the position of the nonzero elements in the vector u -- A zero location is indicated by -1
263 
264  // Apply the fill-reducing permutation
265  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
267  mat = amat.twistedBy(m_Pinv);
268 
269  // Initialization
270  jr.fill(-1);
271  ju.fill(0);
272  u.fill(0);
273 
274  // number of largest elements to keep in each row:
275  int fill_in = static_cast<int> (amat.nonZeros()*m_fillfactor)/n+1;
276  if (fill_in > n) fill_in = n;
277 
278  // number of largest nonzero elements to keep in the L and the U part of the current row:
279  int nnzL = fill_in/2;
280  int nnzU = nnzL;
281  m_lu.reserve(n * (nnzL + nnzU + 1));
282 
283  // global loop over the rows of the sparse matrix
284  for (int ii = 0; ii < n; ii++)
285  {
286  // 1 - copy the lower and the upper part of the row i of mat in the working vector u
287 
288  int sizeu = 1; // number of nonzero elements in the upper part of the current row
289  int sizel = 0; // number of nonzero elements in the lower part of the current row
290  ju(ii) = ii;
291  u(ii) = 0;
292  jr(ii) = ii;
293  RealScalar rownorm = 0;
294 
295  typename FactorType::InnerIterator j_it(mat, ii); // Iterate through the current row ii
296  for (; j_it; ++j_it)
297  {
298  int k = j_it.index();
299  if (k < ii)
300  {
301  // copy the lower part
302  ju(sizel) = k;
303  u(sizel) = j_it.value();
304  jr(k) = sizel;
305  ++sizel;
306  }
307  else if (k == ii)
308  {
309  u(ii) = j_it.value();
310  }
311  else
312  {
313  // copy the upper part
314  int jpos = ii + sizeu;
315  ju(jpos) = k;
316  u(jpos) = j_it.value();
317  jr(k) = jpos;
318  ++sizeu;
319  }
320  rownorm += internal::abs2(j_it.value());
321  }
322 
323  // 2 - detect possible zero row
324  if(rownorm==0)
325  {
326  m_info = NumericalIssue;
327  return;
328  }
329  // Take the 2-norm of the current row as a relative tolerance
330  rownorm = sqrt(rownorm);
331 
332  // 3 - eliminate the previous nonzero rows
333  int jj = 0;
334  int len = 0;
335  while (jj < sizel)
336  {
337  // In order to eliminate in the correct order,
338  // we must select first the smallest column index among ju(jj:sizel)
339  int k;
340  int minrow = ju.segment(jj,sizel-jj).minCoeff(&k); // k is relative to the segment
341  k += jj;
342  if (minrow != ju(jj))
343  {
344  // swap the two locations
345  int j = ju(jj);
346  swap(ju(jj), ju(k));
347  jr(minrow) = jj; jr(j) = k;
348  swap(u(jj), u(k));
349  }
350  // Reset this location
351  jr(minrow) = -1;
352 
353  // Start elimination
354  typename FactorType::InnerIterator ki_it(m_lu, minrow);
355  while (ki_it && ki_it.index() < minrow) ++ki_it;
356  eigen_internal_assert(ki_it && ki_it.col()==minrow);
357  Scalar fact = u(jj) / ki_it.value();
358 
359  // drop too small elements
360  if(abs(fact) <= m_droptol)
361  {
362  jj++;
363  continue;
364  }
365 
366  // linear combination of the current row ii and the row minrow
367  ++ki_it;
368  for (; ki_it; ++ki_it)
369  {
370  Scalar prod = fact * ki_it.value();
371  int j = ki_it.index();
372  int jpos = jr(j);
373  if (jpos == -1) // fill-in element
374  {
375  int newpos;
376  if (j >= ii) // dealing with the upper part
377  {
378  newpos = ii + sizeu;
379  sizeu++;
380  eigen_internal_assert(sizeu<=n);
381  }
382  else // dealing with the lower part
383  {
384  newpos = sizel;
385  sizel++;
386  eigen_internal_assert(sizel<=ii);
387  }
388  ju(newpos) = j;
389  u(newpos) = -prod;
390  jr(j) = newpos;
391  }
392  else
393  u(jpos) -= prod;
394  }
395  // store the pivot element
396  u(len) = fact;
397  ju(len) = minrow;
398  ++len;
399 
400  jj++;
401  } // end of the elimination on the row ii
402 
403  // reset the upper part of the pointer jr to zero
404  for(int k = 0; k <sizeu; k++) jr(ju(ii+k)) = -1;
405 
406  // 4 - partially sort and insert the elements in the m_lu matrix
407 
408  // sort the L-part of the row
409  sizel = len;
410  len = (std::min)(sizel, nnzL);
411  typename Vector::SegmentReturnType ul(u.segment(0, sizel));
412  typename VectorXi::SegmentReturnType jul(ju.segment(0, sizel));
413  QuickSplit(ul, jul, len);
414 
415  // store the largest m_fill elements of the L part
416  m_lu.startVec(ii);
417  for(int k = 0; k < len; k++)
418  m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
419 
420  // store the diagonal element
421  // apply a shifting rule to avoid zero pivots (we are doing an incomplete factorization)
422  if (u(ii) == Scalar(0))
423  u(ii) = sqrt(m_droptol) * rownorm;
424  m_lu.insertBackByOuterInnerUnordered(ii, ii) = u(ii);
425 
426  // sort the U-part of the row
427  // apply the dropping rule first
428  len = 0;
429  for(int k = 1; k < sizeu; k++)
430  {
431  if(abs(u(ii+k)) > m_droptol * rownorm )
432  {
433  ++len;
434  u(ii + len) = u(ii + k);
435  ju(ii + len) = ju(ii + k);
436  }
437  }
438  sizeu = len + 1; // +1 to take into account the diagonal element
439  len = (std::min)(sizeu, nnzU);
440  typename Vector::SegmentReturnType uu(u.segment(ii+1, sizeu-1));
441  typename VectorXi::SegmentReturnType juu(ju.segment(ii+1, sizeu-1));
442  QuickSplit(uu, juu, len);
443 
444  // store the largest elements of the U part
445  for(int k = ii + 1; k < ii + len; k++)
446  m_lu.insertBackByOuterInnerUnordered(ii,ju(k)) = u(k);
447  }
448 
449  m_lu.finalize();
450  m_lu.makeCompressed();
451 
452  m_factorizationIsOk = true;
453  m_info = Success;
454 }
455 
456 namespace internal {
457 
458 template<typename _MatrixType, typename Rhs>
459 struct solve_retval<IncompleteLUT<_MatrixType>, Rhs>
460  : solve_retval_base<IncompleteLUT<_MatrixType>, Rhs>
461 {
462  typedef IncompleteLUT<_MatrixType> Dec;
463  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
464 
465  template<typename Dest> void evalTo(Dest& dst) const
466  {
467  dec()._solve(rhs(),dst);
468  }
469 };
470 
471 } // end namespace internal
472 
473 } // end namespace Eigen
474 
475 #endif // EIGEN_INCOMPLETE_LUT_H
476