Line data Source code
1 : /** \file
2 : * \brief Definitions of reusable thin-SVD row and column deletion updates.
3 : *
4 : * \ingroup gen_math_files
5 : */
6 :
7 : //***********************************************************************//
8 : // Copyright 2026 Jared R. Males (jaredmales@gmail.com)
9 : //
10 : // This file is part of mxlib.
11 : //
12 : // mxlib is free software: you can redistribute it and/or modify
13 : // it under the terms of the GNU General Public License as published by
14 : // the Free Software Foundation, either version 3 of the License, or
15 : // (at your option) any later version.
16 : //
17 : // mxlib is distributed in the hope that it will be useful,
18 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 : // GNU General Public License for more details.
21 : //
22 : // You should have received a copy of the GNU General Public License
23 : // along with mxlib. If not, see <http://www.gnu.org/licenses/>.
24 : //***********************************************************************//
25 :
26 : #include "math/svdDowndate.hpp"
27 :
28 : #include "math/floatUtils.hpp"
29 :
30 : #include <algorithm>
31 : #include <cmath>
32 : #include <cstddef>
33 : #include <cstring>
34 : #include <limits>
35 : #include <new>
36 : #include <stdexcept>
37 : #include <utility>
38 : #include <vector>
39 :
40 : namespace mx
41 : {
42 : namespace math
43 : {
44 :
45 : #ifdef MXLIB_SVD_DELETION_TEST_PRODUCER_ALIGNMENT
46 : static_assert( EIGEN_MAX_ALIGN_BYTES == MXLIB_SVD_DELETION_TEST_PRODUCER_ALIGNMENT,
47 : "The SVD-deletion test producer must use its configured dynamic alignment." );
48 : static_assert( EIGEN_MAX_STATIC_ALIGN_BYTES == MXLIB_SVD_DELETION_TEST_PRODUCER_ALIGNMENT,
49 : "The SVD-deletion test producer must use its configured static alignment." );
50 : #endif
51 :
52 : namespace
53 : {
54 :
55 : // mxlib-owned storage is deliberately unaligned so allocation and release do not depend on consumer Eigen settings.
56 : template <typename realT>
57 : using svdDeletionStorageMatrix =
58 : Eigen::Array<realT, Eigen::Dynamic, Eigen::Dynamic, Eigen::ColMajor | Eigen::DontAlign>;
59 :
60 : // mxlib-owned vector storage follows the same allocator-independent rule.
61 : template <typename realT>
62 : using svdDeletionStorageVector = Eigen::Array<realT, Eigen::Dynamic, 1, Eigen::DontAlign>;
63 :
64 : } // namespace
65 :
66 : template <typename realT, typename abiT>
67 : class svdDeletionResult<realT, abiT>::storage
68 : {
69 : public:
70 : /// Complete descending updated singular spectrum.
71 : svdDeletionStorageVector<realT> m_singularValues;
72 :
73 : /// Complete descending squared singular spectrum.
74 : svdDeletionStorageVector<realT> m_squaredSingularValues;
75 :
76 : /// Preserved-side rotation storage; only the leading output-rank columns are active.
77 : svdDeletionStorageMatrix<realT> m_rotation;
78 :
79 : /// Most recent operation status.
80 : svdDeletionStatus m_status{ svdDeletionStatus::notComputed };
81 :
82 : /// Backend that produced the current outputs.
83 : svdDeletionBackend m_backend{ svdDeletionBackend::stableCore };
84 :
85 : /// Prepared base rank.
86 : Eigen::Index m_baseRank{ 0 };
87 :
88 : /// Prepared output rank.
89 : Eigen::Index m_outputRank{ 0 };
90 :
91 : /// Allocated maximum output rank for the current base rank.
92 : Eigen::Index m_maximumOutputRank{ 0 };
93 :
94 : /// Number of clamped roundoff-scale negative eigenvalues.
95 : Eigen::Index m_clampedEigenvalues{ 0 };
96 :
97 : /// Smallest pre-clamp eigenvalue from the normalized backend PSD core.
98 : realT m_minimumPSDValue{ 0 };
99 :
100 : /// Underlying LAPACK failure code, or zero.
101 : MXLAPACK_INT m_lapackInfo{ 0 };
102 : };
103 :
104 : template <typename realT, typename abiT>
105 : class svdDeletionWorkspace<realT, abiT>::storage
106 : {
107 : public:
108 : /// Gathered deleted rows, with the active operation occupying the leading rows.
109 : svdDeletionStorageMatrix<realT> m_deletedRows;
110 :
111 : /// Base singular values normalized by their maximum to prevent overflow.
112 : svdDeletionStorageVector<realT> m_normalizedSingularValues;
113 :
114 : /// Scaled deleted-row transpose used by the symmetric core.
115 : svdDeletionStorageMatrix<realT> m_scaledDeletedTranspose;
116 :
117 : /// Symmetric covariance or complement matrix overwritten by SYEVR.
118 : svdDeletionStorageMatrix<realT> m_psdCore;
119 :
120 : /// Full eigenvectors returned by SYEVR.
121 : svdDeletionStorageMatrix<realT> m_eigenvectors;
122 :
123 : /// Full ascending eigenvalues returned by SYEVR.
124 : svdDeletionStorageVector<realT> m_eigenvalues;
125 :
126 : /// Complement square root used by the stable backend.
127 : svdDeletionStorageMatrix<realT> m_complementRoot;
128 :
129 : /// Complement-preserving rectangular core overwritten by GESVD.
130 : svdDeletionStorageMatrix<realT> m_stableCore;
131 :
132 : /// Right singular vectors returned row-wise by GESVD.
133 : svdDeletionStorageMatrix<realT> m_rightTranspose;
134 :
135 : /// Singular values returned by GESVD.
136 : svdDeletionStorageVector<realT> m_stableSingularValues;
137 :
138 : /// Mutable diagonal poles used by the structured rank-one deflation pass.
139 : svdDeletionStorageVector<realT> m_secularPoles;
140 :
141 : /// Mutable normalized update components used by structured deflation.
142 : svdDeletionStorageVector<realT> m_secularUpdate;
143 :
144 : /// Compact active poles supplied to LAED9.
145 : svdDeletionStorageVector<realT> m_secularActivePoles;
146 :
147 : /// Compact active update supplied to and overwritten by LAED9.
148 : svdDeletionStorageVector<realT> m_secularActiveUpdate;
149 :
150 : /// Merged active and deflated eigenvalue candidates before sign reversal.
151 : svdDeletionStorageVector<realT> m_secularCandidateValues;
152 :
153 : /// One embedded structured eigenvector transformed back through deflation rotations.
154 : svdDeletionStorageVector<realT> m_secularVector;
155 :
156 : /// Active coordinate indices after structured deflation.
157 : std::vector<Eigen::Index> m_secularActiveIndices;
158 :
159 : /// Deflated coordinate indices after structured deflation.
160 : std::vector<Eigen::Index> m_secularDeflatedIndices;
161 :
162 : /// Merge order for active secular roots and deflated diagonal values.
163 : std::vector<Eigen::Index> m_secularOrder;
164 :
165 : /// First coordinate of each recorded structured-deflation Givens rotation.
166 : std::vector<Eigen::Index> m_secularRotationFirst;
167 :
168 : /// Second coordinate of each recorded structured-deflation Givens rotation.
169 : std::vector<Eigen::Index> m_secularRotationSecond;
170 :
171 : /// Cosine of each recorded structured-deflation Givens rotation.
172 : std::vector<realT> m_secularRotationCosine;
173 :
174 : /// Sine of each recorded structured-deflation Givens rotation.
175 : std::vector<realT> m_secularRotationSine;
176 :
177 : /// Eigenvector support indices required by SYEVR.
178 : std::vector<MXLAPACK_INT> m_support;
179 :
180 : /// Floating-point workspace required by SYEVR.
181 : std::vector<realT> m_syevrWork;
182 :
183 : /// Integer workspace required by SYEVR.
184 : std::vector<MXLAPACK_INT> m_syevrIWork;
185 :
186 : /// Floating-point workspace required by GESVD.
187 : std::vector<realT> m_gesvdWork;
188 :
189 : /// Prepared base rank.
190 : Eigen::Index m_baseRank{ 0 };
191 :
192 : /// Prepared maximum deletion count.
193 : Eigen::Index m_maximumDeleted{ 0 };
194 :
195 : /// Prepared numerical backend.
196 : svdDeletionBackend m_backend{ svdDeletionBackend::stableCore };
197 :
198 : /// True after successful dimensioning and LAPACK workspace queries.
199 : bool m_prepared{ false };
200 :
201 : /// Underlying LAPACK workspace-query failure code, or zero.
202 : MXLAPACK_INT m_lapackInfo{ 0 };
203 : };
204 :
205 : namespace detail
206 : {
207 :
208 : namespace
209 : {
210 :
211 : // Return true when an Eigen dimension can be passed to the configured LAPACK integer ABI.
212 582 : bool validLapackDimension( Eigen::Index dimension /* [in] candidate dimension */ )
213 : {
214 1164 : return dimension >= 0 && static_cast<unsigned long long>( dimension ) <=
215 1164 : static_cast<unsigned long long>( std::numeric_limits<MXLAPACK_INT>::max() );
216 : }
217 :
218 : // Return true when an ABI-v2 signed dimension fits Eigen's configured index type.
219 1384 : bool validAbiDimension( std::int64_t dimension /* [in] ABI-v2 dimension */ )
220 : {
221 2764 : return dimension >= 0 && static_cast<std::uint64_t>( dimension ) <=
222 2764 : static_cast<std::uint64_t>( std::numeric_limits<Eigen::Index>::max() );
223 : }
224 :
225 : // Return true when a contiguous ABI-v2 vector descriptor can be mapped safely.
226 : template <typename realT>
227 135 : bool validAbiVectorView( svdDeletionConstVectorViewV2<realT> view /* [in] borrowed vector descriptor */ )
228 : {
229 270 : return validAbiDimension( view.size ) && ( view.size == 0 || view.data != nullptr ) &&
230 404 : ( view.size == 0 || reinterpret_cast<std::uintptr_t>( view.data ) % alignof( realT ) == 0 ) &&
231 134 : static_cast<std::uint64_t>( view.size ) <=
232 269 : static_cast<std::uint64_t>( std::numeric_limits<std::ptrdiff_t>::max() ) / sizeof( realT );
233 : }
234 :
235 : // Return true when a column-major ABI-v2 matrix descriptor can be mapped safely.
236 : template <typename realT>
237 153 : bool validAbiMatrixView( svdDeletionConstMatrixViewV2<realT> view /* [in] borrowed matrix descriptor */ )
238 : {
239 305 : if( !validAbiDimension( view.rows ) || !validAbiDimension( view.columns ) ||
240 152 : !validAbiDimension( view.outerStride ) )
241 : {
242 1 : return false;
243 : }
244 152 : if( view.rows > 0 && view.columns > 0 && view.data == nullptr )
245 : {
246 2 : return false;
247 : }
248 150 : if( view.rows > 0 && view.columns > 0 && reinterpret_cast<std::uintptr_t>( view.data ) % alignof( realT ) != 0 )
249 : {
250 1 : return false;
251 : }
252 149 : if( view.columns > 1 && view.outerStride < view.rows )
253 : {
254 1 : return false;
255 : }
256 148 : if( view.rows == 0 || view.columns == 0 )
257 : {
258 7 : return true;
259 : }
260 :
261 141 : const std::uint64_t rows = static_cast<std::uint64_t>( view.rows );
262 141 : const std::uint64_t columns = static_cast<std::uint64_t>( view.columns );
263 141 : const std::uint64_t outerStride = static_cast<std::uint64_t>( view.outerStride );
264 141 : const std::uint64_t maximumScalarOffset =
265 : static_cast<std::uint64_t>( std::numeric_limits<std::ptrdiff_t>::max() ) / sizeof( realT );
266 141 : if( rows - 1 > maximumScalarOffset )
267 : {
268 1 : return false;
269 : }
270 140 : if( columns == 1 )
271 : {
272 10 : return true;
273 : }
274 130 : return columns - 1 <= ( maximumScalarOffset - ( rows - 1 ) ) / outerStride;
275 : }
276 :
277 : // Return true when an ABI-v2 signed-index descriptor can be decoded safely.
278 139 : bool validAbiIndexView( svdDeletionConstIndexViewV2 view /* [in] borrowed index descriptor */ )
279 : {
280 139 : if( !validAbiDimension( view.size ) || ( view.size > 0 && view.data == nullptr ) )
281 : {
282 2 : return false;
283 : }
284 137 : if( view.size == 0 )
285 : {
286 11 : return true;
287 : }
288 126 : const bool supportedWidth =
289 126 : view.elementBytes == 1 || view.elementBytes == 2 || view.elementBytes == 4 || view.elementBytes == 8;
290 126 : if( !supportedWidth )
291 : {
292 1 : return false;
293 : }
294 125 : return static_cast<std::uint64_t>( view.size - 1 ) <=
295 125 : static_cast<std::uint64_t>( std::numeric_limits<std::ptrdiff_t>::max() ) /
296 125 : static_cast<std::uint64_t>( view.elementBytes );
297 : }
298 :
299 : // Decode one signed integer from an ABI-v2 index descriptor without type-punning across the boundary.
300 167 : std::int64_t abiIndexAt( svdDeletionConstIndexViewV2 view, /* [in] borrowed index descriptor */
301 : std::int64_t offset /* [in] zero-based element offset */ )
302 : {
303 167 : const std::ptrdiff_t byteOffset = static_cast<std::ptrdiff_t>( offset * view.elementBytes );
304 167 : const auto *bytes = static_cast<const unsigned char *>( view.data ) + byteOffset;
305 167 : switch( view.elementBytes )
306 : {
307 2 : case 1:
308 : {
309 2 : std::int8_t value{ 0 };
310 2 : std::memcpy( &value, bytes, sizeof( value ) );
311 2 : return value;
312 : }
313 2 : case 2:
314 : {
315 2 : std::int16_t value{ 0 };
316 2 : std::memcpy( &value, bytes, sizeof( value ) );
317 2 : return value;
318 : }
319 2 : case 4:
320 : {
321 2 : std::int32_t value{ 0 };
322 2 : std::memcpy( &value, bytes, sizeof( value ) );
323 2 : return value;
324 : }
325 161 : case 8:
326 : {
327 161 : std::int64_t value{ 0 };
328 161 : std::memcpy( &value, bytes, sizeof( value ) );
329 161 : return value;
330 : }
331 : }
332 : return -1; // LCOV_EXCL_LINE -- callers validate the element width before decoding.
333 : }
334 :
335 : // Unaligned Eigen map used only inside the mxlib implementation after ABI-v2 descriptor validation.
336 : template <typename realT>
337 : using svdDeletionAbiMatrixMap =
338 : Eigen::Map<const svdDeletionMatrix<realT>, Eigen::Unaligned, Eigen::OuterStride<Eigen::Dynamic>>;
339 :
340 : // Construct a matrix map after validAbiMatrixView has accepted the descriptor.
341 : template <typename realT>
342 : svdDeletionAbiMatrixMap<realT>
343 141 : mapAbiMatrix( svdDeletionConstMatrixViewV2<realT> view /* [in] validated matrix descriptor */ )
344 : {
345 : return svdDeletionAbiMatrixMap<realT>(
346 141 : view.data,
347 141 : static_cast<Eigen::Index>( view.rows ),
348 141 : static_cast<Eigen::Index>( view.columns ),
349 141 : Eigen::OuterStride<Eigen::Dynamic>( static_cast<Eigen::Index>( view.outerStride ) ) );
350 : }
351 :
352 : // Construct a contiguous vector map after validAbiVectorView has accepted the descriptor.
353 : template <typename realT>
354 : Eigen::Map<const svdDeletionVector<realT>, Eigen::Unaligned>
355 127 : mapAbiVector( svdDeletionConstVectorViewV2<realT> view /* [in] validated vector descriptor */ )
356 : {
357 254 : return Eigen::Map<const svdDeletionVector<realT>, Eigen::Unaligned>( view.data,
358 254 : static_cast<Eigen::Index>( view.size ) );
359 : }
360 :
361 : // Return true when a dense array shape fits Eigen indexing and addressable storage.
362 : template <typename valueT>
363 802 : bool validArrayShape( Eigen::Index rows, /* [in] candidate row count */
364 : Eigen::Index cols /* [in] candidate column count */ )
365 : {
366 802 : if( rows < 0 || cols < 0 || ( rows > 0 && cols > std::numeric_limits<Eigen::Index>::max() / rows ) )
367 : {
368 : return false; // LCOV_EXCL_LINE -- requires dimensions wider than the configured LAPACK integer ABI.
369 : }
370 802 : const Eigen::Index elements = rows * cols;
371 802 : return static_cast<unsigned long long>( elements ) <=
372 802 : static_cast<unsigned long long>( std::numeric_limits<std::size_t>::max() / sizeof( valueT ) );
373 : }
374 :
375 : // Return true when an element count fits addressable `std::vector` storage.
376 : template <typename valueT>
377 344 : bool validVectorSize( Eigen::Index elements /* [in] candidate element count */ )
378 : {
379 688 : return elements >= 0 &&
380 344 : static_cast<unsigned long long>( elements ) <=
381 688 : static_cast<unsigned long long>( std::numeric_limits<std::size_t>::max() / sizeof( valueT ) );
382 : }
383 :
384 : // Return true when every dense Eigen element is finite under mxlib's fast-math-safe predicate.
385 : template <typename derivedT>
386 315 : bool allFinite( const Eigen::DenseBase<derivedT> &values /* [in] dense object to inspect */ )
387 : {
388 1729 : for( Eigen::Index index = 0; index < values.size(); ++index )
389 : {
390 1422 : if( !math::isFinite( values.derived().data()[index] ) )
391 : {
392 8 : return false;
393 : }
394 : }
395 307 : return true;
396 : }
397 :
398 : // Return true when the leading values are finite, nonnegative, and nonincreasing.
399 : template <typename derivedT>
400 30 : bool validDescendingSpectrum( const Eigen::DenseBase<derivedT> &values /* [in] spectrum to inspect */ )
401 : {
402 118 : for( Eigen::Index index = 0; index < values.size(); ++index )
403 : {
404 90 : const auto value = values.derived()( index );
405 90 : if( !math::isFinite( value ) || value < 0 || ( index > 0 && value > values.derived()( index - 1 ) ) )
406 : {
407 2 : return false;
408 : }
409 : }
410 28 : return true;
411 : }
412 :
413 : // Return true when the leading values are finite and nondecreasing.
414 : template <typename derivedT>
415 78 : bool validAscendingSpectrum( const Eigen::DenseBase<derivedT> &values /* [in] spectrum to inspect */ )
416 : {
417 240 : for( Eigen::Index index = 0; index < values.size(); ++index )
418 : {
419 165 : const auto value = values.derived()( index );
420 165 : if( !math::isFinite( value ) || ( index > 0 && value < values.derived()( index - 1 ) ) )
421 : {
422 3 : return false;
423 : }
424 : }
425 75 : return true;
426 : }
427 :
428 : // Invoke an optional deterministic operation hook.
429 : template <typename realT>
430 205 : void callOperationHook( svdDeletionTestOperation operation /* [in] operation about to execute */ )
431 : {
432 205 : auto &hooks = svdDeletionHooks<realT>();
433 205 : if( hooks.operation )
434 : {
435 15 : hooks.operation( operation );
436 : }
437 190 : }
438 :
439 : // Invoke the injected or production SYEVR implementation.
440 : template <typename realT>
441 138 : MXLAPACK_INT callSyevr( char jobz, /* [in] eigenvector request */
442 : char range, /* [in] eigenvalue selection mode */
443 : char uplo, /* [in] populated input triangle */
444 : MXLAPACK_INT n, /* [in] matrix order */
445 : realT *matrix, /* [in,out] symmetric input matrix */
446 : MXLAPACK_INT lda, /* [in] input leading dimension */
447 : realT valueLower, /* [in] lower value bound */
448 : realT valueUpper, /* [in] upper value bound */
449 : MXLAPACK_INT indexLow, /* [in] one-based lower index */
450 : MXLAPACK_INT indexHigh, /* [in] one-based upper index */
451 : realT tolerance, /* [in] convergence tolerance */
452 : MXLAPACK_INT *found, /* [out] selected eigenvalue count */
453 : realT *eigenvalues, /* [out] ascending eigenvalues */
454 : realT *eigenvectors, /* [out] eigenvectors in columns */
455 : MXLAPACK_INT ldz, /* [in] eigenvector leading dimension */
456 : MXLAPACK_INT *support, /* [out] eigenvector support */
457 : realT *work, /* [in,out] floating workspace */
458 : MXLAPACK_INT lwork, /* [in] floating workspace size */
459 : MXLAPACK_INT *iwork, /* [in,out] integer workspace */
460 : MXLAPACK_INT liwork /* [in] integer workspace size */ )
461 : {
462 138 : auto &hooks = svdDeletionHooks<realT>();
463 138 : if( hooks.syevr )
464 : {
465 18 : return hooks.syevr( jobz,
466 : range,
467 : uplo,
468 : n,
469 : matrix,
470 : lda,
471 : valueLower,
472 : valueUpper,
473 : indexLow,
474 : indexHigh,
475 : tolerance,
476 : found,
477 : eigenvalues,
478 : eigenvectors,
479 : ldz,
480 : support,
481 : work,
482 : lwork,
483 : iwork,
484 18 : liwork );
485 : }
486 :
487 120 : return math::syevr<realT>( jobz,
488 : range,
489 : uplo,
490 : n,
491 : matrix,
492 : lda,
493 : valueLower,
494 : valueUpper,
495 : indexLow,
496 : indexHigh,
497 : tolerance,
498 : found,
499 : eigenvalues,
500 : eigenvectors,
501 : ldz,
502 : support,
503 : work,
504 : lwork,
505 : iwork,
506 120 : liwork );
507 : }
508 :
509 : // Invoke the injected or production GESVD implementation.
510 : template <typename realT>
511 68 : MXLAPACK_INT callGesvd( char jobu, /* [in] left singular-vector request */
512 : char jobvt, /* [in] right singular-vector request */
513 : MXLAPACK_INT rows, /* [in] matrix row count */
514 : MXLAPACK_INT cols, /* [in] matrix column count */
515 : realT *matrix, /* [in,out] input matrix */
516 : MXLAPACK_INT lda, /* [in] input leading dimension */
517 : realT *singular, /* [out] descending singular values */
518 : realT *left, /* [out] left singular vectors when requested */
519 : MXLAPACK_INT ldu, /* [in] left-vector leading dimension */
520 : realT *rightT, /* [out] transposed right singular vectors */
521 : MXLAPACK_INT ldvt, /* [in] right-vector leading dimension */
522 : realT *work, /* [in,out] floating workspace */
523 : MXLAPACK_INT lwork /* [in] floating workspace size */ )
524 : {
525 68 : auto &hooks = svdDeletionHooks<realT>();
526 68 : if( hooks.gesvd )
527 : {
528 8 : return hooks.gesvd( jobu, jobvt, rows, cols, matrix, lda, singular, left, ldu, rightT, ldvt, work, lwork );
529 : }
530 :
531 60 : return math::gesvd<realT>( jobu, jobvt, rows, cols, matrix, lda, singular, left, ldu, rightT, ldvt, work, lwork );
532 : }
533 :
534 : // Invoke the injected or production LAED9 structured rank-one eigensolver.
535 : template <typename realT>
536 19 : MXLAPACK_INT callLaed9( realT *eigenvalues, /* [out] ascending updated eigenvalues */
537 : realT *delta, /* [out] LAED9 secular workspace */
538 : realT *eigenvectors, /* [out] updated eigenvectors in columns */
539 : MXLAPACK_INT rank, /* [in] active secular-system dimension */
540 : MXLAPACK_INT leadingDimension, /* [in] matrix leading dimension */
541 : realT rho, /* [in] positive rank-one update weight */
542 : realT *poles, /* [in,out] strictly ascending diagonal poles */
543 : realT *update /* [in,out] normalized update components */ )
544 : {
545 19 : auto &hooks = svdDeletionHooks<realT>();
546 19 : if( hooks.laed9 )
547 : {
548 8 : return hooks.laed9( eigenvalues, delta, eigenvectors, rank, leadingDimension, rho, poles, update );
549 : }
550 :
551 11 : return math::laed9<realT>( eigenvalues,
552 : delta,
553 : eigenvectors,
554 : rank,
555 : 1,
556 : rank,
557 : rank,
558 : leadingDimension,
559 : rho,
560 : poles,
561 : update,
562 11 : leadingDimension );
563 : }
564 :
565 : // Convert a finite LAPACK floating workspace query to an integer size.
566 : template <typename realT>
567 101 : bool querySize( MXLAPACK_INT &size, /* [out] validated integer workspace size */
568 : realT query /* [in] floating-point LAPACK query result */ )
569 : {
570 200 : if( !math::isFinite( query ) || query < realT( 1 ) ||
571 99 : static_cast<long double>( query ) > static_cast<long double>( std::numeric_limits<MXLAPACK_INT>::max() ) )
572 : {
573 2 : return false;
574 : }
575 :
576 99 : size = static_cast<MXLAPACK_INT>( std::ceil( query ) );
577 99 : return size >= 1;
578 : }
579 :
580 : // Return the construction-scale tolerance for a theoretically PSD core.
581 : template <typename realT>
582 74 : realT psdTolerance( Eigen::Index dimension, /* [in] core dimension */
583 : realT scale /* [in] nonnegative core construction scale */ )
584 : {
585 74 : return realT( 64 ) * std::numeric_limits<realT>::epsilon() *
586 74 : static_cast<realT>( std::max<Eigen::Index>( 1, dimension ) ) * scale;
587 : }
588 :
589 : } // namespace
590 :
591 : template <typename realT>
592 612 : svdDeletionTestHooks<realT> &svdDeletionHooks()
593 : {
594 : static svdDeletionTestHooks<realT> hooks;
595 612 : return hooks;
596 : }
597 :
598 : // Implementation access shared by the public deletion entry points.
599 : template <typename realT>
600 : struct svdDeletionImplementation
601 : {
602 : using resultT = svdDeletionResult<realT>;
603 : using workspaceT = svdDeletionWorkspace<realT>;
604 :
605 : // Publish a failure and its optional underlying LAPACK status.
606 59 : static svdDeletionStatus fail( resultT &result, /* [out] result status to invalidate */
607 : svdDeletionStatus status, /* [in] failure status */
608 : MXLAPACK_INT lapackInfo = 0 /* [in] underlying LAPACK status */ )
609 : {
610 59 : if( !result.ensureStorage() )
611 : {
612 1 : return svdDeletionStatus::allocationFailure;
613 : }
614 58 : result.m_storage->m_status = status;
615 58 : result.m_storage->m_lapackInfo = lapackInfo;
616 58 : result.m_storage->m_clampedEigenvalues = 0;
617 58 : if( status != svdDeletionStatus::nonPositiveSemidefinite )
618 : {
619 55 : result.m_storage->m_minimumPSDValue = realT( 0 );
620 : }
621 58 : return status;
622 : }
623 :
624 : // Return true for a supported numerical backend value.
625 198 : static bool validBackend( svdDeletionBackend backend /* [in] backend value to validate */ )
626 : {
627 198 : return backend == svdDeletionBackend::leadingCovariance || backend == svdDeletionBackend::stableCore ||
628 198 : backend == svdDeletionBackend::rankOneSecular;
629 : }
630 :
631 : // Validate singular values and requested rank without checking factor rows.
632 127 : static bool validSingularInputs( svdDeletionConstVectorRef<realT> singularValues, /* [in] base singular values */
633 : Eigen::Index outputRank /* [in] requested output rank */ )
634 : {
635 127 : const Eigen::Index rank = singularValues.size();
636 251 : if( rank <= 0 || !validLapackDimension( rank ) || outputRank <= 0 || outputRank > rank ||
637 124 : !allFinite( singularValues ) )
638 : {
639 4 : return false;
640 : }
641 :
642 462 : for( Eigen::Index index = 0; index < rank; ++index )
643 : {
644 558 : if( singularValues( index ) < realT( 0 ) ||
645 217 : ( index > 0 && singularValues( index ) > singularValues( index - 1 ) ) )
646 : {
647 2 : return false;
648 : }
649 : }
650 121 : return true;
651 : }
652 :
653 : // Validate common core inputs without checking global factor orthonormality.
654 57 : static bool validCoreInputs( svdDeletionConstVectorRef<realT> singularValues, /* [in] base singular values */
655 : svdDeletionConstMatrixRef<realT> deletedRows, /* [in] deleted factor rows */
656 : Eigen::Index outputRank /* [in] requested output rank */ )
657 : {
658 57 : const Eigen::Index rank = singularValues.size();
659 164 : return validSingularInputs( singularValues, outputRank ) && deletedRows.cols() == rank &&
660 52 : validLapackDimension( deletedRows.rows() ) &&
661 169 : deletedRows.rows() <= std::numeric_limits<MXLAPACK_INT>::max() - rank && allFinite( deletedRows );
662 : }
663 :
664 : // Normalize base singular values and return their common scale.
665 94 : static realT normalizeSingularValues( workspaceT &workspace, /* [out] normalized values */
666 : svdDeletionConstVectorRef<realT> singularValues /* [in] base values */ )
667 : {
668 94 : const realT scale = singularValues( 0 );
669 94 : if( scale == realT( 0 ) )
670 : {
671 3 : workspace.m_storage->m_normalizedSingularValues.setZero();
672 3 : return scale;
673 : }
674 :
675 91 : int scaleExponent{ 0 };
676 91 : const realT scaleMantissa = std::frexp( scale, &scaleExponent );
677 365 : for( Eigen::Index index = 0; index < singularValues.size(); ++index )
678 : {
679 274 : int valueExponent{ 0 };
680 274 : const realT valueMantissa = std::frexp( singularValues( index ), &valueExponent );
681 274 : workspace.m_storage->m_normalizedSingularValues( index ) =
682 274 : std::ldexp( valueMantissa / scaleMantissa, valueExponent - scaleExponent );
683 : }
684 91 : return scale;
685 : }
686 :
687 : // Return a ratio without forming an underflowed reciprocal under fast-math builds.
688 5 : static realT normalizedRatio( realT value, /* [in] nonnegative numerator */
689 : realT scale /* [in] positive denominator */ )
690 : {
691 5 : int scaleExponent{ 0 };
692 5 : int valueExponent{ 0 };
693 5 : const realT scaleMantissa = std::frexp( scale, &scaleExponent );
694 5 : const realT valueMantissa = std::frexp( value, &valueExponent );
695 5 : return std::ldexp( valueMantissa / scaleMantissa, valueExponent - scaleExponent );
696 : }
697 :
698 : // Rescale one normalized singular value without squaring before restoring input scale.
699 196 : static bool rescaleSingularValue( realT &singularValue, /* [out] rescaled singular value */
700 : realT &squaredSingularValue, /* [out] rescaled squared value */
701 : realT normalizedValue, /* [in] normalized singular value */
702 : realT scale /* [in] base singular-value scale */ )
703 : {
704 196 : const realT maximum = std::numeric_limits<realT>::max();
705 196 : singularValue = normalizedValue * scale;
706 196 : if( !math::isFinite( singularValue ) || singularValue > std::sqrt( maximum ) )
707 : {
708 3 : return false;
709 : }
710 193 : squaredSingularValue = singularValue * singularValue;
711 193 : return math::isFinite( squaredSingularValue );
712 : }
713 :
714 : // Rescale one normalized squared singular value without overflowing either public representation.
715 111 : static bool rescaleValue( realT &singularValue, /* [out] rescaled singular value */
716 : realT &squaredSingularValue, /* [out] rescaled squared singular value */
717 : realT normalizedSquaredValue, /* [in] normalized squared value */
718 : realT scale /* [in] base singular-value scale */ )
719 : {
720 111 : return rescaleSingularValue( singularValue, squaredSingularValue, std::sqrt( normalizedSquaredValue ), scale );
721 : }
722 :
723 : // Ensure result and worker storage can execute one core operation.
724 100 : static svdDeletionStatus prepare( resultT &result, /* [in,out] output storage */
725 : workspaceT &workspace, /* [in,out] reusable scratch */
726 : Eigen::Index rank, /* [in] base rank */
727 : Eigen::Index deletedCount, /* [in] actual deleted-row count */
728 : Eigen::Index outputRank, /* [in] requested output rank */
729 : svdDeletionBackend backend /* [in] selected backend */ )
730 : {
731 100 : svdDeletionStatus status = result.prepare( rank, outputRank );
732 100 : if( status != svdDeletionStatus::success )
733 : {
734 3 : return status;
735 : }
736 :
737 129 : if( !workspace.prepared() || workspace.baseRank() != rank || workspace.maximumDeleted() < deletedCount ||
738 32 : workspace.backend() != backend )
739 : {
740 69 : status = workspace.prepare( rank, deletedCount, backend );
741 69 : if( status != svdDeletionStatus::success )
742 : {
743 2 : return fail( result, status, workspace.lapackInfo() );
744 : }
745 : }
746 :
747 95 : return svdDeletionStatus::success;
748 : }
749 :
750 : // Copy deleted factor rows into the active portion of a prepared workspace.
751 42 : static void loadDeletedRows( workspaceT &workspace, /* [out] gathered workspace rows */
752 : svdDeletionConstMatrixRef<realT> deletedRows /* [in] rows to copy */ )
753 : {
754 42 : if( deletedRows.rows() > 0 )
755 : {
756 42 : workspace.m_storage->m_deletedRows.matrix().topRows( deletedRows.rows() ) = deletedRows.matrix();
757 : }
758 42 : }
759 :
760 : // Publish the unchanged represented singular system for an empty deletion.
761 11 : static svdDeletionStatus identity( resultT &result, /* [out] identity update */
762 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base singular values */
763 : Eigen::Index outputRank, /* [in] requested output rank */
764 : svdDeletionBackend backend /* [in] reported backend */ )
765 : {
766 11 : const svdDeletionStatus status = result.prepare( singularValues.size(), outputRank );
767 11 : if( status != svdDeletionStatus::success )
768 : {
769 1 : return status;
770 : }
771 :
772 10 : result.m_storage->m_rotation.setZero();
773 30 : for( Eigen::Index index = 0; index < singularValues.size(); ++index )
774 : {
775 22 : result.m_storage->m_singularValues( index ) = singularValues( index );
776 44 : if( singularValues( index ) > realT( 0 ) &&
777 22 : singularValues( index ) > std::sqrt( std::numeric_limits<realT>::max() ) )
778 : {
779 2 : return fail( result, svdDeletionStatus::rescalingOverflow );
780 : }
781 20 : result.m_storage->m_squaredSingularValues( index ) = singularValues( index ) * singularValues( index );
782 20 : if( index < outputRank )
783 : {
784 18 : result.m_storage->m_rotation( index, index ) = realT( 1 );
785 : }
786 : }
787 8 : result.m_storage->m_backend = backend;
788 8 : result.m_storage->m_status = svdDeletionStatus::success;
789 8 : result.m_storage->m_clampedEigenvalues = 0;
790 8 : if( backend == svdDeletionBackend::leadingCovariance || backend == svdDeletionBackend::rankOneSecular )
791 : {
792 5 : const realT scale = singularValues( 0 );
793 5 : realT minimumNormalized{ 0 };
794 5 : if( scale != realT( 0 ) )
795 : {
796 5 : minimumNormalized = normalizedRatio( singularValues( singularValues.size() - 1 ), scale );
797 : }
798 5 : result.m_storage->m_minimumPSDValue = minimumNormalized * minimumNormalized;
799 5 : }
800 : else
801 : {
802 3 : result.m_storage->m_minimumPSDValue = realT( 1 );
803 : }
804 8 : result.m_storage->m_lapackInfo = 0;
805 8 : return result.m_storage->m_status;
806 : }
807 :
808 : // Publish an arbitrary orthonormal system for an identically zero represented matrix.
809 3 : static svdDeletionStatus zeroSystem( resultT &result, /* [out] zero-valued result */
810 : Eigen::Index outputRank, /* [in] requested output rank */
811 : svdDeletionBackend backend /* [in] reported backend */ )
812 : {
813 3 : result.m_storage->m_rotation.setZero();
814 3 : result.m_storage->m_singularValues.setZero();
815 3 : result.m_storage->m_squaredSingularValues.setZero();
816 9 : for( Eigen::Index index = 0; index < outputRank; ++index )
817 : {
818 6 : result.m_storage->m_rotation( index, index ) = realT( 1 );
819 : }
820 3 : result.m_storage->m_backend = backend;
821 3 : result.m_storage->m_status = svdDeletionStatus::success;
822 3 : result.m_storage->m_clampedEigenvalues = 0;
823 3 : result.m_storage->m_minimumPSDValue = realT( 0 );
824 3 : result.m_storage->m_lapackInfo = 0;
825 3 : return result.m_storage->m_status;
826 : }
827 :
828 : // Execute a normalized symmetric covariance solve with deleted rows already loaded.
829 32 : static svdDeletionStatus leadingLoaded( resultT &result, /* [out] updated system */
830 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
831 : Eigen::Index deletedCount, /* [in] active rows */
832 : Eigen::Index outputRank, /* [in] requested rank */
833 : workspaceT &workspace /* [in,out] scratch */ )
834 : {
835 32 : const Eigen::Index rank = singularValues.size();
836 32 : const realT scale = normalizeSingularValues( workspace, singularValues );
837 32 : if( scale == realT( 0 ) )
838 : {
839 1 : return zeroSystem( result, outputRank, svdDeletionBackend::leadingCovariance );
840 : }
841 :
842 77 : for( Eigen::Index column = 0; column < deletedCount; ++column )
843 : {
844 191 : for( Eigen::Index row = 0; row < rank; ++row )
845 : {
846 145 : workspace.m_storage->m_scaledDeletedTranspose( row, column ) =
847 145 : workspace.m_storage->m_normalizedSingularValues( row ) *
848 145 : workspace.m_storage->m_deletedRows( column, row );
849 : }
850 : }
851 :
852 31 : const auto scaledDeleted = workspace.m_storage->m_scaledDeletedTranspose.matrix().leftCols( deletedCount );
853 31 : workspace.m_storage->m_psdCore.matrix().noalias() = -scaledDeleted * scaledDeleted.transpose();
854 124 : for( Eigen::Index index = 0; index < rank; ++index )
855 : {
856 186 : workspace.m_storage->m_psdCore( index, index ) += workspace.m_storage->m_normalizedSingularValues( index ) *
857 93 : workspace.m_storage->m_normalizedSingularValues( index );
858 : }
859 31 : MXLAPACK_INT found{ 0 };
860 31 : const MXLAPACK_INT lapackRank = static_cast<MXLAPACK_INT>( rank );
861 : const MXLAPACK_INT info =
862 186 : callSyevr<realT>( 'V',
863 : 'A',
864 : 'L',
865 : lapackRank,
866 31 : workspace.m_storage->m_psdCore.data(),
867 : lapackRank,
868 : realT( 0 ),
869 : realT( 0 ),
870 : 1,
871 : lapackRank,
872 : math::lamch<realT>( 'S' ),
873 : &found,
874 31 : workspace.m_storage->m_eigenvalues.data(),
875 31 : workspace.m_storage->m_eigenvectors.data(),
876 : lapackRank,
877 31 : workspace.m_storage->m_support.data(),
878 31 : workspace.m_storage->m_syevrWork.data(),
879 31 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_syevrWork.size() ),
880 31 : workspace.m_storage->m_syevrIWork.data(),
881 31 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_syevrIWork.size() ) );
882 31 : if( info != 0 || found != lapackRank )
883 : {
884 2 : return fail( result, svdDeletionStatus::solverFailure, info );
885 : }
886 29 : if( !allFinite( workspace.m_storage->m_eigenvalues ) || !allFinite( workspace.m_storage->m_eigenvectors ) )
887 : {
888 2 : return fail( result, svdDeletionStatus::nonFiniteOutput );
889 : }
890 27 : if( !validAscendingSpectrum( workspace.m_storage->m_eigenvalues ) )
891 : {
892 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
893 : }
894 :
895 26 : result.m_storage->m_minimumPSDValue = workspace.m_storage->m_eigenvalues( 0 );
896 : const realT tolerance =
897 26 : psdTolerance<realT>( rank, realT( 1 ) + static_cast<realT>( scaledDeleted.squaredNorm() ) );
898 26 : if( result.m_storage->m_minimumPSDValue < -tolerance )
899 : {
900 1 : return fail( result, svdDeletionStatus::nonPositiveSemidefinite );
901 : }
902 :
903 25 : result.m_storage->m_clampedEigenvalues = 0;
904 100 : for( Eigen::Index index = 0; index < rank; ++index )
905 : {
906 75 : if( workspace.m_storage->m_eigenvalues( index ) < realT( 0 ) )
907 : {
908 1 : workspace.m_storage->m_eigenvalues( index ) = realT( 0 );
909 1 : ++result.m_storage->m_clampedEigenvalues;
910 : }
911 : }
912 :
913 99 : for( Eigen::Index output = 0; output < rank; ++output )
914 : {
915 75 : const Eigen::Index source = rank - 1 - output;
916 75 : if( !rescaleValue( result.m_storage->m_singularValues( output ),
917 75 : result.m_storage->m_squaredSingularValues( output ),
918 75 : workspace.m_storage->m_eigenvalues( source ),
919 : scale ) )
920 : {
921 1 : return fail( result, svdDeletionStatus::rescalingOverflow );
922 : }
923 74 : if( output < outputRank )
924 : {
925 74 : result.m_storage->m_rotation.matrix().col( output ) =
926 148 : workspace.m_storage->m_eigenvectors.matrix().col( source );
927 : }
928 : }
929 :
930 24 : result.m_storage->m_backend = svdDeletionBackend::leadingCovariance;
931 24 : result.m_storage->m_lapackInfo = 0;
932 24 : result.m_storage->m_status = result.m_storage->m_clampedEigenvalues > 0 ? svdDeletionStatus::successWithClamping
933 : : svdDeletionStatus::success;
934 24 : return result.m_storage->m_status;
935 : }
936 :
937 : // Execute the symmetric covariance deletion core after workspace loading.
938 15 : static svdDeletionStatus leading( resultT &result, /* [out] updated system */
939 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
940 : svdDeletionConstMatrixRef<realT> deletedRows, /* [in] deleted rows */
941 : Eigen::Index outputRank, /* [in] requested rank */
942 : workspaceT &workspace /* [in,out] scratch */ )
943 : {
944 15 : if( !validCoreInputs( singularValues, deletedRows, outputRank ) )
945 : {
946 1 : return fail( result, svdDeletionStatus::invalidInput );
947 : }
948 14 : if( deletedRows.rows() == 0 )
949 : {
950 2 : return identity( result, singularValues, outputRank, svdDeletionBackend::leadingCovariance );
951 : }
952 :
953 12 : const Eigen::Index rank = singularValues.size();
954 : svdDeletionStatus status =
955 12 : prepare( result, workspace, rank, deletedRows.rows(), outputRank, svdDeletionBackend::leadingCovariance );
956 12 : if( status != svdDeletionStatus::success )
957 : {
958 2 : return status;
959 : }
960 10 : loadDeletedRows( workspace, deletedRows );
961 10 : return leadingLoaded( result, singularValues, deletedRows.rows(), outputRank, workspace );
962 : }
963 :
964 : // Execute the normalized diagonal-minus-rank-one secular solve with one deleted row already loaded.
965 23 : static svdDeletionStatus rankOneLoaded( resultT &result, /* [out] updated system */
966 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
967 : Eigen::Index outputRank, /* [in] requested rank */
968 : workspaceT &workspace /* [in,out] scratch */ )
969 : {
970 23 : const Eigen::Index rank = singularValues.size();
971 23 : auto &stored = *workspace.m_storage;
972 23 : const realT scale = normalizeSingularValues( workspace, singularValues );
973 23 : if( scale == realT( 0 ) )
974 : {
975 1 : return zeroSystem( result, outputRank, svdDeletionBackend::rankOneSecular );
976 : }
977 :
978 22 : realT updateNorm{ 0 };
979 22 : realT maximumPole{ 0 };
980 88 : for( Eigen::Index index = 0; index < rank; ++index )
981 : {
982 66 : const realT normalized = stored.m_normalizedSingularValues( index );
983 66 : const realT update = normalized * stored.m_deletedRows( 0, index );
984 66 : const realT pole = -( normalized * normalized );
985 66 : if( !math::isFinite( update ) || !math::isFinite( pole ) )
986 : {
987 : return fail( result, svdDeletionStatus::invalidInput ); // LCOV_EXCL_LINE -- finite normalized products.
988 : }
989 66 : stored.m_scaledDeletedTranspose( index, 0 ) = update;
990 66 : stored.m_secularPoles( index ) = pole;
991 66 : updateNorm = std::hypot( updateNorm, update );
992 66 : maximumPole = std::max( maximumPole, std::abs( pole ) );
993 : }
994 22 : if( updateNorm == realT( 0 ) )
995 : {
996 1 : return identity( result, singularValues, outputRank, svdDeletionBackend::rankOneSecular );
997 : }
998 21 : if( !math::isFinite( updateNorm ) || updateNorm > std::sqrt( std::numeric_limits<realT>::max() ) )
999 : {
1000 1 : return fail( result, svdDeletionStatus::invalidInput );
1001 : }
1002 :
1003 20 : const realT rho = updateNorm * updateNorm;
1004 20 : if( !math::isFinite( rho ) )
1005 : {
1006 : // LCOV_EXCL_START -- updateNorm is bounded before squaring.
1007 : return fail( result, svdDeletionStatus::invalidInput );
1008 : // LCOV_EXCL_STOP
1009 : }
1010 20 : realT maximumUpdate{ 0 };
1011 80 : for( Eigen::Index index = 0; index < rank; ++index )
1012 : {
1013 60 : stored.m_secularUpdate( index ) = stored.m_scaledDeletedTranspose( index, 0 ) / updateNorm;
1014 60 : maximumUpdate = std::max( maximumUpdate, std::abs( stored.m_secularUpdate( index ) ) );
1015 : }
1016 20 : const realT deflationTolerance =
1017 20 : realT( 8 ) * std::numeric_limits<realT>::epsilon() * std::max( maximumPole, maximumUpdate );
1018 20 : if( rho * maximumUpdate <= deflationTolerance )
1019 : {
1020 1 : return identity( result, singularValues, outputRank, svdDeletionBackend::rankOneSecular );
1021 : }
1022 :
1023 19 : Eigen::Index activeCount{ 0 };
1024 19 : Eigen::Index deflatedCount{ 0 };
1025 19 : Eigen::Index rotationCount{ 0 };
1026 19 : Eigen::Index previousActive{ -1 };
1027 77 : for( Eigen::Index next = 0; next < rank; ++next )
1028 : {
1029 58 : if( rho * std::abs( stored.m_secularUpdate( next ) ) <= deflationTolerance )
1030 : {
1031 11 : stored.m_secularDeflatedIndices[deflatedCount++] = next;
1032 11 : continue;
1033 : }
1034 47 : if( previousActive < 0 )
1035 : {
1036 19 : previousActive = next;
1037 19 : continue;
1038 : }
1039 :
1040 28 : const realT previousUpdate = stored.m_secularUpdate( previousActive );
1041 28 : const realT nextUpdate = stored.m_secularUpdate( next );
1042 28 : const realT rotationNorm = std::hypot( nextUpdate, previousUpdate );
1043 28 : const realT cosine = nextUpdate / rotationNorm;
1044 28 : const realT sine = -previousUpdate / rotationNorm;
1045 28 : const realT poleDifference = stored.m_secularPoles( next ) - stored.m_secularPoles( previousActive );
1046 28 : if( std::abs( poleDifference * cosine * sine ) <= deflationTolerance )
1047 : {
1048 4 : stored.m_secularRotationFirst[rotationCount] = previousActive;
1049 4 : stored.m_secularRotationSecond[rotationCount] = next;
1050 4 : stored.m_secularRotationCosine[rotationCount] = cosine;
1051 4 : stored.m_secularRotationSine[rotationCount] = sine;
1052 4 : ++rotationCount;
1053 :
1054 4 : stored.m_secularUpdate( next ) = rotationNorm;
1055 4 : stored.m_secularUpdate( previousActive ) = realT( 0 );
1056 4 : const realT previousPole = stored.m_secularPoles( previousActive );
1057 4 : const realT nextPole = stored.m_secularPoles( next );
1058 4 : stored.m_secularPoles( previousActive ) = previousPole * cosine * cosine + nextPole * sine * sine;
1059 4 : stored.m_secularPoles( next ) = previousPole * sine * sine + nextPole * cosine * cosine;
1060 4 : stored.m_secularDeflatedIndices[deflatedCount++] = previousActive;
1061 4 : previousActive = next;
1062 : }
1063 : else
1064 : {
1065 24 : stored.m_secularActiveIndices[activeCount++] = previousActive;
1066 24 : previousActive = next;
1067 : }
1068 : }
1069 19 : if( previousActive >= 0 )
1070 : {
1071 19 : stored.m_secularActiveIndices[activeCount++] = previousActive;
1072 : }
1073 19 : if( activeCount <= 0 || activeCount + deflatedCount != rank )
1074 : {
1075 : // LCOV_EXCL_START -- the deflation loop partitions every index and retains an active component.
1076 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1077 : // LCOV_EXCL_STOP
1078 : }
1079 :
1080 79 : const auto poleIndexLess = [&]( Eigen::Index left, Eigen::Index right )
1081 : {
1082 60 : const realT leftPole = stored.m_secularPoles( left );
1083 60 : const realT rightPole = stored.m_secularPoles( right );
1084 60 : return leftPole < rightPole || ( leftPole == rightPole && left < right );
1085 : };
1086 19 : std::sort( stored.m_secularActiveIndices.begin(),
1087 19 : stored.m_secularActiveIndices.begin() + activeCount,
1088 : poleIndexLess );
1089 19 : std::sort( stored.m_secularDeflatedIndices.begin(),
1090 19 : stored.m_secularDeflatedIndices.begin() + deflatedCount,
1091 : poleIndexLess );
1092 :
1093 19 : realT activeUpdateNorm{ 0 };
1094 62 : for( Eigen::Index active = 0; active < activeCount; ++active )
1095 : {
1096 43 : const Eigen::Index source = stored.m_secularActiveIndices[active];
1097 43 : stored.m_secularActivePoles( active ) = stored.m_secularPoles( source );
1098 43 : activeUpdateNorm = std::hypot( activeUpdateNorm, stored.m_secularUpdate( source ) );
1099 43 : if( active > 0 && !( stored.m_secularActivePoles( active - 1 ) < stored.m_secularActivePoles( active ) ) )
1100 : {
1101 : // LCOV_EXCL_START -- the preceding pass deflates clustered poles.
1102 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1103 : // LCOV_EXCL_STOP
1104 : }
1105 : }
1106 19 : if( activeUpdateNorm == realT( 0 ) || !math::isFinite( activeUpdateNorm ) )
1107 : {
1108 : // LCOV_EXCL_START -- active updates have already passed the finite nonzero threshold.
1109 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1110 : // LCOV_EXCL_STOP
1111 : }
1112 :
1113 19 : const realT activeRho = rho * activeUpdateNorm * activeUpdateNorm;
1114 19 : if( !math::isFinite( activeRho ) )
1115 : {
1116 : // LCOV_EXCL_START -- orthogonal deflation preserves an active update norm no greater than one.
1117 : return fail( result, svdDeletionStatus::invalidInput );
1118 : // LCOV_EXCL_STOP
1119 : }
1120 19 : if( activeRho == realT( 0 ) )
1121 : {
1122 : // LCOV_EXCL_START -- the whole-update threshold handles tiny rho before deflation.
1123 : return identity( result, singularValues, outputRank, svdDeletionBackend::rankOneSecular );
1124 : // LCOV_EXCL_STOP
1125 : }
1126 62 : for( Eigen::Index active = 0; active < activeCount; ++active )
1127 : {
1128 43 : const Eigen::Index source = stored.m_secularActiveIndices[active];
1129 43 : stored.m_secularActiveUpdate( active ) = stored.m_secularUpdate( source ) / activeUpdateNorm;
1130 : }
1131 :
1132 19 : const MXLAPACK_INT lapackActive = static_cast<MXLAPACK_INT>( activeCount );
1133 19 : const MXLAPACK_INT leadingDimension = static_cast<MXLAPACK_INT>( rank );
1134 19 : const MXLAPACK_INT info = callLaed9<realT>( stored.m_eigenvalues.data(),
1135 : stored.m_psdCore.data(),
1136 : stored.m_eigenvectors.data(),
1137 : lapackActive,
1138 : leadingDimension,
1139 : activeRho,
1140 : stored.m_secularActivePoles.data(),
1141 : stored.m_secularActiveUpdate.data() );
1142 19 : if( info != 0 )
1143 : {
1144 1 : return fail( result, svdDeletionStatus::solverFailure, info );
1145 : }
1146 52 : for( Eigen::Index active = 0; active < activeCount; ++active )
1147 : {
1148 36 : if( !math::isFinite( stored.m_eigenvalues( active ) ) )
1149 : {
1150 1 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1151 : }
1152 121 : for( Eigen::Index row = 0; row < activeCount; ++row )
1153 : {
1154 87 : if( !math::isFinite( stored.m_eigenvectors( row, active ) ) )
1155 : {
1156 1 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1157 : }
1158 : }
1159 : }
1160 16 : if( !validAscendingSpectrum( stored.m_eigenvalues.head( activeCount ) ) )
1161 : {
1162 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1163 : }
1164 :
1165 15 : const realT secularTolerance = realT( 64 ) * std::numeric_limits<realT>::epsilon() *
1166 15 : static_cast<realT>( std::max<Eigen::Index>( 1, rank ) ) *
1167 15 : ( realT( 1 ) + activeRho );
1168 43 : for( Eigen::Index active = 0; active < activeCount; ++active )
1169 : {
1170 29 : const realT lower = stored.m_secularActivePoles( active );
1171 29 : const realT upper = active + 1 < activeCount ? stored.m_secularActivePoles( active + 1 )
1172 14 : : stored.m_secularActivePoles( active ) + activeRho;
1173 57 : if( stored.m_eigenvalues( active ) < lower - secularTolerance ||
1174 28 : stored.m_eigenvalues( active ) > upper + secularTolerance )
1175 : {
1176 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1177 : }
1178 28 : stored.m_secularCandidateValues( active ) = stored.m_eigenvalues( active );
1179 : }
1180 29 : for( Eigen::Index deflated = 0; deflated < deflatedCount; ++deflated )
1181 : {
1182 15 : stored.m_secularCandidateValues( activeCount + deflated ) =
1183 15 : stored.m_secularPoles( stored.m_secularDeflatedIndices[deflated] );
1184 : }
1185 :
1186 14 : Eigen::Index activePosition{ 0 };
1187 14 : Eigen::Index deflatedPosition{ 0 };
1188 14 : Eigen::Index mergedPosition{ 0 };
1189 57 : while( activePosition < activeCount || deflatedPosition < deflatedCount )
1190 : {
1191 62 : const bool takeActive = deflatedPosition >= deflatedCount ||
1192 19 : ( activePosition < activeCount &&
1193 19 : stored.m_secularCandidateValues( activePosition ) <=
1194 19 : stored.m_secularCandidateValues( activeCount + deflatedPosition ) );
1195 43 : stored.m_secularOrder[mergedPosition++] = takeActive ? activePosition++ : activeCount + deflatedPosition++;
1196 : }
1197 14 : if( mergedPosition != rank )
1198 : {
1199 : // LCOV_EXCL_START -- the merge consumes each member of the checked partition once.
1200 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1201 : // LCOV_EXCL_STOP
1202 : }
1203 :
1204 14 : result.m_storage->m_minimumPSDValue = -stored.m_secularCandidateValues( stored.m_secularOrder[rank - 1] );
1205 14 : const realT tolerance = psdTolerance<realT>( rank, realT( 1 ) + rho );
1206 14 : if( result.m_storage->m_minimumPSDValue < -tolerance )
1207 : {
1208 1 : return fail( result, svdDeletionStatus::nonPositiveSemidefinite );
1209 : }
1210 :
1211 13 : result.m_storage->m_clampedEigenvalues = 0;
1212 13 : const realT vectorTolerance = realT( 256 ) * std::numeric_limits<realT>::epsilon() *
1213 13 : static_cast<realT>( std::max<Eigen::Index>( 1, rank ) );
1214 13 : const realT residualTolerance = vectorTolerance * ( realT( 1 ) + rho );
1215 46 : for( Eigen::Index output = 0; output < rank; ++output )
1216 : {
1217 36 : const Eigen::Index candidate = stored.m_secularOrder[output];
1218 36 : realT normalizedSquaredValue = -stored.m_secularCandidateValues( candidate );
1219 36 : if( normalizedSquaredValue < realT( 0 ) )
1220 : {
1221 1 : normalizedSquaredValue = realT( 0 );
1222 1 : ++result.m_storage->m_clampedEigenvalues;
1223 : }
1224 36 : if( !rescaleValue( result.m_storage->m_singularValues( output ),
1225 36 : result.m_storage->m_squaredSingularValues( output ),
1226 : normalizedSquaredValue,
1227 : scale ) )
1228 : {
1229 1 : return fail( result, svdDeletionStatus::rescalingOverflow );
1230 : }
1231 35 : if( output >= outputRank )
1232 : {
1233 3 : continue;
1234 : }
1235 :
1236 32 : stored.m_secularVector.setZero();
1237 32 : if( candidate < activeCount )
1238 : {
1239 71 : for( Eigen::Index active = 0; active < activeCount; ++active )
1240 : {
1241 51 : stored.m_secularVector( stored.m_secularActiveIndices[active] ) =
1242 51 : stored.m_eigenvectors( active, candidate );
1243 : }
1244 : }
1245 : else
1246 : {
1247 12 : const Eigen::Index deflated = candidate - activeCount;
1248 12 : stored.m_secularVector( stored.m_secularDeflatedIndices[deflated] ) = realT( 1 );
1249 : }
1250 48 : for( Eigen::Index rotation = rotationCount; rotation-- > 0; )
1251 : {
1252 16 : const Eigen::Index first = stored.m_secularRotationFirst[rotation];
1253 16 : const Eigen::Index second = stored.m_secularRotationSecond[rotation];
1254 16 : const realT cosine = stored.m_secularRotationCosine[rotation];
1255 16 : const realT sine = stored.m_secularRotationSine[rotation];
1256 16 : const realT firstValue = stored.m_secularVector( first );
1257 16 : const realT secondValue = stored.m_secularVector( second );
1258 16 : stored.m_secularVector( first ) = cosine * firstValue - sine * secondValue;
1259 16 : stored.m_secularVector( second ) = sine * firstValue + cosine * secondValue;
1260 : }
1261 :
1262 32 : realT vectorNorm{ 0 };
1263 32 : Eigen::Index signIndex{ 0 };
1264 32 : realT signMagnitude{ 0 };
1265 147 : for( Eigen::Index row = 0; row < rank; ++row )
1266 : {
1267 115 : const realT value = stored.m_secularVector( row );
1268 115 : if( !math::isFinite( value ) )
1269 : {
1270 : // LCOV_EXCL_START -- finite bounded Givens rotations preserve finite vector entries.
1271 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1272 : // LCOV_EXCL_STOP
1273 : }
1274 115 : vectorNorm = std::hypot( vectorNorm, value );
1275 115 : if( std::abs( value ) > signMagnitude )
1276 : {
1277 49 : signMagnitude = std::abs( value );
1278 49 : signIndex = row;
1279 : }
1280 : }
1281 32 : if( !math::isFinite( vectorNorm ) || std::abs( vectorNorm - realT( 1 ) ) > vectorTolerance )
1282 : {
1283 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1284 : }
1285 31 : if( stored.m_secularVector( signIndex ) < realT( 0 ) )
1286 : {
1287 13 : stored.m_secularVector = -stored.m_secularVector;
1288 : }
1289 :
1290 31 : realT updateDot{ 0 };
1291 143 : for( Eigen::Index row = 0; row < rank; ++row )
1292 : {
1293 112 : updateDot += stored.m_scaledDeletedTranspose( row, 0 ) * stored.m_secularVector( row );
1294 : }
1295 31 : realT residualNorm{ 0 };
1296 143 : for( Eigen::Index row = 0; row < rank; ++row )
1297 : {
1298 112 : const realT normalized = stored.m_normalizedSingularValues( row );
1299 112 : const realT residual = normalized * normalized * stored.m_secularVector( row ) -
1300 112 : stored.m_scaledDeletedTranspose( row, 0 ) * updateDot -
1301 112 : normalizedSquaredValue * stored.m_secularVector( row );
1302 112 : residualNorm = std::hypot( residualNorm, residual );
1303 : }
1304 31 : if( !math::isFinite( residualNorm ) || residualNorm > residualTolerance )
1305 : {
1306 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1307 : }
1308 30 : result.m_storage->m_rotation.matrix().col( output ) = stored.m_secularVector.matrix();
1309 : }
1310 :
1311 10 : result.m_storage->m_backend = svdDeletionBackend::rankOneSecular;
1312 10 : result.m_storage->m_lapackInfo = 0;
1313 10 : result.m_storage->m_status = result.m_storage->m_clampedEigenvalues > 0 ? svdDeletionStatus::successWithClamping
1314 : : svdDeletionStatus::success;
1315 10 : return result.m_storage->m_status;
1316 : }
1317 :
1318 : // Execute the structured diagonal-minus-rank-one deletion core.
1319 17 : static svdDeletionStatus rankOne( resultT &result, /* [out] updated system */
1320 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
1321 : svdDeletionConstMatrixRef<realT> deletedRows, /* [in] deleted rows */
1322 : Eigen::Index outputRank, /* [in] requested rank */
1323 : workspaceT &workspace /* [in,out] scratch */ )
1324 : {
1325 17 : if( !validCoreInputs( singularValues, deletedRows, outputRank ) )
1326 : {
1327 1 : return fail( result, svdDeletionStatus::invalidInput );
1328 : }
1329 16 : if( deletedRows.rows() == 0 )
1330 : {
1331 1 : return identity( result, singularValues, outputRank, svdDeletionBackend::rankOneSecular );
1332 : }
1333 15 : if( deletedRows.rows() != 1 )
1334 : {
1335 1 : return fail( result, svdDeletionStatus::unsupportedDeletionCount );
1336 : }
1337 :
1338 14 : const Eigen::Index rank = singularValues.size();
1339 : const svdDeletionStatus status =
1340 14 : prepare( result, workspace, rank, 1, outputRank, svdDeletionBackend::rankOneSecular );
1341 14 : if( status != svdDeletionStatus::success )
1342 : {
1343 1 : return status;
1344 : }
1345 13 : loadDeletedRows( workspace, deletedRows );
1346 13 : return rankOneLoaded( result, singularValues, outputRank, workspace );
1347 : }
1348 :
1349 : // Execute a normalized complement-preserving SVD with deleted rows already loaded.
1350 39 : static svdDeletionStatus stableLoaded( resultT &result, /* [out] updated system */
1351 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
1352 : Eigen::Index deletedCount, /* [in] active rows */
1353 : Eigen::Index outputRank, /* [in] requested rank */
1354 : workspaceT &workspace /* [in,out] scratch */ )
1355 : {
1356 39 : const Eigen::Index rank = singularValues.size();
1357 39 : const realT scale = normalizeSingularValues( workspace, singularValues );
1358 39 : if( scale == realT( 0 ) )
1359 : {
1360 1 : return zeroSystem( result, outputRank, svdDeletionBackend::stableCore );
1361 : }
1362 :
1363 38 : const auto factorRows = workspace.m_storage->m_deletedRows.matrix().topRows( deletedCount );
1364 38 : auto complement = workspace.m_storage->m_psdCore.matrix().topLeftCorner( deletedCount, deletedCount );
1365 38 : complement.setIdentity();
1366 38 : complement.noalias() -= factorRows * factorRows.transpose();
1367 38 : MXLAPACK_INT found{ 0 };
1368 38 : const MXLAPACK_INT lapackComplement = static_cast<MXLAPACK_INT>( deletedCount );
1369 228 : MXLAPACK_INT info = callSyevr<realT>( 'V',
1370 : 'A',
1371 : 'L',
1372 : lapackComplement,
1373 38 : workspace.m_storage->m_psdCore.data(),
1374 38 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_psdCore.rows() ),
1375 : realT( 0 ),
1376 : realT( 0 ),
1377 : 1,
1378 : lapackComplement,
1379 : math::lamch<realT>( 'S' ),
1380 : &found,
1381 38 : workspace.m_storage->m_eigenvalues.data(),
1382 38 : workspace.m_storage->m_eigenvectors.data(),
1383 38 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_eigenvectors.rows() ),
1384 38 : workspace.m_storage->m_support.data(),
1385 38 : workspace.m_storage->m_syevrWork.data(),
1386 38 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_syevrWork.size() ),
1387 38 : workspace.m_storage->m_syevrIWork.data(),
1388 38 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_syevrIWork.size() ) );
1389 38 : if( info != 0 || found != lapackComplement )
1390 : {
1391 1 : return fail( result, svdDeletionStatus::solverFailure, info );
1392 : }
1393 89 : for( Eigen::Index index = 0; index < deletedCount; ++index )
1394 : {
1395 54 : if( !math::isFinite( workspace.m_storage->m_eigenvalues( index ) ) )
1396 : {
1397 1 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1398 : }
1399 147 : for( Eigen::Index row = 0; row < deletedCount; ++row )
1400 : {
1401 95 : if( !math::isFinite( workspace.m_storage->m_eigenvectors( row, index ) ) )
1402 : {
1403 1 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1404 : }
1405 : }
1406 : }
1407 35 : if( !validAscendingSpectrum( workspace.m_storage->m_eigenvalues.head( deletedCount ) ) )
1408 : {
1409 1 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1410 : }
1411 :
1412 34 : result.m_storage->m_minimumPSDValue = workspace.m_storage->m_eigenvalues( 0 );
1413 : const realT tolerance =
1414 34 : psdTolerance<realT>( deletedCount, realT( 1 ) + static_cast<realT>( factorRows.squaredNorm() ) );
1415 34 : if( result.m_storage->m_minimumPSDValue < -tolerance )
1416 : {
1417 1 : return fail( result, svdDeletionStatus::nonPositiveSemidefinite );
1418 : }
1419 :
1420 33 : result.m_storage->m_clampedEigenvalues = 0;
1421 82 : for( Eigen::Index index = 0; index < deletedCount; ++index )
1422 : {
1423 49 : if( workspace.m_storage->m_eigenvalues( index ) < realT( 0 ) )
1424 : {
1425 3 : workspace.m_storage->m_eigenvalues( index ) = realT( 0 );
1426 3 : ++result.m_storage->m_clampedEigenvalues;
1427 : }
1428 : }
1429 :
1430 : auto complementRoot =
1431 33 : workspace.m_storage->m_complementRoot.matrix().topLeftCorner( deletedCount, deletedCount );
1432 82 : for( Eigen::Index row = 0; row < deletedCount; ++row )
1433 : {
1434 49 : const realT root = std::sqrt( workspace.m_storage->m_eigenvalues( row ) );
1435 138 : for( Eigen::Index column = 0; column < deletedCount; ++column )
1436 : {
1437 89 : complementRoot( row, column ) = root * workspace.m_storage->m_eigenvectors( column, row );
1438 : }
1439 : }
1440 :
1441 33 : auto topCore = workspace.m_storage->m_stableCore.matrix().topRows( rank );
1442 33 : topCore.setIdentity();
1443 33 : topCore.noalias() -= factorRows.transpose() * factorRows;
1444 33 : workspace.m_storage->m_stableCore.matrix().middleRows( rank, deletedCount ).noalias() =
1445 33 : complementRoot * factorRows;
1446 133 : for( Eigen::Index column = 0; column < rank; ++column )
1447 : {
1448 100 : workspace.m_storage->m_stableCore.matrix().col( column ).head( rank + deletedCount ) *=
1449 100 : workspace.m_storage->m_normalizedSingularValues( column );
1450 : }
1451 :
1452 33 : realT unusedLeft{ 0 };
1453 33 : const MXLAPACK_INT coreRows = static_cast<MXLAPACK_INT>( rank + deletedCount );
1454 33 : const MXLAPACK_INT lapackRank = static_cast<MXLAPACK_INT>( rank );
1455 132 : info = callGesvd<realT>( 'N',
1456 : 'S',
1457 : coreRows,
1458 : lapackRank,
1459 33 : workspace.m_storage->m_stableCore.data(),
1460 33 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_stableCore.rows() ),
1461 33 : workspace.m_storage->m_stableSingularValues.data(),
1462 : &unusedLeft,
1463 : 1,
1464 33 : workspace.m_storage->m_rightTranspose.data(),
1465 : lapackRank,
1466 33 : workspace.m_storage->m_gesvdWork.data(),
1467 33 : static_cast<MXLAPACK_INT>( workspace.m_storage->m_gesvdWork.size() ) );
1468 33 : if( info != 0 )
1469 : {
1470 1 : return fail( result, svdDeletionStatus::solverFailure, info );
1471 : }
1472 63 : if( !allFinite( workspace.m_storage->m_stableSingularValues ) ||
1473 31 : !allFinite( workspace.m_storage->m_rightTranspose ) )
1474 : {
1475 2 : return fail( result, svdDeletionStatus::nonFiniteOutput );
1476 : }
1477 30 : if( !validDescendingSpectrum( workspace.m_storage->m_stableSingularValues ) )
1478 : {
1479 2 : return fail( result, svdDeletionStatus::invalidSolverOutput );
1480 : }
1481 :
1482 112 : for( Eigen::Index output = 0; output < rank; ++output )
1483 : {
1484 85 : const realT normalizedValue = workspace.m_storage->m_stableSingularValues( output );
1485 85 : if( !rescaleSingularValue( result.m_storage->m_singularValues( output ),
1486 85 : result.m_storage->m_squaredSingularValues( output ),
1487 : normalizedValue,
1488 : scale ) )
1489 : {
1490 1 : return fail( result, svdDeletionStatus::rescalingOverflow );
1491 : }
1492 84 : if( output < outputRank )
1493 : {
1494 82 : result.m_storage->m_rotation.matrix().col( output ) =
1495 164 : workspace.m_storage->m_rightTranspose.matrix().row( output ).transpose();
1496 : }
1497 : }
1498 :
1499 27 : result.m_storage->m_backend = svdDeletionBackend::stableCore;
1500 27 : result.m_storage->m_lapackInfo = 0;
1501 27 : result.m_storage->m_status = result.m_storage->m_clampedEigenvalues > 0 ? svdDeletionStatus::successWithClamping
1502 : : svdDeletionStatus::success;
1503 27 : return result.m_storage->m_status;
1504 : }
1505 :
1506 : // Execute the complement-preserving small-SVD deletion core after workspace loading.
1507 25 : static svdDeletionStatus stable( resultT &result, /* [out] updated system */
1508 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
1509 : svdDeletionConstMatrixRef<realT> deletedRows, /* [in] deleted rows */
1510 : Eigen::Index outputRank, /* [in] requested rank */
1511 : workspaceT &workspace /* [in,out] scratch */ )
1512 : {
1513 25 : if( !validCoreInputs( singularValues, deletedRows, outputRank ) )
1514 : {
1515 4 : return fail( result, svdDeletionStatus::invalidInput );
1516 : }
1517 21 : if( deletedRows.rows() == 0 )
1518 : {
1519 1 : return identity( result, singularValues, outputRank, svdDeletionBackend::stableCore );
1520 : }
1521 :
1522 20 : const Eigen::Index rank = singularValues.size();
1523 : svdDeletionStatus status =
1524 20 : prepare( result, workspace, rank, deletedRows.rows(), outputRank, svdDeletionBackend::stableCore );
1525 20 : if( status != svdDeletionStatus::success )
1526 : {
1527 1 : return status;
1528 : }
1529 19 : loadDeletedRows( workspace, deletedRows );
1530 19 : return stableLoaded( result, singularValues, deletedRows.rows(), outputRank, workspace );
1531 : }
1532 :
1533 : // Gather arbitrary factor rows and dispatch to the selected core.
1534 70 : static svdDeletionStatus remove( resultT &result, /* [out] updated system */
1535 : svdDeletionConstVectorRef<realT> singularValues, /* [in] base values */
1536 : svdDeletionConstMatrixRef<realT> factor, /* [in] deleted-side factor */
1537 : svdDeletionConstIndexViewV2 indices, /* [in] rows to remove */
1538 : Eigen::Index outputRank, /* [in] requested rank */
1539 : workspaceT &workspace, /* [in,out] scratch */
1540 : svdDeletionBackend backend /* [in] selected backend */ )
1541 : {
1542 70 : const Eigen::Index rank = singularValues.size();
1543 202 : if( !validSingularInputs( singularValues, outputRank ) || !validBackend( backend ) || factor.rows() <= 0 ||
1544 131 : factor.rows() < rank || factor.cols() != rank || !validLapackDimension( factor.rows() ) ||
1545 198 : !validAbiIndexView( indices ) || indices.size >= static_cast<std::int64_t>( factor.rows() ) ||
1546 62 : indices.size > static_cast<std::int64_t>( std::numeric_limits<MXLAPACK_INT>::max() - rank ) )
1547 : {
1548 8 : return fail( result, svdDeletionStatus::invalidInput );
1549 : }
1550 :
1551 62 : Eigen::Index previous = -1;
1552 146 : for( std::int64_t offset = 0; offset < indices.size; ++offset )
1553 : {
1554 86 : const std::int64_t rawIndex = abiIndexAt( indices, offset );
1555 86 : if( !validAbiDimension( rawIndex ) )
1556 : {
1557 1 : return fail( result, svdDeletionStatus::invalidInput );
1558 : }
1559 85 : const Eigen::Index index = static_cast<Eigen::Index>( rawIndex );
1560 85 : if( index < 0 || index >= factor.rows() || index <= previous )
1561 : {
1562 1 : return fail( result, svdDeletionStatus::invalidInput );
1563 : }
1564 84 : previous = index;
1565 : }
1566 :
1567 60 : if( indices.size == 0 )
1568 : {
1569 5 : return identity( result, singularValues, outputRank, backend );
1570 : }
1571 :
1572 55 : if( backend == svdDeletionBackend::rankOneSecular && indices.size != 1 )
1573 : {
1574 1 : return fail( result, svdDeletionStatus::unsupportedDeletionCount );
1575 : }
1576 :
1577 54 : const Eigen::Index deletedCount = static_cast<Eigen::Index>( indices.size );
1578 54 : const svdDeletionStatus status = prepare( result, workspace, rank, deletedCount, outputRank, backend );
1579 54 : if( status != svdDeletionStatus::success )
1580 : {
1581 1 : return status;
1582 : }
1583 :
1584 133 : for( Eigen::Index row = 0; row < deletedCount; ++row )
1585 : {
1586 81 : const Eigen::Index source = static_cast<Eigen::Index>( abiIndexAt( indices, row ) );
1587 339 : for( Eigen::Index column = 0; column < rank; ++column )
1588 : {
1589 259 : const realT value = factor( source, column );
1590 259 : if( !math::isFinite( value ) )
1591 : {
1592 1 : return fail( result, svdDeletionStatus::invalidInput );
1593 : }
1594 258 : workspace.m_storage->m_deletedRows( row, column ) = value;
1595 : }
1596 : }
1597 :
1598 52 : switch( backend )
1599 : {
1600 22 : case svdDeletionBackend::leadingCovariance:
1601 22 : return leadingLoaded( result, singularValues, deletedCount, outputRank, workspace );
1602 20 : case svdDeletionBackend::stableCore:
1603 20 : return stableLoaded( result, singularValues, deletedCount, outputRank, workspace );
1604 10 : case svdDeletionBackend::rankOneSecular:
1605 10 : return rankOneLoaded( result, singularValues, outputRank, workspace );
1606 : }
1607 : // LCOV_EXCL_START -- the backend is validated before the exhaustive switch.
1608 : return fail( result, svdDeletionStatus::invalidInput );
1609 : // LCOV_EXCL_STOP
1610 : }
1611 : };
1612 :
1613 : } // namespace detail
1614 :
1615 36 : const char *svdDeletionBackendName( svdDeletionBackend backend )
1616 : {
1617 36 : switch( backend )
1618 : {
1619 12 : case svdDeletionBackend::leadingCovariance:
1620 12 : return "leadingCovariance";
1621 12 : case svdDeletionBackend::stableCore:
1622 12 : return "stableCore";
1623 11 : case svdDeletionBackend::rankOneSecular:
1624 11 : return "rankOneSecular";
1625 : }
1626 1 : return "unknown";
1627 : }
1628 :
1629 24 : const char *svdDeletionStatusName( svdDeletionStatus status )
1630 : {
1631 24 : switch( status )
1632 : {
1633 1 : case svdDeletionStatus::notComputed:
1634 1 : return "notComputed";
1635 1 : case svdDeletionStatus::success:
1636 1 : return "success";
1637 1 : case svdDeletionStatus::successWithClamping:
1638 1 : return "successWithClamping";
1639 1 : case svdDeletionStatus::invalidInput:
1640 1 : return "invalidInput";
1641 1 : case svdDeletionStatus::allocationFailure:
1642 1 : return "allocationFailure";
1643 1 : case svdDeletionStatus::workspaceQueryFailure:
1644 1 : return "workspaceQueryFailure";
1645 1 : case svdDeletionStatus::solverFailure:
1646 1 : return "solverFailure";
1647 1 : case svdDeletionStatus::nonFiniteOutput:
1648 1 : return "nonFiniteOutput";
1649 1 : case svdDeletionStatus::invalidSolverOutput:
1650 1 : return "invalidSolverOutput";
1651 1 : case svdDeletionStatus::rescalingOverflow:
1652 1 : return "rescalingOverflow";
1653 1 : case svdDeletionStatus::nonPositiveSemidefinite:
1654 1 : return "nonPositiveSemidefinite";
1655 1 : case svdDeletionStatus::factorNotOrthonormal:
1656 1 : return "factorNotOrthonormal";
1657 11 : case svdDeletionStatus::unsupportedDeletionCount:
1658 11 : return "unsupportedDeletionCount";
1659 : }
1660 1 : return "unknown";
1661 : }
1662 :
1663 57 : bool svdDeletionSucceeded( svdDeletionStatus status ) noexcept
1664 : {
1665 57 : return status == svdDeletionStatus::success || status == svdDeletionStatus::successWithClamping;
1666 : }
1667 :
1668 : template <typename realT, typename abiT>
1669 92 : svdDeletionResult<realT, abiT>::svdDeletionResult() : m_storage( std::make_unique<storage>() )
1670 : {
1671 92 : }
1672 :
1673 : template <typename realT, typename abiT>
1674 94 : svdDeletionResult<realT, abiT>::~svdDeletionResult() = default;
1675 :
1676 : template <typename realT, typename abiT>
1677 2 : svdDeletionResult<realT, abiT>::svdDeletionResult( svdDeletionResult &&other ) noexcept = default;
1678 :
1679 : template <typename realT, typename abiT>
1680 : svdDeletionResult<realT, abiT> &
1681 1 : svdDeletionResult<realT, abiT>::operator=( svdDeletionResult &&other ) noexcept = default;
1682 :
1683 : template <typename realT, typename abiT>
1684 179 : bool svdDeletionResult<realT, abiT>::ensureStorage() noexcept
1685 : {
1686 179 : if( m_storage )
1687 : {
1688 175 : return true;
1689 : }
1690 :
1691 : try
1692 : {
1693 4 : detail::callOperationHook<realT>( detail::svdDeletionTestOperation::prepareResult );
1694 2 : m_storage = std::make_unique<storage>();
1695 : }
1696 2 : catch( const std::bad_alloc & )
1697 : {
1698 1 : return false;
1699 : }
1700 1 : catch( const std::length_error & )
1701 : {
1702 1 : return false;
1703 : }
1704 2 : return true;
1705 : }
1706 :
1707 : template <typename realT, typename abiT>
1708 120 : svdDeletionStatus svdDeletionResult<realT, abiT>::prepareAbiV2( std::int64_t abiBaseRank, std::int64_t abiOutputRank )
1709 : {
1710 120 : if( !ensureStorage() )
1711 : {
1712 1 : return svdDeletionStatus::allocationFailure;
1713 : }
1714 119 : auto &stored = *m_storage;
1715 :
1716 119 : if( !detail::validAbiDimension( abiBaseRank ) || !detail::validAbiDimension( abiOutputRank ) )
1717 : {
1718 1 : stored.m_status = svdDeletionStatus::invalidInput;
1719 1 : return stored.m_status;
1720 : }
1721 118 : const Eigen::Index baseRank = static_cast<Eigen::Index>( abiBaseRank );
1722 118 : const Eigen::Index outputRank = static_cast<Eigen::Index>( abiOutputRank );
1723 :
1724 117 : if( baseRank <= 0 || outputRank <= 0 || outputRank > baseRank || !detail::validLapackDimension( baseRank ) ||
1725 235 : !detail::validArrayShape<realT>( baseRank, outputRank ) || !detail::validArrayShape<realT>( outputRank, 1 ) )
1726 : {
1727 2 : stored.m_status = svdDeletionStatus::invalidInput;
1728 2 : return stored.m_status;
1729 : }
1730 :
1731 : try
1732 : {
1733 116 : if( stored.m_baseRank != baseRank || stored.m_maximumOutputRank < outputRank )
1734 : {
1735 92 : detail::callOperationHook<realT>( detail::svdDeletionTestOperation::prepareResult );
1736 86 : stored.m_singularValues.resize( baseRank );
1737 86 : stored.m_squaredSingularValues.resize( baseRank );
1738 86 : stored.m_rotation.resize( baseRank, outputRank );
1739 86 : stored.m_maximumOutputRank = outputRank;
1740 : }
1741 : }
1742 6 : catch( const std::bad_alloc & )
1743 : {
1744 5 : stored.m_status = svdDeletionStatus::allocationFailure;
1745 5 : return stored.m_status;
1746 : }
1747 1 : catch( const std::length_error & )
1748 : {
1749 1 : stored.m_status = svdDeletionStatus::allocationFailure;
1750 1 : return stored.m_status;
1751 : }
1752 :
1753 110 : stored.m_baseRank = baseRank;
1754 110 : stored.m_outputRank = outputRank;
1755 110 : stored.m_status = svdDeletionStatus::notComputed;
1756 110 : stored.m_lapackInfo = 0;
1757 110 : stored.m_clampedEigenvalues = 0;
1758 110 : stored.m_minimumPSDValue = realT( 0 );
1759 110 : return svdDeletionStatus::success;
1760 : }
1761 :
1762 : template <typename realT, typename abiT>
1763 10 : svdDeletionStatus svdDeletionResult<realT, abiT>::status() const noexcept
1764 : {
1765 10 : return m_storage ? m_storage->m_status : svdDeletionStatus::notComputed;
1766 : }
1767 :
1768 : template <typename realT, typename abiT>
1769 14 : svdDeletionBackend svdDeletionResult<realT, abiT>::backend() const noexcept
1770 : {
1771 14 : return m_storage ? m_storage->m_backend : svdDeletionBackend::stableCore;
1772 : }
1773 :
1774 : template <typename realT, typename abiT>
1775 81 : svdDeletionConstVectorViewV2<realT> svdDeletionResult<realT, abiT>::singularValuesViewAbiV2() const noexcept
1776 : {
1777 81 : return m_storage
1778 81 : ? svdDeletionConstVectorViewV2<realT>{ m_storage->m_singularValues.data(),
1779 80 : static_cast<std::int64_t>( m_storage->m_singularValues.size() ) }
1780 161 : : svdDeletionConstVectorViewV2<realT>{};
1781 : }
1782 :
1783 : template <typename realT, typename abiT>
1784 102 : svdDeletionConstVectorViewV2<realT> svdDeletionResult<realT, abiT>::squaredSingularValuesViewAbiV2() const noexcept
1785 : {
1786 102 : return m_storage ? svdDeletionConstVectorViewV2<realT>{ m_storage->m_squaredSingularValues.data(),
1787 : static_cast<std::int64_t>(
1788 101 : m_storage->m_squaredSingularValues.size() ) }
1789 203 : : svdDeletionConstVectorViewV2<realT>{};
1790 : }
1791 :
1792 : template <typename realT, typename abiT>
1793 50 : svdDeletionConstMatrixViewV2<realT> svdDeletionResult<realT, abiT>::rotationViewAbiV2() const noexcept
1794 : {
1795 50 : if( !m_storage )
1796 : {
1797 1 : return {};
1798 : }
1799 49 : return { m_storage->m_rotation.data(),
1800 49 : static_cast<std::int64_t>( m_storage->m_rotation.rows() ),
1801 98 : static_cast<std::int64_t>( m_storage->m_outputRank ),
1802 98 : static_cast<std::int64_t>( m_storage->m_rotation.outerStride() ) };
1803 : }
1804 :
1805 : template <typename realT, typename abiT>
1806 36 : std::int64_t svdDeletionResult<realT, abiT>::baseRank() const noexcept
1807 : {
1808 36 : return m_storage ? m_storage->m_baseRank : 0;
1809 : }
1810 :
1811 : template <typename realT, typename abiT>
1812 38 : std::int64_t svdDeletionResult<realT, abiT>::outputRank() const noexcept
1813 : {
1814 38 : return m_storage ? m_storage->m_outputRank : 0;
1815 : }
1816 :
1817 : template <typename realT, typename abiT>
1818 2 : std::int64_t svdDeletionResult<realT, abiT>::maximumOutputRank() const noexcept
1819 : {
1820 2 : return m_storage ? m_storage->m_maximumOutputRank : 0;
1821 : }
1822 :
1823 : template <typename realT, typename abiT>
1824 4 : std::int64_t svdDeletionResult<realT, abiT>::clampedEigenvalues() const noexcept
1825 : {
1826 4 : return m_storage ? m_storage->m_clampedEigenvalues : 0;
1827 : }
1828 :
1829 : template <typename realT, typename abiT>
1830 11 : realT svdDeletionResult<realT, abiT>::minimumPSDValue() const noexcept
1831 : {
1832 11 : return m_storage ? m_storage->m_minimumPSDValue : realT( 0 );
1833 : }
1834 :
1835 : template <typename realT, typename abiT>
1836 7 : MXLAPACK_INT svdDeletionResult<realT, abiT>::lapackInfo() const noexcept
1837 : {
1838 7 : return m_storage ? m_storage->m_lapackInfo : 0;
1839 : }
1840 :
1841 : template <typename realT, typename abiT>
1842 93 : svdDeletionWorkspace<realT, abiT>::svdDeletionWorkspace() : m_storage( std::make_unique<storage>() )
1843 : {
1844 93 : }
1845 :
1846 : template <typename realT, typename abiT>
1847 95 : svdDeletionWorkspace<realT, abiT>::~svdDeletionWorkspace()
1848 : {
1849 95 : }
1850 :
1851 : template <typename realT, typename abiT>
1852 2 : svdDeletionWorkspace<realT, abiT>::svdDeletionWorkspace( svdDeletionWorkspace &&other ) noexcept = default;
1853 :
1854 : template <typename realT, typename abiT>
1855 : svdDeletionWorkspace<realT, abiT> &
1856 1 : svdDeletionWorkspace<realT, abiT>::operator=( svdDeletionWorkspace &&other ) noexcept = default;
1857 :
1858 : template <typename realT, typename abiT>
1859 103 : bool svdDeletionWorkspace<realT, abiT>::ensureStorage() noexcept
1860 : {
1861 103 : if( m_storage )
1862 : {
1863 99 : return true;
1864 : }
1865 :
1866 : try
1867 : {
1868 4 : detail::callOperationHook<realT>( detail::svdDeletionTestOperation::prepareWorkspace );
1869 2 : m_storage = std::make_unique<storage>();
1870 : }
1871 2 : catch( const std::bad_alloc & )
1872 : {
1873 1 : return false;
1874 : }
1875 1 : catch( const std::length_error & )
1876 : {
1877 1 : return false;
1878 : }
1879 2 : return true;
1880 : }
1881 :
1882 : template <typename realT, typename abiT>
1883 103 : svdDeletionStatus svdDeletionWorkspace<realT, abiT>::prepareAbiV2( std::int64_t abiBaseRank,
1884 : std::int64_t abiMaximumDeleted,
1885 : svdDeletionBackend backend )
1886 : {
1887 103 : if( !ensureStorage() )
1888 : {
1889 2 : return svdDeletionStatus::allocationFailure;
1890 : }
1891 101 : auto &stored = *m_storage;
1892 101 : stored.m_lapackInfo = 0;
1893 :
1894 101 : if( !detail::validAbiDimension( abiBaseRank ) || !detail::validAbiDimension( abiMaximumDeleted ) )
1895 : {
1896 1 : return svdDeletionStatus::invalidInput;
1897 : }
1898 100 : const Eigen::Index baseRank = static_cast<Eigen::Index>( abiBaseRank );
1899 100 : const Eigen::Index maximumDeleted = static_cast<Eigen::Index>( abiMaximumDeleted );
1900 100 : if( baseRank <= 0 || maximumDeleted < 0 || !detail::svdDeletionImplementation<realT>::validBackend( backend ) ||
1901 299 : !detail::validLapackDimension( baseRank ) || !detail::validLapackDimension( maximumDeleted ) ||
1902 99 : maximumDeleted > std::numeric_limits<MXLAPACK_INT>::max() - baseRank )
1903 : {
1904 2 : return svdDeletionStatus::invalidInput;
1905 : }
1906 98 : if( backend == svdDeletionBackend::rankOneSecular && maximumDeleted > 1 )
1907 : {
1908 1 : return svdDeletionStatus::unsupportedDeletionCount;
1909 : }
1910 97 : if( stored.m_prepared && stored.m_baseRank == baseRank && stored.m_maximumDeleted >= maximumDeleted &&
1911 6 : stored.m_backend == backend )
1912 : {
1913 1 : return svdDeletionStatus::success;
1914 : }
1915 :
1916 96 : const Eigen::Index eigenSize = backend == svdDeletionBackend::stableCore ? maximumDeleted : baseRank;
1917 96 : const Eigen::Index syevrSize = backend == svdDeletionBackend::rankOneSecular ? 0 : eigenSize;
1918 96 : if( syevrSize > std::numeric_limits<Eigen::Index>::max() / 2 )
1919 : {
1920 : return svdDeletionStatus::invalidInput; // LCOV_EXCL_LINE -- requires an Eigen index narrower than LAPACK.
1921 : }
1922 96 : const Eigen::Index supportSize = syevrSize > 0 ? 2 * std::max<Eigen::Index>( 1, syevrSize ) : 0;
1923 96 : Eigen::Index minimumGesvdWork{ 0 };
1924 96 : bool validShapes =
1925 288 : detail::validArrayShape<realT>( maximumDeleted, baseRank ) && detail::validArrayShape<realT>( baseRank, 1 ) &&
1926 384 : detail::validArrayShape<realT>( eigenSize, eigenSize ) && detail::validArrayShape<realT>( eigenSize, 1 ) &&
1927 96 : detail::validVectorSize<MXLAPACK_INT>( supportSize );
1928 96 : if( backend == svdDeletionBackend::leadingCovariance )
1929 : {
1930 33 : validShapes = validShapes && detail::validArrayShape<realT>( baseRank, maximumDeleted );
1931 : }
1932 63 : else if( backend == svdDeletionBackend::stableCore )
1933 : {
1934 80 : validShapes = validShapes && detail::validArrayShape<realT>( maximumDeleted, maximumDeleted ) &&
1935 120 : detail::validArrayShape<realT>( baseRank + maximumDeleted, baseRank ) &&
1936 39 : detail::validArrayShape<realT>( baseRank, baseRank );
1937 40 : const unsigned long long rank = static_cast<unsigned long long>( baseRank );
1938 40 : const unsigned long long deleted = static_cast<unsigned long long>( maximumDeleted );
1939 40 : const unsigned long long lapackMaximum =
1940 : static_cast<unsigned long long>( std::numeric_limits<MXLAPACK_INT>::max() );
1941 40 : validShapes = validShapes && rank <= lapackMaximum / 5 && rank <= ( lapackMaximum - deleted ) / 4;
1942 40 : if( validShapes )
1943 : {
1944 39 : minimumGesvdWork = std::max<Eigen::Index>( 1, std::max( 5 * baseRank, 4 * baseRank + maximumDeleted ) );
1945 39 : validShapes = detail::validVectorSize<realT>( minimumGesvdWork );
1946 : }
1947 : }
1948 : else
1949 : {
1950 69 : validShapes = validShapes && detail::validArrayShape<realT>( baseRank, baseRank ) &&
1951 69 : detail::validVectorSize<realT>( baseRank ) && detail::validVectorSize<Eigen::Index>( baseRank );
1952 : }
1953 96 : if( !validShapes )
1954 : {
1955 1 : return svdDeletionStatus::invalidInput;
1956 : }
1957 :
1958 : try
1959 : {
1960 95 : detail::callOperationHook<realT>( detail::svdDeletionTestOperation::prepareWorkspace );
1961 :
1962 92 : stored.m_deletedRows.resize( maximumDeleted, baseRank );
1963 92 : stored.m_normalizedSingularValues.resize( baseRank );
1964 92 : stored.m_scaledDeletedTranspose.resize( 0, 0 );
1965 92 : stored.m_complementRoot.resize( 0, 0 );
1966 92 : stored.m_stableCore.resize( 0, 0 );
1967 92 : stored.m_rightTranspose.resize( 0, 0 );
1968 92 : stored.m_stableSingularValues.resize( 0 );
1969 92 : stored.m_secularPoles.resize( 0 );
1970 92 : stored.m_secularUpdate.resize( 0 );
1971 92 : stored.m_secularActivePoles.resize( 0 );
1972 92 : stored.m_secularActiveUpdate.resize( 0 );
1973 92 : stored.m_secularCandidateValues.resize( 0 );
1974 92 : stored.m_secularVector.resize( 0 );
1975 92 : stored.m_secularActiveIndices.clear();
1976 92 : stored.m_secularDeflatedIndices.clear();
1977 92 : stored.m_secularOrder.clear();
1978 92 : stored.m_secularRotationFirst.clear();
1979 92 : stored.m_secularRotationSecond.clear();
1980 92 : stored.m_secularRotationCosine.clear();
1981 92 : stored.m_secularRotationSine.clear();
1982 :
1983 92 : if( backend == svdDeletionBackend::leadingCovariance )
1984 : {
1985 33 : stored.m_scaledDeletedTranspose.resize( baseRank, maximumDeleted );
1986 33 : stored.m_psdCore.resize( baseRank, baseRank );
1987 : }
1988 59 : else if( backend == svdDeletionBackend::stableCore )
1989 : {
1990 36 : stored.m_psdCore.resize( maximumDeleted, maximumDeleted );
1991 36 : stored.m_complementRoot.resize( maximumDeleted, maximumDeleted );
1992 36 : stored.m_stableCore.resize( baseRank + maximumDeleted, baseRank );
1993 36 : stored.m_rightTranspose.resize( baseRank, baseRank );
1994 36 : stored.m_stableSingularValues.resize( baseRank );
1995 : }
1996 : else
1997 : {
1998 23 : stored.m_scaledDeletedTranspose.resize( baseRank, 1 );
1999 23 : stored.m_psdCore.resize( baseRank, baseRank );
2000 23 : stored.m_secularPoles.resize( baseRank );
2001 23 : stored.m_secularUpdate.resize( baseRank );
2002 23 : stored.m_secularActivePoles.resize( baseRank );
2003 23 : stored.m_secularActiveUpdate.resize( baseRank );
2004 23 : stored.m_secularCandidateValues.resize( baseRank );
2005 23 : stored.m_secularVector.resize( baseRank );
2006 23 : stored.m_secularActiveIndices.resize( static_cast<std::size_t>( baseRank ) );
2007 23 : stored.m_secularDeflatedIndices.resize( static_cast<std::size_t>( baseRank ) );
2008 23 : stored.m_secularOrder.resize( static_cast<std::size_t>( baseRank ) );
2009 23 : stored.m_secularRotationFirst.resize( static_cast<std::size_t>( baseRank ) );
2010 23 : stored.m_secularRotationSecond.resize( static_cast<std::size_t>( baseRank ) );
2011 23 : stored.m_secularRotationCosine.resize( static_cast<std::size_t>( baseRank ) );
2012 23 : stored.m_secularRotationSine.resize( static_cast<std::size_t>( baseRank ) );
2013 : }
2014 :
2015 92 : stored.m_eigenvectors.resize( eigenSize, eigenSize );
2016 92 : stored.m_eigenvalues.resize( eigenSize );
2017 92 : stored.m_support.resize( static_cast<std::size_t>( supportSize ) );
2018 92 : stored.m_syevrWork.clear();
2019 92 : stored.m_syevrIWork.clear();
2020 92 : stored.m_gesvdWork.clear();
2021 :
2022 92 : if( syevrSize > 0 )
2023 : {
2024 69 : stored.m_psdCore.matrix().setIdentity();
2025 69 : stored.m_eigenvectors.setZero();
2026 69 : stored.m_eigenvalues.setZero();
2027 69 : realT workQuery{ 0 };
2028 69 : MXLAPACK_INT integerWorkQuery{ 0 };
2029 69 : MXLAPACK_INT found{ 0 };
2030 69 : const MXLAPACK_INT lapackEigenSize = static_cast<MXLAPACK_INT>( syevrSize );
2031 69 : const MXLAPACK_INT info = detail::callSyevr<realT>( 'V',
2032 : 'A',
2033 : 'L',
2034 : lapackEigenSize,
2035 : stored.m_psdCore.data(),
2036 : lapackEigenSize,
2037 : realT( 0 ),
2038 : realT( 0 ),
2039 : 1,
2040 : lapackEigenSize,
2041 : math::lamch<realT>( 'S' ),
2042 : &found,
2043 : stored.m_eigenvalues.data(),
2044 : stored.m_eigenvectors.data(),
2045 : lapackEigenSize,
2046 : stored.m_support.data(),
2047 : &workQuery,
2048 : -1,
2049 : &integerWorkQuery,
2050 : -1 );
2051 69 : MXLAPACK_INT realWorkSize{ 0 };
2052 67 : if( info != 0 || !detail::querySize( realWorkSize, workQuery ) || integerWorkQuery < 1 ||
2053 201 : !detail::validVectorSize<realT>( realWorkSize ) ||
2054 65 : !detail::validVectorSize<MXLAPACK_INT>( integerWorkQuery ) )
2055 : {
2056 4 : clear();
2057 4 : stored.m_lapackInfo = info;
2058 4 : return svdDeletionStatus::workspaceQueryFailure;
2059 : }
2060 65 : stored.m_syevrWork.resize( static_cast<std::size_t>( realWorkSize ) );
2061 65 : stored.m_syevrIWork.resize( static_cast<std::size_t>( integerWorkQuery ) );
2062 : }
2063 :
2064 88 : if( backend == svdDeletionBackend::stableCore )
2065 : {
2066 35 : stored.m_stableCore.setZero();
2067 35 : stored.m_rightTranspose.setZero();
2068 35 : stored.m_stableSingularValues.setZero();
2069 35 : realT unusedLeft{ 0 };
2070 35 : realT workQuery{ 0 };
2071 35 : const MXLAPACK_INT rows = static_cast<MXLAPACK_INT>( baseRank + maximumDeleted );
2072 35 : const MXLAPACK_INT cols = static_cast<MXLAPACK_INT>( baseRank );
2073 35 : const MXLAPACK_INT info = detail::callGesvd<realT>( 'N',
2074 : 'S',
2075 : rows,
2076 : cols,
2077 : stored.m_stableCore.data(),
2078 : rows,
2079 : stored.m_stableSingularValues.data(),
2080 : &unusedLeft,
2081 : 1,
2082 : stored.m_rightTranspose.data(),
2083 : cols,
2084 : &workQuery,
2085 : -1 );
2086 35 : MXLAPACK_INT workSize{ 0 };
2087 35 : if( info != 0 || !detail::querySize( workSize, workQuery ) )
2088 : {
2089 2 : clear();
2090 2 : stored.m_lapackInfo = info;
2091 2 : return svdDeletionStatus::workspaceQueryFailure;
2092 : }
2093 33 : workSize = std::max( workSize, static_cast<MXLAPACK_INT>( minimumGesvdWork ) );
2094 33 : if( !detail::validVectorSize<realT>( workSize ) )
2095 : {
2096 : clear(); // LCOV_EXCL_LINE -- requires size_t narrower than the configured LAPACK integer ABI.
2097 : return svdDeletionStatus::workspaceQueryFailure; // LCOV_EXCL_LINE
2098 : }
2099 33 : stored.m_gesvdWork.resize( static_cast<std::size_t>( workSize ) );
2100 : }
2101 : }
2102 3 : catch( const std::bad_alloc & )
2103 : {
2104 2 : clear();
2105 2 : return svdDeletionStatus::allocationFailure;
2106 : }
2107 1 : catch( const std::length_error & )
2108 : {
2109 1 : clear();
2110 1 : return svdDeletionStatus::allocationFailure;
2111 : }
2112 :
2113 86 : stored.m_baseRank = baseRank;
2114 86 : stored.m_maximumDeleted = maximumDeleted;
2115 86 : stored.m_backend = backend;
2116 86 : stored.m_prepared = true;
2117 86 : return svdDeletionStatus::success;
2118 : }
2119 :
2120 : template <typename realT, typename abiT>
2121 12 : void svdDeletionWorkspace<realT, abiT>::clear() noexcept
2122 : {
2123 12 : if( !m_storage )
2124 : {
2125 1 : return;
2126 : }
2127 :
2128 11 : m_storage->m_deletedRows.resize( 0, 0 );
2129 11 : m_storage->m_normalizedSingularValues.resize( 0 );
2130 11 : m_storage->m_scaledDeletedTranspose.resize( 0, 0 );
2131 11 : m_storage->m_psdCore.resize( 0, 0 );
2132 11 : m_storage->m_eigenvectors.resize( 0, 0 );
2133 11 : m_storage->m_eigenvalues.resize( 0 );
2134 11 : m_storage->m_complementRoot.resize( 0, 0 );
2135 11 : m_storage->m_stableCore.resize( 0, 0 );
2136 11 : m_storage->m_rightTranspose.resize( 0, 0 );
2137 11 : m_storage->m_stableSingularValues.resize( 0 );
2138 11 : m_storage->m_secularPoles.resize( 0 );
2139 11 : m_storage->m_secularUpdate.resize( 0 );
2140 11 : m_storage->m_secularActivePoles.resize( 0 );
2141 11 : m_storage->m_secularActiveUpdate.resize( 0 );
2142 11 : m_storage->m_secularCandidateValues.resize( 0 );
2143 11 : m_storage->m_secularVector.resize( 0 );
2144 11 : m_storage->m_secularActiveIndices.clear();
2145 11 : m_storage->m_secularDeflatedIndices.clear();
2146 11 : m_storage->m_secularOrder.clear();
2147 11 : m_storage->m_secularRotationFirst.clear();
2148 11 : m_storage->m_secularRotationSecond.clear();
2149 11 : m_storage->m_secularRotationCosine.clear();
2150 11 : m_storage->m_secularRotationSine.clear();
2151 11 : m_storage->m_support.clear();
2152 11 : m_storage->m_syevrWork.clear();
2153 11 : m_storage->m_syevrIWork.clear();
2154 11 : m_storage->m_gesvdWork.clear();
2155 11 : m_storage->m_baseRank = 0;
2156 11 : m_storage->m_maximumDeleted = 0;
2157 11 : m_storage->m_backend = svdDeletionBackend::stableCore;
2158 11 : m_storage->m_prepared = false;
2159 11 : m_storage->m_lapackInfo = 0;
2160 : }
2161 :
2162 : template <typename realT, typename abiT>
2163 107 : bool svdDeletionWorkspace<realT, abiT>::prepared() const noexcept
2164 : {
2165 107 : return m_storage && m_storage->m_prepared;
2166 : }
2167 :
2168 : template <typename realT, typename abiT>
2169 36 : std::int64_t svdDeletionWorkspace<realT, abiT>::baseRank() const noexcept
2170 : {
2171 36 : return m_storage ? m_storage->m_baseRank : 0;
2172 : }
2173 :
2174 : template <typename realT, typename abiT>
2175 38 : std::int64_t svdDeletionWorkspace<realT, abiT>::maximumDeleted() const noexcept
2176 : {
2177 38 : return m_storage ? m_storage->m_maximumDeleted : 0;
2178 : }
2179 :
2180 : template <typename realT, typename abiT>
2181 34 : svdDeletionBackend svdDeletionWorkspace<realT, abiT>::backend() const noexcept
2182 : {
2183 34 : return m_storage ? m_storage->m_backend : svdDeletionBackend::stableCore;
2184 : }
2185 :
2186 : template <typename realT, typename abiT>
2187 8 : MXLAPACK_INT svdDeletionWorkspace<realT, abiT>::lapackInfo() const noexcept
2188 : {
2189 8 : return m_storage ? m_storage->m_lapackInfo : 0;
2190 : }
2191 :
2192 : namespace
2193 : {
2194 :
2195 : template <typename realT>
2196 14 : svdDeletionStatus validateSvdDeletionFactorImpl( svdDeletionConstMatrixRef<realT> factor, realT tolerance )
2197 : {
2198 27 : if( factor.rows() < factor.cols() || factor.cols() <= 0 || !detail::validLapackDimension( factor.rows() ) ||
2199 12 : !detail::validLapackDimension( factor.cols() ) || !math::isFinite( tolerance ) || tolerance < realT( 0 ) ||
2200 27 : !detail::validArrayShape<realT>( factor.cols(), factor.cols() ) || !detail::allFinite( factor ) )
2201 : {
2202 4 : return svdDeletionStatus::invalidInput;
2203 : }
2204 :
2205 : try
2206 : {
2207 10 : detail::callOperationHook<realT>( detail::svdDeletionTestOperation::validateFactor );
2208 8 : svdDeletionStorageMatrix<realT> gram( factor.cols(), factor.cols() );
2209 8 : gram.matrix().noalias() = factor.matrix().transpose() * factor.matrix();
2210 8 : if( !detail::allFinite( gram ) )
2211 : {
2212 1 : return svdDeletionStatus::invalidInput;
2213 : }
2214 26 : for( Eigen::Index index = 0; index < gram.rows(); ++index )
2215 : {
2216 19 : gram( index, index ) -= realT( 1 );
2217 : }
2218 :
2219 7 : if( tolerance == realT( 0 ) )
2220 : {
2221 7 : tolerance = realT( 64 ) * std::numeric_limits<realT>::epsilon() *
2222 7 : static_cast<realT>( std::max( factor.rows(), factor.cols() ) );
2223 : }
2224 7 : if( gram.abs().maxCoeff() > tolerance )
2225 : {
2226 1 : return svdDeletionStatus::factorNotOrthonormal;
2227 : }
2228 8 : }
2229 2 : catch( const std::bad_alloc & )
2230 : {
2231 1 : return svdDeletionStatus::allocationFailure;
2232 : }
2233 1 : catch( const std::length_error & )
2234 : {
2235 1 : return svdDeletionStatus::allocationFailure;
2236 : }
2237 :
2238 6 : return svdDeletionStatus::success;
2239 : }
2240 :
2241 : } // namespace
2242 :
2243 2 : svdDeletionStatus detail::validateSvdDeletionFactorAbiV2( svdDeletionConstMatrixViewV2<float> factor, float tolerance )
2244 : {
2245 2 : if( !detail::validAbiMatrixView( factor ) )
2246 : {
2247 1 : return svdDeletionStatus::invalidInput;
2248 : }
2249 1 : return validateSvdDeletionFactorImpl<float>( detail::mapAbiMatrix( factor ), tolerance );
2250 : }
2251 :
2252 19 : svdDeletionStatus detail::validateSvdDeletionFactorAbiV2( svdDeletionConstMatrixViewV2<double> factor,
2253 : double tolerance )
2254 : {
2255 19 : if( !detail::validAbiMatrixView( factor ) )
2256 : {
2257 6 : return svdDeletionStatus::invalidInput;
2258 : }
2259 13 : return validateSvdDeletionFactorImpl<double>( detail::mapAbiMatrix( factor ), tolerance );
2260 : }
2261 :
2262 : template <typename realT>
2263 14 : svdDeletionStatus detail::svdDeletionLeadingCoreAbiV2( svdDeletionResult<realT> &result,
2264 : svdDeletionConstVectorViewV2<realT> singularValues,
2265 : svdDeletionConstMatrixViewV2<realT> deletedRows,
2266 : std::int64_t outputRank,
2267 : svdDeletionWorkspace<realT> &workspace )
2268 : {
2269 26 : if( !detail::validAbiVectorView( singularValues ) || !detail::validAbiMatrixView( deletedRows ) ||
2270 12 : !detail::validAbiDimension( outputRank ) )
2271 : {
2272 2 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2273 : }
2274 12 : return detail::svdDeletionImplementation<realT>::leading( result,
2275 24 : detail::mapAbiVector( singularValues ),
2276 24 : detail::mapAbiMatrix( deletedRows ),
2277 : static_cast<Eigen::Index>( outputRank ),
2278 12 : workspace );
2279 : }
2280 :
2281 : template <typename realT>
2282 15 : svdDeletionStatus detail::svdDeletionStableCoreAbiV2( svdDeletionResult<realT> &result,
2283 : svdDeletionConstVectorViewV2<realT> singularValues,
2284 : svdDeletionConstMatrixViewV2<realT> deletedRows,
2285 : std::int64_t outputRank,
2286 : svdDeletionWorkspace<realT> &workspace )
2287 : {
2288 29 : if( !detail::validAbiVectorView( singularValues ) || !detail::validAbiMatrixView( deletedRows ) ||
2289 14 : !detail::validAbiDimension( outputRank ) )
2290 : {
2291 1 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2292 : }
2293 14 : return detail::svdDeletionImplementation<realT>::stable( result,
2294 28 : detail::mapAbiVector( singularValues ),
2295 28 : detail::mapAbiMatrix( deletedRows ),
2296 : static_cast<Eigen::Index>( outputRank ),
2297 14 : workspace );
2298 : }
2299 :
2300 : template <typename realT>
2301 32 : svdDeletionStatus detail::svdDeletionCoreAbiV2( svdDeletionResult<realT> &result,
2302 : svdDeletionConstVectorViewV2<realT> singularValues,
2303 : svdDeletionConstMatrixViewV2<realT> deletedRows,
2304 : std::int64_t outputRank,
2305 : svdDeletionWorkspace<realT> &workspace,
2306 : svdDeletionBackend backend )
2307 : {
2308 64 : if( !detail::validAbiVectorView( singularValues ) || !detail::validAbiMatrixView( deletedRows ) ||
2309 64 : !detail::validAbiDimension( outputRank ) || !detail::svdDeletionImplementation<realT>::validBackend( backend ) )
2310 : {
2311 1 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2312 : }
2313 31 : switch( backend )
2314 : {
2315 3 : case svdDeletionBackend::leadingCovariance:
2316 3 : return detail::svdDeletionImplementation<realT>::leading( result,
2317 6 : detail::mapAbiVector( singularValues ),
2318 6 : detail::mapAbiMatrix( deletedRows ),
2319 : static_cast<Eigen::Index>( outputRank ),
2320 3 : workspace );
2321 11 : case svdDeletionBackend::stableCore:
2322 11 : return detail::svdDeletionImplementation<realT>::stable( result,
2323 22 : detail::mapAbiVector( singularValues ),
2324 22 : detail::mapAbiMatrix( deletedRows ),
2325 : static_cast<Eigen::Index>( outputRank ),
2326 11 : workspace );
2327 17 : case svdDeletionBackend::rankOneSecular:
2328 17 : return detail::svdDeletionImplementation<realT>::rankOne( result,
2329 34 : detail::mapAbiVector( singularValues ),
2330 34 : detail::mapAbiMatrix( deletedRows ),
2331 : static_cast<Eigen::Index>( outputRank ),
2332 17 : workspace );
2333 : }
2334 : // LCOV_EXCL_START -- the backend is validated before the exhaustive switch.
2335 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2336 : // LCOV_EXCL_STOP
2337 : }
2338 :
2339 : template <typename realT>
2340 63 : svdDeletionStatus detail::svdRemoveRowsAbiV2( svdDeletionResult<realT> &result,
2341 : svdDeletionConstVectorViewV2<realT> singularValues,
2342 : svdDeletionConstMatrixViewV2<realT> leftFactor,
2343 : svdDeletionConstIndexViewV2 deletedIndices,
2344 : std::int64_t outputRank,
2345 : svdDeletionWorkspace<realT> &workspace,
2346 : svdDeletionBackend backend )
2347 : {
2348 126 : if( !detail::validAbiVectorView( singularValues ) || !detail::validAbiMatrixView( leftFactor ) ||
2349 126 : !detail::validAbiIndexView( deletedIndices ) || !detail::validAbiDimension( outputRank ) )
2350 : {
2351 3 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2352 : }
2353 60 : return detail::svdDeletionImplementation<realT>::remove( result,
2354 120 : detail::mapAbiVector( singularValues ),
2355 120 : detail::mapAbiMatrix( leftFactor ),
2356 : deletedIndices,
2357 : static_cast<Eigen::Index>( outputRank ),
2358 : workspace,
2359 60 : backend );
2360 : }
2361 :
2362 : template <typename realT>
2363 11 : svdDeletionStatus detail::svdRemoveColumnsAbiV2( svdDeletionResult<realT> &result,
2364 : svdDeletionConstVectorViewV2<realT> singularValues,
2365 : svdDeletionConstMatrixViewV2<realT> rightFactor,
2366 : svdDeletionConstIndexViewV2 deletedIndices,
2367 : std::int64_t outputRank,
2368 : svdDeletionWorkspace<realT> &workspace,
2369 : svdDeletionBackend backend )
2370 : {
2371 22 : if( !detail::validAbiVectorView( singularValues ) || !detail::validAbiMatrixView( rightFactor ) ||
2372 22 : !detail::validAbiIndexView( deletedIndices ) || !detail::validAbiDimension( outputRank ) )
2373 : {
2374 1 : return detail::svdDeletionImplementation<realT>::fail( result, svdDeletionStatus::invalidInput );
2375 : }
2376 10 : return detail::svdDeletionImplementation<realT>::remove( result,
2377 20 : detail::mapAbiVector( singularValues ),
2378 20 : detail::mapAbiMatrix( rightFactor ),
2379 : deletedIndices,
2380 : static_cast<Eigen::Index>( outputRank ),
2381 : workspace,
2382 10 : backend );
2383 : }
2384 :
2385 : template class svdDeletionResult<float>;
2386 : template class svdDeletionResult<double>;
2387 : template class svdDeletionWorkspace<float>;
2388 : template class svdDeletionWorkspace<double>;
2389 :
2390 : template svdDeletionStatus detail::svdDeletionLeadingCoreAbiV2<float>( svdDeletionResult<float> &,
2391 : svdDeletionConstVectorViewV2<float>,
2392 : svdDeletionConstMatrixViewV2<float>,
2393 : std::int64_t,
2394 : svdDeletionWorkspace<float> & );
2395 : template svdDeletionStatus detail::svdDeletionLeadingCoreAbiV2<double>( svdDeletionResult<double> &,
2396 : svdDeletionConstVectorViewV2<double>,
2397 : svdDeletionConstMatrixViewV2<double>,
2398 : std::int64_t,
2399 : svdDeletionWorkspace<double> & );
2400 :
2401 : template svdDeletionStatus detail::svdDeletionStableCoreAbiV2<float>( svdDeletionResult<float> &,
2402 : svdDeletionConstVectorViewV2<float>,
2403 : svdDeletionConstMatrixViewV2<float>,
2404 : std::int64_t,
2405 : svdDeletionWorkspace<float> & );
2406 : template svdDeletionStatus detail::svdDeletionStableCoreAbiV2<double>( svdDeletionResult<double> &,
2407 : svdDeletionConstVectorViewV2<double>,
2408 : svdDeletionConstMatrixViewV2<double>,
2409 : std::int64_t,
2410 : svdDeletionWorkspace<double> & );
2411 :
2412 : template svdDeletionStatus detail::svdDeletionCoreAbiV2<float>( svdDeletionResult<float> &,
2413 : svdDeletionConstVectorViewV2<float>,
2414 : svdDeletionConstMatrixViewV2<float>,
2415 : std::int64_t,
2416 : svdDeletionWorkspace<float> &,
2417 : svdDeletionBackend );
2418 : template svdDeletionStatus detail::svdDeletionCoreAbiV2<double>( svdDeletionResult<double> &,
2419 : svdDeletionConstVectorViewV2<double>,
2420 : svdDeletionConstMatrixViewV2<double>,
2421 : std::int64_t,
2422 : svdDeletionWorkspace<double> &,
2423 : svdDeletionBackend );
2424 :
2425 : template svdDeletionStatus detail::svdRemoveRowsAbiV2<float>( svdDeletionResult<float> &,
2426 : svdDeletionConstVectorViewV2<float>,
2427 : svdDeletionConstMatrixViewV2<float>,
2428 : svdDeletionConstIndexViewV2,
2429 : std::int64_t,
2430 : svdDeletionWorkspace<float> &,
2431 : svdDeletionBackend );
2432 : template svdDeletionStatus detail::svdRemoveRowsAbiV2<double>( svdDeletionResult<double> &,
2433 : svdDeletionConstVectorViewV2<double>,
2434 : svdDeletionConstMatrixViewV2<double>,
2435 : svdDeletionConstIndexViewV2,
2436 : std::int64_t,
2437 : svdDeletionWorkspace<double> &,
2438 : svdDeletionBackend );
2439 :
2440 : template svdDeletionStatus detail::svdRemoveColumnsAbiV2<float>( svdDeletionResult<float> &,
2441 : svdDeletionConstVectorViewV2<float>,
2442 : svdDeletionConstMatrixViewV2<float>,
2443 : svdDeletionConstIndexViewV2,
2444 : std::int64_t,
2445 : svdDeletionWorkspace<float> &,
2446 : svdDeletionBackend );
2447 : template svdDeletionStatus detail::svdRemoveColumnsAbiV2<double>( svdDeletionResult<double> &,
2448 : svdDeletionConstVectorViewV2<double>,
2449 : svdDeletionConstMatrixViewV2<double>,
2450 : svdDeletionConstIndexViewV2,
2451 : std::int64_t,
2452 : svdDeletionWorkspace<double> &,
2453 : svdDeletionBackend );
2454 :
2455 : template detail::svdDeletionTestHooks<float> &detail::svdDeletionHooks<float>();
2456 : template detail::svdDeletionTestHooks<double> &detail::svdDeletionHooks<double>();
2457 :
2458 : } // namespace math
2459 : } // namespace mx
|
