Lines 61-69
Link Here
|
61 |
> struct qr_preconditioner_impl {}; |
61 |
> struct qr_preconditioner_impl {}; |
62 |
|
62 |
|
63 |
template<typename MatrixType, int QRPreconditioner, int Case> |
63 |
template<typename MatrixType, int QRPreconditioner, int Case> |
64 |
struct qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> |
64 |
class qr_preconditioner_impl<MatrixType, QRPreconditioner, Case, false> |
65 |
{ |
65 |
{ |
66 |
static bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) |
66 |
public: |
|
|
67 |
typedef typename MatrixType::Index Index; |
68 |
void resize(Index, Index) {} |
69 |
bool run(JacobiSVD<MatrixType, QRPreconditioner>&, const MatrixType&) |
67 |
{ |
70 |
{ |
68 |
return false; |
71 |
return false; |
69 |
} |
72 |
} |
Lines 72-205
Link Here
|
72 |
/*** preconditioner using FullPivHouseholderQR ***/ |
75 |
/*** preconditioner using FullPivHouseholderQR ***/ |
73 |
|
76 |
|
74 |
template<typename MatrixType> |
77 |
template<typename MatrixType> |
75 |
struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
78 |
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
76 |
{ |
79 |
{ |
77 |
static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
80 |
public: |
|
|
81 |
typedef typename MatrixType::Index Index; |
82 |
|
83 |
void resize(Index rows, Index cols) |
84 |
{ |
85 |
if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = FullPivHouseholderQR<MatrixType>(rows, cols); |
86 |
} |
87 |
|
88 |
bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
78 |
{ |
89 |
{ |
79 |
if(matrix.rows() > matrix.cols()) |
90 |
if(matrix.rows() > matrix.cols()) |
80 |
{ |
91 |
{ |
81 |
FullPivHouseholderQR<MatrixType> qr(matrix); |
92 |
m_qr.compute(matrix); |
82 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
93 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
83 |
if(svd.m_computeFullU) svd.m_matrixU = qr.matrixQ(); |
94 |
if(svd.m_computeFullU) svd.m_matrixU = m_qr.matrixQ(); |
84 |
if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); |
95 |
if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); |
85 |
return true; |
96 |
return true; |
86 |
} |
97 |
} |
87 |
return false; |
98 |
return false; |
88 |
} |
99 |
} |
|
|
100 |
private: |
101 |
FullPivHouseholderQR<MatrixType> m_qr; |
89 |
}; |
102 |
}; |
90 |
|
103 |
|
91 |
template<typename MatrixType> |
104 |
template<typename MatrixType> |
92 |
struct qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
105 |
class qr_preconditioner_impl<MatrixType, FullPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
93 |
{ |
106 |
{ |
94 |
static bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
107 |
public: |
|
|
108 |
typedef typename MatrixType::Index Index; |
109 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
110 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
111 |
TransposeTypeWithSameStorageOrder; |
112 |
|
113 |
void resize(Index rows, Index cols) |
114 |
{ |
115 |
if (cols != m_qr.rows() || rows != m_qr.cols()) |
116 |
{ |
117 |
m_qr = FullPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows); |
118 |
} |
119 |
m_adjoint.resize(cols, rows); |
120 |
} |
121 |
|
122 |
bool run(JacobiSVD<MatrixType, FullPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
95 |
{ |
123 |
{ |
96 |
if(matrix.cols() > matrix.rows()) |
124 |
if(matrix.cols() > matrix.rows()) |
97 |
{ |
125 |
{ |
98 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
126 |
m_adjoint = matrix.adjoint(); |
99 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
127 |
m_qr.compute(m_adjoint); |
100 |
TransposeTypeWithSameStorageOrder; |
128 |
|
101 |
FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); |
129 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
102 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
130 |
if(svd.m_computeFullV) svd.m_matrixV = m_qr.matrixQ(); |
103 |
if(svd.m_computeFullV) svd.m_matrixV = qr.matrixQ(); |
131 |
if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); |
104 |
if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); |
|
|
105 |
return true; |
132 |
return true; |
106 |
} |
133 |
} |
107 |
else return false; |
134 |
else return false; |
108 |
} |
135 |
} |
|
|
136 |
private: |
137 |
FullPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; |
138 |
TransposeTypeWithSameStorageOrder m_adjoint; |
109 |
}; |
139 |
}; |
110 |
|
140 |
|
111 |
/*** preconditioner using ColPivHouseholderQR ***/ |
141 |
/*** preconditioner using ColPivHouseholderQR ***/ |
112 |
|
142 |
|
113 |
template<typename MatrixType> |
143 |
template<typename MatrixType> |
114 |
struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
144 |
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
115 |
{ |
145 |
{ |
116 |
static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
146 |
public: |
|
|
147 |
typedef typename MatrixType::Index Index; |
148 |
|
149 |
void resize(Index rows, Index cols) |
150 |
{ |
151 |
if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = ColPivHouseholderQR<MatrixType>(rows, cols); |
152 |
} |
153 |
|
154 |
bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
117 |
{ |
155 |
{ |
118 |
if(matrix.rows() > matrix.cols()) |
156 |
if(matrix.rows() > matrix.cols()) |
119 |
{ |
157 |
{ |
120 |
ColPivHouseholderQR<MatrixType> qr(matrix); |
158 |
m_qr.compute(matrix); |
121 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
159 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
122 |
if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); |
160 |
if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ(); |
123 |
else if(svd.m_computeThinU) { |
161 |
else if(svd.m_computeThinU) { |
124 |
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); |
162 |
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); |
125 |
qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); |
163 |
m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); |
126 |
} |
164 |
} |
127 |
if(svd.computeV()) svd.m_matrixV = qr.colsPermutation(); |
165 |
if(svd.computeV()) svd.m_matrixV = m_qr.colsPermutation(); |
128 |
return true; |
166 |
return true; |
129 |
} |
167 |
} |
130 |
return false; |
168 |
return false; |
131 |
} |
169 |
} |
|
|
170 |
|
171 |
private: |
172 |
ColPivHouseholderQR<MatrixType> m_qr; |
132 |
}; |
173 |
}; |
133 |
|
174 |
|
134 |
template<typename MatrixType> |
175 |
template<typename MatrixType> |
135 |
struct qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
176 |
class qr_preconditioner_impl<MatrixType, ColPivHouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
136 |
{ |
177 |
{ |
137 |
static bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
178 |
public: |
|
|
179 |
typedef typename MatrixType::Index Index; |
180 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
181 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
182 |
TransposeTypeWithSameStorageOrder; |
183 |
|
184 |
void resize(Index rows, Index cols) |
185 |
{ |
186 |
if (cols != m_qr.rows() || rows != m_qr.cols()) |
187 |
{ |
188 |
m_qr = ColPivHouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows); |
189 |
} |
190 |
m_adjoint.resize(cols, rows); |
191 |
} |
192 |
|
193 |
bool run(JacobiSVD<MatrixType, ColPivHouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
138 |
{ |
194 |
{ |
139 |
if(matrix.cols() > matrix.rows()) |
195 |
if(matrix.cols() > matrix.rows()) |
140 |
{ |
196 |
{ |
141 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
197 |
m_adjoint = matrix.adjoint(); |
142 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
198 |
m_qr.compute(m_adjoint); |
143 |
TransposeTypeWithSameStorageOrder; |
199 |
|
144 |
ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); |
200 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
145 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
201 |
if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ(); |
146 |
if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); |
|
|
147 |
else if(svd.m_computeThinV) { |
202 |
else if(svd.m_computeThinV) { |
148 |
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); |
203 |
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); |
149 |
qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); |
204 |
m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); |
150 |
} |
205 |
} |
151 |
if(svd.computeU()) svd.m_matrixU = qr.colsPermutation(); |
206 |
if(svd.computeU()) svd.m_matrixU = m_qr.colsPermutation(); |
152 |
return true; |
207 |
return true; |
153 |
} |
208 |
} |
154 |
else return false; |
209 |
else return false; |
155 |
} |
210 |
} |
|
|
211 |
|
212 |
private: |
213 |
ColPivHouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; |
214 |
TransposeTypeWithSameStorageOrder m_adjoint; |
156 |
}; |
215 |
}; |
157 |
|
216 |
|
158 |
/*** preconditioner using HouseholderQR ***/ |
217 |
/*** preconditioner using HouseholderQR ***/ |
159 |
|
218 |
|
160 |
template<typename MatrixType> |
219 |
template<typename MatrixType> |
161 |
struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
220 |
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreRowsThanCols, true> |
162 |
{ |
221 |
{ |
163 |
static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
222 |
public: |
|
|
223 |
typedef typename MatrixType::Index Index; |
224 |
|
225 |
void resize(Index rows, Index cols) |
226 |
{ |
227 |
if (rows != m_qr.rows() || cols != m_qr.cols()) m_qr = HouseholderQR<MatrixType>(rows, cols); |
228 |
} |
229 |
|
230 |
bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
164 |
{ |
231 |
{ |
165 |
if(matrix.rows() > matrix.cols()) |
232 |
if(matrix.rows() > matrix.cols()) |
166 |
{ |
233 |
{ |
167 |
HouseholderQR<MatrixType> qr(matrix); |
234 |
m_qr.compute(matrix); |
168 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
235 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.cols(),matrix.cols()).template triangularView<Upper>(); |
169 |
if(svd.m_computeFullU) svd.m_matrixU = qr.householderQ(); |
236 |
if(svd.m_computeFullU) svd.m_matrixU = m_qr.householderQ(); |
170 |
else if(svd.m_computeThinU) { |
237 |
else if(svd.m_computeThinU) { |
171 |
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); |
238 |
svd.m_matrixU.setIdentity(matrix.rows(), matrix.cols()); |
172 |
qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); |
239 |
m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixU); |
173 |
} |
240 |
} |
174 |
if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); |
241 |
if(svd.computeV()) svd.m_matrixV.setIdentity(matrix.cols(), matrix.cols()); |
175 |
return true; |
242 |
return true; |
176 |
} |
243 |
} |
177 |
return false; |
244 |
return false; |
178 |
} |
245 |
} |
|
|
246 |
private: |
247 |
HouseholderQR<MatrixType> m_qr; |
179 |
}; |
248 |
}; |
180 |
|
249 |
|
181 |
template<typename MatrixType> |
250 |
template<typename MatrixType> |
182 |
struct qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
251 |
class qr_preconditioner_impl<MatrixType, HouseholderQRPreconditioner, PreconditionIfMoreColsThanRows, true> |
183 |
{ |
252 |
{ |
184 |
static bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
253 |
public: |
|
|
254 |
typedef typename MatrixType::Index Index; |
255 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
256 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
257 |
TransposeTypeWithSameStorageOrder; |
258 |
|
259 |
void resize(Index rows, Index cols) |
260 |
{ |
261 |
if (cols != m_qr.rows() || rows != m_qr.cols()) |
262 |
{ |
263 |
m_qr = HouseholderQR<TransposeTypeWithSameStorageOrder>(cols, rows); |
264 |
} |
265 |
m_adjoint.resize(cols, rows); |
266 |
} |
267 |
|
268 |
bool run(JacobiSVD<MatrixType, HouseholderQRPreconditioner>& svd, const MatrixType& matrix) |
185 |
{ |
269 |
{ |
186 |
if(matrix.cols() > matrix.rows()) |
270 |
if(matrix.cols() > matrix.rows()) |
187 |
{ |
271 |
{ |
188 |
typedef Matrix<typename MatrixType::Scalar, MatrixType::ColsAtCompileTime, MatrixType::RowsAtCompileTime, |
272 |
m_adjoint = matrix.adjoint(); |
189 |
MatrixType::Options, MatrixType::MaxColsAtCompileTime, MatrixType::MaxRowsAtCompileTime> |
273 |
m_qr.compute(m_adjoint); |
190 |
TransposeTypeWithSameStorageOrder; |
274 |
|
191 |
HouseholderQR<TransposeTypeWithSameStorageOrder> qr(matrix.adjoint()); |
275 |
svd.m_workMatrix = m_qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
192 |
svd.m_workMatrix = qr.matrixQR().block(0,0,matrix.rows(),matrix.rows()).template triangularView<Upper>().adjoint(); |
276 |
if(svd.m_computeFullV) svd.m_matrixV = m_qr.householderQ(); |
193 |
if(svd.m_computeFullV) svd.m_matrixV = qr.householderQ(); |
|
|
194 |
else if(svd.m_computeThinV) { |
277 |
else if(svd.m_computeThinV) { |
195 |
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); |
278 |
svd.m_matrixV.setIdentity(matrix.cols(), matrix.rows()); |
196 |
qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); |
279 |
m_qr.householderQ().applyThisOnTheLeft(svd.m_matrixV); |
197 |
} |
280 |
} |
198 |
if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); |
281 |
if(svd.computeU()) svd.m_matrixU.setIdentity(matrix.rows(), matrix.rows()); |
199 |
return true; |
282 |
return true; |
200 |
} |
283 |
} |
201 |
else return false; |
284 |
else return false; |
202 |
} |
285 |
} |
|
|
286 |
|
287 |
private: |
288 |
HouseholderQR<TransposeTypeWithSameStorageOrder> m_qr; |
289 |
TransposeTypeWithSameStorageOrder m_adjoint; |
203 |
}; |
290 |
}; |
204 |
|
291 |
|
205 |
/*** 2x2 SVD implementation |
292 |
/*** 2x2 SVD implementation |
Lines 510-515
Link Here
|
510 |
friend struct internal::svd_precondition_2x2_block_to_be_real; |
597 |
friend struct internal::svd_precondition_2x2_block_to_be_real; |
511 |
template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> |
598 |
template<typename __MatrixType, int _QRPreconditioner, int _Case, bool _DoAnything> |
512 |
friend struct internal::qr_preconditioner_impl; |
599 |
friend struct internal::qr_preconditioner_impl; |
|
|
600 |
|
601 |
internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows> m_qr_precond_morecols; |
602 |
internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols> m_qr_precond_morerows; |
513 |
}; |
603 |
}; |
514 |
|
604 |
|
515 |
template<typename MatrixType, int QRPreconditioner> |
605 |
template<typename MatrixType, int QRPreconditioner> |
Lines 541-546
Link Here
|
541 |
: m_computeThinV ? m_diagSize |
631 |
: m_computeThinV ? m_diagSize |
542 |
: 0); |
632 |
: 0); |
543 |
m_workMatrix.resize(m_diagSize, m_diagSize); |
633 |
m_workMatrix.resize(m_diagSize, m_diagSize); |
|
|
634 |
m_qr_precond_morecols.resize(rows, cols); |
635 |
m_qr_precond_morerows.resize(rows, cols); |
544 |
} |
636 |
} |
545 |
|
637 |
|
546 |
template<typename MatrixType, int QRPreconditioner> |
638 |
template<typename MatrixType, int QRPreconditioner> |
Lines 554-562
Link Here
|
554 |
const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); |
646 |
const RealScalar precision = RealScalar(2) * NumTraits<Scalar>::epsilon(); |
555 |
|
647 |
|
556 |
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ |
648 |
/*** step 1. The R-SVD step: we use a QR decomposition to reduce to the case of a square matrix */ |
557 |
|
649 |
if(!m_qr_precond_morecols.run(*this, matrix) && !m_qr_precond_morerows.run(*this, matrix)) |
558 |
if(!internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreColsThanRows>::run(*this, matrix) |
|
|
559 |
&& !internal::qr_preconditioner_impl<MatrixType, QRPreconditioner, internal::PreconditionIfMoreRowsThanCols>::run(*this, matrix)) |
560 |
{ |
650 |
{ |
561 |
m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); |
651 |
m_workMatrix = matrix.block(0,0,m_diagSize,m_diagSize); |
562 |
if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); |
652 |
if(m_computeFullU) m_matrixU.setIdentity(m_rows,m_rows); |