BKLDLT.h
Go to the documentation of this file.
1 // Copyright (C) 2019-2025 Yixuan Qiu <yixuan.qiu@cos.name>
2 //
3 // This Source Code Form is subject to the terms of the Mozilla
4 // Public License v. 2.0. If a copy of the MPL was not distributed
5 // with this file, You can obtain one at https://mozilla.org/MPL/2.0/.
6 
7 #ifndef SPECTRA_BK_LDLT_H
8 #define SPECTRA_BK_LDLT_H
9 
10 #include <Eigen/Core>
11 #include <vector>
12 #include <stdexcept>
13 #include <type_traits> // std::is_same
14 
15 #include "../Util/CompInfo.h"
16 
17 namespace Spectra {
18 
19 // We need a generic conj() function for both real and complex values,
20 // and hope that conj(x) == x if x is real-valued. However, in STL,
21 // conj(x) == std::complex(x, 0) for such cases, meaning that the
22 // return value type is not necessarily the same as x. To avoid this
23 // inconvenience, we define a simple class that does this task
24 //
25 // Similarly, define a real(x) function that returns x itself if
26 // x is real-valued, and returns std::complex(x, 0) if x is complex-valued
27 template <typename Scalar>
28 struct ScalarOp
29 {
30  static Scalar conj(const Scalar& x)
31  {
32  return x;
33  }
34 
35  static Scalar real(const Scalar& x)
36  {
37  return x;
38  }
39 };
40 // Specialization for complex values
41 template <typename RealScalar>
42 struct ScalarOp<std::complex<RealScalar>>
43 {
44  static std::complex<RealScalar> conj(const std::complex<RealScalar>& x)
45  {
46  using std::conj;
47  return conj(x);
48  }
49 
50  static std::complex<RealScalar> real(const std::complex<RealScalar>& x)
51  {
52  return std::complex<RealScalar>(x.real(), RealScalar(0));
53  }
54 };
55 
56 // Bunch-Kaufman LDLT decomposition
57 // References:
58 // 1. Bunch, J. R., & Kaufman, L. (1977). Some stable methods for calculating inertia and solving symmetric linear systems.
59 // Mathematics of computation, 31(137), 163-179.
60 // 2. Golub, G. H., & Van Loan, C. F. (2012). Matrix computations (Vol. 3). JHU press. Section 4.4.
61 // 3. Bunch-Parlett diagonal pivoting <http://oz.nthu.edu.tw/~d947207/Chap13_GE3.ppt>
62 // 4. Ashcraft, C., Grimes, R. G., & Lewis, J. G. (1998). Accurate symmetric indefinite linear equation solvers.
63 // SIAM Journal on Matrix Analysis and Applications, 20(2), 513-561.
64 template <typename Scalar = double>
65 class BKLDLT
66 {
67 private:
68  // The real part type of the matrix element
77 
79  Vector m_data; // storage for a lower-triangular matrix
80  std::vector<Scalar*> m_colptr; // pointers to columns
81  IntVector m_perm; // [-2, -1, 3, 1, 4, 5]: 0 <-> 2, 1 <-> 1, 2 <-> 3, 3 <-> 1, 4 <-> 4, 5 <-> 5
82  std::vector<std::pair<Index, Index>> m_permc; // compressed version of m_perm: [(0, 2), (2, 3), (3, 1)]
83 
84  bool m_computed;
86 
87  // Access to elements
88  // Pointer to the k-th column
89  Scalar* col_pointer(Index k) { return m_colptr[k]; }
90  // A[i, j] -> m_colptr[j][i - j], i >= j
91  Scalar& coeff(Index i, Index j) { return m_colptr[j][i - j]; }
92  const Scalar& coeff(Index i, Index j) const { return m_colptr[j][i - j]; }
93  // A[i, i] -> m_colptr[i][0]
94  Scalar& diag_coeff(Index i) { return m_colptr[i][0]; }
95  const Scalar& diag_coeff(Index i) const { return m_colptr[i][0]; }
96 
97  // Compute column pointers
99  {
100  m_colptr.clear();
101  m_colptr.reserve(m_n);
102  Scalar* head = m_data.data();
103 
104  for (Index i = 0; i < m_n; i++)
105  {
106  m_colptr.push_back(head);
107  head += (m_n - i);
108  }
109  }
110 
111  // Copy mat - shift * I to m_data
112  template <typename Derived>
113  void copy_data(const Eigen::MatrixBase<Derived>& mat, int uplo, const RealScalar& shift)
114  {
115  // If mat is an expression, first evaluate it into a temporary object
116  // This can be achieved by assigning mat to a const Eigen::Ref<const Matrix>&
117  // If mat is a plain object, no temporary object is created
119 
120  // Efficient copying for column-major matrices with lower triangular part
121  if ((!Derived::PlainObject::IsRowMajor) && uplo == Eigen::Lower)
122  {
123  for (Index j = 0; j < m_n; j++)
124  {
125  const Scalar* begin = &src.coeffRef(j, j);
126  const Index len = m_n - j;
127  std::copy(begin, begin + len, col_pointer(j));
128  diag_coeff(j) -= Scalar(shift);
129  }
130  }
131  else
132  {
133  Scalar* dest = m_data.data();
134  for (Index j = 0; j < m_n; j++)
135  {
136  for (Index i = j; i < m_n; i++, dest++)
137  {
138  if (uplo == Eigen::Lower)
139  *dest = src.coeff(i, j);
140  else
141  *dest = ScalarOp<Scalar>::conj(src.coeff(j, i));
142  }
143  diag_coeff(j) -= Scalar(shift);
144  }
145  }
146  }
147 
148  // Compute compressed permutations
150  {
151  for (Index i = 0; i < m_n; i++)
152  {
153  // Recover the permutation action
154  const Index perm = (m_perm[i] >= 0) ? (m_perm[i]) : (-m_perm[i] - 1);
155  if (perm != i)
156  m_permc.push_back(std::make_pair(i, perm));
157  }
158  }
159 
160  // Working on the A[k:end, k:end] submatrix
161  // Exchange k <-> r
162  // Assume r >= k
164  {
165  m_perm[k] = r;
166 
167  // No permutation
168  if (k == r)
169  return;
170 
171  // A[k, k] <-> A[r, r]
173 
174  // A[(r+1):end, k] <-> A[(r+1):end, r]
175  std::swap_ranges(&coeff(r + 1, k), col_pointer(k + 1), &coeff(r + 1, r));
176 
177  // A[(k+1):(r-1), k] <-> A[r, (k+1):(r-1)]
178  // Note: for Hermitian matrices, also need to do conjugate
179  Scalar* src = &coeff(k + 1, k);
181  {
182  // Simple swapping for real values
183  for (Index j = k + 1; j < r; j++, src++)
184  {
185  std::swap(*src, coeff(r, j));
186  }
187  }
188  else
189  {
190  // For complex values
191  for (Index j = k + 1; j < r; j++, src++)
192  {
193  const Scalar src_conj = ScalarOp<Scalar>::conj(*src);
194  *src = ScalarOp<Scalar>::conj(coeff(r, j));
195  coeff(r, j) = src_conj;
196  }
197  }
198 
199  // A[r, k] <- Conj(A[r, k])
201  {
202  coeff(r, k) = ScalarOp<Scalar>::conj(coeff(r, k));
203  }
204  }
205 
206  // Working on the A[k:end, k:end] submatrix
207  // Exchange [k+1, k] <-> [r, p]
208  // Assume p >= k, r >= k+1
210  {
211  pivoting_1x1(k, p);
212  pivoting_1x1(k + 1, r);
213 
214  // A[k+1, k] <-> A[r, k]
215  std::swap(coeff(k + 1, k), coeff(r, k));
216 
217  // Use negative signs to indicate a 2x2 block
218  // Also minus one to distinguish a negative zero from a positive zero
219  m_perm[k] = -m_perm[k] - 1;
220  m_perm[k + 1] = -m_perm[k + 1] - 1;
221  }
222 
223  // A[r1, c1:c2] <-> A[r2, c1:c2]
224  // Assume r2 >= r1 > c2 >= c1
226  {
227  if (r1 == r2)
228  return;
229 
230  for (Index j = c1; j <= c2; j++)
231  {
232  std::swap(coeff(r1, j), coeff(r2, j));
233  }
234  }
235 
236  // lambda = |A[r, k]| = max{|A[k+1, k]|, ..., |A[end, k]|}
237  // Largest (in magnitude) off-diagonal element in the first column of the current reduced matrix
238  // r is the row index
239  // Assume k < end
241  {
242  using std::abs;
243 
244  const Scalar* head = col_pointer(k); // => A[k, k]
245  const Scalar* end = col_pointer(k + 1);
246  // Start with r=k+1, lambda=A[k+1, k]
247  r = k + 1;
248  RealScalar lambda = abs(head[1]);
249  // Scan remaining elements
250  for (const Scalar* ptr = head + 2; ptr < end; ptr++)
251  {
252  const RealScalar abs_elem = abs(*ptr);
253  if (lambda < abs_elem)
254  {
255  lambda = abs_elem;
256  r = k + (ptr - head);
257  }
258  }
259 
260  return lambda;
261  }
262 
263  // sigma = |A[p, r]| = max {|A[k, r]|, ..., |A[end, r]|} \ {A[r, r]}
264  // Largest (in magnitude) off-diagonal element in the r-th column of the current reduced matrix
265  // p is the row index
266  // Assume k < r < end
268  {
269  using std::abs;
270 
271  // First search A[r+1, r], ..., A[end, r], which has the same task as find_lambda()
272  // If r == end, we skip this search
274  if (r < m_n - 1)
275  sigma = find_lambda(r, p);
276 
277  // Then search A[k, r], ..., A[r-1, r], which maps to A[r, k], ..., A[r, r-1]
278  for (Index j = k; j < r; j++)
279  {
280  const RealScalar abs_elem = abs(coeff(r, j));
281  if (sigma < abs_elem)
282  {
283  sigma = abs_elem;
284  p = j;
285  }
286  }
287 
288  return sigma;
289  }
290 
291  // Generate permutations and apply to A
292  // Return true if the resulting pivoting is 1x1, and false if 2x2
294  {
295  using std::abs;
296 
297  Index r = k, p = k;
298  const RealScalar lambda = find_lambda(k, r);
299 
300  // If lambda=0, no need to interchange
301  if (lambda > RealScalar(0))
302  {
303  const RealScalar abs_akk = abs(diag_coeff(k));
304  // If |A[k, k]| >= alpha * lambda, no need to interchange
305  if (abs_akk < alpha * lambda)
306  {
307  const RealScalar sigma = find_sigma(k, r, p);
308 
309  // If sigma * |A[k, k]| >= alpha * lambda^2, no need to interchange
310  if (sigma * abs_akk < alpha * lambda * lambda)
311  {
312  if (abs_akk >= alpha * sigma)
313  {
314  // Permutation on A
315  pivoting_1x1(k, r);
316 
317  // Permutation on L
318  interchange_rows(k, r, 0, k - 1);
319  return true;
320  }
321  else
322  {
323  // There are two versions of permutation here
324  // 1. A[k+1, k] <-> A[r, k]
325  // 2. A[k+1, k] <-> A[r, p], where p >= k and r >= k+1
326  //
327  // Version 1 and 2 are used by Ref[1] and Ref[2], respectively
328 
329  // Version 1 implementation
330  p = k;
331 
332  // Version 2 implementation
333  // [r, p] and [p, r] are symmetric, but we need to make sure
334  // p >= k and r >= k+1, so it is safe to always make r > p
335  // One exception is when min{r,p} == k+1, in which case we make
336  // r = k+1, so that only one permutation needs to be performed
337  /* const Index rp_min = std::min(r, p);
338  const Index rp_max = std::max(r, p);
339  if (rp_min == k + 1)
340  {
341  r = rp_min; p = rp_max;
342  } else {
343  r = rp_max; p = rp_min;
344  } */
345 
346  // Right now we use Version 1 since it reduces the overhead of interchange
347 
348  // Permutation on A
349  pivoting_2x2(k, r, p);
350 
351  // Permutation on L
352  interchange_rows(k, p, 0, k - 1);
353  interchange_rows(k + 1, r, 0, k - 1);
354  return false;
355  }
356  }
357  }
358  }
359 
360  return true;
361  }
362 
363  // E = [e11, e12]
364  // [e21, e22]
365  // Overwrite E with inv(E)
366  void inverse_inplace_2x2(Scalar& e11, Scalar& e21, Scalar& e22) const
367  {
368  // inv(E) = [d11, d12], d11 = e22/delta, d21 = -e21/delta, d22 = e11/delta
369  // [d21, d22]
370  // delta = e11 * e22 - e12 * e21
371  const Scalar e12 = ScalarOp<Scalar>::conj(e21);
372  const Scalar delta = e11 * e22 - e12 * e21;
373  std::swap(e11, e22);
374  e11 /= delta;
375  e22 /= delta;
376  e21 = -e21 / delta;
377  }
378 
379  // E = [e11, e12]
380  // [e21, e22]
381  // Overwrite b with x = inv(E) * b, which is equivalent to solving E * x = b
383  const Scalar& e11, const Scalar& e21, const Scalar& e22,
384  Scalar& b1, Scalar& b2) const
385  {
386  using std::abs;
387 
388  const Scalar e12 = ScalarOp<Scalar>::conj(e21);
389  const RealScalar e11_abs = abs(e11);
390  const RealScalar e21_abs = abs(e21);
391  // If |e11| >= |e21|, no need to exchange rows
392  if (e11_abs >= e21_abs)
393  {
394  const Scalar fac = e21 / e11;
395  const Scalar x2 = (b2 - fac * b1) / (e22 - fac * e12);
396  const Scalar x1 = (b1 - e12 * x2) / e11;
397  b1 = x1;
398  b2 = x2;
399  }
400  else
401  {
402  // Exchange row 1 and row 2, so the system becomes
403  // E* = [e21, e22], b* = [b2], x* = [x1]
404  // [e11, e12] [b1] [x2]
405  const Scalar fac = e11 / e21;
406  const Scalar x2 = (b1 - fac * b2) / (e12 - fac * e22);
407  const Scalar x1 = (b2 - e22 * x2) / e21;
408  b1 = x1;
409  b2 = x2;
410  }
411  }
412 
413  // Compute C * inv(E), which is equivalent to solving X * E = C
414  // X [n x 2], E [2 x 2], C [n x 2]
415  // X = [x1, x2], E = [e11, e12], C = [c1 c2]
416  // [e21, e22]
418  const Scalar& e11, const Scalar& e21, const Scalar& e22,
419  const MapVec& c1, const MapVec& c2,
421  {
422  using std::abs;
423 
424  const Scalar e12 = ScalarOp<Scalar>::conj(e21);
425  const RealScalar e11_abs = abs(e11);
426  const RealScalar e12_abs = abs(e12);
427  // If |e11| >= |e12|, no need to exchange rows
428  if (e11_abs >= e12_abs)
429  {
430  const Scalar fac = e12 / e11;
431  // const Scalar x2 = (c2 - fac * c1) / (e22 - fac * e21);
432  // const Scalar x1 = (c1 - e21 * x2) / e11;
433  x.col(1).array() = (c2 - fac * c1).array() / (e22 - fac * e21);
434  x.col(0).array() = (c1 - e21 * x.col(1)).array() / e11;
435  }
436  else
437  {
438  // Exchange column 1 and column 2, so the system becomes
439  // X* = [x1, x2], E* = [e12, e11], C* = [c2 c1]
440  // [e22, e21]
441  const Scalar fac = e11 / e12;
442  // const Scalar x2 = (c1 - fac * c2) / (e21 - fac * e22);
443  // const Scalar x1 = (c2 - e22 * x2) / e12;
444  x.col(1).array() = (c1 - fac * c2).array() / (e21 - fac * e22);
445  x.col(0).array() = (c2 - e22 * x.col(1)).array() / e12;
446  }
447  }
448 
449  // Return value is the status, CompInfo::Successful/NumericalIssue
451  {
452  // A[k, k] is known to be real-valued, so we force its imaginary
453  // part to be zero when Scalar is a complex type
454  // Interestingly, this has a significant effect on the accuracy
455  // and numerical stability of the final solution
456  const Scalar akk = ScalarOp<Scalar>::real(diag_coeff(k));
457  diag_coeff(k) = akk;
458  // Return CompInfo::NumericalIssue if not invertible
459  if (akk == Scalar(0))
461 
462  // [inverse]
463  // diag_coeff(k) = Scalar(1) / akk;
464 
465  // B -= l * l^H / A[k, k], B := A[(k+1):end, (k+1):end], l := L[(k+1):end, k]
466  Scalar* lptr = col_pointer(k) + 1;
467  const Index ldim = m_n - k - 1;
468  MapVec l(lptr, ldim);
469  for (Index j = 0; j < ldim; j++)
470  {
471  Scalar l_conj = ScalarOp<Scalar>::conj(lptr[j]);
472  MapVec(col_pointer(j + k + 1), ldim - j).noalias() -= (l_conj / akk) * l.tail(ldim - j);
473  }
474 
475  // l /= A[k, k]
476  l /= akk;
477 
478  return CompInfo::Successful;
479  }
480 
481  // Return value is the status, CompInfo::Successful/NumericalIssue
483  {
484  Scalar& e11 = diag_coeff(k);
485  Scalar& e21 = coeff(k + 1, k);
486  Scalar& e22 = diag_coeff(k + 1);
487 
488  // A[k, k] and A[k+1, k+1] are known to be real-valued,
489  // so we force their imaginary parts to be zero when Scalar
490  // is a complex type
491  // Interestingly, this has a significant effect on the accuracy
492  // and numerical stability of the final solution
493  e11 = ScalarOp<Scalar>::real(e11);
494  e22 = ScalarOp<Scalar>::real(e22);
495  Scalar e12 = ScalarOp<Scalar>::conj(e21);
496  // Return CompInfo::NumericalIssue if not invertible
497  if (e11 * e22 - e12 * e21 == Scalar(0))
499 
500  // [inverse]
501  // inverse_inplace_2x2(e11, e21, e22);
502 
503  // X = l * inv(E), l := L[(k+2):end, k:(k+1)]
504  Scalar* l1ptr = &coeff(k + 2, k);
505  Scalar* l2ptr = &coeff(k + 2, k + 1);
506  const Index ldim = m_n - k - 2;
507  MapVec l1(l1ptr, ldim), l2(l2ptr, ldim);
508 
510  // [inverse]
511  // e12 = ScalarOp<Scalar>::conj(e21);
512  // X.col(0).noalias() = l1 * e11 + l2 * e21;
513  // X.col(1).noalias() = l1 * e12 + l2 * e22;
514  // [solve]
515  solve_left_2x2(e11, e21, e22, l1, l2, X);
516 
517  // B -= l * inv(E) * l^H = X * l^H, B = A[(k+2):end, (k+2):end]
518  for (Index j = 0; j < ldim; j++)
519  {
520  const Scalar l1j_conj = ScalarOp<Scalar>::conj(l1ptr[j]);
521  const Scalar l2j_conj = ScalarOp<Scalar>::conj(l2ptr[j]);
522  MapVec(col_pointer(j + k + 2), ldim - j).noalias() -= (X.col(0).tail(ldim - j) * l1j_conj + X.col(1).tail(ldim - j) * l2j_conj);
523  }
524 
525  // l = X
526  l1.noalias() = X.col(0);
527  l2.noalias() = X.col(1);
528 
529  return CompInfo::Successful;
530  }
531 
532 public:
533  BKLDLT() :
535  {}
536 
537  // Factorize mat - shift * I
538  template <typename Derived>
539  BKLDLT(const Eigen::MatrixBase<Derived>& mat, int uplo = Eigen::Lower, const RealScalar& shift = RealScalar(0)) :
541  {
542  compute(mat, uplo, shift);
543  }
544 
545  template <typename Derived>
546  void compute(const Eigen::MatrixBase<Derived>& mat, int uplo = Eigen::Lower, const RealScalar& shift = RealScalar(0))
547  {
548  using std::abs;
549 
550  m_n = mat.rows();
551  if (m_n != mat.cols())
552  throw std::invalid_argument("BKLDLT: matrix must be square");
553 
554  m_perm.setLinSpaced(m_n, 0, m_n - 1);
555  m_permc.clear();
556 
557  // Copy data
558  m_data.resize((m_n * (m_n + 1)) / 2);
559  compute_pointer();
560  copy_data(mat, uplo, shift);
561 
562  const RealScalar alpha = (1.0 + std::sqrt(17.0)) / 8.0;
563  Index k = 0;
564  for (k = 0; k < m_n - 1; k++)
565  {
566  // 1. Interchange rows and columns of A, and save the result to m_perm
567  bool is_1x1 = permutate_mat(k, alpha);
568 
569  // 2. Gaussian elimination
570  if (is_1x1)
571  {
573  }
574  else
575  {
577  k++;
578  }
579 
580  // 3. Check status
582  break;
583  }
584  // Invert the last 1x1 block if it exists
585  if (k == m_n - 1)
586  {
587  const Scalar akk = ScalarOp<Scalar>::real(diag_coeff(k));
588  diag_coeff(k) = akk;
589  if (akk == Scalar(0))
591 
592  // [inverse]
593  // diag_coeff(k) = Scalar(1) / diag_coeff(k);
594  }
595 
597 
598  m_computed = true;
599  }
600 
601  // Solve Ax=b
603  {
604  if (!m_computed)
605  throw std::logic_error("BKLDLT: need to call compute() first");
606 
607  // PAP' = LD(L^H), A = P'LD(L^H)P
608  // Ax = b ==> P'LD(L^H)Px = b ==> LD(L^H)Px = Pb
609  // 1. b -> Pb
610  Scalar* x = b.data();
611  MapVec res(x, m_n);
612  Index npermc = m_permc.size();
613  for (Index i = 0; i < npermc; i++)
614  {
615  std::swap(x[m_permc[i].first], x[m_permc[i].second]);
616  }
617 
618  // z = D(L^H)Px
619  // 2. Lz = Pb
620  // If m_perm[end] < 0, then end with m_n - 3, otherwise end with m_n - 2
621  const Index end = (m_perm[m_n - 1] < 0) ? (m_n - 3) : (m_n - 2);
622  for (Index i = 0; i <= end; i++)
623  {
624  const Index b1size = m_n - i - 1;
625  const Index b2size = b1size - 1;
626  if (m_perm[i] >= 0)
627  {
628  MapConstVec l(&coeff(i + 1, i), b1size);
629  res.segment(i + 1, b1size).noalias() -= l * x[i];
630  }
631  else
632  {
633  MapConstVec l1(&coeff(i + 2, i), b2size);
634  MapConstVec l2(&coeff(i + 2, i + 1), b2size);
635  res.segment(i + 2, b2size).noalias() -= (l1 * x[i] + l2 * x[i + 1]);
636  i++;
637  }
638  }
639 
640  // w = (L^H)Px
641  // 3. Dw = z
642  for (Index i = 0; i < m_n; i++)
643  {
644  const Scalar e11 = diag_coeff(i);
645  if (m_perm[i] >= 0)
646  {
647  // [inverse]
648  // x[i] *= e11;
649  // [solve]
650  x[i] /= e11;
651  }
652  else
653  {
654  const Scalar e21 = coeff(i + 1, i), e22 = diag_coeff(i + 1);
655  // [inverse]
656  // const Scalar e12 = ScalarOp<Scalar>::conj(e21);
657  // const Scalar wi = x[i] * e11 + x[i + 1] * e12;
658  // x[i + 1] = x[i] * e21 + x[i + 1] * e22;
659  // x[i] = wi;
660  // [solve]
661  solve_inplace_2x2(e11, e21, e22, x[i], x[i + 1]);
662 
663  i++;
664  }
665  }
666 
667  // y = Px
668  // 4. (L^H)y = w
669  // If m_perm[end] < 0, then start with m_n - 3, otherwise start with m_n - 2
670  Index i = (m_perm[m_n - 1] < 0) ? (m_n - 3) : (m_n - 2);
671  for (; i >= 0; i--)
672  {
673  const Index ldim = m_n - i - 1;
674  MapConstVec l(&coeff(i + 1, i), ldim);
675  x[i] -= l.dot(res.segment(i + 1, ldim));
676 
677  if (m_perm[i] < 0)
678  {
679  MapConstVec l2(&coeff(i + 1, i - 1), ldim);
680  x[i - 1] -= l2.dot(res.segment(i + 1, ldim));
681  i--;
682  }
683  }
684 
685  // 5. x = P'y
686  for (Index i = npermc - 1; i >= 0; i--)
687  {
688  std::swap(x[m_permc[i].first], x[m_permc[i].second]);
689  }
690  }
691 
693  {
694  Vector res = b;
696  return res;
697  }
698 
699  CompInfo info() const { return m_info; }
700 };
701 
702 } // namespace Spectra
703 
704 #endif // SPECTRA_BK_LDLT_H
Spectra::BKLDLT::BKLDLT
BKLDLT(const Eigen::MatrixBase< Derived > &mat, int uplo=Eigen::Lower, const RealScalar &shift=RealScalar(0))
Definition: BKLDLT.h:539
perm
idx_t idx_t idx_t idx_t idx_t * perm
Definition: include/metis.h:223
simple_graph::b1
Vector2 b1(2, -1)
Spectra::BKLDLT::m_data
Vector m_data
Definition: BKLDLT.h:79
array
int array[24]
Definition: Map_general_stride.cpp:1
Spectra::BKLDLT::inverse_inplace_2x2
void inverse_inplace_2x2(Scalar &e11, Scalar &e21, Scalar &e22) const
Definition: BKLDLT.h:366
Spectra::CompInfo::Successful
@ Successful
Computation was successful.
alpha
RealScalar alpha
Definition: level1_cplx_impl.h:147
Spectra::BKLDLT::solve_inplace
void solve_inplace(GenericVector b) const
Definition: BKLDLT.h:602
l2
gtsam::Key l2
Definition: testLinearContainerFactor.cpp:24
Spectra::BKLDLT::copy_data
void copy_data(const Eigen::MatrixBase< Derived > &mat, int uplo, const RealScalar &shift)
Definition: BKLDLT.h:113
r2
static const double r2
Definition: testSmartRangeFactor.cpp:32
b
Scalar * b
Definition: benchVecAdd.cpp:17
simple_graph::b2
Vector2 b2(4, -5)
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
Spectra::BKLDLT::coeff
const Scalar & coeff(Index i, Index j) const
Definition: BKLDLT.h:92
Spectra::BKLDLT::col_pointer
Scalar * col_pointer(Index k)
Definition: BKLDLT.h:89
Spectra::BKLDLT::m_permc
std::vector< std::pair< Index, Index > > m_permc
Definition: BKLDLT.h:82
Spectra::BKLDLT::BKLDLT
BKLDLT()
Definition: BKLDLT.h:533
Spectra::BKLDLT::coeff
Scalar & coeff(Index i, Index j)
Definition: BKLDLT.h:91
X
#define X
Definition: icosphere.cpp:20
Spectra::BKLDLT::compute
void compute(const Eigen::MatrixBase< Derived > &mat, int uplo=Eigen::Lower, const RealScalar &shift=RealScalar(0))
Definition: BKLDLT.h:546
mat
MatrixXf mat
Definition: Tutorial_AdvancedInitialization_CommaTemporary.cpp:1
res
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Definition: PartialRedux_count.cpp:3
Eigen::PlainObjectBase::resize
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Spectra::BKLDLT::MapVec
Eigen::Map< Vector > MapVec
Definition: BKLDLT.h:72
Spectra::BKLDLT::m_n
Index m_n
Definition: BKLDLT.h:78
r1
static const double r1
Definition: testSmartRangeFactor.cpp:32
rows
int rows
Definition: Tutorial_commainit_02.cpp:1
sampling::sigma
static const double sigma
Definition: testGaussianBayesNet.cpp:170
c1
static double c1
Definition: airy.c:54
Spectra::BKLDLT::diag_coeff
const Scalar & diag_coeff(Index i) const
Definition: BKLDLT.h:95
Spectra::BKLDLT::interchange_rows
void interchange_rows(Index r1, Index r2, Index c1, Index c2)
Definition: BKLDLT.h:225
Spectra::BKLDLT::Index
Eigen::Index Index
Definition: BKLDLT.h:70
gtsam.examples.PlanarManipulatorExample.delta
def delta(g0, g1)
Definition: PlanarManipulatorExample.py:45
Spectra::BKLDLT::diag_coeff
Scalar & diag_coeff(Index i)
Definition: BKLDLT.h:94
Spectra::BKLDLT::m_info
CompInfo m_info
Definition: BKLDLT.h:85
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
Spectra::CompInfo
CompInfo
Definition: CompInfo.h:17
x1
Pose3 x1
Definition: testPose3.cpp:663
Spectra::BKLDLT::find_sigma
RealScalar find_sigma(Index k, Index r, Index &p)
Definition: BKLDLT.h:267
Spectra::BKLDLT::solve_inplace_2x2
void solve_inplace_2x2(const Scalar &e11, const Scalar &e21, const Scalar &e22, Scalar &b1, Scalar &b2) const
Definition: BKLDLT.h:382
Spectra::BKLDLT::pivoting_2x2
void pivoting_2x2(Index k, Index r, Index p)
Definition: BKLDLT.h:209
l
static const Line3 l(Rot3(), 1, 1)
Spectra::BKLDLT::RealScalar
typename Eigen::NumTraits< Scalar >::Real RealScalar
Definition: BKLDLT.h:69
std::swap
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Definition: NearestNeighbor.hpp:827
Spectra::ScalarOp::real
static Scalar real(const Scalar &x)
Definition: BKLDLT.h:35
Spectra::BKLDLT::permutate_mat
bool permutate_mat(Index k, const RealScalar &alpha)
Definition: BKLDLT.h:293
head
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE FixedSegmentReturnType< internal::get_fixed_value< NType >::value >::Type head(NType n)
Definition: BlockMethods.h:1208
conj
AnnoyingScalar conj(const AnnoyingScalar &x)
Definition: AnnoyingScalar.h:104
Eigen::Lower
@ Lower
Definition: Constants.h:209
Spectra::BKLDLT::solve_left_2x2
void solve_left_2x2(const Scalar &e11, const Scalar &e21, const Scalar &e22, const MapVec &c1, const MapVec &c2, Eigen::Matrix< Scalar, Eigen::Dynamic, 2 > &x) const
Definition: BKLDLT.h:417
Eigen::Map
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
lambda
static double lambda[]
Definition: jv.c:524
Spectra::BKLDLT::solve
Vector solve(ConstGenericVector &b) const
Definition: BKLDLT.h:692
Spectra::BKLDLT::compress_permutation
void compress_permutation()
Definition: BKLDLT.h:149
Spectra::ScalarOp
Definition: BKLDLT.h:28
Spectra::CompInfo::NumericalIssue
@ NumericalIssue
Spectra::BKLDLT::compute_pointer
void compute_pointer()
Definition: BKLDLT.h:98
RealScalar
NumTraits< Scalar >::Real RealScalar
Definition: bench_gemm.cpp:47
Eigen::Ref
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:281
l1
gtsam::Key l1
Definition: testLinearContainerFactor.cpp:24
Spectra::BKLDLT::gaussian_elimination_1x1
CompInfo gaussian_elimination_1x1(Index k)
Definition: BKLDLT.h:450
std
Definition: BFloat16.h:88
p
float * p
Definition: Tutorial_Map_using.cpp:9
Spectra::ScalarOp< std::complex< RealScalar > >::real
static std::complex< RealScalar > real(const std::complex< RealScalar > &x)
Definition: BKLDLT.h:50
Spectra::BKLDLT::m_perm
IntVector m_perm
Definition: BKLDLT.h:81
Spectra
Definition: LOBPCGSolver.h:19
c2
static double c2
Definition: airy.c:55
Spectra::BKLDLT::m_computed
bool m_computed
Definition: BKLDLT.h:84
Eigen::PlainObjectBase::data
EIGEN_DEVICE_FUNC const EIGEN_STRONG_INLINE Scalar * data() const
Definition: PlainObjectBase.h:247
Eigen::Matrix< Scalar, Eigen::Dynamic, 1 >
Spectra::BKLDLT::pivoting_1x1
void pivoting_1x1(Index k, Index r)
Definition: BKLDLT.h:163
abs
#define abs(x)
Definition: datatypes.h:17
len
size_t len(handle h)
Get the length of a Python object.
Definition: pytypes.h:2446
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
Spectra::ScalarOp::conj
static Scalar conj(const Scalar &x)
Definition: BKLDLT.h:30
Spectra::BKLDLT::find_lambda
RealScalar find_lambda(Index k, Index &r)
Definition: BKLDLT.h:240
Spectra::BKLDLT::m_colptr
std::vector< Scalar * > m_colptr
Definition: BKLDLT.h:80
Eigen::placeholders::end
static const EIGEN_DEPRECATED end_t end
Definition: IndexedViewHelper.h:181
Spectra::CompInfo::NotComputed
@ NotComputed
Spectra::BKLDLT::gaussian_elimination_2x2
CompInfo gaussian_elimination_2x2(Index k)
Definition: BKLDLT.h:482
Spectra::ScalarOp< std::complex< RealScalar > >::conj
static std::complex< RealScalar > conj(const std::complex< RealScalar > &x)
Definition: BKLDLT.h:44
Spectra::BKLDLT::info
CompInfo info() const
Definition: BKLDLT.h:699
test_callbacks.value
value
Definition: test_callbacks.py:160
x2
Pose3 x2(Rot3::Ypr(0.0, 0.0, 0.0), l2)
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
Spectra::BKLDLT
Definition: BKLDLT.h:65
Eigen::GenericNumTraits::Real
T Real
Definition: NumTraits.h:164
complex
Definition: datatypes.h:12
Scalar
SCALAR Scalar
Definition: bench_gemm.cpp:46
Eigen::Index
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74


gtsam
Author(s):
autogenerated on Sun Feb 16 2025 04:01:02