surge_sparse/
klu.rs

1// SPDX-License-Identifier: LicenseRef-PolyForm-Noncommercial-1.0.0
2//! SuiteSparse KLU sparse LU factorization wrapper.
3
4use std::ffi::c_int;
5use std::os::raw::c_void;
6use std::ptr;
7
8use crate::csc::{CscMatrix, try_usize_to_i32, validate_csc_pattern};
9use crate::error::{SparseError, SparseResult};
10
11const STRICT_RCOND_THRESHOLD: f64 = 1e-12;
12
13#[repr(C)]
14struct KluCommon {
15    tol: f64,
16    memgrow: f64,
17    initmem_amd: f64,
18    initmem: f64,
19    maxwork: f64,
20    btf: c_int,
21    ordering: c_int,
22    scale: c_int,
23    user_order:
24        Option<unsafe extern "C" fn(i32, *mut i32, *mut i32, *mut i32, *mut KluCommon) -> i32>,
25    user_data: *mut c_void,
26    halt_if_singular: c_int,
27    status: c_int,
28    nrealloc: c_int,
29    structural_rank: i32,
30    numerical_rank: i32,
31    singular_col: i32,
32    noffdiag: i32,
33    flops: f64,
34    rcond: f64,
35    condest: f64,
36    rgrowth: f64,
37    work: f64,
38    memusage: usize,
39    mempeak: usize,
40}
41
42enum KluSymbolic {}
43enum KluNumeric {}
44
45unsafe extern "C" {
46    fn klu_defaults(common: *mut KluCommon) -> c_int;
47    fn klu_analyze(
48        n: i32,
49        ap: *const i32,
50        ai: *const i32,
51        common: *mut KluCommon,
52    ) -> *mut KluSymbolic;
53    fn klu_factor(
54        ap: *const i32,
55        ai: *const i32,
56        ax: *const f64,
57        symbolic: *mut KluSymbolic,
58        common: *mut KluCommon,
59    ) -> *mut KluNumeric;
60    fn klu_refactor(
61        ap: *const i32,
62        ai: *const i32,
63        ax: *const f64,
64        symbolic: *mut KluSymbolic,
65        numeric: *mut KluNumeric,
66        common: *mut KluCommon,
67    ) -> c_int;
68    fn klu_solve(
69        symbolic: *mut KluSymbolic,
70        numeric: *mut KluNumeric,
71        ldim: i32,
72        nrhs: i32,
73        b: *mut f64,
74        common: *mut KluCommon,
75    ) -> c_int;
76    fn klu_tsolve(
77        symbolic: *mut KluSymbolic,
78        numeric: *mut KluNumeric,
79        ldim: i32,
80        nrhs: i32,
81        b: *mut f64,
82        common: *mut KluCommon,
83    ) -> c_int;
84    fn klu_free_symbolic(symbolic: *mut *mut KluSymbolic, common: *mut KluCommon) -> c_int;
85    fn klu_free_numeric(numeric: *mut *mut KluNumeric, common: *mut KluCommon) -> c_int;
86    fn klu_rcond(
87        symbolic: *mut KluSymbolic,
88        numeric: *mut KluNumeric,
89        common: *mut KluCommon,
90    ) -> c_int;
91}
92
93/// SuiteSparse KLU factorization for a fixed square CSC sparsity pattern.
94pub struct KluSolver {
95    common: KluCommon,
96    symbolic: *mut KluSymbolic,
97    numeric: *mut KluNumeric,
98    dim: i32,
99    nnz: usize,
100    col_ptrs: Vec<i32>,
101    row_indices: Vec<i32>,
102}
103
104unsafe impl Send for KluSolver {}
105
106impl KluSolver {
107    /// Analyze a square CSC sparsity pattern.
108    pub fn new(dim: usize, col_ptrs: &[usize], row_indices: &[usize]) -> SparseResult<Self> {
109        validate_csc_pattern(dim, dim, col_ptrs, row_indices)?;
110        Self::analyze(dim, col_ptrs, row_indices)
111    }
112
113    /// Analyze the sparsity pattern from a validated CSC matrix.
114    pub fn from_csc<T>(matrix: &CscMatrix<T>) -> SparseResult<Self> {
115        if !matrix.is_square() {
116            return Err(SparseError::MatrixNotSquare {
117                nrows: matrix.nrows(),
118                ncols: matrix.ncols(),
119            });
120        }
121        Self::analyze(matrix.nrows(), matrix.col_ptrs(), matrix.row_indices())
122    }
123
124    fn analyze(dim: usize, col_ptrs: &[usize], row_indices: &[usize]) -> SparseResult<Self> {
125        if dim == 0 {
126            return Err(SparseError::EmptyMatrix);
127        }
128
129        let dim = try_usize_to_i32("matrix dimension", dim)?;
130        let col_ptrs = col_ptrs
131            .iter()
132            .copied()
133            .map(|value| try_usize_to_i32("column pointer", value))
134            .collect::<SparseResult<Vec<_>>>()?;
135        let row_indices = row_indices
136            .iter()
137            .copied()
138            .map(|value| try_usize_to_i32("row index", value))
139            .collect::<SparseResult<Vec<_>>>()?;
140
141        let mut common = unsafe { std::mem::zeroed() };
142        unsafe {
143            klu_defaults(&mut common);
144        }
145
146        let symbolic =
147            unsafe { klu_analyze(dim, col_ptrs.as_ptr(), row_indices.as_ptr(), &mut common) };
148        if symbolic.is_null() {
149            return Err(SparseError::KluAnalyzeFailed);
150        }
151
152        Ok(Self {
153            common,
154            symbolic,
155            numeric: ptr::null_mut(),
156            dim,
157            nnz: row_indices.len(),
158            col_ptrs,
159            row_indices,
160        })
161    }
162
163    fn clear_numeric(&mut self) {
164        if !self.numeric.is_null() {
165            unsafe {
166                klu_free_numeric(&mut self.numeric, &mut self.common);
167            }
168        }
169    }
170
171    pub fn factor(&mut self, values: &[f64]) -> SparseResult<()> {
172        self.validate_values(values)?;
173
174        self.clear_numeric();
175
176        self.numeric = unsafe {
177            klu_factor(
178                self.col_ptrs.as_ptr(),
179                self.row_indices.as_ptr(),
180                values.as_ptr(),
181                self.symbolic,
182                &mut self.common,
183            )
184        };
185        if self.numeric.is_null() {
186            return Err(SparseError::KluFactorFailed);
187        }
188
189        if let Err(error) = self.finish_factorization(SparseError::KluFactorFailed) {
190            self.clear_numeric();
191            return Err(error);
192        }
193        Ok(())
194    }
195
196    pub fn refactor(&mut self, values: &[f64]) -> SparseResult<()> {
197        self.validate_values(values)?;
198        if self.numeric.is_null() {
199            return Err(SparseError::NotFactorized);
200        }
201
202        let ok = unsafe {
203            klu_refactor(
204                self.col_ptrs.as_ptr(),
205                self.row_indices.as_ptr(),
206                values.as_ptr(),
207                self.symbolic,
208                self.numeric,
209                &mut self.common,
210            )
211        };
212        if ok == 0 {
213            self.clear_numeric();
214            return Err(SparseError::KluRefactorFailed);
215        }
216
217        if let Err(error) = self.finish_factorization(SparseError::KluRefactorFailed) {
218            self.clear_numeric();
219            return Err(error);
220        }
221        Ok(())
222    }
223
224    pub fn solve(&mut self, rhs: &mut [f64]) -> SparseResult<()> {
225        self.validate_rhs(rhs)?;
226
227        let ok = unsafe {
228            klu_solve(
229                self.symbolic,
230                self.numeric,
231                self.dim,
232                1,
233                rhs.as_mut_ptr(),
234                &mut self.common,
235            )
236        };
237        if ok == 0 {
238            return Err(SparseError::KluSolveFailed);
239        }
240        Ok(())
241    }
242
243    pub fn solve_transpose(&mut self, rhs: &mut [f64]) -> SparseResult<()> {
244        self.validate_rhs(rhs)?;
245
246        let ok = unsafe {
247            klu_tsolve(
248                self.symbolic,
249                self.numeric,
250                self.dim,
251                1,
252                rhs.as_mut_ptr(),
253                &mut self.common,
254            )
255        };
256        if ok == 0 {
257            return Err(SparseError::KluSolveFailed);
258        }
259        Ok(())
260    }
261
262    pub fn rcond(&self) -> f64 {
263        self.common.rcond
264    }
265
266    fn validate_values(&self, values: &[f64]) -> SparseResult<()> {
267        if values.len() != self.nnz {
268            return Err(SparseError::ValueCountMismatch {
269                expected: self.nnz,
270                found: values.len(),
271            });
272        }
273        Ok(())
274    }
275
276    fn validate_rhs(&self, rhs: &[f64]) -> SparseResult<()> {
277        if self.numeric.is_null() {
278            return Err(SparseError::NotFactorized);
279        }
280        if rhs.len() != self.dim as usize {
281            return Err(SparseError::RhsLengthMismatch {
282                expected: self.dim as usize,
283                found: rhs.len(),
284            });
285        }
286        Ok(())
287    }
288
289    fn finish_factorization(&mut self, ill_conditioned_error: SparseError) -> SparseResult<()> {
290        let ok = unsafe { klu_rcond(self.symbolic, self.numeric, &mut self.common) };
291        if ok == 0 {
292            return Err(SparseError::KluRcondFailed);
293        }
294        let rcond = self.common.rcond;
295        if !rcond.is_finite() || rcond < STRICT_RCOND_THRESHOLD {
296            return Err(match ill_conditioned_error {
297                SparseError::KluFactorFailed => SparseError::KluIllConditioned {
298                    rcond,
299                    threshold: STRICT_RCOND_THRESHOLD,
300                },
301                SparseError::KluRefactorFailed => SparseError::KluIllConditioned {
302                    rcond,
303                    threshold: STRICT_RCOND_THRESHOLD,
304                },
305                other => other,
306            });
307        }
308        Ok(())
309    }
310
311    /// Solve AX = B in-place for multiple RHS columns.
312    ///
313    /// `rhs` must contain `dim * nrhs` values laid out column-major, so
314    /// column `j` occupies `rhs[j * dim .. (j + 1) * dim]`.
315    pub fn solve_many(&mut self, rhs: &mut [f64], nrhs: usize) -> SparseResult<()> {
316        if self.numeric.is_null() {
317            return Err(SparseError::NotFactorized);
318        }
319        let expected = self.dim as usize * nrhs;
320        if rhs.len() != expected {
321            return Err(SparseError::RhsLengthMismatch {
322                expected,
323                found: rhs.len(),
324            });
325        }
326        if nrhs == 0 {
327            return Ok(());
328        }
329        let nrhs_i32 = try_usize_to_i32("nrhs", nrhs)?;
330        let ok = unsafe {
331            klu_solve(
332                self.symbolic,
333                self.numeric,
334                self.dim,
335                nrhs_i32,
336                rhs.as_mut_ptr(),
337                &mut self.common,
338            )
339        };
340        if ok == 0 {
341            return Err(SparseError::KluSolveFailed);
342        }
343        Ok(())
344    }
345}
346
347impl Drop for KluSolver {
348    fn drop(&mut self) {
349        unsafe {
350            if !self.numeric.is_null() {
351                klu_free_numeric(&mut self.numeric, &mut self.common);
352            }
353            if !self.symbolic.is_null() {
354                klu_free_symbolic(&mut self.symbolic, &mut self.common);
355            }
356        }
357    }
358}
359
360#[cfg(test)]
361mod tests {
362    use super::*;
363    use crate::{CscMatrix, Triplet};
364
365    /// 2x2 identity: [[1,0],[0,1]] * x = [3,7] → x = [3,7]
366    #[test]
367    fn klu_identity_2x2() {
368        let mat = CscMatrix::try_from_triplets(
369            2,
370            2,
371            &[
372                Triplet {
373                    row: 0,
374                    col: 0,
375                    val: 1.0,
376                },
377                Triplet {
378                    row: 1,
379                    col: 1,
380                    val: 1.0,
381                },
382            ],
383        )
384        .unwrap();
385
386        let mut solver = KluSolver::from_csc(&mat).unwrap();
387        solver.factor(mat.values()).unwrap();
388
389        let mut rhs = vec![3.0, 7.0];
390        solver.solve(&mut rhs).unwrap();
391        assert!((rhs[0] - 3.0).abs() < 1e-14);
392        assert!((rhs[1] - 7.0).abs() < 1e-14);
393    }
394
395    /// 3x3 lower-triangular: [[2,0,0],[1,3,0],[0,4,5]] * x = [6,10,29] → x = [3,7/3,...]
396    #[test]
397    fn klu_lower_triangular() {
398        let triplets = vec![
399            Triplet {
400                row: 0,
401                col: 0,
402                val: 2.0,
403            },
404            Triplet {
405                row: 1,
406                col: 0,
407                val: 1.0,
408            },
409            Triplet {
410                row: 1,
411                col: 1,
412                val: 3.0,
413            },
414            Triplet {
415                row: 2,
416                col: 1,
417                val: 4.0,
418            },
419            Triplet {
420                row: 2,
421                col: 2,
422                val: 5.0,
423            },
424        ];
425        let mat = CscMatrix::try_from_triplets(3, 3, &triplets).unwrap();
426        let mut solver = KluSolver::from_csc(&mat).unwrap();
427        solver.factor(mat.values()).unwrap();
428
429        // Solve Ax = b where b = A * [1, 2, 3]
430        let x_exact = [1.0, 2.0, 3.0];
431        let mut rhs = vec![2.0 * 1.0, 1.0 * 1.0 + 3.0 * 2.0, 4.0 * 2.0 + 5.0 * 3.0];
432        solver.solve(&mut rhs).unwrap();
433        for i in 0..3 {
434            assert!(
435                (rhs[i] - x_exact[i]).abs() < 1e-12,
436                "x[{i}] = {} expected {}",
437                rhs[i],
438                x_exact[i],
439            );
440        }
441    }
442
443    /// Sparse 4x4 SPD matrix solve.
444    #[test]
445    fn klu_sparse_4x4() {
446        // Tridiagonal: [[4,-1,0,0],[-1,4,-1,0],[0,-1,4,-1],[0,0,-1,4]]
447        let triplets = vec![
448            Triplet {
449                row: 0,
450                col: 0,
451                val: 4.0,
452            },
453            Triplet {
454                row: 1,
455                col: 0,
456                val: -1.0,
457            },
458            Triplet {
459                row: 0,
460                col: 1,
461                val: -1.0,
462            },
463            Triplet {
464                row: 1,
465                col: 1,
466                val: 4.0,
467            },
468            Triplet {
469                row: 2,
470                col: 1,
471                val: -1.0,
472            },
473            Triplet {
474                row: 1,
475                col: 2,
476                val: -1.0,
477            },
478            Triplet {
479                row: 2,
480                col: 2,
481                val: 4.0,
482            },
483            Triplet {
484                row: 3,
485                col: 2,
486                val: -1.0,
487            },
488            Triplet {
489                row: 2,
490                col: 3,
491                val: -1.0,
492            },
493            Triplet {
494                row: 3,
495                col: 3,
496                val: 4.0,
497            },
498        ];
499        let mat = CscMatrix::try_from_triplets(4, 4, &triplets).unwrap();
500        let mut solver = KluSolver::from_csc(&mat).unwrap();
501        solver.factor(mat.values()).unwrap();
502
503        let mut rhs = vec![1.0, 0.0, 0.0, 1.0];
504        solver.solve(&mut rhs).unwrap();
505
506        // Verify A * x ≈ b (original rhs)
507        let b_orig = [1.0, 0.0, 0.0, 1.0];
508        let ax = [
509            4.0 * rhs[0] - rhs[1],
510            -rhs[0] + 4.0 * rhs[1] - rhs[2],
511            -rhs[1] + 4.0 * rhs[2] - rhs[3],
512            -rhs[2] + 4.0 * rhs[3],
513        ];
514        for i in 0..4 {
515            assert!(
516                (ax[i] - b_orig[i]).abs() < 1e-12,
517                "residual[{i}] = {:.2e}",
518                (ax[i] - b_orig[i]).abs(),
519            );
520        }
521    }
522
523    /// Refactoring with new values on the same pattern.
524    #[test]
525    fn klu_refactor_same_pattern() {
526        let triplets = vec![
527            Triplet {
528                row: 0,
529                col: 0,
530                val: 2.0,
531            },
532            Triplet {
533                row: 1,
534                col: 1,
535                val: 3.0,
536            },
537        ];
538        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
539        let mut solver = KluSolver::from_csc(&mat).unwrap();
540        solver.factor(mat.values()).unwrap();
541
542        let mut rhs = vec![4.0, 9.0];
543        solver.solve(&mut rhs).unwrap();
544        assert!((rhs[0] - 2.0).abs() < 1e-14);
545        assert!((rhs[1] - 3.0).abs() < 1e-14);
546
547        // Refactor with different values: [[5,0],[0,10]]
548        solver.refactor(&[5.0, 10.0]).unwrap();
549        let mut rhs2 = vec![15.0, 30.0];
550        solver.solve(&mut rhs2).unwrap();
551        assert!((rhs2[0] - 3.0).abs() < 1e-14);
552        assert!((rhs2[1] - 3.0).abs() < 1e-14);
553    }
554
555    /// Transpose solve: A^T x = b.
556    #[test]
557    fn klu_transpose_solve() {
558        let triplets = vec![
559            Triplet {
560                row: 0,
561                col: 0,
562                val: 1.0,
563            },
564            Triplet {
565                row: 1,
566                col: 0,
567                val: 2.0,
568            },
569            Triplet {
570                row: 1,
571                col: 1,
572                val: 3.0,
573            },
574        ];
575        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
576        let mut solver = KluSolver::from_csc(&mat).unwrap();
577        solver.factor(mat.values()).unwrap();
578
579        // A^T = [[1,2],[0,3]], A^T * x = [5, 6] → x = [1, 2]
580        let mut rhs = vec![5.0, 6.0];
581        solver.solve_transpose(&mut rhs).unwrap();
582        assert!((rhs[0] - 1.0).abs() < 1e-14);
583        assert!((rhs[1] - 2.0).abs() < 1e-14);
584    }
585
586    /// Multiple RHS solve.
587    #[test]
588    fn klu_solve_many() {
589        let triplets = vec![
590            Triplet {
591                row: 0,
592                col: 0,
593                val: 2.0,
594            },
595            Triplet {
596                row: 1,
597                col: 1,
598                val: 5.0,
599            },
600        ];
601        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
602        let mut solver = KluSolver::from_csc(&mat).unwrap();
603        solver.factor(mat.values()).unwrap();
604
605        // Two RHS stored column-major: [4,10, 6,15]
606        let mut rhs = vec![4.0, 10.0, 6.0, 15.0];
607        solver.solve_many(&mut rhs, 2).unwrap();
608        assert!((rhs[0] - 2.0).abs() < 1e-14); // 4/2
609        assert!((rhs[1] - 2.0).abs() < 1e-14); // 10/5
610        assert!((rhs[2] - 3.0).abs() < 1e-14); // 6/2
611        assert!((rhs[3] - 3.0).abs() < 1e-14); // 15/5
612    }
613
614    /// Error: solve before factoring.
615    #[test]
616    fn klu_solve_before_factor_errors() {
617        let triplets = vec![
618            Triplet {
619                row: 0,
620                col: 0,
621                val: 1.0,
622            },
623            Triplet {
624                row: 1,
625                col: 1,
626                val: 1.0,
627            },
628        ];
629        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
630        let mut solver = KluSolver::from_csc(&mat).unwrap();
631
632        let mut rhs = vec![1.0, 2.0];
633        let err = solver.solve(&mut rhs).unwrap_err();
634        assert!(matches!(err, SparseError::NotFactorized));
635    }
636
637    /// Error: wrong RHS length.
638    #[test]
639    fn klu_rhs_length_mismatch() {
640        let triplets = vec![
641            Triplet {
642                row: 0,
643                col: 0,
644                val: 1.0,
645            },
646            Triplet {
647                row: 1,
648                col: 1,
649                val: 1.0,
650            },
651        ];
652        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
653        let mut solver = KluSolver::from_csc(&mat).unwrap();
654        solver.factor(mat.values()).unwrap();
655
656        let mut rhs = vec![1.0, 2.0, 3.0]; // wrong length
657        let err = solver.solve(&mut rhs).unwrap_err();
658        assert!(matches!(err, SparseError::RhsLengthMismatch { .. }));
659    }
660
661    /// COO to CSC with duplicate entries are summed.
662    #[test]
663    fn csc_triplet_duplicates_summed() {
664        let triplets = vec![
665            Triplet {
666                row: 0,
667                col: 0,
668                val: 1.0,
669            },
670            Triplet {
671                row: 0,
672                col: 0,
673                val: 2.0,
674            },
675            Triplet {
676                row: 1,
677                col: 1,
678                val: 5.0,
679            },
680        ];
681        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
682        assert_eq!(mat.nnz(), 2); // duplicates merged
683        assert!((mat.values()[0] - 3.0_f64).abs() < 1e-14); // 1 + 2
684    }
685
686    /// Large sparse identity: factor and solve.
687    #[test]
688    fn klu_large_identity() {
689        let n = 500;
690        let triplets: Vec<Triplet<f64>> = (0..n)
691            .map(|i| Triplet {
692                row: i,
693                col: i,
694                val: 1.0,
695            })
696            .collect();
697        let mat = CscMatrix::try_from_triplets(n, n, &triplets).unwrap();
698        let mut solver = KluSolver::from_csc(&mat).unwrap();
699        solver.factor(mat.values()).unwrap();
700
701        let mut rhs: Vec<f64> = (0..n).map(|i| i as f64).collect();
702        solver.solve(&mut rhs).unwrap();
703        for (i, val) in rhs.iter().enumerate() {
704            assert!((val - i as f64).abs() < 1e-12);
705        }
706    }
707
708    /// Structurally singular matrix: column 1 has no entries.
709    /// The pattern [[1,0],[0,0]] has an empty column so factor should fail.
710    #[test]
711    fn klu_structurally_singular_empty_column() {
712        // Only one entry at (0,0), column 1 is entirely empty.
713        let triplets = vec![Triplet {
714            row: 0,
715            col: 0,
716            val: 1.0,
717        }];
718        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
719        let mut solver = KluSolver::from_csc(&mat).unwrap();
720        let result = solver.factor(mat.values());
721        assert!(
722            result.is_err(),
723            "factoring a structurally singular matrix should fail"
724        );
725    }
726
727    /// Numerically singular matrix: [[1,1],[1,1]] has rank 1.
728    /// Strict factorization should reject the low reciprocal condition number.
729    #[test]
730    fn klu_numerically_singular_matrix() {
731        let triplets = vec![
732            Triplet {
733                row: 0,
734                col: 0,
735                val: 1.0,
736            },
737            Triplet {
738                row: 1,
739                col: 0,
740                val: 1.0,
741            },
742            Triplet {
743                row: 0,
744                col: 1,
745                val: 1.0,
746            },
747            Triplet {
748                row: 1,
749                col: 1,
750                val: 1.0,
751            },
752        ];
753        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
754        let mut solver = KluSolver::from_csc(&mat).unwrap();
755        let err = solver.factor(mat.values()).unwrap_err();
756        assert!(
757            matches!(
758                err,
759                SparseError::KluIllConditioned { .. } | SparseError::KluFactorFailed
760            ),
761            "singular matrix should be rejected, got {err:?}"
762        );
763    }
764
765    /// Ill-conditioned matrix should produce a very small rcond.
766    #[test]
767    fn klu_ill_conditioned_rcond() {
768        // Hilbert-like 3x3 matrix is ill-conditioned.
769        // [[1, 1/2, 1/3], [1/2, 1/3, 1/4], [1/3, 1/4, 1/5]]
770        let triplets = vec![
771            Triplet {
772                row: 0,
773                col: 0,
774                val: 1.0,
775            },
776            Triplet {
777                row: 1,
778                col: 0,
779                val: 0.5,
780            },
781            Triplet {
782                row: 2,
783                col: 0,
784                val: 1.0 / 3.0,
785            },
786            Triplet {
787                row: 0,
788                col: 1,
789                val: 0.5,
790            },
791            Triplet {
792                row: 1,
793                col: 1,
794                val: 1.0 / 3.0,
795            },
796            Triplet {
797                row: 2,
798                col: 1,
799                val: 0.25,
800            },
801            Triplet {
802                row: 0,
803                col: 2,
804                val: 1.0 / 3.0,
805            },
806            Triplet {
807                row: 1,
808                col: 2,
809                val: 0.25,
810            },
811            Triplet {
812                row: 2,
813                col: 2,
814                val: 0.2,
815            },
816        ];
817        let mat = CscMatrix::try_from_triplets(3, 3, &triplets).unwrap();
818        let mut solver = KluSolver::from_csc(&mat).unwrap();
819        solver.factor(mat.values()).unwrap();
820
821        // Hilbert matrices are notoriously ill-conditioned; rcond should be small.
822        assert!(
823            solver.rcond() < 0.1,
824            "rcond for a 3x3 Hilbert matrix should be small, got {}",
825            solver.rcond()
826        );
827        assert!(solver.rcond() > 0.0, "rcond should still be positive");
828    }
829
830    /// Refactor with values that become zero (diagonal zeros).
831    #[test]
832    fn klu_refactor_with_zero_values() {
833        // Start with a well-conditioned diagonal.
834        let triplets = vec![
835            Triplet {
836                row: 0,
837                col: 0,
838                val: 2.0,
839            },
840            Triplet {
841                row: 0,
842                col: 1,
843                val: 1.0,
844            },
845            Triplet {
846                row: 1,
847                col: 0,
848                val: 1.0,
849            },
850            Triplet {
851                row: 1,
852                col: 1,
853                val: 3.0,
854            },
855        ];
856        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
857        let mut solver = KluSolver::from_csc(&mat).unwrap();
858        solver.factor(mat.values()).unwrap();
859
860        // Verify the initial factorization works.
861        let mut rhs = vec![3.0, 4.0];
862        solver.solve(&mut rhs).unwrap();
863
864        // Refactor with values making the matrix singular: [[0,1],[1,0]] has det = -1 but
865        // refactor with [[1,1],[1,1]] to make it singular.
866        let singular_vals: Vec<f64> = vec![1.0, 1.0, 1.0, 1.0]; // [[1,1],[1,1]]
867        let err = solver.refactor(&singular_vals).unwrap_err();
868        assert!(
869            matches!(
870                err,
871                SparseError::KluIllConditioned { .. } | SparseError::KluRefactorFailed
872            ),
873            "refactoring to singular values should be rejected, got {err:?}"
874        );
875
876        let mut stale_rhs = vec![3.0, 4.0];
877        let solve_err = solver.solve(&mut stale_rhs).unwrap_err();
878        assert!(
879            matches!(solve_err, SparseError::NotFactorized),
880            "failed refactor must invalidate numeric factors, got {solve_err:?}"
881        );
882    }
883
884    /// solve_many with multiple RHS columns on a non-diagonal matrix.
885    #[test]
886    fn klu_solve_many_non_diagonal() {
887        // A = [[2, 1], [1, 3]]
888        let triplets = vec![
889            Triplet {
890                row: 0,
891                col: 0,
892                val: 2.0,
893            },
894            Triplet {
895                row: 1,
896                col: 0,
897                val: 1.0,
898            },
899            Triplet {
900                row: 0,
901                col: 1,
902                val: 1.0,
903            },
904            Triplet {
905                row: 1,
906                col: 1,
907                val: 3.0,
908            },
909        ];
910        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
911        let mut solver = KluSolver::from_csc(&mat).unwrap();
912        solver.factor(mat.values()).unwrap();
913
914        // Three RHS columns: b1 = A*[1,0], b2 = A*[0,1], b3 = A*[1,1]
915        // b1 = [2,1], b2 = [1,3], b3 = [3,4]
916        let mut rhs = vec![2.0, 1.0, 1.0, 3.0, 3.0, 4.0];
917        solver.solve_many(&mut rhs, 3).unwrap();
918        assert!((rhs[0] - 1.0).abs() < 1e-12, "col1[0]={}", rhs[0]);
919        assert!((rhs[1] - 0.0).abs() < 1e-12, "col1[1]={}", rhs[1]);
920        assert!((rhs[2] - 0.0).abs() < 1e-12, "col2[0]={}", rhs[2]);
921        assert!((rhs[3] - 1.0).abs() < 1e-12, "col2[1]={}", rhs[3]);
922        assert!((rhs[4] - 1.0).abs() < 1e-12, "col3[0]={}", rhs[4]);
923        assert!((rhs[5] - 1.0).abs() < 1e-12, "col3[1]={}", rhs[5]);
924    }
925
926    /// solve_many with zero RHS columns is a no-op.
927    #[test]
928    fn klu_solve_many_zero_columns() {
929        let triplets = vec![
930            Triplet {
931                row: 0,
932                col: 0,
933                val: 1.0,
934            },
935            Triplet {
936                row: 1,
937                col: 1,
938                val: 1.0,
939            },
940        ];
941        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
942        let mut solver = KluSolver::from_csc(&mat).unwrap();
943        solver.factor(mat.values()).unwrap();
944
945        let mut rhs = vec![];
946        solver.solve_many(&mut rhs, 0).unwrap();
947        assert!(rhs.is_empty());
948    }
949
950    /// solve_many rejects incorrect buffer length.
951    #[test]
952    fn klu_solve_many_wrong_length() {
953        let triplets = vec![
954            Triplet {
955                row: 0,
956                col: 0,
957                val: 1.0,
958            },
959            Triplet {
960                row: 1,
961                col: 1,
962                val: 1.0,
963            },
964        ];
965        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
966        let mut solver = KluSolver::from_csc(&mat).unwrap();
967        solver.factor(mat.values()).unwrap();
968
969        let mut rhs = vec![1.0, 2.0, 3.0]; // dim=2, nrhs=2 needs 4 elements
970        let err = solver.solve_many(&mut rhs, 2).unwrap_err();
971        assert!(matches!(err, SparseError::RhsLengthMismatch { .. }));
972    }
973
974    /// Wrong value count for factor should fail.
975    #[test]
976    fn klu_factor_wrong_value_count() {
977        let triplets = vec![
978            Triplet {
979                row: 0,
980                col: 0,
981                val: 1.0,
982            },
983            Triplet {
984                row: 1,
985                col: 1,
986                val: 1.0,
987            },
988        ];
989        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
990        let mut solver = KluSolver::from_csc(&mat).unwrap();
991
992        let err = solver.factor(&[1.0]).unwrap_err(); // expects 2 values
993        assert!(matches!(err, SparseError::ValueCountMismatch { .. }));
994    }
995
996    /// Refactor before any factor should fail.
997    #[test]
998    fn klu_refactor_before_factor() {
999        let triplets = vec![
1000            Triplet {
1001                row: 0,
1002                col: 0,
1003                val: 1.0,
1004            },
1005            Triplet {
1006                row: 1,
1007                col: 1,
1008                val: 1.0,
1009            },
1010        ];
1011        let mat = CscMatrix::try_from_triplets(2, 2, &triplets).unwrap();
1012        let mut solver = KluSolver::from_csc(&mat).unwrap();
1013
1014        let err = solver.refactor(&[2.0, 3.0]).unwrap_err();
1015        assert!(matches!(err, SparseError::NotFactorized));
1016    }
1017
1018    /// Constructing a KluSolver with a non-square matrix should fail.
1019    #[test]
1020    fn klu_rejects_non_square_from_csc() {
1021        let mat =
1022            CscMatrix::try_new(2, 3, vec![0, 1, 2, 3], vec![0, 1, 0], vec![1.0, 2.0, 3.0]).unwrap();
1023        let result = KluSolver::from_csc(&mat);
1024        assert!(matches!(result, Err(SparseError::MatrixNotSquare { .. })));
1025    }
1026
1027    /// Constructing a KluSolver with a 0x0 matrix should fail.
1028    #[test]
1029    fn klu_rejects_empty_matrix() {
1030        let mat = CscMatrix::<f64>::try_new(0, 0, vec![0], vec![], vec![]).unwrap();
1031        let result = KluSolver::from_csc(&mat);
1032        assert!(matches!(result, Err(SparseError::EmptyMatrix)));
1033    }
1034}
surge_sparse/klu.rs

surge_sparse/
klu.rs
