Molassembler  3.0.0
Molecule graph and conformer library
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
SylvestersCriterion.h
Go to the documentation of this file.
1 
9 #ifndef INCLUDE_TEMPLE_SYLVESTERS_CRITERION_H
10 #define INCLUDE_TEMPLE_SYLVESTERS_CRITERION_H
11 
12 #include <Eigen/Core>
13 
14 namespace Scine {
15 namespace Molassembler {
16 namespace Temple {
17 
25 template<typename Derived>
26 bool positiveDefinite(const Eigen::MatrixBase<Derived>& matrix) {
27  using Scalar = typename Eigen::MatrixBase<Derived>::Scalar;
28  assert(matrix.rows() == matrix.cols());
29 
30  const unsigned N = matrix.rows();
31 
32  // 1st leading principal minor
33  if(matrix(0, 0) <= Scalar {0.0}) {
34  return false;
35  }
36 
37  if(N == 1) {
38  return true;
39  }
40 
41  /* Test leading principal minors of sizes 2-4 using fixed-size matrix
42  * subexpressions, those should be significantly faster
43  */
44  // 2nd leading principal minor
45  if(matrix.template block<2, 2>(0, 0).determinant() <= Scalar {0.0}) {
46  return false;
47  }
48 
49  if(N == 2) {
50  return true;
51  }
52 
53  // 3rd leading principal minor
54  if(matrix.template block<3, 3>(0, 0).determinant() <= Scalar {0.0}) {
55  return false;
56  }
57 
58  if(N == 3) {
59  return true;
60  }
61 
62  // 4th leading principal minor
63  if(matrix.template block<4, 4>(0, 0).determinant() <= Scalar {0.0}) {
64  return false;
65  }
66 
67  if(N == 4) {
68  return true;
69  }
70 
71  /* There is limited speed advantage to fixed-size matrix sub-expressions from
72  * here, so we can now loop over the rest:
73  */
74  for(unsigned I = 5; I <= N; ++I) {
75  if(matrix.block(0, 0, I, I).determinant() <= Scalar {0.0}) {
76  return false;
77  }
78  }
79 
80  return true;
81 }
82 
90 template<typename Derived>
91 bool positiveSemidefinite(const Eigen::MatrixBase<Derived>& matrix) {
92  using Scalar = typename Eigen::MatrixBase<Derived>::Scalar;
93  assert(matrix.rows() == matrix.cols());
94 
95  const unsigned N = matrix.rows();
96 
97  // Principal minors of size I = 1
98  if((matrix.array() < Scalar {0.0}).any()) {
99  return false;
100  }
101 
102  if(N == 1) {
103  return true;
104  }
105 
106  // Principal minors of size I = 2
107  for(unsigned i = 0; i < N - 1; ++i) {
108  for(unsigned j = 0; j < N - 1; ++j) {
109  if(matrix.template block<2, 2>(i, j).determinant() < Scalar {0.0}) {
110  return false;
111  }
112  }
113  }
114 
115  if(N == 2) {
116  return true;
117  }
118 
119  // Principal minors of size I = 3
120  for(unsigned i = 0; i < N - 2; ++i) {
121  for(unsigned j = 0; j < N - 2; ++j) {
122  if(matrix.template block<3, 3>(i, j).determinant() < Scalar {0.0}) {
123  return false;
124  }
125  }
126  }
127 
128  if(N == 3) {
129  return true;
130  }
131 
132  // Principal minors of size I = 4
133  for(unsigned i = 0; i < N - 3; ++i) {
134  for(unsigned j = 0; j < N - 3; ++j) {
135  if(matrix.template block<4, 4>(i, j).determinant() < Scalar {0.0}) {
136  return false;
137  }
138  }
139  }
140 
141  if(N == 4) {
142  return true;
143  }
144 
145  // Remaining principal minors
146  for(unsigned I = 5; I <= N; ++I) {
147  for(unsigned i = 0; i < N - I; ++i) {
148  for(unsigned j = 0; j < N - I; ++j) {
149  if(matrix.block(i, j, I, I).determinant() < Scalar {0.0}) {
150  return false;
151  }
152  }
153  }
154  }
155 
156  return true;
157 }
158 
159 } // namespace Temple
160 } // namespace Molassembler
161 } // namespace Scine
162 
163 #endif
bool positiveSemidefinite(const Eigen::MatrixBase< Derived > &matrix)
Determine whether a matrix is positive semidefinite.
Definition: SylvestersCriterion.h:91
bool positiveDefinite(const Eigen::MatrixBase< Derived > &matrix)
Determine whether a matrix is positive definite.
Definition: SylvestersCriterion.h:26