mxlib
c++ tools for analyzing astronomical data and other tasks by Jared R. Males. [git repo]
Loading...
Searching...
No Matches
eigenLapack_test.cpp
Go to the documentation of this file.
1/** \file eigenLapack_test.cpp
2 * \brief Tests Eigen-compatible BLAS/LAPACK helpers.
3 */
4#include "../../catch2/catch.hpp"
5
6#define MX_NO_ERROR_REPORTS
7
8#include "../../../include/math/eigenLapack.hpp"
9
10#include <algorithm>
11#include <cstdlib>
12#include <limits>
13#include <type_traits>
14#include <vector>
15
16/** \cond */
17/// Compile the symmetric-eigensolver workspace implementation for a representative scalar type.
18template struct mx::math::syevrMem<double>;
19
20/** Verify that the owning mx::math::syevrMem workspace cannot be copied. */
21TEST_CASE( "syevrMem workspaces are non-copyable", "[math::syevrMem]" )
22{
23 REQUIRE_FALSE( std::is_copy_constructible_v<mx::math::syevrMem<double>> );
24 REQUIRE_FALSE( std::is_copy_assignable_v<mx::math::syevrMem<double>> );
25}
26
27using matrixT = Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic>;
28
29int allocationCall{ 0 };
30int failAllocationOnCall{ 0 };
31int queryCall{ 0 };
32int solveCall{ 0 };
33MXLAPACK_INT queryInfo{ 0 };
34MXLAPACK_INT solveInfo{ 0 };
35bool overrideQuerySizes{ false };
36double injectedWorkSize{ 0.0 };
37MXLAPACK_INT injectedIntegerWorkSize{ 0 };
38std::vector<double> injectedEigenvalues;
39int nonfiniteEigenvectorRow{ -1 };
40int nonfiniteEigenvectorColumn{ -1 };
41
42/// Reset all deterministic eigenLapack failure controls.
43void resetControls()
44{
45 allocationCall = 0;
46 failAllocationOnCall = 0;
47 queryCall = 0;
48 solveCall = 0;
49 queryInfo = 0;
50 solveInfo = 0;
51 overrideQuerySizes = false;
52 injectedWorkSize = 0.0;
53 injectedIntegerWorkSize = 0;
54 injectedEigenvalues.clear();
55 nonfiniteEigenvectorRow = -1;
56 nonfiniteEigenvectorColumn = -1;
57 mx::math::detail::eigenLapackTestHooks<double>::allocator = nullptr;
58 mx::math::detail::eigenLapackTestHooks<double>::solver = nullptr;
59}
60
61/// Restore production allocation and solver functions at scope exit.
62struct hookReset
63{
64 /// Start a test scope from the production hook state.
65 hookReset()
66 {
67 resetControls();
68 }
69
70 /// Restore the production hook state even when a Catch2 assertion aborts the scope.
71 ~hookReset()
72 {
73 resetControls();
74 }
75};
76
77/// Return a diagonal square matrix with the supplied diagonal entries.
78matrixT diagonalMatrix( const std::vector<double> &diagonal /**< [in] diagonal entries */ )
79{
80 matrixT matrix( diagonal.size(), diagonal.size() );
81 matrix.setZero();
82 for( std::size_t index = 0; index < diagonal.size(); ++index )
83 {
84 matrix( index, index ) = diagonal[index];
85 }
86 return matrix;
87}
88
89/// Allocate workspace storage unless the configured allocation call must fail.
90void *controlledAllocation( std::size_t bytes /**< [in] requested byte count */ )
91{
92 ++allocationCall;
93 if( allocationCall == failAllocationOnCall )
94 {
95 return nullptr;
96 }
97 return std::malloc( bytes );
98}
99
100/// Provide controlled SYEVR query, failure, and result behavior.
101MXLAPACK_INT controlledSyevr( char jobz, /**< [in] eigenvector request */
102 char range, /**< [in] eigenvalue selection mode */
103 char uplo, /**< [in] populated input triangle */
104 MXLAPACK_INT n, /**< [in] matrix order */
105 double *matrix, /**< [in,out] input matrix */
106 MXLAPACK_INT lda, /**< [in] input leading dimension */
107 double valueLower, /**< [in] lower value bound */
108 double valueUpper, /**< [in] upper value bound */
109 MXLAPACK_INT indexLow, /**< [in] one-based lower index */
110 MXLAPACK_INT indexHigh, /**< [in] one-based upper index */
111 double tolerance, /**< [in] convergence tolerance */
112 MXLAPACK_INT *found, /**< [out] selected eigenvalue count */
113 double *eigenvalues, /**< [out] eigenvalues */
114 double *eigenvectors, /**< [out] eigenvectors */
115 MXLAPACK_INT ldz, /**< [in] eigenvector leading dimension */
116 MXLAPACK_INT *support, /**< [out] eigenvector support */
117 double *work, /**< [in,out] floating workspace */
118 MXLAPACK_INT lwork, /**< [in] floating workspace size */
119 MXLAPACK_INT *iwork, /**< [in,out] integer workspace */
120 MXLAPACK_INT liwork /**< [in] integer workspace size */ )
121{
122 static_cast<void>( jobz );
123 static_cast<void>( uplo );
124 static_cast<void>( matrix );
125 static_cast<void>( lda );
126 static_cast<void>( valueLower );
127 static_cast<void>( valueUpper );
128 static_cast<void>( tolerance );
129 static_cast<void>( support );
130
131 if( lwork == -1 || liwork == -1 )
132 {
133 ++queryCall;
134 if( queryInfo != 0 )
135 {
136 return queryInfo;
137 }
138
139 work[0] = overrideQuerySizes ? injectedWorkSize : std::max<MXLAPACK_INT>( 1, 26 * n );
140 iwork[0] = overrideQuerySizes ? injectedIntegerWorkSize : std::max<MXLAPACK_INT>( 1, 10 * n );
141 return 0;
142 }
143
144 ++solveCall;
145 if( solveInfo != 0 )
146 {
147 return solveInfo;
148 }
149
150 const MXLAPACK_INT first = range == 'I' ? indexLow - 1 : 0;
151 const MXLAPACK_INT count = range == 'I' ? indexHigh - indexLow + 1 : n;
152 *found = count;
153
154 for( MXLAPACK_INT column = 0; column < count; ++column )
155 {
156 const MXLAPACK_INT sourceColumn = first + column;
157 eigenvalues[column] = injectedEigenvalues.empty() ? sourceColumn + 1 : injectedEigenvalues[sourceColumn];
158
159 for( MXLAPACK_INT row = 0; row < n; ++row )
160 {
161 eigenvectors[column * ldz + row] = row == sourceColumn ? 1.0 : 0.0;
162 }
163 }
164
165 if( nonfiniteEigenvectorRow >= 0 && nonfiniteEigenvectorColumn >= 0 )
166 {
167 eigenvectors[nonfiniteEigenvectorColumn * ldz + nonfiniteEigenvectorRow] =
168 std::numeric_limits<double>::infinity();
169 }
170
171 return 0;
172}
173/** \endcond */
174
175/** \defgroup eigenLapack_unit_tests eigenLapack Unit Tests
176 * \ingroup math_unit_tests
177 */
178
179namespace unitTest::math_eigenLapack_test
180{
181
182/** \brief Verifies eigenSYEVR full/ranged solves, triangle selection, workspace reuse, and explicit release.
183 *
184 * \ingroup eigenLapack_unit_tests
185 */
186TEST_CASE( "eigenSYEVR solves selected ranges and reuses workspaces", "[math::eigenLapack][eigenSYEVR]" )
187{
188 hookReset reset;
189
190 matrixT eigenvectors;
191 matrixT eigenvalues;
192
193 matrixT covariance = diagonalMatrix( { 1.0, 4.0, 9.0 } );
194 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance ) == 0 );
195 REQUIRE( eigenvectors.rows() == 3 );
196 REQUIRE( eigenvectors.cols() == 3 );
197 REQUIRE( eigenvalues.rows() == 3 );
198 REQUIRE( eigenvalues( 0 ) == Approx( 1.0 ) );
199 REQUIRE( eigenvalues( 1 ) == Approx( 4.0 ) );
200 REQUIRE( eigenvalues( 2 ) == Approx( 9.0 ) );
201
202 mx::math::syevrMem<double> workspace;
203 covariance = diagonalMatrix( { 1.0, 4.0, 9.0 } );
204 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 1, 3, 'L', &workspace ) == 0 );
205 REQUIRE( eigenvectors.rows() == 3 );
206 REQUIRE( eigenvectors.cols() == 2 );
207 REQUIRE( eigenvalues( 0 ) == Approx( 4.0 ) );
208 REQUIRE( eigenvalues( 1 ) == Approx( 9.0 ) );
209
210 const auto *support = workspace.iSuppZ;
211 const auto *minimumWork = workspace.minWork;
212 const auto *work = workspace.work;
213 const auto *minimumIntegerWork = workspace.minIWork;
214 const auto *integerWork = workspace.iWork;
215
216 covariance = diagonalMatrix( { 1.0, 4.0, 9.0 } );
217 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 1, 3, 'L', &workspace ) == 0 );
218 REQUIRE( workspace.iSuppZ == support );
219 REQUIRE( workspace.minWork == minimumWork );
220 REQUIRE( workspace.work == work );
221 REQUIRE( workspace.minIWork == minimumIntegerWork );
222 REQUIRE( workspace.iWork == integerWork );
223
224 covariance = diagonalMatrix( { 1.0, 2.0, 3.0, 4.0, 5.0, 6.0 } );
225 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -1, 'L', &workspace ) == 0 );
226 REQUIRE( workspace.sizeISuppZ >= 12 );
227 REQUIRE( workspace.sizeMinWork >= 156 );
228 REQUIRE( workspace.sizeMinIWork >= 60 );
229
230 matrixT upperCovariance( 2, 2 );
231 upperCovariance( 0, 0 ) = 2.0;
232 upperCovariance( 0, 1 ) = 1.0;
233 upperCovariance( 1, 0 ) = 99.0;
234 upperCovariance( 1, 1 ) = 2.0;
235 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, upperCovariance, 0, -1, 'U', &workspace ) == 0 );
236 REQUIRE( eigenvalues( 0 ) == Approx( 1.0 ) );
237 REQUIRE( eigenvalues( 1 ) == Approx( 3.0 ) );
238
239 workspace.free();
240 REQUIRE( workspace.iSuppZ == nullptr );
241 REQUIRE( workspace.minWork == nullptr );
242 REQUIRE( workspace.work == nullptr );
243 REQUIRE( workspace.minIWork == nullptr );
244 REQUIRE( workspace.iWork == nullptr );
245 REQUIRE( workspace.sizeISuppZ == 0 );
246 REQUIRE( workspace.sizeMinWork == 0 );
247 REQUIRE( workspace.sizeWork == 0 );
248 REQUIRE( workspace.sizeMinIWork == 0 );
249 REQUIRE( workspace.sizeIWork == 0 );
250 REQUIRE( workspace.n == 0 );
251
252 covariance = diagonalMatrix( { 2.0, 5.0 } );
253 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -1, 'L', &workspace ) == 0 );
254 REQUIRE( eigenvalues( 0 ) == Approx( 2.0 ) );
255 REQUIRE( eigenvalues( 1 ) == Approx( 5.0 ) );
256
257 workspace.free();
258 workspace.free();
259}
260
261/** \brief Verifies eigenSYEVR rejects empty, non-square, and invalid half-open eigenvalue ranges before LAPACK.
262 *
263 * \ingroup eigenLapack_unit_tests
264 */
265TEST_CASE( "eigenSYEVR rejects invalid matrix and range geometry", "[math::eigenLapack][eigenSYEVR]" )
266{
267 hookReset reset;
268
269 matrixT eigenvectors;
270 matrixT eigenvalues;
271
272 matrixT empty( 0, 0 );
273 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, empty ) == -1 );
274
275 matrixT rectangular( 2, 3 );
276 rectangular.setZero();
277 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, rectangular ) == -1 );
278
279 matrixT covariance = diagonalMatrix( { 1.0, 2.0 } );
280 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, 0 ) == -1 );
281 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, -1, 1 ) == -1 );
282 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, 3 ) == -1 );
283 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -2 ) == -1 );
284}
285
286/** \brief Verifies eigenSYEVR propagates every workspace-allocation failure and both LAPACK solver failures.
287 *
288 * \ingroup eigenLapack_unit_tests
289 */
290TEST_CASE( "eigenSYEVR reports allocation and solver failures", "[math::eigenLapack][eigenSYEVR]" )
291{
292 hookReset reset;
293 matrixT eigenvectors;
294 matrixT eigenvalues;
295
296 mx::math::detail::eigenLapackTestHooks<double>::allocator = &controlledAllocation;
297 failAllocationOnCall = 1;
298 matrixT covariance = diagonalMatrix( { 1.0, 2.0 } );
299 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance ) == -1000 );
300
301 resetControls();
302 mx::math::detail::eigenLapackTestHooks<double>::allocator = &controlledAllocation;
303 failAllocationOnCall = 5;
304 covariance = diagonalMatrix( { 1.0, 2.0 } );
305 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance ) == -1000 );
306
307 for( int failedCall = 1; failedCall <= 5; ++failedCall )
308 {
309 resetControls();
310 mx::math::detail::eigenLapackTestHooks<double>::allocator = &controlledAllocation;
311 failAllocationOnCall = failedCall;
312
313 mx::math::syevrMem<double> workspace;
314 covariance = diagonalMatrix( { 1.0, 2.0 } );
315 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -1, 'L', &workspace ) == -1000 );
316 REQUIRE( allocationCall == failedCall );
317 }
318
319 resetControls();
320 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
321 queryInfo = 17;
322 covariance = diagonalMatrix( { 1.0, 2.0 } );
323 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance ) == 17 );
324 REQUIRE( queryCall == 1 );
325 REQUIRE( solveCall == 0 );
326
327 resetControls();
328 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
329 solveInfo = 23;
330 mx::math::syevrMem<double> workspace;
331 covariance = diagonalMatrix( { 1.0, 2.0 } );
332 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -1, 'L', &workspace ) == 23 );
333 REQUIRE( queryCall == 1 );
334 REQUIRE( solveCall == 1 );
335}
336
337/** \brief Verifies eigenSYEVR rejects nonpositive, nonfinite, and unrepresentable queried workspace sizes.
338 *
339 * \ingroup eigenLapack_unit_tests
340 */
341TEST_CASE( "eigenSYEVR validates queried workspace sizes", "[math::eigenLapack][eigenSYEVR]" )
342{
343 hookReset reset;
344 matrixT eigenvectors;
345 matrixT eigenvalues;
346 matrixT covariance;
347
348 const std::vector<double> invalidFloatingSizes{ 0.0,
349 -1.0,
350 std::numeric_limits<double>::quiet_NaN(),
351 std::numeric_limits<double>::infinity(),
352 static_cast<double>( std::numeric_limits<MXLAPACK_INT>::max() ) +
353 1.0 };
354
355 for( const double invalidSize : invalidFloatingSizes )
356 {
357 resetControls();
358 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
359 overrideQuerySizes = true;
360 injectedWorkSize = invalidSize;
361 injectedIntegerWorkSize = 10;
362 covariance = diagonalMatrix( { 1.0, 2.0 } );
363 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance ) == -1 );
364 REQUIRE( solveCall == 0 );
365 }
366
367 for( const MXLAPACK_INT invalidSize : { MXLAPACK_INT{ 0 }, MXLAPACK_INT{ -1 } } )
368 {
369 resetControls();
370 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
371 overrideQuerySizes = true;
372 injectedWorkSize = 10;
373 injectedIntegerWorkSize = invalidSize;
374 mx::math::syevrMem<double> workspace;
375 covariance = diagonalMatrix( { 1.0, 2.0 } );
376 REQUIRE( mx::math::eigenSYEVR( eigenvectors, eigenvalues, covariance, 0, -1, 'L', &workspace ) == -1 );
377 REQUIRE( solveCall == 0 );
378 }
379}
380
381/** \brief Verifies calcEigenVecs mode clamping, raw eigenvectors, scaling, and caller-owned workspace reuse.
382 *
383 * \ingroup eigenLapack_unit_tests
384 */
385TEST_CASE( "calcEigenVecs clamps mode requests and controls normalization", "[math::eigenLapack][calcEigenVecs]" )
386{
387 hookReset reset;
388 matrixT eigenvectors;
389 matrixT eigenvalues;
390 matrixT covariance = diagonalMatrix( { 1.0, 4.0, 9.0 } );
391
392 double eigenTime{ -1.0 };
393 REQUIRE( mx::math::calcEigenVecs<
394 double>( eigenvectors, eigenvalues, covariance, 1, false, false, nullptr, &eigenTime ) == 0 );
395 REQUIRE( eigenvectors.rows() == 3 );
396 REQUIRE( eigenvectors.cols() == 1 );
397 REQUIRE( eigenvalues( 0 ) == Approx( 9.0 ) );
398 REQUIRE( eigenvectors.matrix().col( 0 ).norm() == Approx( 1.0 ).margin( 1e-12 ) );
399 REQUIRE( eigenTime >= 0.0 );
400
401 mx::math::syevrMem<double> workspace;
402 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, 0, false, false, &workspace ) == 0 );
403 REQUIRE( eigenvectors.cols() == 3 );
404 REQUIRE( eigenvalues( 0 ) == Approx( 1.0 ) );
405 REQUIRE( eigenvalues( 1 ) == Approx( 4.0 ) );
406 REQUIRE( eigenvalues( 2 ) == Approx( 9.0 ) );
407
408 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, -2, false, true, &workspace ) == 0 );
409 REQUIRE( eigenvectors.cols() == 3 );
410
411 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, 8, false, false, &workspace ) == 0 );
412 REQUIRE( eigenvectors.cols() == 3 );
413
414 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, 2, true, false, &workspace ) == 0 );
415 REQUIRE( eigenvectors.cols() == 2 );
416 REQUIRE( eigenvalues( 0 ) == Approx( 4.0 ) );
417 REQUIRE( eigenvalues( 1 ) == Approx( 9.0 ) );
418 REQUIRE( eigenvectors.matrix().col( 0 ).norm() == Approx( 0.5 ).margin( 1e-12 ) );
419 REQUIRE( eigenvectors.matrix().col( 1 ).norm() == Approx( 1.0 / 3.0 ).margin( 1e-12 ) );
420}
421
422/** \brief Verifies calcEigenVecs check-mode handling of zero, negative, tiny, and nonfinite solver results.
423 *
424 * \ingroup eigenLapack_unit_tests
425 */
426TEST_CASE( "calcEigenVecs validates exceptional normalized results", "[math::eigenLapack][calcEigenVecs]" )
427{
428 hookReset reset;
429 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
430 injectedEigenvalues = { 0.0, -1.0, std::numeric_limits<double>::quiet_NaN(), 1e-20, 4.0 };
431 nonfiniteEigenvectorRow = 0;
432 nonfiniteEigenvectorColumn = 4;
433
434 matrixT covariance = diagonalMatrix( { 1.0, 2.0, 3.0, 4.0, 5.0 } );
435 matrixT eigenvectors;
436 matrixT eigenvalues;
437 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, 0, true, true ) == 0 );
438
439 REQUIRE( eigenvectors.rows() == 5 );
440 REQUIRE( eigenvectors.cols() == 5 );
441 REQUIRE( eigenvectors.matrix().col( 0 ).norm() == Approx( 0.0 ) );
442 REQUIRE( eigenvectors.matrix().col( 1 ).norm() == Approx( 0.0 ) );
443 REQUIRE( eigenvectors.matrix().col( 2 ).norm() == Approx( 0.0 ) );
444 REQUIRE( eigenvectors.matrix().col( 3 ).norm() == Approx( 1e10 ).epsilon( 1e-12 ) );
445 REQUIRE( eigenvectors.matrix().col( 4 ).norm() == Approx( 0.0 ) );
446 REQUIRE( eigenvectors.isFinite().all() );
447 REQUIRE( eigenvalues( 0 ) == Approx( 0.0 ) );
448 REQUIRE( eigenvalues( 1 ) == Approx( -1.0 ) );
449 REQUIRE( mx::math::isNan( eigenvalues( 2 ) ) );
450 REQUIRE( eigenvalues( 3 ) == Approx( 1e-20 ) );
451 REQUIRE( eigenvalues( 4 ) == Approx( 4.0 ) );
452}
453
454/** \brief Verifies calcEigenVecs rejects invalid covariance geometry and maps query/solve failures to -1.
455 *
456 * \ingroup eigenLapack_unit_tests
457 */
458TEST_CASE( "calcEigenVecs rejects invalid inputs and solver failures", "[math::eigenLapack][calcEigenVecs]" )
459{
460 hookReset reset;
461 matrixT eigenvectors;
462 matrixT eigenvalues;
463
464 matrixT rectangular( 2, 3 );
465 rectangular.setZero();
466 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, rectangular ) == -1 );
467
468 matrixT empty( 0, 0 );
469 mx::math::syevrMem<double> workspace;
470 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, empty, 0, false, false, &workspace ) == -1 );
471
472 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
473 queryInfo = 29;
474 matrixT covariance = diagonalMatrix( { 1.0, 2.0 } );
475 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance ) == -1 );
476
477 resetControls();
478 mx::math::detail::eigenLapackTestHooks<double>::solver = &controlledSyevr;
479 solveInfo = 31;
480 REQUIRE( mx::math::calcEigenVecs( eigenvectors, eigenvalues, covariance, 0, false, false, &workspace ) == -1 );
481}
482
483/** \brief Verifies covariance construction and KL-mode calculation for a small reference matrix.
484 *
485 * \ingroup eigenLapack_unit_tests
486 */
487TEST_CASE( "eigenSYRK and calcKLModes form normalized KL modes", "[math::eigenLapack]" )
488{
489 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> references( 3, 2 );
490 references << 1.0, 0.0, 0.0, 2.0, 0.0, 0.0;
491
492 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> covariance;
493 mx::math::eigenSYRK( covariance, references );
494
495 REQUIRE( covariance.rows() == 2 );
496 REQUIRE( covariance.cols() == 2 );
497 REQUIRE( covariance( 0, 0 ) == Approx( 1.0 ) );
498 REQUIRE( covariance( 1, 0 ) == Approx( 0.0 ) );
499 REQUIRE( covariance( 1, 1 ) == Approx( 4.0 ) );
500
501 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> modes;
502 double eigenTime{ -1 };
503 double modeTime{ -1 };
504 REQUIRE( mx::math::calcKLModes<double>( modes, covariance, references, 0, nullptr, &eigenTime, &modeTime ) == 0 );
505 REQUIRE( modes.rows() == 2 );
506 REQUIRE( modes.cols() == 3 );
507 REQUIRE( modes.matrix().row( 0 ).norm() == Approx( 1.0 ).margin( 1e-12 ) );
508 REQUIRE( modes.matrix().row( 1 ).norm() == Approx( 1.0 ).margin( 1e-12 ) );
509 REQUIRE( modes.matrix().row( 0 ).dot( modes.matrix().row( 1 ) ) == Approx( 0.0 ).margin( 1e-12 ) );
510 REQUIRE( eigenTime >= 0 );
511 REQUIRE( modeTime >= 0 );
512
513 mx::math::syevrMem<double> workspace;
514 REQUIRE( mx::math::calcKLModes( modes, covariance, references, 1, &workspace ) == 0 );
515 REQUIRE( modes.rows() == 1 );
516 REQUIRE( modes.cols() == 3 );
517 REQUIRE( modes.matrix().row( 0 ).norm() == Approx( 1.0 ).margin( 1e-12 ) );
518
519 const auto *work = workspace.work;
520 const auto *iWork = workspace.iWork;
521 REQUIRE( mx::math::calcKLModes( modes, covariance, references, 1, &workspace ) == 0 );
522 REQUIRE( workspace.work == work );
523 REQUIRE( workspace.iWork == iWork );
524}
525
526/** \brief Verifies that calcKLModes rejects incompatible and non-square covariance matrices without leaking workspaces.
527 *
528 * \ingroup eigenLapack_unit_tests
529 */
530TEST_CASE( "calcKLModes rejects invalid covariance geometry", "[math::eigenLapack]" )
531{
532 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> modes;
533
534 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> covariance( 2, 2 );
535 covariance.matrix().setIdentity();
536 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> wrongReferenceCount( 3, 3 );
537 wrongReferenceCount.setZero();
538 REQUIRE( mx::math::calcKLModes( modes, covariance, wrongReferenceCount ) == -1 );
539
540 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> nonSquareCovariance( 2, 3 );
541 nonSquareCovariance.setZero();
542 Eigen::Array<double, Eigen::Dynamic, Eigen::Dynamic> references( 3, 2 );
543 references.setZero();
544 REQUIRE( mx::math::calcKLModes( modes, nonSquareCovariance, references ) == -1 );
545}
546
547} // namespace unitTest::math_eigenLapack_test
eigenLapack.hpp
Interfaces to Lapack and BLAS for Eigen-like arrays.
TEST_CASE
TEST_CASE("Loading aoAtmosphere config settings", "[ao::analysis::aoAtmosphere]")
Verify parsing and validation of atmosphere configuration settings.
Definition aoAtmosphere_test.cpp:23
mx::math::calcKLModes
MXLAPACK_INT calcKLModes(eigenT &klModes, eigenT &cv, const eigenT1 &Rims, int n_modes=0, syevrMem< _evCalcT > *mem=0, double *t_eigenv=nullptr, double *t_klim=nullptr)
Calculate the K-L modes, or principle components, given a covariance matrix.
Definition eigenLapack.hpp:679
mx::math::calcEigenVecs
MXLAPACK_INT calcEigenVecs(eigenT &evecs, eigenT &evals, eigenT &cv, int nVecs=0, bool normalize=false, bool check=false, syevrMem< _evCalcT > *mem=0, double *t_eigenv=nullptr)
Calculate the eigenvectors and eigenvalues given a triangular matrix.
Definition eigenLapack.hpp:537
mx::math::eigenSYRK
void eigenSYRK(eigenT1 &cv, const eigenT2 &ims)
Calculates the lower triangular part of the covariance matrix of ims.
Definition eigenLapack.hpp:230
mx::math::eigenSYEVR
MXLAPACK_INT eigenSYEVR(arrT &eigvec, arrT &eigval, arrT &X, int ev0=0, int ev1=-1, char UPLO='L', syevrMem< typename arrT::Scalar > *mem=0)
Calculate select eigenvalues and eigenvectors of an Eigen Array.
Definition eigenLapack.hpp:360
mx::math::isNan
bool isNan(realT value)
Test whether a floating-point value is NaN, including under finite-math-only optimization.
Definition floatUtils.hpp:39
mx::math::syevrMem
A struct to hold the working memory for eigenSYEVR and maintain it between calls if desired.
Definition eigenLapack.hpp:255
mx::math::syevrMem::iSuppZ
MXLAPACK_INT * iSuppZ
LAPACK eigenvector-support workspace owned by this object.
Definition eigenLapack.hpp:290
mx::math::syevrMem::work
floatT * work
Optimized floating-point LAPACK workspace owned by this object.
Definition eigenLapack.hpp:296
mx::math::syevrMem::sizeISuppZ
MXLAPACK_INT sizeISuppZ
Capacity of the eigenvector-support workspace.
Definition eigenLapack.hpp:275
mx::math::syevrMem::free
void free()
Release all workspace allocations and reset the cached LAPACK configuration.
Definition eigenLapack.hpp:331
mx::math::syevrMem::n
MXLAPACK_INT n
Matrix order associated with the cached workspace query.
Definition eigenLapack.hpp:260
mx::math::syevrMem::sizeMinWork
MXLAPACK_INT sizeMinWork
Capacity of the minimum floating-point query workspace.
Definition eigenLapack.hpp:278
mx::math::syevrMem::sizeMinIWork
MXLAPACK_INT sizeMinIWork
Capacity of the minimum integer query workspace.
Definition eigenLapack.hpp:284
mx::math::syevrMem::iWork
MXLAPACK_INT * iWork
Optimized integer LAPACK workspace owned by this object.
Definition eigenLapack.hpp:302
mx::math::syevrMem::minIWork
MXLAPACK_INT * minIWork
Minimum integer workspace used for LAPACK size queries.
Definition eigenLapack.hpp:299
mx::math::syevrMem::sizeWork
MXLAPACK_INT sizeWork
Capacity of the optimized floating-point workspace.
Definition eigenLapack.hpp:281
mx::math::syevrMem::minWork
floatT * minWork
Minimum floating-point workspace used for LAPACK size queries.
Definition eigenLapack.hpp:293
mx::math::syevrMem::sizeIWork
MXLAPACK_INT sizeIWork
Capacity of the optimized integer workspace.
Definition eigenLapack.hpp:287