1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
use sprs::errors::LinalgError;
use sprs::{CsStructureI, CsStructureViewI, PermOwnedI, SpIndex};
use suitesparse_camd_sys::*;
const MAX_INT32: usize = std::u32::MAX as usize;
pub fn try_camd<I, Iptr>(
mat: CsStructureViewI<I, Iptr>,
) -> Result<PermOwnedI<I>, LinalgError>
where
I: SpIndex,
Iptr: SpIndex,
{
let n = mat.rows();
if n != mat.cols() {
return Err(LinalgError::NonSquareMatrix);
}
let mut control = [0.; CAMD_CONTROL];
let mut info = [0.; CAMD_INFO];
let (camd_res, perm) = if n <= MAX_INT32 {
let constraint: *const SuiteSparseInt = std::ptr::null();
let mat: CsStructureI<SuiteSparseInt, SuiteSparseInt> =
mat.to_other_types();
let mut perm: Vec<SuiteSparseInt> = vec![0; n];
let camd_res = unsafe {
let indptr_proper = mat.proper_indptr();
camd_order(
n as SuiteSparseInt,
indptr_proper.as_ptr(),
mat.indices().as_ptr(),
perm.as_mut_ptr(),
control.as_mut_ptr(),
info.as_mut_ptr(),
constraint as *mut SuiteSparseInt,
) as isize
};
let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect();
(camd_res, perm)
} else {
let constraint: *const SuiteSparseLong = std::ptr::null();
let mat: CsStructureI<SuiteSparseLong, SuiteSparseLong> =
mat.to_other_types();
let mut perm: Vec<SuiteSparseLong> = vec![0; n];
let camd_res = unsafe {
let indptr_proper = mat.proper_indptr();
camd_l_order(
n as SuiteSparseLong,
indptr_proper.as_ptr(),
mat.indices().as_ptr(),
perm.as_mut_ptr(),
control.as_mut_ptr(),
info.as_mut_ptr(),
constraint as *mut SuiteSparseLong,
)
} as isize;
let perm = perm.iter().map(|&i| I::from_usize(i as usize)).collect();
(camd_res, perm)
};
if camd_res != CAMD_OK {
panic!("CsMat invariants have been broken");
}
Ok(PermOwnedI::new(perm))
}
pub fn camd<I, Iptr>(mat: CsStructureViewI<I, Iptr>) -> PermOwnedI<I>
where
I: SpIndex,
Iptr: SpIndex,
{
try_camd(mat).unwrap()
}
#[cfg(test)]
mod tests {
use sprs::CsMatI;
#[test]
fn try_camd() {
let mat = CsMatI::new_csc(
(4, 4),
vec![0, 2, 4, 6, 8],
vec![0, 3, 1, 2, 1, 2, 0, 3],
vec![1., 2., 21., 6., 6., 2., 2., 8.],
);
let res = super::try_camd(mat.structure_view());
assert!(res.is_ok());
}
}