1#![forbid(unsafe_code)]
2
3use core::fmt;
6
7#[derive(Clone, Copy, Debug, PartialEq, Eq)]
16#[non_exhaustive]
17pub enum ArithmeticOperation {
18 MatrixInfinityNorm,
20 SymmetryCheck,
22 LuFactorization,
24 LdltFactorization,
26 LuSolve,
28 LdltSolve,
30 Determinant,
32 DeterminantErrorBound,
34 VectorDotProduct,
36 VectorSquaredNorm,
38}
39
40impl fmt::Display for ArithmeticOperation {
41 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
42 f.write_str(match self {
43 Self::MatrixInfinityNorm => "matrix infinity norm",
44 Self::SymmetryCheck => "symmetry check",
45 Self::LuFactorization => "LU factorization",
46 Self::LdltFactorization => "LDLT factorization",
47 Self::LuSolve => "LU solve",
48 Self::LdltSolve => "LDLT solve",
49 Self::Determinant => "determinant",
50 Self::DeterminantErrorBound => "determinant error bound",
51 Self::VectorDotProduct => "vector dot product",
52 Self::VectorSquaredNorm => "vector squared norm",
53 })
54 }
55}
56
57#[derive(Clone, Copy, Debug, PartialEq, Eq)]
66#[non_exhaustive]
67pub enum FactorizationKind {
68 Lu,
70 Ldlt,
72}
73
74impl fmt::Display for FactorizationKind {
75 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
76 f.write_str(match self {
77 Self::Lu => "LU",
78 Self::Ldlt => "LDLT",
79 })
80 }
81}
82
83#[derive(Clone, Copy, Debug, PartialEq, Eq)]
98#[non_exhaustive]
99pub enum InvalidToleranceReason {
100 Negative,
102 NotFinite,
104}
105
106#[derive(Clone, Copy, Debug, PartialEq, Eq)]
121#[non_exhaustive]
122pub enum NonFiniteLocation {
123 #[non_exhaustive]
125 MatrixCell {
126 row: usize,
128 col: usize,
130 },
131 #[non_exhaustive]
133 VectorEntry {
134 index: usize,
136 },
137 #[non_exhaustive]
139 Step {
140 index: usize,
142 },
143 Scalar,
145}
146
147#[derive(Clone, Copy, Debug, PartialEq, Eq)]
163#[non_exhaustive]
164pub enum NonFiniteOrigin {
165 Input,
167 #[non_exhaustive]
169 Computation {
170 operation: ArithmeticOperation,
172 },
173}
174
175#[derive(Clone, Copy, Debug, PartialEq)]
194#[non_exhaustive]
195pub enum PositiveSemidefiniteViolation {
196 #[non_exhaustive]
198 NegativePivot {
199 value: f64,
201 },
202 #[non_exhaustive]
204 ZeroPivotCoupling {
205 row: usize,
207 value: f64,
209 },
210}
211
212#[derive(Clone, Copy, Debug, PartialEq)]
226#[non_exhaustive]
227pub enum SingularityReason {
228 Exact,
230 #[non_exhaustive]
232 Numerical {
233 factorization: FactorizationKind,
235 pivot_magnitude: f64,
237 tolerance: f64,
239 },
240}
241
242#[derive(Clone, Copy, Debug, PartialEq, Eq)]
262#[non_exhaustive]
263pub enum UnrepresentableReason {
264 RequiresRounding,
267 NotFinite,
269}
270
271#[derive(Clone, Copy, Debug, PartialEq)]
290#[non_exhaustive]
291pub enum LaError {
292 #[non_exhaustive]
294 Singular {
295 pivot_col: usize,
297 reason: SingularityReason,
299 },
300 #[non_exhaustive]
302 NonFinite {
303 location: NonFiniteLocation,
305 origin: NonFiniteOrigin,
307 },
308 #[non_exhaustive]
310 Unrepresentable {
311 index: Option<usize>,
313 reason: UnrepresentableReason,
315 },
316 #[non_exhaustive]
318 DeterminantScaleOverflow {
319 dim: usize,
321 min_exponent: i32,
323 },
324 #[non_exhaustive]
326 UnsupportedDimension {
327 requested: usize,
329 max: usize,
331 },
332 #[non_exhaustive]
334 IndexOutOfBounds {
335 row: usize,
337 col: usize,
339 dim: usize,
341 },
342 #[non_exhaustive]
344 InvalidTolerance {
345 value: f64,
347 reason: InvalidToleranceReason,
349 },
350 #[non_exhaustive]
352 Asymmetric {
353 row: usize,
355 col: usize,
357 dim: usize,
359 upper: f64,
361 lower: f64,
363 allowed_abs_diff: f64,
365 },
366 #[non_exhaustive]
372 NotPositiveSemidefinite {
373 pivot_col: usize,
375 violation: PositiveSemidefiniteViolation,
377 },
378}
379
380impl LaError {
381 #[inline]
384 #[must_use]
385 pub const fn singular_exact(pivot_col: usize) -> Self {
386 Self::Singular {
387 pivot_col,
388 reason: SingularityReason::Exact,
389 }
390 }
391
392 #[inline]
396 #[must_use]
397 pub const fn singular_numerical(
398 pivot_col: usize,
399 factorization: FactorizationKind,
400 pivot_magnitude: f64,
401 tolerance: f64,
402 ) -> Self {
403 Self::Singular {
404 pivot_col,
405 reason: SingularityReason::Numerical {
406 factorization,
407 pivot_magnitude,
408 tolerance,
409 },
410 }
411 }
412
413 #[inline]
416 #[must_use]
417 pub const fn non_finite_input_matrix(row: usize, col: usize) -> Self {
418 Self::NonFinite {
419 location: NonFiniteLocation::MatrixCell { row, col },
420 origin: NonFiniteOrigin::Input,
421 }
422 }
423
424 #[inline]
427 #[must_use]
428 pub const fn non_finite_input_vector(index: usize) -> Self {
429 Self::NonFinite {
430 location: NonFiniteLocation::VectorEntry { index },
431 origin: NonFiniteOrigin::Input,
432 }
433 }
434
435 #[inline]
438 #[must_use]
439 pub const fn non_finite_input_scalar() -> Self {
440 Self::NonFinite {
441 location: NonFiniteLocation::Scalar,
442 origin: NonFiniteOrigin::Input,
443 }
444 }
445
446 #[inline]
449 #[must_use]
450 pub const fn non_finite_computation_matrix(
451 operation: ArithmeticOperation,
452 row: usize,
453 col: usize,
454 ) -> Self {
455 Self::NonFinite {
456 location: NonFiniteLocation::MatrixCell { row, col },
457 origin: NonFiniteOrigin::Computation { operation },
458 }
459 }
460
461 #[inline]
464 #[must_use]
465 pub const fn non_finite_computation_step(operation: ArithmeticOperation, index: usize) -> Self {
466 Self::NonFinite {
467 location: NonFiniteLocation::Step { index },
468 origin: NonFiniteOrigin::Computation { operation },
469 }
470 }
471
472 #[inline]
475 #[must_use]
476 pub const fn non_finite_computation_scalar(operation: ArithmeticOperation) -> Self {
477 Self::NonFinite {
478 location: NonFiniteLocation::Scalar,
479 origin: NonFiniteOrigin::Computation { operation },
480 }
481 }
482
483 #[inline]
486 #[must_use]
487 pub const fn unrepresentable(index: Option<usize>, reason: UnrepresentableReason) -> Self {
488 Self::Unrepresentable { index, reason }
489 }
490
491 #[inline]
494 #[must_use]
495 pub const fn unrepresentable_reason(&self) -> Option<UnrepresentableReason> {
496 match self {
497 Self::Unrepresentable { reason, .. } => Some(*reason),
498 _ => None,
499 }
500 }
501
502 #[inline]
505 #[must_use]
506 pub const fn requires_rounding(&self) -> bool {
507 matches!(
508 self,
509 Self::Unrepresentable {
510 reason: UnrepresentableReason::RequiresRounding,
511 ..
512 }
513 )
514 }
515
516 #[inline]
519 #[must_use]
520 pub const fn determinant_scale_overflow(dim: usize, min_exponent: i32) -> Self {
521 Self::DeterminantScaleOverflow { dim, min_exponent }
522 }
523
524 #[inline]
527 #[must_use]
528 pub const fn unsupported_dimension(requested: usize, max: usize) -> Self {
529 Self::UnsupportedDimension { requested, max }
530 }
531
532 #[inline]
535 #[must_use]
536 pub const fn index_out_of_bounds(row: usize, col: usize, dim: usize) -> Self {
537 Self::IndexOutOfBounds { row, col, dim }
538 }
539
540 #[inline]
548 #[must_use]
549 pub const fn invalid_tolerance(value: f64) -> Self {
550 let reason = if value.is_finite() {
551 InvalidToleranceReason::Negative
552 } else {
553 InvalidToleranceReason::NotFinite
554 };
555 Self::InvalidTolerance { value, reason }
556 }
557
558 #[inline]
562 #[must_use]
563 pub const fn asymmetric(
564 row: usize,
565 col: usize,
566 dim: usize,
567 upper: f64,
568 lower: f64,
569 allowed_abs_diff: f64,
570 ) -> Self {
571 Self::Asymmetric {
572 row,
573 col,
574 dim,
575 upper,
576 lower,
577 allowed_abs_diff,
578 }
579 }
580
581 #[inline]
584 #[must_use]
585 pub const fn not_positive_semidefinite_negative(pivot_col: usize, value: f64) -> Self {
586 Self::NotPositiveSemidefinite {
587 pivot_col,
588 violation: PositiveSemidefiniteViolation::NegativePivot { value },
589 }
590 }
591
592 #[inline]
596 #[must_use]
597 pub const fn not_positive_semidefinite_zero_coupling(
598 pivot_col: usize,
599 row: usize,
600 value: f64,
601 ) -> Self {
602 Self::NotPositiveSemidefinite {
603 pivot_col,
604 violation: PositiveSemidefiniteViolation::ZeroPivotCoupling { row, value },
605 }
606 }
607}
608
609fn write_non_finite_location(
613 f: &mut fmt::Formatter<'_>,
614 location: NonFiniteLocation,
615) -> fmt::Result {
616 match location {
617 NonFiniteLocation::MatrixCell { row, col } => {
618 write!(f, "matrix cell ({row}, {col})")
619 }
620 NonFiniteLocation::VectorEntry { index } => write!(f, "vector entry {index}"),
621 NonFiniteLocation::Step { index } => write!(f, "step {index}"),
622 NonFiniteLocation::Scalar => f.write_str("scalar value"),
623 }
624}
625
626fn write_non_finite(
629 f: &mut fmt::Formatter<'_>,
630 location: NonFiniteLocation,
631 origin: NonFiniteOrigin,
632) -> fmt::Result {
633 match (location, origin) {
634 (NonFiniteLocation::Scalar, NonFiniteOrigin::Input) => {
635 f.write_str("non-finite scalar input")
636 }
637 (NonFiniteLocation::Scalar, NonFiniteOrigin::Computation { operation }) => {
638 write!(f, "non-finite scalar result computed during {operation}")
639 }
640 (location, NonFiniteOrigin::Input) => {
641 f.write_str("non-finite input value at ")?;
642 write_non_finite_location(f, location)
643 }
644 (location, NonFiniteOrigin::Computation { operation }) => {
645 write!(f, "non-finite value computed during {operation} at ")?;
646 write_non_finite_location(f, location)
647 }
648 }
649}
650
651impl fmt::Display for LaError {
652 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
653 match *self {
654 Self::Singular {
655 pivot_col,
656 reason: SingularityReason::Exact,
657 } => write!(f, "matrix is exactly singular at pivot column {pivot_col}"),
658 Self::Singular {
659 pivot_col,
660 reason:
661 SingularityReason::Numerical {
662 factorization,
663 pivot_magnitude,
664 tolerance,
665 },
666 } => write!(
667 f,
668 "matrix is numerically singular during {factorization} factorization at pivot column {pivot_col}: pivot magnitude {pivot_magnitude} <= tolerance {tolerance}"
669 ),
670 Self::NonFinite { location, origin } => write_non_finite(f, location, origin),
671 Self::Unrepresentable {
672 index: Some(index),
673 reason: UnrepresentableReason::RequiresRounding,
674 } => write!(
675 f,
676 "exact result requires rounding to fit finite f64 at index {index}"
677 ),
678 Self::Unrepresentable {
679 index: None,
680 reason: UnrepresentableReason::RequiresRounding,
681 } => f.write_str("exact result requires rounding to fit finite f64"),
682 Self::Unrepresentable {
683 index: Some(index),
684 reason: UnrepresentableReason::NotFinite,
685 } => write!(
686 f,
687 "exact result has no finite f64 representation after rounding at index {index}"
688 ),
689 Self::Unrepresentable {
690 index: None,
691 reason: UnrepresentableReason::NotFinite,
692 } => f.write_str("exact result has no finite f64 representation after rounding"),
693 Self::DeterminantScaleOverflow { dim, min_exponent } => write!(
694 f,
695 "exact determinant scale exponent overflows for dimension {dim} with minimum entry exponent {min_exponent}"
696 ),
697 Self::UnsupportedDimension { requested, max } => write!(
698 f,
699 "unsupported matrix dimension {requested}; maximum stack-dispatch dimension is {max}"
700 ),
701 Self::IndexOutOfBounds { row, col, dim } => write!(
702 f,
703 "matrix index ({row}, {col}) is out of bounds for dimension {dim}"
704 ),
705 Self::InvalidTolerance {
706 value,
707 reason: InvalidToleranceReason::Negative,
708 } => write!(f, "invalid tolerance {value}; expected value >= 0"),
709 Self::InvalidTolerance {
710 value,
711 reason: InvalidToleranceReason::NotFinite,
712 } => write!(f, "invalid tolerance {value}; expected a finite value"),
713 Self::Asymmetric {
714 row,
715 col,
716 dim,
717 upper,
718 lower,
719 allowed_abs_diff,
720 } => write!(
721 f,
722 "matrix is not symmetric for dimension {dim}: entry ({row}, {col}) = {upper} and entry ({col}, {row}) = {lower} differ by more than allowed absolute difference {allowed_abs_diff}"
723 ),
724 Self::NotPositiveSemidefinite {
725 pivot_col,
726 violation: PositiveSemidefiniteViolation::NegativePivot { value },
727 } => write!(
728 f,
729 "LDLT rejected the matrix at pivot column {pivot_col}: computed diagonal value {value} < 0"
730 ),
731 Self::NotPositiveSemidefinite {
732 pivot_col,
733 violation: PositiveSemidefiniteViolation::ZeroPivotCoupling { row, value },
734 } => write!(
735 f,
736 "LDLT rejected the matrix at pivot column {pivot_col}: computed zero diagonal has non-zero coupling at row {row} with value {value}"
737 ),
738 }
739 }
740}
741
742impl std::error::Error for LaError {}
743
744#[cfg(test)]
745mod tests {
746 use std::error::Error;
747
748 use super::*;
749 use crate::MAX_STACK_MATRIX_DISPATCH_DIM;
750
751 #[test]
752 fn category_displays_are_concise() {
753 assert_eq!(FactorizationKind::Lu.to_string(), "LU");
754 assert_eq!(FactorizationKind::Ldlt.to_string(), "LDLT");
755 assert_eq!(
756 ArithmeticOperation::MatrixInfinityNorm.to_string(),
757 "matrix infinity norm"
758 );
759 assert_eq!(
760 ArithmeticOperation::SymmetryCheck.to_string(),
761 "symmetry check"
762 );
763 assert_eq!(
764 ArithmeticOperation::LuFactorization.to_string(),
765 "LU factorization"
766 );
767 assert_eq!(
768 ArithmeticOperation::LdltFactorization.to_string(),
769 "LDLT factorization"
770 );
771 assert_eq!(ArithmeticOperation::LuSolve.to_string(), "LU solve");
772 assert_eq!(ArithmeticOperation::LdltSolve.to_string(), "LDLT solve");
773 assert_eq!(ArithmeticOperation::Determinant.to_string(), "determinant");
774 assert_eq!(
775 ArithmeticOperation::DeterminantErrorBound.to_string(),
776 "determinant error bound"
777 );
778 assert_eq!(
779 ArithmeticOperation::VectorDotProduct.to_string(),
780 "vector dot product"
781 );
782 assert_eq!(
783 ArithmeticOperation::VectorSquaredNorm.to_string(),
784 "vector squared norm"
785 );
786 }
787
788 #[test]
789 fn singular_constructors_and_displays_preserve_reason() {
790 let exact = LaError::singular_exact(3);
791 assert_eq!(
792 exact,
793 LaError::Singular {
794 pivot_col: 3,
795 reason: SingularityReason::Exact,
796 }
797 );
798 assert_eq!(
799 exact.to_string(),
800 "matrix is exactly singular at pivot column 3"
801 );
802
803 let numerical = LaError::singular_numerical(2, FactorizationKind::Lu, 1e-14, 1e-12);
804 assert_eq!(
805 numerical,
806 LaError::Singular {
807 pivot_col: 2,
808 reason: SingularityReason::Numerical {
809 factorization: FactorizationKind::Lu,
810 pivot_magnitude: 1e-14,
811 tolerance: 1e-12,
812 },
813 }
814 );
815 assert_eq!(
816 numerical.to_string(),
817 "matrix is numerically singular during LU factorization at pivot column 2: pivot magnitude 0.00000000000001 <= tolerance 0.000000000001"
818 );
819 }
820
821 #[test]
822 fn non_finite_constructors_preserve_location_and_origin() {
823 assert_eq!(
824 LaError::non_finite_input_matrix(1, 2),
825 LaError::NonFinite {
826 location: NonFiniteLocation::MatrixCell { row: 1, col: 2 },
827 origin: NonFiniteOrigin::Input,
828 }
829 );
830 assert_eq!(
831 LaError::non_finite_input_vector(3).to_string(),
832 "non-finite input value at vector entry 3"
833 );
834 assert_eq!(
835 LaError::non_finite_input_scalar().to_string(),
836 "non-finite scalar input"
837 );
838 assert_eq!(
839 LaError::non_finite_computation_matrix(ArithmeticOperation::LuFactorization, 2, 1)
840 .to_string(),
841 "non-finite value computed during LU factorization at matrix cell (2, 1)"
842 );
843 assert_eq!(
844 LaError::non_finite_computation_step(ArithmeticOperation::LuSolve, 1).to_string(),
845 "non-finite value computed during LU solve at step 1"
846 );
847 assert_eq!(
848 LaError::non_finite_computation_scalar(ArithmeticOperation::Determinant).to_string(),
849 "non-finite scalar result computed during determinant"
850 );
851 }
852
853 #[test]
854 fn unrepresentable_helpers_preserve_recovery_reason() {
855 let rounding = LaError::unrepresentable(Some(2), UnrepresentableReason::RequiresRounding);
856 let scalar_rounding =
857 LaError::unrepresentable(None, UnrepresentableReason::RequiresRounding);
858 let indexed_not_finite =
859 LaError::unrepresentable(Some(2), UnrepresentableReason::NotFinite);
860 let not_finite = LaError::unrepresentable(None, UnrepresentableReason::NotFinite);
861 assert_eq!(
862 rounding.unrepresentable_reason(),
863 Some(UnrepresentableReason::RequiresRounding)
864 );
865 assert!(rounding.requires_rounding());
866 assert_eq!(
867 rounding.to_string(),
868 "exact result requires rounding to fit finite f64 at index 2"
869 );
870 assert_eq!(
871 scalar_rounding.to_string(),
872 "exact result requires rounding to fit finite f64"
873 );
874 assert_eq!(
875 indexed_not_finite.to_string(),
876 "exact result has no finite f64 representation after rounding at index 2"
877 );
878 assert_eq!(
879 not_finite.to_string(),
880 "exact result has no finite f64 representation after rounding"
881 );
882 assert!(!not_finite.requires_rounding());
883 assert_eq!(LaError::singular_exact(0).unrepresentable_reason(), None);
884 }
885
886 #[test]
887 fn invalid_tolerance_classifies_non_finite_before_negative() {
888 assert_eq!(
889 LaError::invalid_tolerance(-1.0),
890 LaError::InvalidTolerance {
891 value: -1.0,
892 reason: InvalidToleranceReason::Negative,
893 }
894 );
895 assert_eq!(
896 LaError::invalid_tolerance(f64::NEG_INFINITY),
897 LaError::InvalidTolerance {
898 value: f64::NEG_INFINITY,
899 reason: InvalidToleranceReason::NotFinite,
900 }
901 );
902 assert_eq!(
903 LaError::invalid_tolerance(-1.0).to_string(),
904 "invalid tolerance -1; expected value >= 0"
905 );
906 assert_eq!(
907 LaError::invalid_tolerance(f64::NEG_INFINITY).to_string(),
908 "invalid tolerance -inf; expected a finite value"
909 );
910 }
911
912 #[test]
913 fn asymmetric_error_retains_observed_values_and_bound() {
914 let err = LaError::asymmetric(0, 2, 3, 1.0, 1.5, 1e-12);
915 assert_eq!(
916 err,
917 LaError::Asymmetric {
918 row: 0,
919 col: 2,
920 dim: 3,
921 upper: 1.0,
922 lower: 1.5,
923 allowed_abs_diff: 1e-12,
924 }
925 );
926 assert_eq!(
927 err.to_string(),
928 "matrix is not symmetric for dimension 3: entry (0, 2) = 1 and entry (2, 0) = 1.5 differ by more than allowed absolute difference 0.000000000001"
929 );
930 }
931
932 #[test]
933 fn positive_semidefinite_errors_preserve_distinct_violations() {
934 assert_eq!(
935 LaError::not_positive_semidefinite_negative(1, -3.0).to_string(),
936 "LDLT rejected the matrix at pivot column 1: computed diagonal value -3 < 0"
937 );
938 assert_eq!(
939 LaError::not_positive_semidefinite_zero_coupling(0, 1, 2.0).to_string(),
940 "LDLT rejected the matrix at pivot column 0: computed zero diagonal has non-zero coupling at row 1 with value 2"
941 );
942 }
943
944 #[test]
945 fn remaining_helpers_and_displays_preserve_fields() {
946 assert_eq!(
947 LaError::determinant_scale_overflow(3, -1074).to_string(),
948 "exact determinant scale exponent overflows for dimension 3 with minimum entry exponent -1074"
949 );
950 assert_eq!(
951 LaError::unsupported_dimension(8, MAX_STACK_MATRIX_DISPATCH_DIM).to_string(),
952 "unsupported matrix dimension 8; maximum stack-dispatch dimension is 7"
953 );
954 assert_eq!(
955 LaError::index_out_of_bounds(3, 0, 3).to_string(),
956 "matrix index (3, 0) is out of bounds for dimension 3"
957 );
958 }
959
960 #[test]
961 fn is_std_error_with_no_source() {
962 let err = LaError::singular_exact(0);
963 let error: &dyn Error = &err;
964 assert!(error.source().is_none());
965 }
966}